notebooks/CVX_MVEE_algorithm.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Minimum Volume Ellipsoids を解くアルゴリズム\n",
    "\n",
    "参考\n",
    "* [Minimum Volume Ellipsoids](https://www.amazon.co.jp/-/en/Michael-J-Todd/dp/1611974372)\n",
    "* [Derivative of ln (det X)](https://statisticaloddsandends.wordpress.com/2018/05/24/derivative-of-ln-det-x/)\n",
    "* [ベクトルの内積とは？](https://risalc.info/src/inner-product.html#defR)\n",
    "* [安易に逆行列を数値計算するのはやめよう](https://opqrstuvcut.github.io/blog/posts/%E5%AE%89%E6%98%93%E3%81%AB%E9%80%86%E8%A1%8C%E5%88%97%E3%82%92%E6%95%B0%E5%80%A4%E8%A8%88%E7%AE%97%E3%81%99%E3%82%8B%E3%81%AE%E3%81%AF%E3%82%84%E3%82%81%E3%82%88%E3%81%86/)\n",
    "* [TFのrank-one Cholesky更新](https://www.tensorflow.org/probability/api_docs/python/tfp/math/cholesky_update)\n",
    "* [Wikipedia (Cholesky decomposition)](https://en.wikipedia.org/wiki/Cholesky_decomposition)\n",
    "\n",
    "表記\n",
    "* 行列の内積（[ベクトルの内積とは？](https://risalc.info/src/inner-product.html#defR)参照）：$\\langle A, B \\rangle := \\operatorname{Trace}(A^T B)$\n",
    "\n",
    "今回は[前回](CVX_MVEE_algorithm.ipynb)の続きです。\n",
    "次の主問題と双対問題を解くことを考えます。\n",
    "\n",
    "$$\n",
    "\\text{ 主問題 }\\;  \n",
    "\\begin{cases}\n",
    "&\\min_{H\\in S^n} f(H) := -\\ln \\det (H)\\\\\n",
    "&x_i^T H x_i \\leq n\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\text{ 双対問題 }\\; \n",
    "\\begin{cases}\n",
    "&\\max_{u\\in \\mathbb{R}^m} g(u) := \\ln \\det (XUX^T) \\\\\n",
    "& e^T u = 1, \\\\\n",
    "& u \\geq 0\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "この二つはどうやったら効率よく解けそうでしょうか？\n",
    "とりあえず二つの勾配について見てみましょう：\n",
    "\n",
    "**主問題**\n",
    "* $\\ln\\det$の微分は次の結果が有名です（[Derivative of ln (det X)](https://statisticaloddsandends.wordpress.com/2018/05/24/derivative-of-ln-det-x/)）：$\\nabla(-\\ln\\det)(H)=-H^{-1}$\n",
    "* $\\ln \\det(H)$の二回方向微分（導出は教科書のAppendix A.4を参照してください）：$E_1$の方向に方向微分した後、$E_2$の方向に方向微分すると、次の$E_1$と$E_2$についての関数が得られます\n",
    "    * $\\nabla^2 \\ln \\det(H)[E_1, E_2]=\\langle -(H^{-1}E_1 H^{-1}), E_2\\rangle$\n",
    "\n",
    "**双対問題**\n",
    "\n",
    "* $w(u):=\\nabla g(u) = (x_i^T H(u) x_i)^m_{i=1}$\n",
    "* $(\\nabla^2 g(u))_{i, j} = -(x_i^T H(u) x_j)^2, \\; i,j=1, \\dots, m$\n",
    "    * ここで、$H(u):=(XUX^T)^{-1}$です\n",
    "\n",
    "このことから、主問題でも双対問題でも二回微分を考えることはそう難しくないことがわかります。\n",
    "しかし、$n$が特に大きいときは、$n\\times n$の行列$H$の逆行列を各イテレーションで考える作業はコストが大きいですね。\n",
    "\n",
    "この問題は、双対問題の1次微分を使った座標上昇法を考えれば回避できます。$H(u)$を計算する際に逆行列が必要になりますが、座標上昇法みたいに変数を一つ変えるだけであれば、逆行列は効率よく計算できます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import jax.numpy as jnp\n",
    "import jax\n",
    "\n",
    "\n",
    "def log_det(X):\n",
    "    return -jnp.log(jnp.linalg.det(X))\n",
    "\n",
    "\n",
    "def mat_prod(A, B):\n",
    "    return jnp.trace(A.T @ B)\n",
    "\n",
    "\n",
    "n = 2\n",
    "L = jnp.array([[4., 0.], [2., 3.]])\n",
    "H = L @ L.T\n",
    "H_inv = jnp.linalg.inv(H)\n",
    "\n",
    "\n",
    "# X = Hでの一回微分\n",
    "jax_grad = jax.grad(log_det)(H)  # jaxによる勾配\n",
    "analytical_grad = -H_inv  # 解析解\n",
    "assert jax_grad.shape == H.shape"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 座標上昇法を使ったアルゴリズム\n",
    "\n",
    "表記\n",
    "* $i$：座標上昇法で選択した座標\n",
    "* $b\\in [m]$：適当なインデックス\n",
    "\n",
    "$u\\in \\mathbb{R}^m$を双対問題の実行可能解とします。つまり、$u$は非負かつ$e^Tu=1$、そして$g(u)=\\ln\\det(XUX^T)$は有限であり、$XUX^T$は正則であるとします。\n",
    "この$u$を使って、\n",
    "\n",
    "$$w:=w(u)=\\nabla g(u) = (x_b^T H(u) x_b)^m_{b=1}$$\n",
    "\n",
    "を持っているとしましょう。ここで、$H:=H(u):=(XUX^T)^{-1}$です。\n",
    "\n",
    "さて、この$u$について\n",
    "\n",
    "$$u_+:=(1-\\tau)u + \\tau e_i$$\n",
    "\n",
    "で更新することを考えます。ここで$\\tau\\neq 1$であり、$e_i$は$i$軸に沿った単位ベクトルです。$u=e_i$を代入すると、これも双対問題の実行可能解であることがわかります。\n",
    "これは次のように書き換えることも可能です：\n",
    "\n",
    "$$\\hat{u}:=u + \\lambda e_i$$\n",
    "\n",
    "ここで、$\\lambda:=\\tau / (1-\\tau)$です。つまり、$u_+=(1-\\tau)\\hat{u}$になっています。\n",
    "$-u_i \\leq \\lambda < \\infty$である限り、$u_+$は非負かつ$e_i$と一致しないことがわかります。\n",
    "\n",
    "このような更新を行うと、更新後の$XUX^T$に良い性質が生まれます。更新後の行列は\n",
    "\n",
    "$$XU_+X^T=(1+\\lambda)^{-1}(XUX^T + \\lambda x_i x_i^T)$$\n",
    "\n",
    "の形になります。この行列式と逆行列は次のように効率よく計算することができます。\n",
    "\n",
    "---\n",
    "\n",
    "$$\\det(XU_+X^T)=(1+\\lambda)^{-n}(1+\\lambda w_i)\\det (XUX^T)$$\n",
    "\n",
    "また、$1+\\lambda w_i$が正であれば、$XU_+X^T$は正則であり、\n",
    "$$(XU_+X^T)^{-1}=(1+\\lambda)\\left(H-\\frac{\\lambda}{1+\\lambda w_i}Hx_ix_i^TH\\right)$$\n",
    "\n",
    "---\n",
    "\n",
    "これより、更新後の目的関数もきれいに書くことができます。\n",
    "\n",
    "$$\n",
    "g(u_+) - g(u)=-n \\ln (1+\\lambda) + \\ln (1 + \\lambda w_i)\n",
    "$$\n",
    "\n",
    "です。\n",
    "\n",
    "この差分を最大化させるように$\\lambda$を決定すると、最適な更新ステップがわかりそうですね。実際、\n",
    "$$\\lambda^*:= \\frac{w_i-n}{(n-1)w_i}$$\n",
    "でこの更新量は最大になり、$\\lambda^*\\geq -u_i$である限り、$(n-1)\\ln\\left(\\frac{(n-1)w_i}{n(w_i-1)}\\right)+\\ln\\left(\\frac{w_i}{n}\\right)$だけ増えます。\n",
    "$\\lambda^*\\geq -u_i$が満たされないときは最適な$\\lambda$は$-u_i$になります。\n",
    "\n",
    "また、更新後の$w = (x_b^T H(u) x_b)^m_{b=1}$は\n",
    "\n",
    "$$w_b(u_+)=(1+\\lambda)\\left(w_b - \\frac{\\lambda}{1+\\lambda w_i}(\\hat{x}^T x_b)^2\\right)$$\n",
    "が成立します。ここで、$\\hat{x}:=Hx_i$としました。また、とくに\n",
    "\n",
    "$$w_i(u_+)=(1+\\lambda)\\left(w_i - \\frac{\\lambda w_i^2}{1+\\lambda w_i}\\right)=\\frac{(1+\\lambda)w_i}{1+\\lambda w_i}$$\n",
    "\n",
    "が成立します。"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### インデックス$i$の選び方\n",
    "\n",
    "教科書では選び方のモチベーションが３つ説明されていますが、ここでは簡単なやつを紹介します（詳細はp.28参照）。\n",
    "\n",
    "目的関数$g(u)$のテイラー展開を考えると、\n",
    "\n",
    "$$\n",
    "g(u+\\Delta u)\\approx \\tilde{g}(u+\\Delta u):= g(u)+w^T \\Delta u\n",
    "$$\n",
    "\n",
    "であるため、$e_i$の方向としては$w_i$が一番大きいインデックスを選べば良いことがわかります。これをするのが**Frank-Wolfeアルゴリズム**です。\n",
    "\n",
    "また、Frank-Wolfeに追加して、最も悪い$w_j$から逃げるように更新することも考えられます。これはAtwoodによって提案されているっぽいので、こちらはWolfe-Atwoodアルゴリズムと呼ばれています。"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Frank-Wolfeアルゴリズム\n",
    "\n",
    "以上をまとめて、Frank-Wolfeアルゴリズムを見てみましょう。\n",
    "Frank-Wolfeアルゴリズムの一番大変な部分は$w=(x_i^T H x_i)_i^m$であり、この計算の効率化がアルゴリズムの高速化に繋がります。\n",
    "上で説明したように、更新後の$w$は\n",
    "$$w_b(u_+)=(1+\\lambda)\\left(w_b - \\frac{\\lambda}{1+\\lambda w_i}(\\hat{x}_i^T x_b)^2\\right)$$\n",
    "であり、$\\hat{x}_i:=Hx_i$が効率的に計算できればよいわけですね。\n",
    "\n",
    "ここで、$XUX^T$のscaled Cholesky分解を\n",
    "$$XUX^T=LL^T$$\n",
    "とします。ここで$L$は下三角行列です。\n",
    "\n",
    "このとき、$\\hat{x}_i=L^{-T}L^{-1}x_i$の形になっています。\n",
    "実は、単純に$L^{-T}L^{-1}$を逆行列で計算すると計算量は$O(n^3)$ですが、この式を変形して連立一次方程式を解くと計算量は$O(n^2)$でいけます。\n",
    "まず、$z_i=L^{-1}x_i$を考えましょう。変形すると$Lz_i = x_i$であり、これを解いて$z_i$を求めます。\n",
    "次に$L^T \\hat{x}_i = z_i$を解けば$\\hat{x}_i$が求まります。\n",
    "（参考：[安易に逆行列を数値計算するのはやめよう](https://opqrstuvcut.github.io/blog/posts/%E5%AE%89%E6%98%93%E3%81%AB%E9%80%86%E8%A1%8C%E5%88%97%E3%82%92%E6%95%B0%E5%80%A4%E8%A8%88%E7%AE%97%E3%81%99%E3%82%8B%E3%81%AE%E3%81%AF%E3%82%84%E3%82%81%E3%82%88%E3%81%86/)）\n",
    "\n",
    "つまりCholesky分解があれば$O(n^2)$で$w$が求まるので、$L$も$L_+$に更新しましょう。\n",
    "ここで、\n",
    "$$XU_+X^T=\\frac{1}{(1+\\lambda)} LL^T + \\frac{\\lambda}{1 + \\lambda} x_i x_i^T$$\n",
    "であることに注意しましょう。\n",
    "つまり、$\\bar{L}:= L / \\sqrt{1 + \\lambda}$と$\\bar{x}:=\\sqrt{\\frac{\\lambda}{1 + \\lambda}}x_i$として、\n",
    "$$XU_+X^T=\\bar{L}\\bar{L}^T + \\bar{x}_i \\bar{x}_i^T$$\n",
    "になります。\n",
    "このようにCholesky分解をrank-oneで更新する場合、Cholesky分解は$O(n^2)$でできることが知られています（参考：[TFのrank-one Cholesky更新](https://www.tensorflow.org/probability/api_docs/python/tfp/math/cholesky_update)）。\n",
    "\n",
    "[Wikipedia (Cholesky decomposition)](https://en.wikipedia.org/wiki/Cholesky_decomposition)を参考に、rank-one Cholesky updateを実装してみましょう。\n",
    "また、$\\lambda$の符号が負の場合にも対応できるようにしておきます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import jax.numpy as jnp\n",
    "import numpy as np\n",
    "import jax\n",
    "\n",
    "\n",
    "@jax.jit\n",
    "def cholupdate(L, x, lam):\n",
    "    sign = jnp.sign(lam)\n",
    "    L = L / jnp.sqrt(1 + lam)\n",
    "    x = x * jnp.sqrt(sign * lam / (1 + lam))\n",
    "\n",
    "    n = jnp.size(x)\n",
    "    for k in range(n):\n",
    "        r = jnp.sqrt(L[k,k]**2 + sign * x[k]**2)\n",
    "        c = r / L[k,k]\n",
    "        s = x[k] / L[k,k]\n",
    "        L = L.at[k,k].set(r)\n",
    "        if k < n:\n",
    "            L = L.at[k+1:n, k].set((L[k+1:n, k] + sign * s*x[k+1:n]) / c)\n",
    "            x = x.at[k+1:n].set(c*x[k+1:n] - s*L[k+1:n, k])\n",
    "    return L\n",
    "\n",
    "\n",
    "# テスト\n",
    "n = 3\n",
    "m = 200\n",
    "X = jnp.array(np.random.randn(n, m))\n",
    "XUX = X @ jnp.diag(jnp.ones(m) / m) @ X.T\n",
    "L = jnp.linalg.cholesky(XUX)\n",
    "\n",
    "i = 3\n",
    "x = X[:, i]\n",
    "lam = 0.1\n",
    "L_updated = cholupdate(L, x, lam)\n",
    "XUX_updated = L @ L.T / (1 + lam) + x[:, None] @ x[None, :] * lam / (1 + lam)\n",
    "L_true = jnp.linalg.cholesky(XUX_updated)\n",
    "np.testing.assert_allclose(L_updated, L_true, atol=1e-5)\n",
    "\n",
    "lam_neg = -0.1\n",
    "L_updated = cholupdate(L, x, lam_neg)\n",
    "XUX_updated = L @ L.T / (1 + lam_neg) + x[:, None] @ x[None, :] * lam_neg / (1 + lam_neg)\n",
    "L_true = jnp.linalg.cholesky(XUX_updated)\n",
    "np.testing.assert_allclose(L_updated, L_true, atol=1e-5)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以上の内容をまとめてみましょう：\n",
    "\n",
    "**Frank-Wolfe アルゴリズム**\n",
    "\n",
    "0. 初期化です。\n",
    "   1. 双対問題の実行可能解$u$を一つ選択します。$u=e/m$とします。（他にも選び方は色々あります。教科書p.30参照。）\n",
    "   2. $\\varepsilon > 0$を選択します。\n",
    "   3. 初期のCholesky 分解$XUX^T=LL^T$を計算します。\n",
    "   4. 初期の$w=(x_b^T H x_b)_b^m$を計算します。\n",
    "1. $\\varepsilon_+ := \\max_b (w_b - n) / n$を計算し、最大値を与える$b$を$i$とします。計算量は$O(m)$です。\n",
    "2. もし$\\varepsilon_+ \\leq \\varepsilon$であれば、$u$は$\\varepsilon$最適な解なので終了します。（[CVX_minimum_volume_ellipsoids.ipynb](CVX_minimum_volume_ellipsoids.ipynb)を参考にしてください。）\n",
    "3. $\\varepsilon_+ > \\varepsilon$ならば$\\lambda^*:= (w_i-n) / ((n-1)w_i)$を計算し、$u \\leftarrow (1+\\lambda^*)^{-1}(u + \\lambda^* e_i)$で更新します。\n",
    "6. $\\bar{L}:= L / \\sqrt{1 + \\lambda}$と$\\bar{x}:=\\sqrt{\\frac{\\lambda}{1 + \\lambda}}x_i$として、$L\\leftarrow\\operatorname{cholupdate}(\\bar{L}, \\bar{x})$で$L$を更新します。\n",
    "7. 連立一次方程式$L z_i = x_i$を解いて$z_i$を計算します。\n",
    "8. 連立一次方程式$L^T \\hat{x}_i = z_i$を解いて$\\hat{x}_i$を計算します。\n",
    "9. $w_b=(1+\\lambda)\\left(w_b - \\frac{\\lambda}{1+\\lambda w_i}(\\hat{x}_i^T x_b)^2\\right)$によって$w$を更新します。\n",
    "10. ステップ１に戻ります。"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "コレスキー分解を使ったやつはあまり安定しなかったので、先に逆行列を丁寧に計算する方を実装してみます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGdCAYAAAA8F1jjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAYW0lEQVR4nO3db2xV9f3A8U8LchFtLyKlpbPyRxSHU/wXGowmIxJKw4huxqhRV4xxhrgZrX8mi4I4sYgJM/5C5Mkc+kTdE3m0mSgRzULFqCHLjC7AcMVAi2LoBRbLKOf34Bdv1p+AOLk935bXKzkJ95xz7/lwuOG+c++5bVWWZVkAACSoOu8BAACORagAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQrJF5D/B9HTlyJHbt2hU1NTVRVVWV9zgAwAnIsiz2798fjY2NUV197PdNhnyo7Nq1K5qamvIeAwD4L+zcuTPOOeecY24f8qFSU1MTEf/3F62trc15GgDgRJRKpWhqaiq/jh/LkA+Vrz/uqa2tFSoAMMR822UbLqYFAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGRVNFTeeeedWLhwYTQ2NkZVVVWsX79+wPZFixZFVVXVgGX+/PmVHAkAGEIqGioHDx6MmTNnxpo1a465z/z582P37t3l5eWXX67kSADAEDKykg/e2toara2tx92nUChEQ0NDJccAAIao3K9R2bhxY0yYMCGmT58eixcvjr179x53/76+viiVSgMWAGB4yjVU5s+fHy+99FJs2LAhnn766Xj77bejtbU1+vv7j3mfjo6OKBaL5aWpqWkQJwYABlNVlmXZoByoqipee+21uP7664+5zz/+8Y8477zz4s0334xrr732qPv09fVFX19f+XapVIqmpqbo7e2N2trakz02AFABpVIpisXit75+5/7Rz3+aOnVqjB8/PrZt23bMfQqFQtTW1g5YAIDhKalQ+eyzz2Lv3r0xceLEvEcBABJQ0W/9HDhwYMC7Izt27IgtW7bEuHHjYty4cbF8+fK44YYboqGhIbZv3x4PP/xwTJs2LVpaWio5FgAwRFQ0VN5///2YM2dO+XZ7e3tERLS1tcXzzz8ff/3rX+PFF1+Mffv2RWNjY8ybNy9++9vfRqFQqORYAMAQMWgX01bKiV6MAwCkY0heTAsA8J+ECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyapoqLzzzjuxcOHCaGxsjKqqqli/fv2A7VmWxdKlS2PixIlx+umnx9y5c2Pr1q2VHAkAGEIqGioHDx6MmTNnxpo1a466fdWqVfHcc8/F2rVrY/PmzXHGGWdES0tLfPXVV5UcCwAYIkZW8sFbW1ujtbX1qNuyLItnn302Hn300bjuuusiIuKll16K+vr6WL9+fdx8882VHA0AGAJyu0Zlx44d0d3dHXPnzi2vKxaL0dzcHJ2dnce8X19fX5RKpQELADA85RYq3d3dERFRX18/YH19fX1529F0dHREsVgsL01NTRWdEwDIz5D71s+SJUuit7e3vOzcuTPvkQCACsktVBoaGiIioqenZ8D6np6e8rajKRQKUVtbO2ABAIan3EJlypQp0dDQEBs2bCivK5VKsXnz5pg9e3ZeYwEACanot34OHDgQ27ZtK9/esWNHbNmyJcaNGxfnnntu3HffffHkk0/G+eefH1OmTInHHnssGhsb4/rrr6/kWADAEFHRUHn//fdjzpw55dvt7e0REdHW1hbr1q2Lhx9+OA4ePBi/+MUvYt++fXH11VfH66+/HqNHj67kWADAEFGVZVmW9xDfR6lUimKxGL29va5XAYAh4kRfv4fct34AgFOHUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBk5R4qjz/+eFRVVQ1YLrzwwrzHAgASMDLvASIiLrroonjzzTfLt0eOTGIsACBnSRTByJEjo6GhIe8xAIDE5P7RT0TE1q1bo7GxMaZOnRq33nprdHV1HXPfvr6+KJVKAxYAYHjKPVSam5tj3bp18frrr8fzzz8fO3bsiGuuuSb2799/1P07OjqiWCyWl6ampkGeGAAYLFVZlmV5D/Gf9u3bF5MmTYrVq1fHnXfe+Y3tfX190dfXV75dKpWiqakpent7o7a2djBHBQD+S6VSKYrF4re+fidxjcp/Gjt2bFxwwQWxbdu2o24vFApRKBQGeSoAIA+5f/Tz/x04cCC2b98eEydOzHsUACBnuYfKgw8+GG+//XZ8+umnsWnTpvjpT38aI0aMiFtuuSXv0QCAnOX+0c9nn30Wt9xyS+zduzfq6uri6quvjnfffTfq6uryHg0AyFnuofLKK6/kPQIAkKjcP/oBADgWoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKSCJU1a9bE5MmTY/To0dHc3Bzvvfde3iMBAAnIPVReffXVaG9vj2XLlsWHH34YM2fOjJaWltizZ0/eowEAOcs9VFavXh133XVX3HHHHTFjxoxYu3ZtjBkzJl544YW8RwMAcpZrqBw6dCg++OCDmDt3bnlddXV1zJ07Nzo7O496n76+viiVSgMWAGB4yjVUvvjii+jv74/6+voB6+vr66O7u/uo9+no6IhisVhempqaBmNUACAHuX/0810tWbIkent7y8vOnTvzHgkAqJCReR58/PjxMWLEiOjp6RmwvqenJxoaGo56n0KhEIVCYTDGAwBylus7KqNGjYorrrgiNmzYUF535MiR2LBhQ8yePTvHyQCAFOT6jkpERHt7e7S1tcWVV14Zs2bNimeffTYOHjwYd9xxR96jAQA5yz1Ubrrppvj8889j6dKl0d3dHZdeemm8/vrr37jAFgA49VRlWZblPcT3USqVolgsRm9vb9TW1uY9DgBwAk709XvIfesHADh1CBUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJKVa6hMnjw5qqqqBiwrV67McyQAICEj8x7giSeeiLvuuqt8u6amJsdpAICU5B4qNTU10dDQkPcYAECCcr9GZeXKlXH22WfHZZddFs8880wcPnz4uPv39fVFqVQasAAAw1Ou76jce++9cfnll8e4ceNi06ZNsWTJkti9e3esXr36mPfp6OiI5cuXD+KUAEBeqrIsy07mAz7yyCPx9NNPH3efjz/+OC688MJvrH/hhRfi7rvvjgMHDkShUDjqffv6+qKvr698u1QqRVNTU/T29kZtbe33Gx4AGBSlUimKxeK3vn6f9FD5/PPPY+/evcfdZ+rUqTFq1KhvrP/oo4/iRz/6UXzyyScxffr0Ezreif5FAYB0nOjr90n/6Keuri7q6ur+q/tu2bIlqqurY8KECSd5KgBgKMrtGpXOzs7YvHlzzJkzJ2pqaqKzszPuv//+uO222+Kss87KaywAICG5hUqhUIhXXnklHn/88ejr64spU6bE/fffH+3t7XmNBAAkJrdQufzyy+Pdd9/N6/AAwBCQ+89RAQA4FqECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyKhYqK1asiKuuuirGjBkTY8eOPeo+XV1dsWDBghgzZkxMmDAhHnrooTh8+HClRgIAhpiRlXrgQ4cOxY033hizZ8+O3//+99/Y3t/fHwsWLIiGhobYtGlT7N69O37+85/HaaedFk899VSlxgIAhpCqLMuySh5g3bp1cd9998W+ffsGrP/zn/8cP/nJT2LXrl1RX18fERFr166NX//61/H555/HqFGjTujxS6VSFIvF6O3tjdra2pM9PgBQASf6+p3bNSqdnZ1x8cUXlyMlIqKlpSVKpVJ89NFHx7xfX19flEqlAQsAMDzlFird3d0DIiUiyre7u7uPeb+Ojo4oFovlpampqaJzAgD5+U6h8sgjj0RVVdVxl08++aRSs0ZExJIlS6K3t7e87Ny5s6LHAwDy850upn3ggQdi0aJFx91n6tSpJ/RYDQ0N8d577w1Y19PTU952LIVCIQqFwgkdAwAY2r5TqNTV1UVdXd1JOfDs2bNjxYoVsWfPnpgwYUJERLzxxhtRW1sbM2bMOCnHAACGtop9Pbmrqyu+/PLL6Orqiv7+/tiyZUtEREybNi3OPPPMmDdvXsyYMSNuv/32WLVqVXR3d8ejjz4a99xzj3dMAICIqODXkxctWhQvvvjiN9a/9dZb8eMf/zgiIv75z3/G4sWLY+PGjXHGGWdEW1tbrFy5MkaOPPF+8vVkABh6TvT1u+I/R6XShAoADD3J/xwVAIBvI1QAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZFXslxIOlq9/A0CpVMp5EgDgRH39uv1tv8lnyIfK/v37IyKiqakp50kAgO9q//79USwWj7l9yP9SwiNHjsSuXbuipqYmqqqqTug+pVIpmpqaYufOnX6R4SBxzgefcz64nO/B55wPvpN5zrMsi/3790djY2NUVx/7SpQh/45KdXV1nHPOOf/VfWtraz25B5lzPvic88HlfA8+53zwnaxzfrx3Ur7mYloAIFlCBQBI1ikZKoVCIZYtWxaFQiHvUU4Zzvngc84Hl/M9+JzzwZfHOR/yF9MCAMPXKfmOCgAwNAgVACBZQgUASJZQAQCSdcqFyooVK+Kqq66KMWPGxNixY4+6T1dXVyxYsCDGjBkTEyZMiIceeigOHz48uIMOY5MnT46qqqoBy8qVK/Mea1hZs2ZNTJ48OUaPHh3Nzc3x3nvv5T3SsPX4449/4/l84YUX5j3WsPLOO+/EwoULo7GxMaqqqmL9+vUDtmdZFkuXLo2JEyfG6aefHnPnzo2tW7fmM+ww8W3nfNGiRd943s+fP78is5xyoXLo0KG48cYbY/HixUfd3t/fHwsWLIhDhw7Fpk2b4sUXX4x169bF0qVLB3nS4e2JJ56I3bt3l5df/epXeY80bLz66qvR3t4ey5Ytiw8//DBmzpwZLS0tsWfPnrxHG7YuuuiiAc/nv/zlL3mPNKwcPHgwZs6cGWvWrDnq9lWrVsVzzz0Xa9eujc2bN8cZZ5wRLS0t8dVXXw3ypMPHt53ziIj58+cPeN6//PLLlRkmO0X94Q9/yIrF4jfW/+lPf8qqq6uz7u7u8rrnn38+q62tzfr6+gZxwuFr0qRJ2e9+97u8xxi2Zs2ald1zzz3l2/39/VljY2PW0dGR41TD17Jly7KZM2fmPcYpIyKy1157rXz7yJEjWUNDQ/bMM8+U1+3bty8rFArZyy+/nMOEw8//P+dZlmVtbW3ZddddNyjHP+XeUfk2nZ2dcfHFF0d9fX15XUtLS5RKpfjoo49ynGx4WblyZZx99tlx2WWXxTPPPOOjtZPk0KFD8cEHH8TcuXPL66qrq2Pu3LnR2dmZ42TD29atW6OxsTGmTp0at956a3R1deU90iljx44d0d3dPeA5XywWo7m52XO+wjZu3BgTJkyI6dOnx+LFi2Pv3r0VOc6Q/6WEJ1t3d/eASImI8u3u7u48Rhp27r333rj88stj3LhxsWnTpliyZEns3r07Vq9enfdoQ94XX3wR/f39R30Of/LJJzlNNbw1NzfHunXrYvr06bF79+5Yvnx5XHPNNfG3v/0tampq8h5v2Pv6/+WjPef9n1058+fPj5/97GcxZcqU2L59e/zmN7+J1tbW6OzsjBEjRpzUYw2LUHnkkUfi6aefPu4+H3/8sQvcKui7/Bu0t7eX111yySUxatSouPvuu6Ojo8OPwmbIaW1tLf/5kksuiebm5pg0aVL88Y9/jDvvvDPHyaBybr755vKfL7744rjkkkvivPPOi40bN8a11157Uo81LELlgQceiEWLFh13n6lTp57QYzU0NHzjGxI9PT3lbRzd9/k3aG5ujsOHD8enn34a06dPr8B0p47x48fHiBEjys/Zr/X09Hj+DpKxY8fGBRdcENu2bct7lFPC18/rnp6emDhxYnl9T09PXHrppTlNdeqZOnVqjB8/PrZt2yZUjqauri7q6upOymPNnj07VqxYEXv27IkJEyZERMQbb7wRtbW1MWPGjJNyjOHo+/wbbNmyJaqrq8vnm//eqFGj4oorrogNGzbE9ddfHxERR44ciQ0bNsQvf/nLfIc7RRw4cCC2b98et99+e96jnBKmTJkSDQ0NsWHDhnKYlEql2Lx58zG/3cnJ99lnn8XevXsHxOLJMixC5bvo6uqKL7/8Mrq6uqK/vz+2bNkSERHTpk2LM888M+bNmxczZsyI22+/PVatWhXd3d3x6KOPxj333ONjiZOgs7MzNm/eHHPmzImampro7OyM+++/P2677bY466yz8h5vWGhvb4+2tra48sorY9asWfHss8/GwYMH44477sh7tGHpwQcfjIULF8akSZNi165dsWzZshgxYkTccssteY82bBw4cGDAO1Q7duyILVu2xLhx4+Lcc8+N++67L5588sk4//zzY8qUKfHYY49FY2NjOdb57o53zseNGxfLly+PG264IRoaGmL79u3x8MMPx7Rp06KlpeXkDzMo3y1KSFtbWxYR31jeeuut8j6ffvpp1tramp1++unZ+PHjswceeCD797//nd/Qw8gHH3yQNTc3Z8ViMRs9enT2wx/+MHvqqaeyr776Ku/RhpX/+Z//yc4999xs1KhR2axZs7J3330375GGrZtuuimbOHFiNmrUqOwHP/hBdtNNN2Xbtm3Le6xh5a233jrq/9ttbW1Zlv3fV5Qfe+yxrL6+PisUCtm1116b/f3vf8936CHueOf8X//6VzZv3rysrq4uO+2007JJkyZld91114Af63EyVWVZlp38/AEA+P78HBUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBk/S/i7BmbbD95swAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import jax.numpy as jnp\n",
    "import jax\n",
    "\n",
    "m = 100\n",
    "prngkey = jax.random.PRNGKey(5)\n",
    "Y = jax.random.normal(prngkey, (2, m)) * 3\n",
    "Y = Y + jax.random.uniform(prngkey, (2, m)) * 5\n",
    "Y = jnp.array(Y)  # 元データ\n",
    "X = jnp.vstack((Y, jnp.array([1,] * m)))  # 中心を計算する用\n",
    "n = 2 + 1\n",
    "assert X.shape == (n, m)\n",
    "\n",
    "# Frank-Wolfe アルゴリズム\n",
    "vmap_dot = jax.vmap(jnp.dot)\n",
    "\n",
    "# >>>> 初期化 >>>>\n",
    "# 0.1\n",
    "u = jnp.ones(m) / m\n",
    "\n",
    "# 0.2\n",
    "epsilon = 0.01  # 精度\n",
    "\n",
    "# 0.4\n",
    "XUX = X @ jnp.diag(u) @ X.T\n",
    "H = jnp.linalg.inv(XUX)\n",
    "vmap_xHx = jax.vmap(lambda x: x @ H @ x)\n",
    "w = vmap_xHx(X.T)\n",
    "assert w.shape == (m,)\n",
    "\n",
    "# <<<< 初期化 <<<<\n",
    "\n",
    "@jax.jit\n",
    "def one_step_exact_FW(u, L, w):\n",
    "    # １ステップFrank-Wolfeアルゴリズム\n",
    "\n",
    "    # 1. iの計算\n",
    "    i = w.argmax()\n",
    "    cur_epsilon = ((w - n) / n)[i]\n",
    "    w_i = w[i]\n",
    "    u_i = u[i]\n",
    "\n",
    "    # 3. 更新幅の計算とuの更新\n",
    "    lam = (w_i - n) / ((n - 1) * w_i)\n",
    "    lam = jnp.maximum(lam, -u_i)\n",
    "    u = u.at[i].set(u_i + lam)\n",
    "    u = u / (1 + lam)\n",
    "\n",
    "    # 0.3：コレスキー分解\n",
    "    XUX = X @ jnp.diag(u) @ X.T\n",
    "    H = jnp.linalg.inv(XUX)\n",
    "    vmap_xHx = jax.vmap(lambda x: x @ H @ x)\n",
    "    w_updated = vmap_xHx(X.T)\n",
    "    assert w.shape == (m,)\n",
    "    return cur_epsilon, u, L, w_updated\n",
    "\n",
    "\n",
    "# 描画\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from celluloid import Camera\n",
    "\n",
    "def plot_ellipsoid(u):\n",
    "    # データの描画\n",
    "    x, y = Y\n",
    "    plt.scatter(x, y)\n",
    "\n",
    "    # 楕円の描画\n",
    "    y_ = Y @ u\n",
    "    theta = np.linspace(0, 2 * np.pi, 100)\n",
    "    x = np.cos(theta)\n",
    "    y = np.sin(theta)\n",
    "    XUX = X @ jnp.diag(u) @ X.T\n",
    "    H = jnp.linalg.inv(XUX)\n",
    "    H = H[:-1, :-1]\n",
    "    L_plot = np.linalg.cholesky(H)\n",
    "    L_plot_T = np.linalg.inv(L_plot).T\n",
    "    x2, y2 = [], []\n",
    "    for _x, _y in zip(x, y):\n",
    "        v = y_ + np.sqrt(n) * np.dot(L_plot_T, np.array([_x, _y]))\n",
    "        x2.append(v[0])\n",
    "        y2.append(v[1])\n",
    "    plt.scatter(x2, y2, s=30)\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "camera = Camera(fig)\n",
    "\n",
    "cur_epsilon = jnp.inf\n",
    "count = 0\n",
    "while True:\n",
    "    cur_epsilon, new_u, new_L, new_w = one_step_exact_FW(u, L, w)\n",
    "    if cur_epsilon <= epsilon:\n",
    "        break\n",
    "    u, L, w = new_u, new_L, new_w\n",
    "    count += 1\n",
    "    if count < 10:\n",
    "        plot_ellipsoid(u)\n",
    "        camera.snap()\n",
    "\n",
    "\n",
    "animation = camera.animate()\n",
    "animation.save(\"./figs/MVEE_exact.gif\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![fig](figs/MVEE_exact.gif)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "コレスキー分解を使ったやつもやってみます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import jax.numpy as jnp\n",
    "import jax\n",
    "\n",
    "m = 10000\n",
    "prngkey = jax.random.PRNGKey(10)\n",
    "Y = jax.random.normal(prngkey, (2, m)) * 3\n",
    "Y = Y + jax.random.uniform(prngkey, (2, m)) * 5\n",
    "Y = jnp.array(Y)  # 元データ\n",
    "X = jnp.vstack((Y, jnp.array([1,] * m)))  # 中心を計算する用\n",
    "n = 2 + 1\n",
    "assert X.shape == (n, m)\n",
    "\n",
    "# Frank-Wolfe アルゴリズム\n",
    "vmap_dot = jax.vmap(jnp.dot)\n",
    "\n",
    "# >>>> 初期化 >>>>\n",
    "# 0.1\n",
    "u = jnp.ones(m) / m\n",
    "\n",
    "# 0.2\n",
    "epsilon = 0.01  # 精度\n",
    "\n",
    "# 0.3：コレスキー分解\n",
    "XUX = X @ jnp.diag(u) @ X.T\n",
    "L = jnp.linalg.cholesky(XUX)\n",
    "np.testing.assert_allclose(XUX, L @ L.T, atol=1e-5)\n",
    "\n",
    "# 0.4\n",
    "H = jnp.linalg.inv(XUX)\n",
    "vmap_xHx = jax.vmap(lambda x: x @ H @ x)\n",
    "w = vmap_xHx(X.T)\n",
    "assert w.shape == (m,)\n",
    "\n",
    "# <<<< 初期化 <<<<\n",
    "\n",
    "@jax.jit\n",
    "def one_step_FW(u, L, w):\n",
    "    # １ステップFrank-Wolfeアルゴリズム\n",
    "    # 1. iの計算\n",
    "    i = w.argmax()\n",
    "    cur_epsilon = ((w - n) / n)[i]\n",
    "    w_i = w[i]\n",
    "    x_i = X[:, i]\n",
    "    u_i = u[i]\n",
    "\n",
    "    # 3. 更新幅の計算とuの更新\n",
    "    lam = (w_i - n) / ((n - 1) * w_i)\n",
    "    lam = jnp.maximum(lam, -u_i)\n",
    "    u = u.at[i].set(u_i + lam)\n",
    "    u = u / (1 + lam)\n",
    "\n",
    "    # 5. wの更新\n",
    "    z_i = jnp.linalg.solve(L, x_i)\n",
    "    H_x_i = jnp.linalg.solve(L.T, z_i)\n",
    "    xx2 = (H_x_i @ X) ** 2\n",
    "    w_updated = (1 + lam) * (w - (lam / (1 + lam * w_i) * xx2))\n",
    "\n",
    "    # 4. 下三角行列の更新\n",
    "    L = cholupdate(L, x_i, lam)\n",
    "\n",
    "    return cur_epsilon, u, L, w_updated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGdCAYAAAA8F1jjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAexElEQVR4nO3db2wUBf7H8c8W6EKhu6X0z9KjQCsK4p+KKKX+uYg2tA0SUc6oUa8YwhksECiIcBEQFIqQcMYLJ0+E4gPB8wEYvNOIPSi5UDBimotGCHB4xbRbsKa7tHdsoZ3fA+P+7NFCObqd78L7lUzizs7ufocV953ZmdXjOI4jAAAAgxLcHgAAAKA7hAoAADCLUAEAAGYRKgAAwCxCBQAAmEWoAAAAswgVAABgFqECAADM6u/2ANeqo6ND9fX1Sk5OlsfjcXscAADQA47j6Ny5c8rKylJCQvfHTeI+VOrr65Wdne32GAAA4H9w+vRpjRgxotv74z5UkpOTJf20oz6fz+VpAABAT4TDYWVnZ0c/x7sT96Hy89c9Pp+PUAEAIM5c6bQNTqYFAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWTENlYqKCt17771KTk5WRkaGZsyYoWPHjnXa5vz58yorK9OwYcM0ZMgQzZw5U42NjbEcCwAAxImYhkp1dbXKysp06NAh7d27VxcuXNDUqVPV2toa3WbRokXas2ePPvzwQ1VXV6u+vl5PPPFELMcCAABxwuM4jtNXL3b27FllZGSourpav/71rxUKhZSenq73339fv/nNbyRJR48e1a233qqamhpNnjz5is8ZDofl9/sVCoX4ZVoAAOJETz+/+/QclVAoJElKTU2VJB05ckQXLlxQYWFhdJtx48Zp5MiRqqmp6cvRAACAQX32//rp6OjQwoULdf/99+v222+XJAWDQSUmJiolJaXTtpmZmQoGg10+TyQSUSQSid4Oh8MxmxkAALirz46olJWV6euvv9bOnTuv6XkqKirk9/ujS3Z2di9NCAAArOmTUJk3b54+/vhj7du3TyNGjIiuDwQCamtrU3Nzc6ftGxsbFQgEunyu5cuXKxQKRZfTp0/HcnQAAOCimIaK4ziaN2+edu3apb/97W/KycnpdP/EiRM1YMAAVVVVRdcdO3ZMdXV1Kigo6PI5vV6vfD5fpwUAAFyfYnqOSllZmd5//3199NFHSk5Ojp534vf7NWjQIPn9fs2ePVvl5eVKTU2Vz+fT/PnzVVBQ0KMrfgAAwPUtppcnezyeLtdv27ZNs2bNkvTTD74tXrxYO3bsUCQSUVFRkf70pz91+9XPf+PyZAAA4k9PP7/79HdUYoFQAQAg/pj8HRUAAICrQagAAACzCBUAAGAWoQIAAMwiVAAAgFmECgAAMItQAQAAZhEqAADALEIFAACYRagAAACzCBUAAGAWoQIAAMwiVAAAgFmECgAAMItQAQAAZhEqAADALEIFAACYRagAAACzCBUAAGAWoQIAAMwiVAAAgFmECgAAMItQAQAAZhEqAADALEIFAACYRagAAACzCBUAAGAWoQIAAMwiVAAAgFmECgAAMItQAQAAZhEqAADALEIFAACYRagAAACzCBUAAGAWoQIAAMwiVAAAgFmECgAAMItQAQAAZsU0VA4cOKDp06crKytLHo9Hu3fv7nT/rFmz5PF4Oi3FxcWxHAkAAMSRmIZKa2ur8vLytHnz5m63KS4uVkNDQ3TZsWNHLEcCAABxpH8sn7ykpEQlJSWX3cbr9SoQCMRyDAAAEKdcP0dl//79ysjI0NixYzV37lw1NTVddvtIJKJwONxpAQAA1ydXQ6W4uFjvvfeeqqqq9Oabb6q6ulolJSVqb2/v9jEVFRXy+/3RJTs7uw8nBgAAfcnjOI7TJy/k8WjXrl2aMWNGt9v885//1E033aTPP/9cjzzySJfbRCIRRSKR6O1wOKzs7GyFQiH5fL7eHhsAAMRAOByW3++/4ue361/9/FJubq7S0tJ04sSJbrfxer3y+XydFgAAcH0yFSrff/+9mpqaNHz4cLdHAQAABsT0qp+WlpZOR0dOnTql2tpapaamKjU1VatXr9bMmTMVCAR08uRJLV26VGPGjFFRUVEsxwIAAHEipqHy5ZdfasqUKdHb5eXlkqTS0lK98847+sc//qHt27erublZWVlZmjp1ql5//XV5vd5YjgUAAOJEn51MGys9PRkHAADYEZcn0wIAAPwSoQIAAMwiVAAAgFmECgAAMItQAQAAZhEqAADALEIFAACYRagAAACzCBUAAGAWoQIAAMwiVAAAgFmECgAAMItQAQAAZhEqAADALEIFAACYRagAAACzCBUAAGAWoQIAAMwiVAAAgFmECgAAMItQAQAAZhEqAADALEIFAACYRagAAACzCBUAAGAWoQIAAMwiVAAAgFmECgAAMItQAQAAZhEqAADALEIFAACYRagAAACzCBUAAGAWoQIAAMwiVAAAgFmECgAAMItQAQAAZhEqAADALEIFAACYFdNQOXDggKZPn66srCx5PB7t3r270/2O42jlypUaPny4Bg0apMLCQh0/fjyWIwEAgDgS01BpbW1VXl6eNm/e3OX9GzZs0Ntvv60tW7bo8OHDGjx4sIqKinT+/PlYjgUAAOJE/1g+eUlJiUpKSrq8z3EcvfXWW3r11Vf12GOPSZLee+89ZWZmavfu3Xr66adjORoAAIgDrp2jcurUKQWDQRUWFkbX+f1+5efnq6amptvHRSIRhcPhTgsAALg+uRYqwWBQkpSZmdlpfWZmZvS+rlRUVMjv90eX7OzsmM4JAADcE3dX/SxfvlyhUCi6nD592u2RAABAjLgWKoFAQJLU2NjYaX1jY2P0vq54vV75fL5OCwAAuD65Fio5OTkKBAKqqqqKrguHwzp8+LAKCgrcGgsAABgS06t+WlpadOLEiejtU6dOqba2VqmpqRo5cqQWLlyoN954QzfffLNycnK0YsUKZWVlacaMGbEcCwAAxImYhsqXX36pKVOmRG+Xl5dLkkpLS1VZWamlS5eqtbVVv/vd79Tc3KwHHnhAn376qQYOHBjLsQAAQJzwOI7juD3EtQiHw/L7/QqFQpyvAgBAnOjp53fcXfUDAABuHIQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYJbrofLaa6/J4/F0WsaNG+f2WAAAwID+bg8gSbfddps+//zz6O3+/U2MBQAAXGaiCPr3769AIOD2GAAAwBjXv/qRpOPHjysrK0u5ubl69tlnVVdX1+22kUhE4XC40wIAAK5ProdKfn6+Kisr9emnn+qdd97RqVOn9OCDD+rcuXNdbl9RUSG/3x9dsrOz+3hiAADQVzyO4zhuD/FLzc3NGjVqlDZt2qTZs2dfcn8kElEkEoneDofDys7OVigUks/n68tRAQDA/ygcDsvv91/x89vEOSq/lJKSoltuuUUnTpzo8n6v1yuv19vHUwEAADe4/tXPf2tpadHJkyc1fPhwt0cBAAAucz1UlixZourqan333Xc6ePCgHn/8cfXr10/PPPOM26MBAACXuf7Vz/fff69nnnlGTU1NSk9P1wMPPKBDhw4pPT3d7dEAAIDLXA+VnTt3uj0CAAAwyvWvfgAAALpDqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwy0SobN68WaNHj9bAgQOVn5+vL774wu2RAACAAa6HygcffKDy8nKtWrVKX331lfLy8lRUVKQzZ864PRoAAHCZ66GyadMmzZkzRy+88ILGjx+vLVu2KCkpSVu3bnV7NAAA4DJXQ6WtrU1HjhxRYWFhdF1CQoIKCwtVU1PT5WMikYjC4XCnBQAAXJ9cDZUffvhB7e3tyszM7LQ+MzNTwWCwy8dUVFTI7/dHl+zs7L4YFQAAuMD1r36u1vLlyxUKhaLL6dOn3R4JAADESH83XzwtLU39+vVTY2Njp/WNjY0KBAJdPsbr9crr9fbFeAAAwGWuHlFJTEzUxIkTVVVVFV3X0dGhqqoqFRQUuDgZAACwwNUjKpJUXl6u0tJS3XPPPZo0aZLeeusttba26oUXXnB7NAAA4DLXQ+Wpp57S2bNntXLlSgWDQd1111369NNPLznBFgAA3Hg8juM4bg9xLcLhsPx+v0KhkHw+n9vjAACAHujp53fcXfUDAABuHIQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYJaroTJ69Gh5PJ5Oy/r1690cCQAAGNLf7QHWrFmjOXPmRG8nJye7OA0AALDE9VBJTk5WIBBwewwAAGCQ6+eorF+/XsOGDdOECRO0ceNGXbx48bLbRyIRhcPhTgsAALg+uXpEZcGCBbr77ruVmpqqgwcPavny5WpoaNCmTZu6fUxFRYVWr17dh1MCAAC3eBzHcXrzCZctW6Y333zzstt8++23Gjdu3CXrt27dqhdffFEtLS3yer1dPjYSiSgSiURvh8NhZWdnKxQKyefzXdvwAACgT4TDYfn9/it+fvd6qJw9e1ZNTU2X3SY3N1eJiYmXrP/mm290++236+jRoxo7dmyPXq+nOwoAAOzo6ed3r3/1k56ervT09P/psbW1tUpISFBGRkYvTwUAAOKRa+eo1NTU6PDhw5oyZYqSk5NVU1OjRYsW6bnnntPQoUPdGgsAABjiWqh4vV7t3LlTr732miKRiHJycrRo0SKVl5e7NRIAADDGtVC5++67dejQIbdeHgAAxAHXf0cFAACgO4QKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYFbMQmXt2rW67777lJSUpJSUlC63qaur07Rp05SUlKSMjAy9/PLLunjxYqxGAgAAcaZ/rJ64ra1NTz75pAoKCvTuu+9ecn97e7umTZumQCCggwcPqqGhQb/97W81YMAArVu3LlZjAQCAOOJxHMeJ5QtUVlZq4cKFam5u7rT+k08+0aOPPqr6+nplZmZKkrZs2aJXXnlFZ8+eVWJiYo+ePxwOy+/3KxQKyefz9fb4AAAgBnr6+e3aOSo1NTW64447opEiSUVFRQqHw/rmm2+6fVwkElE4HO60AACA65NroRIMBjtFiqTo7WAw2O3jKioq5Pf7o0t2dnZM5wQAAO65qlBZtmyZPB7PZZejR4/GalZJ0vLlyxUKhaLL6dOnY/p6AADAPVd1Mu3ixYs1a9asy26Tm5vbo+cKBAL64osvOq1rbGyM3tcdr9crr9fbo9cAAADx7apCJT09Xenp6b3ywgUFBVq7dq3OnDmjjIwMSdLevXvl8/k0fvz4XnkNAAAQ32J2eXJdXZ1+/PFH1dXVqb29XbW1tZKkMWPGaMiQIZo6darGjx+v559/Xhs2bFAwGNSrr76qsrIyjpgAAABJMbw8edasWdq+ffsl6/ft26eHHnpIkvSvf/1Lc+fO1f79+zV48GCVlpZq/fr16t+/5/3E5ckAAMSfnn5+x/x3VGKNUAEAIP6Y/x0VAACAKyFUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJhFqAAAALMIFQAAYBahAgAAzCJUAACAWYQKAAAwi1ABAABmESoAAMAsQgUAAJgVs1BZu3at7rvvPiUlJSklJaXLbTwezyXLzp07YzUSAACIM/1j9cRtbW168sknVVBQoHfffbfb7bZt26bi4uLo7e6iBgAA3HhiFiqrV6+WJFVWVl52u5SUFAUCgViNAQAA4pjr56iUlZUpLS1NkyZN0tatW+U4jtsjAQAAI2J2RKUn1qxZo4cfflhJSUn67LPP9NJLL6mlpUULFizo9jGRSESRSCR6OxwO98WoAADABVd1RGXZsmVdngD7y+Xo0aM9fr4VK1bo/vvv14QJE/TKK69o6dKl2rhx42UfU1FRIb/fH12ys7OvZhcAAEAc8ThX8V3L2bNn1dTUdNltcnNzlZiYGL1dWVmphQsXqrm5+YrP/5e//EWPPvqozp8/L6/X2+U2XR1Ryc7OVigUks/n69mOAAAAV4XDYfn9/it+fl/VVz/p6elKT0+/5uG6U1tbq6FDh3YbKZLk9Xovez8AALh+xOwclbq6Ov3444+qq6tTe3u7amtrJUljxozRkCFDtGfPHjU2Nmry5MkaOHCg9u7dq3Xr1mnJkiWxGgkAAMSZmIXKypUrtX379ujtCRMmSJL27dunhx56SAMGDNDmzZu1aNEiOY6jMWPGaNOmTZozZ06sRgIAAHHmqs5RsSgUCiklJUWnT5/mHBUAAOLEz+eYNjc3y+/3d7udq5cn94Zz585JElf/AAAQh86dO3fZUIn7IyodHR2qr69XcnKyPB6P2+NI+v9KvJGO8rDPN8Y+SzfmfrPP7PP1ys19dhxH586dU1ZWlhISuv+1lLg/opKQkKARI0a4PUaXfD7fDfMv+8/Y5xvHjbjf7PONgX3uO5c7kvIz139CHwAAoDuECgAAMItQiQGv16tVq1bdUD9Mxz7fOG7E/Wafbwzss01xfzItAAC4fnFEBQAAmEWoAAAAswgVAABgFqECAADMIlR62dq1a3XfffcpKSlJKSkpXW7j8XguWXbu3Nm3g/ainuxzXV2dpk2bpqSkJGVkZOjll1/WxYsX+3bQGBs9evQl7+v69evdHqtXbd68WaNHj9bAgQOVn5+vL774wu2RYua111675P0cN26c22P1qgMHDmj69OnKysqSx+PR7t27O93vOI5Wrlyp4cOHa9CgQSosLNTx48fdGbaXXGmfZ82adcn7Xlxc7M6wvaSiokL33nuvkpOTlZGRoRkzZujYsWOdtjl//rzKyso0bNgwDRkyRDNnzlRjY6NLE3dGqPSytrY2Pfnkk5o7d+5lt9u2bZsaGhqiy4wZM/pmwBi40j63t7dr2rRpamtr08GDB7V9+3ZVVlZq5cqVfTxp7K1Zs6bT+zp//ny3R+o1H3zwgcrLy7Vq1Sp99dVXysvLU1FRkc6cOeP2aDFz2223dXo///73v7s9Uq9qbW1VXl6eNm/e3OX9GzZs0Ntvv60tW7bo8OHDGjx4sIqKinT+/Pk+nrT3XGmfJam4uLjT+75jx44+nLD3VVdXq6ysTIcOHdLevXt14cIFTZ06Va2trdFtFi1apD179ujDDz9UdXW16uvr9cQTT7g49S84iIlt27Y5fr+/y/skObt27erTefpCd/v817/+1UlISHCCwWB03TvvvOP4fD4nEon04YSxNWrUKOcPf/iD22PEzKRJk5yysrLo7fb2dicrK8upqKhwcarYWbVqlZOXl+f2GH3mv/+71NHR4QQCAWfjxo3Rdc3NzY7X63V27NjhwoS9r6v/FpeWljqPPfaYK/P0lTNnzjiSnOrqasdxfnpfBwwY4Hz44YfRbb799ltHklNTU+PWmFEcUXFJWVmZ0tLSNGnSJG3dulXOdfxzNjU1NbrjjjuUmZkZXVdUVKRwOKxvvvnGxcl63/r16zVs2DBNmDBBGzduvG6+3mpra9ORI0dUWFgYXZeQkKDCwkLV1NS4OFlsHT9+XFlZWcrNzdWzzz6ruro6t0fqM6dOnVIwGOz0nvv9fuXn51/X77kk7d+/XxkZGRo7dqzmzp2rpqYmt0fqVaFQSJKUmpoqSTpy5IguXLjQ6b0eN26cRo4caeK9jvv/KWE8WrNmjR5++GElJSXps88+00svvaSWlhYtWLDA7dFiIhgMdooUSdHbwWDQjZFiYsGCBbr77ruVmpqqgwcPavny5WpoaNCmTZvcHu2a/fDDD2pvb+/yfTx69KhLU8VWfn6+KisrNXbsWDU0NGj16tV68MEH9fXXXys5Odnt8WLu57+bXb3n19Pf2/9WXFysJ554Qjk5OTp58qR+//vfq6SkRDU1NerXr5/b412zjo4OLVy4UPfff79uv/12ST+914mJiZecY2jlvSZUemDZsmV68803L7vNt99+2+MT7VasWBH95wkTJqi1tVUbN240FSq9vc/x6mr+HMrLy6Pr7rzzTiUmJurFF19URUWF6Z+nRtdKSkqi/3znnXcqPz9fo0aN0p///GfNnj3bxckQS08//XT0n++44w7deeeduummm7R//3498sgjLk7WO8rKyvT111/H1flWhEoPLF68WLNmzbrsNrm5uf/z8+fn5+v1119XJBIx84HWm/scCAQuuTrk57PJA4HA/zRfX7mWP4f8/HxdvHhR3333ncaOHRuD6fpOWlqa+vXrd8lVAI2Njebfw96SkpKiW265RSdOnHB7lD7x8/va2Nio4cOHR9c3Njbqrrvucmmqvpebm6u0tDSdOHEi7kNl3rx5+vjjj3XgwAGNGDEiuj4QCKitrU3Nzc2djqpY+ftNqPRAenq60tPTY/b8tbW1Gjp0qJlIkXp3nwsKCrR27VqdOXNGGRkZkqS9e/fK5/Np/PjxvfIasXItfw61tbVKSEiI7nM8S0xM1MSJE1VVVRW9Qq2jo0NVVVWaN2+eu8P1kZaWFp08eVLPP/+826P0iZycHAUCAVVVVUXDJBwO6/Dhw1e8qvF68v3336upqalTrMUbx3E0f/587dq1S/v371dOTk6n+ydOnKgBAwaoqqpKM2fOlCQdO3ZMdXV1KigocGPkTgiVXlZXV6cff/xRdXV1am9vV21trSRpzJgxGjJkiPbs2aPGxkZNnjxZAwcO1N69e7Vu3TotWbLE3cGvwZX2eerUqRo/fryef/55bdiwQcFgUK+++qrKyspMxdm1qKmp0eHDhzVlyhQlJyerpqZGixYt0nPPPaehQ4e6PV6vKC8vV2lpqe655x5NmjRJb731llpbW/XCCy+4PVpMLFmyRNOnT9eoUaNUX1+vVatWqV+/fnrmmWfcHq3XtLS0dDpCdOrUKdXW1io1NVUjR47UwoUL9cYbb+jmm29WTk6OVqxYoaysrLj+OYXL7XNqaqpWr16tmTNnKhAI6OTJk1q6dKnGjBmjoqIiF6e+NmVlZXr//ff10UcfKTk5OXreid/v16BBg+T3+zV79myVl5crNTVVPp9P8+fPV0FBgSZPnuzy9OLy5N5WWlrqSLpk2bdvn+M4jvPJJ584d911lzNkyBBn8ODBTl5enrNlyxanvb3d3cGvwZX22XEc57vvvnNKSkqcQYMGOWlpac7ixYudCxcuuDd0Lzty5IiTn5/v+P1+Z+DAgc6tt97qrFu3zjl//rzbo/WqP/7xj87IkSOdxMREZ9KkSc6hQ4fcHilmnnrqKWf48OFOYmKi86tf/cp56qmnnBMnTrg9Vq/at29fl393S0tLHcf56RLlFStWOJmZmY7X63UeeeQR59ixY+4OfY0ut8///ve/nalTpzrp6enOgAEDnFGjRjlz5szp9NMK8air/ZXkbNu2LbrNf/7zH+ell15yhg4d6iQlJTmPP/6409DQ4N7Qv+BxnOv4ulgAABDX+B0VAABgFqECAADMIlQAAIBZhAoAADCLUAEAAGYRKgAAwCxCBQAAmEWoAAAAswgVAABgFqECAADMIlQAAIBZhAoAADDr/wBrKxfLgsxL8QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import HTML\n",
    "from celluloid import Camera\n",
    "\n",
    "\n",
    "def plot_ellipsoid(u):\n",
    "    # データの描画\n",
    "    x, y = Y\n",
    "    plt.scatter(x, y)\n",
    "\n",
    "    # 楕円の描画\n",
    "    y_ = Y @ u\n",
    "    theta = np.linspace(0, 2 * np.pi, 100)\n",
    "    x = np.cos(theta)\n",
    "    y = np.sin(theta)\n",
    "    XUX = X @ jnp.diag(u) @ X.T\n",
    "    H = jnp.linalg.inv(XUX)\n",
    "    H = H[:-1, :-1]\n",
    "    L_plot = np.linalg.cholesky(H)\n",
    "    L_plot_T = np.linalg.inv(L_plot).T\n",
    "    x2, y2 = [], []\n",
    "    for _x, _y in zip(x, y):\n",
    "        v = y_ + np.sqrt(n) * np.dot(L_plot_T, np.array([_x, _y]))\n",
    "        x2.append(v[0])\n",
    "        y2.append(v[1])\n",
    "    plt.scatter(x2, y2, s=30)\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "camera = Camera(fig)\n",
    "\n",
    "cur_epsilon = jnp.inf\n",
    "count = 0\n",
    "while True:\n",
    "    cur_epsilon, new_u, new_L, new_w = one_step_FW(u, L, w)\n",
    "    if cur_epsilon <= epsilon:\n",
    "        break\n",
    "    u, L, w = new_u, new_L, new_w\n",
    "    count += 1\n",
    "    if count < 10:\n",
    "        plot_ellipsoid(u)\n",
    "        camera.snap()\n",
    "\n",
    "\n",
    "animation = camera.animate()\n",
    "animation.save(\"./figs/MVEE.gif\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![figs](figs/MVEE.gif)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Core-setの獲得\n",
    "\n",
    "主問題と双対問題を改めて確認しましょう。\n",
    "\n",
    "---\n",
    "\n",
    "**主問題と双対問題の直感**\n",
    "\n",
    "$$\n",
    "\\text{ 主問題 }\\;  \n",
    "\\begin{cases}\n",
    "&\\min_{H\\in S^n} f(H) := -\\ln \\det (H) \\quad &&\\text{ （楕円の面積を最小にする） }\\\\\n",
    "&x_i^T H x_i \\leq n \\quad &&\\text {（楕円はデータを全て含む）}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "主問題の解釈は簡単ですね。双対問題はどうでしょうか。\n",
    "\n",
    "$$\n",
    "\\text{ 双対問題 }\\; \n",
    "\\begin{cases}\n",
    "&\\max_{u\\in \\mathbb{R}^m} g(u) := \\ln \\det (XUX^T) \\\\\n",
    "& e^T u = 1\\\\\n",
    "& u \\geq 0\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "* TODO: こっからまじで雰囲気で書いています。線形代数の回を作ったら書き直します。一応[Derivative of ln (det X)](https://statisticaloddsandends.wordpress.com/2018/05/24/derivative-of-ln-det-x/)に基づいた直感を書いています。\n",
    "    * $XUX^T$は$\\sum_i^n u_i x_i x_i^T$であるから、$\\det (XUX^T)$は「$U$によって選択されたデータが作る平行六面体の体積」と解釈することができます\n",
    "    * また、$H(u)=(XUX^T)^{-1}$によって作られる楕円は、一番大きい平行六面体を囲むように作られた楕円と解釈できます。\n",
    "\n",
    "つまり、$u_i=0$のデータ$x_i$は最小楕円を作る上では重要ではないデータであり、$\\{x_i: u_i > 0\\}$であるデータだけでも同じ最小楕円が得られます。このような$\\{x_i: u_i > 0\\}$は**core-set**と呼ばれます。\n",
    "\n",
    "---\n",
    "\n",
    "これを踏まえて、$\\varepsilon$-core setを導入します。\n",
    "まず、[CVX_minimum_volume_ellipsoids.ipynb](CVX_minimum_volume_ellipsoids.ipynb)でやった近似最適解の性質を思い出しましょう。\n",
    "\n",
    "* 実行可能解$u$は次を満たすとき$\\varepsilon$-主実行可能であると言う。\n",
    "    * $H(u):=(XU^TX)^{-1}$が $x_i^T H(u)x_i \\leq (1 + \\varepsilon) n, \\quad i=1,\\dots, m$ を満たす。\n",
    "\n",
    "よって、$u$が$\\varepsilon$-主実行可能解であるならば、$H(u)$によって作られる楕円は最適解よりも大きくなり、かつ全てのデータを含みます。\n",
    "また、書き換えると$\\mathcal{E}([1+\\varepsilon]^{-1}H(u))$が$X$を含むこともわかります。\n",
    "$H(u)$を$[1+\\varepsilon]^{-1}$でスケールしているので、その体積は$[1+\\varepsilon]^{-n}$倍されます。\n",
    "\n",
    "このような$u$について$S=\\{x_i; u_i >0\\}$とすると、$g(u)$は$S$を含む最小楕円の近似解でもあるわけです。\n",
    "$\\varepsilon$-主実行可能解は$n\\ln(1+\\varepsilon)$-以内最適であることを思い出すと、$S$が生成する楕円の最小楕円との違いは$[1+\\varepsilon]^{-n}$のオーダー程度になります（多分）。\n",
    "このような$S$は$\\varepsilon$-core-setと呼ばれます。\n",
    "\n",
    "\n",
    "**データの削減**\n",
    "\n",
    "TODO: 詳細書く。\n",
    "\n",
    "実はアルゴリズムが$\\varepsilon$-主実行可能解（もしくは$\\varepsilon$-最適解）で停止すると、FWアルゴリズムは$4n(\\ln \\ln n+4) + 28 n / \\varepsilon$のサイズのコアセットを出力します。\n",
    "（これは特殊な初期条件についてしか成り立たないことに注意しましょう。詳しくは教科書p.43）\n",
    "\n",
    "詳細は省略しますが、$\\delta:=\\delta(u)=(w_i(u) - n) / n > 0$とすると、\n",
    "\n",
    "$$\n",
    "w_i(u) < n\\left(1 + \\frac{\\delta n}{2} - \\sqrt{\\delta n - \\delta + \\frac{\\delta^2 n^2}{4}}\\right)\n",
    "$$\n",
    "を満たす$x_i$は無視することでcore-setを作ることができます（教科書p.44）。\n",
    "\n",
    "このデータの削減をWolfe-Atwoodのアルゴリズム（教科書p.29）と合体させて確認してみましょう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGdCAYAAAA8F1jjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAbiElEQVR4nO3df2xV9f348Vf5dRFtL1YKpVtFKhPHpkzdrPjRZE7CjzgzN2PUOAfGuc0xDVTnZFEQFYuaMONC9B8n7o+p2x+6bNnMlAhmoWLm0iwzaoThioEWh6FX2LxoOd8/9t2NHS2i4/a8C49HchLOj3vvi5Mr95l7z73WZFmWBQBAgkbkPQAAwGCECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkalfcA/6v9+/fH9u3bo7a2NmpqavIeBwA4BFmWxbvvvhtNTU0xYsTg75sM+1DZvn17NDc35z0GAPAJbNu2LT796U8Pun/Yh0ptbW1E/PsvWldXl/M0AMChKJVK0dzcXHkdH8ywD5X/fNxTV1cnVABgmPmoyzZcTAsAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMmqaqi88MILcfHFF0dTU1PU1NTE008/3W//woULo6ampt8yb968ao4EAAwjVQ2VvXv3xsyZM2PNmjWDHjNv3rzYsWNHZXn88cerORIAMIyMquadz58/P+bPn3/QYwqFQjQ2NlZzDABgmMr9GpX169fHxIkTY/r06XH99dfHrl27Dnp8uVyOUqnUbwEAjky5hsq8efPi5z//eaxbty7uvffe2LBhQ8yfPz/6+voGvU17e3sUi8XK0tzcPIQTAwBDqSbLsmxIHqimJp566qm45JJLBj3mb3/7W5x88snx3HPPxYUXXjjgMeVyOcrlcmW9VCpFc3Nz9Pb2Rl1d3eEeGwCoglKpFMVi8SNfv3P/6OfDWlpaYsKECbF58+ZBjykUClFXV9dvAQCOTEmFyltvvRW7du2KyZMn5z0KAJCAqn7rZ8+ePf3eHdm6dWt0dnZGfX191NfXx4oVK+LSSy+NxsbG2LJlS9xyyy0xbdq0mDt3bjXHAgCGiaqGyp/+9Ke44IILKuttbW0REbFgwYJ46KGH4i9/+Us89thjsXv37mhqaoo5c+bEXXfdFYVCoZpjAQDDxJBdTFsth3oxDgCQjmF5MS0AwIcJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZVQ2VF154IS6++OJoamqKmpqaePrpp/vtz7Isli1bFpMnT45jjjkmZs+eHW+88UY1RwIAhpGqhsrevXtj5syZsWbNmgH333ffffHggw/Gww8/HJs2bYpjjz025s6dG++99141xwIAholR1bzz+fPnx/z58wfcl2VZPPDAA3HbbbfF1772tYiI+PnPfx6TJk2Kp59+Oq644opqjgYADAO5XaOydevW6O7ujtmzZ1e2FYvFaG1tjY6OjkFvVy6Xo1Qq9VsAgCNTbqHS3d0dERGTJk3qt33SpEmVfQNpb2+PYrFYWZqbm6s6JwCQn2H3rZ+lS5dGb29vZdm2bVveIwEAVZJbqDQ2NkZERE9PT7/tPT09lX0DKRQKUVdX128BAI5MuYXK1KlTo7GxMdatW1fZViqVYtOmTTFr1qy8xgIAElLVb/3s2bMnNm/eXFnfunVrdHZ2Rn19fZx44omxePHiuPvuu+Mzn/lMTJ06NW6//fZoamqKSy65pJpjAQDDRFVD5U9/+lNccMEFlfW2traIiFiwYEGsXbs2brnllti7d2985zvfid27d8d5550XzzzzTIwdO7aaYwEAw0RNlmVZ3kP8L0qlUhSLxejt7XW9CgAME4f6+j3svvUDABw9hAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLJyD5U77rgjampq+i2nnnpq3mMBAAkYlfcAERGf+9zn4rnnnqusjxqVxFgAQM6SKIJRo0ZFY2Nj3mMAAInJ/aOfiIg33ngjmpqaoqWlJa666qro6uoa9NhyuRylUqnfAgAcmXIPldbW1li7dm0888wz8dBDD8XWrVvj/PPPj3fffXfA49vb26NYLFaW5ubmIZ4YABgqNVmWZXkP8WG7d++OKVOmxOrVq+Paa689YH+5XI5yuVxZL5VK0dzcHL29vVFXVzeUowIAn1CpVIpisfiRr99JXKPyYePHj49TTjklNm/ePOD+QqEQhUJhiKcCAPKQ+0c//23Pnj2xZcuWmDx5ct6jAAA5yz1Ubr755tiwYUO8+eabsXHjxvj6178eI0eOjCuvvDLv0QCAnOX+0c9bb70VV155ZezatSsaGhrivPPOixdffDEaGhryHg0AyFnuofLEE0/kPQIAkKjcP/oBABiMUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSlUSorFmzJk466aQYO3ZstLa2xksvvZT3SABAAnIPlSeffDLa2tpi+fLl8ec//zlmzpwZc+fOjZ07d+Y9GgCQs9xDZfXq1XHdddfFNddcEzNmzIiHH344xo0bFz/72c/yHg0AyFmuobJv3754+eWXY/bs2ZVtI0aMiNmzZ0dHR8eAtymXy1EqlfotAMCRKddQ+cc//hF9fX0xadKkftsnTZoU3d3dA96mvb09isViZWlubh6KUQGAHOT+0c/HtXTp0ujt7a0s27Zty3skAKBKRuX54BMmTIiRI0dGT09Pv+09PT3R2Ng44G0KhUIUCoWhGA8AyFmu76iMGTMmzjrrrFi3bl1l2/79+2PdunUxa9asHCcDAFKQ6zsqERFtbW2xYMGC+OIXvxhnn312PPDAA7F379645ppr8h4NAMhZ7qFy+eWXx9tvvx3Lli2L7u7u+MIXvhDPPPPMARfYAgBHn5osy7K8h/hflEqlKBaL0dvbG3V1dXmPAwAcgkN9/R523/oBAI4eQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIVq6hctJJJ0VNTU2/ZdWqVXmOBAAkZFTeA9x5551x3XXXVdZra2tznAYASEnuoVJbWxuNjY15jwEAJCj3a1RWrVoVJ5xwQpxxxhlx//33xwcffHDQ48vlcpRKpX4LAHBkyvUdlRtvvDHOPPPMqK+vj40bN8bSpUtjx44dsXr16kFv097eHitWrBjCKQGAvNRkWZYdzju89dZb49577z3oMa+++mqceuqpB2z/2c9+Ft/97ndjz549USgUBrxtuVyOcrlcWS+VStHc3By9vb1RV1f3vw0PAAyJUqkUxWLxI1+/D3uovP3227Fr166DHtPS0hJjxow5YPsrr7wSn//85+O1116L6dOnH9LjHepfFABIx6G+fh/2j34aGhqioaHhE922s7MzRowYERMnTjzMUwEAw1Fu16h0dHTEpk2b4oILLoja2tro6OiIJUuWxDe/+c04/vjj8xoLAEhIbqFSKBTiiSeeiDvuuCPK5XJMnTo1lixZEm1tbXmNBAAkJrdQOfPMM+PFF1/M6+EBgGEg999RAQAYjFABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEhW1UJl5cqVce6558a4ceNi/PjxAx7T1dUVF110UYwbNy4mTpwYP/zhD+ODDz6o1kgAwDAzqlp3vG/fvrjsssti1qxZ8cgjjxywv6+vLy666KJobGyMjRs3xo4dO+Jb3/pWjB49Ou65555qjQUADCM1WZZl1XyAtWvXxuLFi2P37t39tv/+97+Pr371q7F9+/aYNGlSREQ8/PDD8aMf/SjefvvtGDNmzCHdf6lUimKxGL29vVFXV3e4xwcAquBQX79zu0alo6MjTjvttEqkRETMnTs3SqVSvPLKK4PerlwuR6lU6rcAAEem3EKlu7u7X6RERGW9u7t70Nu1t7dHsVisLM3NzVWdEwDIz8cKlVtvvTVqamoOurz22mvVmjUiIpYuXRq9vb2VZdu2bVV9PAAgPx/rYtqbbropFi5ceNBjWlpaDum+Ghsb46WXXuq3raenp7JvMIVCIQqFwiE9BgAwvH2sUGloaIiGhobD8sCzZs2KlStXxs6dO2PixIkREfHss89GXV1dzJgx47A8BgAwvFXt68ldXV3xzjvvRFdXV/T19UVnZ2dEREybNi2OO+64mDNnTsyYMSOuvvrquO+++6K7uztuu+22WLRokXdMAICIqOLXkxcuXBiPPfbYAduff/75+PKXvxwREX//+9/j+uuvj/Xr18exxx4bCxYsiFWrVsWoUYfeT76eDADDz6G+flf9d1SqTagAwPCT/O+oAAB8FKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJKtqobJy5co499xzY9y4cTF+/PgBj6mpqTlgeeKJJ6o1EgAwzIyq1h3v27cvLrvsspg1a1Y88sgjgx736KOPxrx58yrrg0UNAHD0qVqorFixIiIi1q5de9Djxo8fH42NjdUaAwAYxqoWKodq0aJF8e1vfztaWlrie9/7XlxzzTVRU1Mz6PHlcjnK5XJlvbe3NyIiSqVS1WcFAA6P/7xuZ1l20ONyDZU777wzvvKVr8S4cePiD3/4Q3z/+9+PPXv2xI033jjobdrb2yvv1nxYc3NzNUcFAKrg3XffjWKxOOj+muyjUuZDbr311rj33nsPesyrr74ap556amV97dq1sXjx4ti9e/dH3v+yZcvi0UcfjW3btg16zH+/o7J///5455134oQTTjjoOzHDSalUiubm5ti2bVvU1dXlPU5SnJuBOS+Dc24G5rwMzrkZ2OE+L1mWxbvvvhtNTU0xYsTg3+35WO+o3HTTTbFw4cKDHtPS0vJx7rKf1tbWuOuuu6JcLkehUBjwmEKhcMC+I/UC3Lq6Ov+RDMK5GZjzMjjnZmDOy+Ccm4EdzvNysHdS/uNjhUpDQ0M0NDR84oE+SmdnZxx//PGDRgoAcHSp2jUqXV1d8c4770RXV1f09fVFZ2dnRERMmzYtjjvuuPjNb34TPT09cc4558TYsWPj2WefjXvuuSduvvnmao0EAAwzVQuVZcuWxWOPPVZZP+OMMyIi4vnnn48vf/nLMXr06FizZk0sWbIksiyLadOmxerVq+O6666r1kjDRqFQiOXLl3tnaQDOzcCcl8E5NwNzXgbn3Awsr/PysS6mBQAYSv5fPwBAsoQKAJAsoQIAJEuoAADJEiqJWblyZZx77rkxbty4QX/IrqurKy666KIYN25cTJw4MX74wx/GBx98MLSDJuCkk06KmpqafsuqVavyHisXa9asiZNOOinGjh0bra2t8dJLL+U9Uq7uuOOOA54bH/7F7KPJCy+8EBdffHE0NTVFTU1NPP300/32Z1kWy5Yti8mTJ8cxxxwTs2fPjjfeeCOfYYfQR52XhQsXHvAcmjdvXj7DDqH29vb40pe+FLW1tTFx4sS45JJL4vXXX+93zHvvvReLFi2KE044IY477ri49NJLo6enp2ozCZXE7Nu3Ly677LK4/vrrB9zf19cXF110Uezbty82btwYjz32WKxduzaWLVs2xJOm4c4774wdO3ZUlhtuuCHvkYbck08+GW1tbbF8+fL485//HDNnzoy5c+fGzp078x4tV5/73Of6PTf++Mc/5j1SLvbu3RszZ86MNWvWDLj/vvvuiwcffDAefvjh2LRpUxx77LExd+7ceO+994Z40qH1UeclImLevHn9nkOPP/74EE6Yjw0bNsSiRYvixRdfjGeffTbef//9mDNnTuzdu7dyzJIlS+I3v/lN/OpXv4oNGzbE9u3b4xvf+Eb1hspI0qOPPpoVi8UDtv/ud7/LRowYkXV3d1e2PfTQQ1ldXV1WLpeHcML8TZkyJfvJT36S9xi5O/vss7NFixZV1vv6+rKmpqasvb09x6nytXz58mzmzJl5j5GciMieeuqpyvr+/fuzxsbG7P77769s2717d1YoFLLHH388hwnz8d/nJcuybMGCBdnXvva1XOZJyc6dO7OIyDZs2JBl2b+fH6NHj85+9atfVY559dVXs4jIOjo6qjKDd1SGmY6OjjjttNNi0qRJlW1z586NUqkUr7zySo6T5WPVqlVxwgknxBlnnBH333//UfcR2L59++Lll1+O2bNnV7aNGDEiZs+eHR0dHTlOlr833ngjmpqaoqWlJa666qro6urKe6TkbN26Nbq7u/s9f4rFYrS2th71z5+IiPXr18fEiRNj+vTpcf3118euXbvyHmnI9fb2RkREfX19RES8/PLL8f777/d7zpx66qlx4oknVu05U7VfpqU6uru7+0VKRFTWu7u78xgpNzfeeGOceeaZUV9fHxs3boylS5fGjh07YvXq1XmPNmT+8Y9/RF9f34DPiddeey2nqfLX2toaa9eujenTp8eOHTtixYoVcf7558df//rXqK2tzXu8ZPzn34yBnj9H278n/23evHnxjW98I6ZOnRpbtmyJH//4xzF//vzo6OiIkSNH5j3ekNi/f38sXrw4/u///i8+//nPR8S/nzNjxow54BrKaj5nhMoQuPXWW+Pee+896DGvvvrqUXux34d9nHPV1tZW2Xb66afHmDFj4rvf/W60t7f76euj3Pz58yt/Pv3006O1tTWmTJkSv/zlL+Paa6/NcTKGiyuuuKLy59NOOy1OP/30OPnkk2P9+vVx4YUX5jjZ0Fm0aFH89a9/zf36LqEyBG666aZYuHDhQY9paWk5pPtqbGw84Bsd/7naurGx8RPNl5L/5Vy1trbGBx98EG+++WZMnz69CtOlZ8KECTFy5MgDrrjv6ek5Ip4Ph8v48ePjlFNOic2bN+c9SlL+8xzp6emJyZMnV7b39PTEF77whZymSlNLS0tMmDAhNm/efFSEyg9+8IP47W9/Gy+88EJ8+tOfrmxvbGyMffv2xe7du/u9q1LNf3OEyhBoaGiIhoaGw3Jfs2bNipUrV8bOnTtj4sSJERHx7LPPRl1dXcyYMeOwPEae/pdz1dnZGSNGjKicl6PBmDFj4qyzzop169bFJZdcEhH/frt23bp18YMf/CDf4RKyZ8+e2LJlS1x99dV5j5KUqVOnRmNjY6xbt64SJqVSKTZt2jToNw+PVm+99Vbs2rWrX9AdibIsixtuuCGeeuqpWL9+fUydOrXf/rPOOitGjx4d69ati0svvTQiIl5//fXo6uqKWbNmVWUmoZKYrq6ueOedd6Krqyv6+vqis7MzIiKmTZsWxx13XMyZMydmzJgRV199ddx3333R3d0dt912WyxatOio+rijo6MjNm3aFBdccEHU1tZGR0dHLFmyJL75zW/G8ccfn/d4Q6qtrS0WLFgQX/ziF+Pss8+OBx54IPbu3RvXXHNN3qPl5uabb46LL744pkyZEtu3b4/ly5fHyJEj48orr8x7tCG3Z8+efu8kbd26NTo7O6O+vj5OPPHEWLx4cdx9993xmc98JqZOnRq33357NDU1VcL3SHWw81JfXx8rVqyISy+9NBobG2PLli1xyy23xLRp02Lu3Lk5Tl19ixYtil/84hfx61//OmprayvXnRSLxTjmmGOiWCzGtddeG21tbVFfXx91dXVxww03xKxZs+Kcc86pzlBV+S4Rn9iCBQuyiDhgef755yvHvPnmm9n8+fOzY445JpswYUJ20003Ze+//35+Q+fg5ZdfzlpbW7NisZiNHTs2++xnP5vdc8892XvvvZf3aLn46U9/mp144onZmDFjsrPPPjt78cUX8x4pV5dffnk2efLkbMyYMdmnPvWp7PLLL882b96c91i5eP755wf8N2XBggVZlv37K8q33357NmnSpKxQKGQXXnhh9vrrr+c79BA42Hn55z//mc2ZMydraGjIRo8enU2ZMiW77rrr+v0sxJFqoHMSEdmjjz5aOeZf//pX9v3vfz87/vjjs3HjxmVf//rXsx07dlRtppr/PxgAQHL8jgoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECy/h9gdOHcK4VzDwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import jax.numpy as jnp\n",
    "import jax\n",
    "\n",
    "m = 1000\n",
    "prngkey = jax.random.PRNGKey(10)\n",
    "Y = jax.random.normal(prngkey, (2, m)) * 3\n",
    "Y = Y + jax.random.uniform(prngkey, (2, m)) * 5\n",
    "Y = jnp.array(Y)  # 元データ\n",
    "X = jnp.vstack((Y, jnp.array([1,] * m)))  # 中心を計算する用\n",
    "n = 2 + 1\n",
    "assert X.shape == (n, m)\n",
    "\n",
    "# Frank-Wolfe アルゴリズム\n",
    "vmap_dot = jax.vmap(jnp.dot)\n",
    "\n",
    "# >>>> 初期化 >>>>\n",
    "# 0.1\n",
    "u = jnp.ones(m) / m\n",
    "\n",
    "# 0.2\n",
    "epsilon = 0.01  # 精度\n",
    "\n",
    "# 0.3：コレスキー分解\n",
    "XUX = X @ jnp.diag(u) @ X.T\n",
    "L = jnp.linalg.cholesky(XUX)\n",
    "np.testing.assert_allclose(XUX, L @ L.T, atol=1e-5)\n",
    "\n",
    "# 0.4\n",
    "H = jnp.linalg.inv(XUX)\n",
    "vmap_xHx = jax.vmap(lambda x: x @ H @ x)\n",
    "w = vmap_xHx(X.T)\n",
    "assert w.shape == (m,)\n",
    "\n",
    "# <<<< 初期化 <<<<\n",
    "\n",
    "# @jax.jit\n",
    "def one_step_WA_pos(u, L, w):\n",
    "    # １ステップWolfe-Atwoodアルゴリズム。最大化方向の更新\n",
    "    # 1. iの計算\n",
    "    i = w.argmax()\n",
    "    cur_epsilon = ((w - n) / n)[i]\n",
    "    w_i = w[i]\n",
    "    x_i = X[:, i]\n",
    "    u_i = u[i]\n",
    "\n",
    "    # 3. 更新幅の計算とuの更新\n",
    "    lam = (w_i - n) / ((n - 1) * w_i)\n",
    "    lam = jnp.maximum(lam, -u_i)\n",
    "    u = u.at[i].set(u_i + lam)\n",
    "    u = u / (1 + lam)\n",
    "\n",
    "    # 5. wの更新\n",
    "    z_i = jnp.linalg.solve(L, x_i)\n",
    "    H_x_i = jnp.linalg.solve(L.T, z_i)\n",
    "    xx2 = (H_x_i @ X) ** 2\n",
    "    w_updated = (1 + lam) * (w - (lam / (1 + lam * w_i) * xx2))\n",
    "\n",
    "    # 4. 下三角行列の更新\n",
    "    L = cholupdate(L, x_i, lam)\n",
    "\n",
    "    return cur_epsilon, u, L, w_updated\n",
    "\n",
    "\n",
    "# @jax.jit\n",
    "def one_step_WA_neg(u, L, w):\n",
    "    # １ステップWolfe-Atwoodアルゴリズム。最悪値の除去方向の更新\n",
    "    # 1. iの計算\n",
    "    w_ignore_u0 = w + (u == 0) * jnp.inf\n",
    "    i = w_ignore_u0.argmin()\n",
    "    cur_epsilon = ((n - w) / n)[i]\n",
    "    w_i = w[i]\n",
    "    x_i = X[:, i]\n",
    "    u_i = u[i]\n",
    "\n",
    "    # 3. 更新幅の計算とuの更新\n",
    "    lam = (w_i - n) / ((n - 1) * w_i)\n",
    "    lam = jnp.maximum(lam, -u_i)\n",
    "    u = u.at[i].set(u_i + lam)\n",
    "    u = u / (1 + lam)\n",
    "\n",
    "    # 5. wの更新\n",
    "    z_i = jnp.linalg.solve(L, x_i)\n",
    "    H_x_i = jnp.linalg.solve(L.T, z_i)\n",
    "    xx2 = (H_x_i @ X) ** 2\n",
    "    w_updated = (1 + lam) * (w - (lam / (1 + lam * w_i) * xx2))\n",
    "\n",
    "    # 4. 下三角行列の更新\n",
    "    L = cholupdate(L, x_i, lam)\n",
    "\n",
    "    return cur_epsilon, u, L, w_updated\n",
    "\n",
    "\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import HTML\n",
    "from celluloid import Camera\n",
    "\n",
    "\n",
    "def plot_ellipsoid(cur_epsilon, w, u):\n",
    "    # データの描画\n",
    "    x, y = Y\n",
    "    p1 = plt.scatter(x, y, s=3, color=\"lightblue\")\n",
    "\n",
    "    delta = cur_epsilon\n",
    "    thresh = ((delta ** 2) * (n ** 2)) / 4\n",
    "    thresh = 1 + delta * n / 2 - jnp.sqrt(delta * n - delta + thresh)\n",
    "    thresh = n * thresh\n",
    "    essential_data = w >= thresh\n",
    "\n",
    "    # Core-setの描画\n",
    "    x, y = Y[:, essential_data]\n",
    "    num_coreset = essential_data.sum()\n",
    "    p2 = plt.scatter(x, y, s=30, color=\"red\")\n",
    "\n",
    "    # 楕円の描画\n",
    "    y_ = Y @ u\n",
    "    theta = np.linspace(0, 2 * np.pi, 100)\n",
    "    x = np.cos(theta)\n",
    "    y = np.sin(theta)\n",
    "    XUX = X @ jnp.diag(u) @ X.T\n",
    "    H = jnp.linalg.inv(XUX)\n",
    "    H = H[:-1, :-1]\n",
    "    L_plot = np.linalg.cholesky(H)\n",
    "    L_plot_T = np.linalg.inv(L_plot).T\n",
    "    x2, y2 = [], []\n",
    "    for _x, _y in zip(x, y):\n",
    "        v = y_ + np.sqrt(n) * np.dot(L_plot_T, np.array([_x, _y]))\n",
    "        x2.append(v[0])\n",
    "        y2.append(v[1])\n",
    "    p3 = plt.scatter(x2, y2, s=30, color=\"green\")\n",
    "\n",
    "    plt.legend([p1, p2, p3], [f\"Dataset: size={m}\", f\"Core set: size={num_coreset}\", \"Minimum Volume Ellipsoid\"])\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "camera = Camera(fig)\n",
    "\n",
    "i = w.argmax()\n",
    "j = w.argmin()\n",
    "cur_epsilon_pos = ((w - n) / n)[i]\n",
    "cur_epsilon_neg = ((n - w) / n)[i]\n",
    "count = 0\n",
    "while True:\n",
    "    if count < 30:\n",
    "        plot_ellipsoid(cur_epsilon_pos, w, u)\n",
    "        camera.snap()\n",
    "\n",
    "    cur_epsilon_pos, u, L, w = one_step_WA_pos(u, L, w)\n",
    "    cur_epsilon_neg, u, L, w = one_step_WA_neg(u, L, w)\n",
    "    if max(cur_epsilon_pos, cur_epsilon_neg) <= epsilon:\n",
    "        break\n",
    "    count += 1\n",
    "\n",
    "animation = camera.animate()\n",
    "animation.save(\"./figs/MVEE-coreset.gif\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![fig](figs/MVEE-coreset.gif)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 初期化について\n",
    "\n",
    "これまでのアルゴリズムは一様に初期化を行っていましたが、計算量はかかりますがもっと良い初期化が存在します。\n",
    "まず、次のBHアルゴリズムを走らせます\n",
    "\n",
    "0. 適当な非ゼロの定数$c_1\\in \\mathbb{R}^n$を選択し、$j=1$とします。\n",
    "1. $\\overline{z}_j=\\arg\\max_{x \\in X} c_j^T x$と$\\underline{z}_j = \\arg\\min_{x \\in X} c_j^T x$とします。$y_j:=\\overline{z}_j - \\underline{z}_j$とします。\n",
    "2. $j<n$なら、$y_1, \\dots, y_j$に対して直角な$c_{j+1}$を選択して、$j$を一つ増やします。そして１に戻ります。\n",
    "\n",
    "これによって得られた集合$Z:=\\{\\overline{z}_j, \\underline{z}_j: j=1, \\dots, n\\}$について一様な重みを割り当てて作る$u_{KY}$で初期化をします（KumarとYildirimが考えたのでKYとつけてます）。\n",
    "\n",
    "これによって初期化されたアルゴリズムはコアセットを少しだけ小さくできます（3.2章参照）。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGdCAYAAAA8F1jjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAbmElEQVR4nO3da4xU9f348c9yW0R2B1cWlq0rgigWL3grK1aTWgmXWOMtRo1aMMYaRY2sVxoFUXFRE2psiDxRwQdV6wNtbFpTJYJpWDFqNo1WjVDsYmAXxbAjtCy6nt8D/47un4tYnT3fxdcr+SbOmTMzH48j887MmaEiy7IsAAAS1C/vAQAA9kSoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkKwBeQ/wfX3xxRexcePGqKqqioqKirzHAQD2QZZl8emnn0Z9fX3067fn9036fKhs3LgxGhoa8h4DAPgfbNiwIQ455JA9Xt/nQ6WqqioivvwXra6uznkaAGBfFIvFaGhoKL2O70mfD5WvPu6prq4WKgDQx3zbaRtOpgUAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFllDZVXXnklzj777Kivr4+Kiop47rnnelw/a9asqKio6LGmT59ezpEAgD6krKGyffv2mDhxYixZsmSP+0yfPj02bdpUWk8++WQ5RwIA+pAB5bzzGTNmxIwZM/a6T2VlZdTV1ZVzDACgj8r9HJWVK1fGiBEjYvz48XHNNdfEli1b9rp/V1dXFIvFHgsA2D/lGirTp0+PJ554IlasWBH3339/rFq1KmbMmBHd3d17vE1zc3MUCoXSamho6MWJAYDeVJFlWdYrD1RREc8++2yce+65e9znX//6Vxx++OHx0ksvxZlnnrnbfbq6uqKrq6t0uVgsRkNDQ3R2dkZ1dfUPPTYAUAbFYjEKhcK3vn7n/tHPN40dOzaGDx8ea9eu3eM+lZWVUV1d3WMBAPunpELlww8/jC1btsSoUaPyHgUASEBZv/Wzbdu2Hu+OrF+/PlpbW6OmpiZqampiwYIFccEFF0RdXV2sW7cubr311hg3blxMmzatnGMBAH1EWUPl9ddfjzPOOKN0uampKSIiZs6cGY888kj84x//iOXLl8fWrVujvr4+pk6dGvfcc09UVlaWcywAoI/otZNpy2VfT8YBANLRJ0+mBQD4JqECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJCssobKK6+8EmeffXbU19dHRUVFPPfccz2uz7Is5s2bF6NGjYoDDjggpkyZEu+//345RwIA+pCyhsr27dtj4sSJsWTJkt1e/8ADD8TDDz8cS5cujTVr1sSBBx4Y06ZNix07dpRzLACgjxhQzjufMWNGzJgxY7fXZVkWDz30UNxxxx1xzjnnRETEE088ESNHjoznnnsuLr744nKOBgD0Abmdo7J+/fpob2+PKVOmlLYVCoVobGyMlpaWPd6uq6srisVijwUA7J9yC5X29vaIiBg5cmSP7SNHjixdtzvNzc1RKBRKq6GhoaxzAgD56XPf+pk7d250dnaW1oYNG/IeCQAok9xCpa6uLiIiOjo6emzv6OgoXbc7lZWVUV1d3WMBAPun3EJlzJgxUVdXFytWrChtKxaLsWbNmpg8eXJeYwEACSnrt362bdsWa9euLV1ev359tLa2Rk1NTRx66KFx4403xr333htHHHFEjBkzJu68886or6+Pc889t5xjAQB9RFlD5fXXX48zzjijdLmpqSkiImbOnBnLli2LW2+9NbZv3x6/+c1vYuvWrXHaaafFCy+8EIMHDy7nWABAH1GRZVmW9xDfR7FYjEKhEJ2dnc5XAYA+Yl9fv/vct34AgB8PoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECycg+Vu+66KyoqKnqso446Ku+xAIAEDMh7gIiIo48+Ol566aXS5QEDkhgLAMhZEkUwYMCAqKury3sMACAxuX/0ExHx/vvvR319fYwdOzYuvfTSaGtry3skACABub+j0tjYGMuWLYvx48fHpk2bYsGCBXH66afHW2+9FVVVVbvs39XVFV1dXaXLxWKxN8cFAHpRRZZlWd5DfNPWrVtj9OjRsXjx4rjyyit3uf6uu+6KBQsW7LK9s7Mzqqure2NEAOB7KhaLUSgUvvX1O4mPfr5p2LBhceSRR8batWt3e/3cuXOjs7OztDZs2NDLEwIAvSW5UNm2bVusW7cuRo0atdvrKysro7q6uscCAPZPuYfKzTffHKtWrYoPPvggVq9eHeedd170798/LrnkkrxHAwBylvvJtB9++GFccsklsWXLlqitrY3TTjstXn311aitrc17NAAgZ7mHylNPPZX3CABAonL/6AcAYE+ECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAspIIlSVLlsRhhx0WgwcPjsbGxnjttdfyHgkASEDuofL0009HU1NTzJ8/P958882YOHFiTJs2LTZv3pz3aABAznIPlcWLF8dVV10VV1xxRUyYMCGWLl0aQ4YMicceeyzv0QCAnOUaKjt37ow33ngjpkyZUtrWr1+/mDJlSrS0tOz2Nl1dXVEsFnssAGD/lGuofPzxx9Hd3R0jR47ssX3kyJHR3t6+29s0NzdHoVAorYaGht4YFQDIQe4f/XxXc+fOjc7OztLasGFD3iMBAGUyIM8HHz58ePTv3z86Ojp6bO/o6Ii6urrd3qaysjIqKyt7YzwAIGe5vqMyaNCgOOmkk2LFihWlbV988UWsWLEiJk+enONkAEAKcn1HJSKiqakpZs6cGSeffHJMmjQpHnroodi+fXtcccUVeY8GAOQs91C56KKL4qOPPop58+ZFe3t7HH/88fHCCy/scoItAPDjU5FlWZb3EN9HsViMQqEQnZ2dUV1dnfc4AMA+2NfX7z73rR8A4MdDqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJCvXUDnssMOioqKix1q0aFGeIwEACRmQ9wB33313XHXVVaXLVVVVOU4DAKQk91CpqqqKurq6vMcAABKU+zkqixYtioMPPjhOOOGEePDBB+Pzzz/f6/5dXV1RLBZ7LABg/5TrOyo33HBDnHjiiVFTUxOrV6+OuXPnxqZNm2Lx4sV7vE1zc3MsWLCgF6cEAPJSkWVZ9kPe4e233x7333//Xvd555134qijjtpl+2OPPRZXX311bNu2LSorK3d7266urujq6ipdLhaL0dDQEJ2dnVFdXf39hgcAekWxWIxCofCtr98/eKh89NFHsWXLlr3uM3bs2Bg0aNAu299+++045phj4t13343x48fv0+Pt678oAJCOfX39/sE/+qmtrY3a2tr/6batra3Rr1+/GDFixA88FQDQF+V2jkpLS0usWbMmzjjjjKiqqoqWlpaYM2dOXHbZZXHQQQflNRYAkJDcQqWysjKeeuqpuOuuu6KrqyvGjBkTc+bMiaamprxGAgASk1uonHjiifHqq6/m9fAAQB+Q+++oAADsiVABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBIllABAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWWULlYULF8app54aQ4YMiWHDhu12n7a2tjjrrLNiyJAhMWLEiLjlllvi888/L9dIAEAfM6Bcd7xz58648MILY/LkyfHoo4/ucn13d3ecddZZUVdXF6tXr45NmzbFr3/96xg4cGDcd9995RoLAOhDKrIsy8r5AMuWLYsbb7wxtm7d2mP7X//61/jVr34VGzdujJEjR0ZExNKlS+O2226Ljz76KAYNGrRP918sFqNQKERnZ2dUV1f/0OMDAGWwr6/fuZ2j0tLSEscee2wpUiIipk2bFsViMd5+++28xgIAElK2j36+TXt7e49IiYjS5fb29j3erqurK7q6ukqXi8VieQYEAHL3nd5Ruf3226OiomKv69133y3XrBER0dzcHIVCobQaGhrK+ngAQH6+0zsqN910U8yaNWuv+4wdO3af7quuri5ee+21Hts6OjpK1+3J3Llzo6mpqXS5WCyKFQDYT32nUKmtrY3a2tof5IEnT54cCxcujM2bN8eIESMiIuLFF1+M6urqmDBhwh5vV1lZGZWVlT/IDABA2sp2jkpbW1t88skn0dbWFt3d3dHa2hoREePGjYuhQ4fG1KlTY8KECXH55ZfHAw88EO3t7XHHHXfE7NmzhQgAEBFl/HryrFmzYvny5btsf/nll+MXv/hFRET8+9//jmuuuSZWrlwZBx54YMycOTMWLVoUAwbsez/5ejIA9D37+vpd9t9RKTehAgB9T/K/owIA8G2ECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAsoQKAJAsoQIAJEuoAADJEioAQLKECgCQLKECACRLqAAAyRIqAECyhAoAkCyhAgAkS6gAAMkSKgBAssoWKgsXLoxTTz01hgwZEsOGDdvtPhUVFbusp556qlwjAQB9zIBy3fHOnTvjwgsvjMmTJ8ejjz66x/0ef/zxmD59eunynqIGAPjxKVuoLFiwICIili1bttf9hg0bFnV1deUaAwDow3I/R2X27NkxfPjwmDRpUjz22GORZdle9+/q6opisdhjAQD7p7K9o7Iv7r777vjlL38ZQ4YMib/97W9x7bXXxrZt2+KGG27Y422am5tL79Z8k2ABgL7jq9ftb3uDIrLv4LbbbssiYq/rnXfe6XGbxx9/PCsUCvt0/3feeWd2yCGH7HWfHTt2ZJ2dnaX1z3/+81tnsizLsiwrzbVhw4a9vu5/p3dUbrrpppg1a9Ze9xk7dux3ucseGhsb45577omurq6orKzc7T6VlZU9rhs6dGhs2LAhqqqqoqKi4n9+7P1BsViMhoaG2LBhQ1RXV+c9Tu4cj685Fj05Hl9zLHpyPL5W7mORZVl8+umnUV9fv9f9vlOo1NbWRm1t7fcabG9aW1vjoIMO2mOk7E6/fv3ikEMOKdtMfVF1dfWP/n+wb3I8vuZY9OR4fM2x6Mnx+Fo5j0WhUPjWfcp2jkpbW1t88skn0dbWFt3d3dHa2hoREePGjYuhQ4fG888/Hx0dHXHKKafE4MGD48UXX4z77rsvbr755nKNBAD0MWULlXnz5sXy5ctLl0844YSIiHj55ZfjF7/4RQwcODCWLFkSc+bMiSzLYty4cbF48eK46qqryjUSANDHlC1Uli1bttffUJk+fXqPH3rj+6usrIz58+d/p4/O9meOx9cci54cj685Fj05Hl9L5VhUZNm3fS8IACAfuf/gGwDAnggVACBZQgUASJZQAQCSJVT2EwsXLoxTTz01hgwZEsOGDdvtPm1tbXHWWWfFkCFDYsSIEXHLLbfE559/3ruD5uSwww6LioqKHmvRokV5j9VrlixZEocddlgMHjw4Ghsb47XXXst7pF5311137fIcOOqoo/Ieq9e88sorcfbZZ0d9fX1UVFTEc8891+P6LMti3rx5MWrUqDjggANiypQp8f777+czbC/4tuMxa9asXZ4v++s3VZubm+NnP/tZVFVVxYgRI+Lcc8+N9957r8c+O3bsiNmzZ8fBBx8cQ4cOjQsuuCA6Ojp6ZT6hsp/YuXNnXHjhhXHNNdfs9vru7u4466yzYufOnbF69epYvnx5LFu2LObNm9fLk+bn7rvvjk2bNpXW9ddfn/dIveLpp5+OpqammD9/frz55psxceLEmDZtWmzevDnv0Xrd0Ucf3eM58Pe//z3vkXrN9u3bY+LEibFkyZLdXv/AAw/Eww8/HEuXLo01a9bEgQceGNOmTYsdO3b08qS949uOR8SXP6PxzefLk08+2YsT9p5Vq1bF7Nmz49VXX40XX3wxPvvss5g6dWps3769tM+cOXPi+eefj2eeeSZWrVoVGzdujPPPP793Btynvy2QPmNPfwnkX/7yl6xfv35Ze3t7adsjjzySVVdXZ11dXb04YT5Gjx6d/e53v8t7jFxMmjQpmz17dulyd3d3Vl9fnzU3N+c4Ve+bP39+NnHixLzHSEJEZM8++2zp8hdffJHV1dVlDz74YGnb1q1bs8rKyuzJJ5/MYcLe9f8fjyzLspkzZ2bnnHNOLvPkbfPmzVlEZKtWrcqy7MvnwsCBA7NnnnmmtM8777yTRUTW0tJS9nm8o/Ij0dLSEscee2yMHDmytG3atGlRLBbj7bffznGy3rNo0aI4+OCD44QTTogHH3zwR/Gx186dO+ONN96IKVOmlLb169cvpkyZEi0tLTlOlo/3338/6uvrY+zYsXHppZdGW1tb3iMlYf369dHe3t7jeVIoFKKxsfFH+Tz5ysqVK2PEiBExfvz4uOaaa2LLli15j9QrOjs7IyKipqYmIiLeeOON+Oyzz3o8P4466qg49NBDe+X5UbZfpiUt7e3tPSIlIkqX29vb8xipV91www1x4oknRk1NTaxevTrmzp0bmzZtisWLF+c9Wll9/PHH0d3dvdv/9u+++25OU+WjsbExli1bFuPHj49NmzbFggUL4vTTT4+33norqqqq8h4vV1/9GbC758mP4c+H3Zk+fXqcf/75MWbMmFi3bl389re/jRkzZkRLS0v0798/7/HK5osvvogbb7wxfv7zn8cxxxwTEV8+PwYNGrTL+Y+99fwQKgm7/fbb4/7779/rPu+8886P6oTAb/oux6epqam07bjjjotBgwbF1VdfHc3Nzbn/PDS9Y8aMGaV/Pu6446KxsTFGjx4df/zjH+PKK6/McTJSdPHFF5f++dhjj43jjjsuDj/88Fi5cmWceeaZOU5WXrNnz4633norqfO3hErCbrrpppg1a9Ze9xk7duw+3VddXd0u3/T46ozturq6/2m+vH2f49PY2Biff/55fPDBBzF+/PgyTJeG4cOHR//+/Xc5O7+jo6PP/nf/oQwbNiyOPPLIWLt2bd6j5O6r50JHR0eMGjWqtL2joyOOP/74nKZKy9ixY2P48OGxdu3a/TZUrrvuuvjzn/8cr7zyShxyyCGl7XV1dbFz587YunVrj3dVeuvPEaGSsNra2qitrf1B7mvy5MmxcOHC2Lx5c4wYMSIiIl588cWorq6OCRMm/CCP0du+z/FpbW2Nfv36lY7F/mrQoEFx0kknxYoVK+Lcc8+NiC/f2l2xYkVcd911+Q6Xs23btsW6devi8ssvz3uU3I0ZMybq6upixYoVpTApFouxZs2aPX6T8Mfmww8/jC1btvQIuf1FlmVx/fXXx7PPPhsrV66MMWPG9Lj+pJNOioEDB8aKFSviggsuiIiI9957L9ra2mLy5Mlln0+o7Cfa2trik08+iba2tuju7o7W1taIiBg3blwMHTo0pk6dGhMmTIjLL788HnjggWhvb4877rgjZs+evd9/9NHS0hJr1qyJM844I6qqqqKlpSXmzJkTl112WRx00EF5j1d2TU1NMXPmzDj55JNj0qRJ8dBDD8X27dvjiiuuyHu0XnXzzTfH2WefHaNHj46NGzfG/Pnzo3///nHJJZfkPVqv2LZtW493j9avXx+tra1RU1MThx56aNx4441x7733xhFHHBFjxoyJO++8M+rr60uBu7/Z2/GoqamJBQsWxAUXXBB1dXWxbt26uPXWW2PcuHExbdq0HKcuj9mzZ8cf/vCH+NOf/hRVVVWl804KhUIccMABUSgU4sorr4ympqaoqamJ6urquP7662Py5MlxyimnlH/Asn+viF4xc+bMLCJ2WS+//HJpnw8++CCbMWNGdsABB2TDhw/Pbrrppuyzzz7Lb+he8sYbb2SNjY1ZoVDIBg8enP30pz/N7rvvvmzHjh15j9Zrfv/732eHHnpoNmjQoGzSpEnZq6++mvdIve6iiy7KRo0alQ0aNCj7yU9+kl100UXZ2rVr8x6r17z88su7/TNi5syZWZZ9+RXlO++8Mxs5cmRWWVmZnXnmmdl7772X79BltLfj8Z///CebOnVqVltbmw0cODAbPXp0dtVVV/X4eYf9ye6OQ0Rkjz/+eGmf//73v9m1116bHXTQQdmQIUOy8847L9u0aVOvzFfx/4YEAEiO31EBAJIlVACAZAkVACBZQgUASJZQAQCSJVQAgGQJFQAgWUIFAEiWUAEAkiVUAIBkCRUAIFlCBQBI1v8BhOBsRD5a3WUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import jax.numpy as jnp\n",
    "import jax\n",
    "\n",
    "m = 1000\n",
    "prngkey = jax.random.PRNGKey(10)\n",
    "Y = jax.random.normal(prngkey, (2, m)) * 3\n",
    "Y = Y + jax.random.uniform(prngkey, (2, m)) * 5\n",
    "Y = jnp.array(Y)  # 元データ\n",
    "X = jnp.vstack((Y, jnp.array([1,] * m)))  # 中心を計算する用\n",
    "n = 2 + 1\n",
    "assert X.shape == (n, m)\n",
    "\n",
    "# Frank-Wolfe アルゴリズム\n",
    "vmap_dot = jax.vmap(jnp.dot)\n",
    "\n",
    "# >>>> 初期化 >>>>\n",
    "# 0.1\n",
    "c = jax.random.normal(key3, (n,))\n",
    "u = jnp.zeros(m)\n",
    "\n",
    "# 一回目\n",
    "cx = c @ X\n",
    "idx_max, idx_min = jnp.argmax(cx), jnp.argmin(cx)\n",
    "u = u.at[idx_max].set(1).at[idx_min].set(1)\n",
    "y1 = X.T[idx_max] - X.T[idx_min]\n",
    "\n",
    "# 二回目\n",
    "c = jnp.cross(y1, c)\n",
    "cx = c @ X\n",
    "idx_max, idx_min = jnp.argmax(cx), jnp.argmin(cx)\n",
    "u = u.at[idx_max].set(1).at[idx_min].set(1)\n",
    "y2 = X.T[idx_max] - X.T[idx_min]\n",
    "\n",
    "# 三回目\n",
    "c = jnp.cross(y1, y2)\n",
    "cx = c @ X\n",
    "idx_max, idx_min = jnp.argmax(cx), jnp.argmin(cx)\n",
    "u = u.at[idx_max].set(1).at[idx_min].set(1)\n",
    "\n",
    "u = u / u.sum()\n",
    "\n",
    "\n",
    "# 0.2\n",
    "epsilon = 0.01  # 精度\n",
    "\n",
    "# 0.3：コレスキー分解\n",
    "XUX = X @ jnp.diag(u) @ X.T\n",
    "L = jnp.linalg.cholesky(XUX)\n",
    "np.testing.assert_allclose(XUX, L @ L.T, atol=1e-5)\n",
    "\n",
    "# 0.4\n",
    "H = jnp.linalg.inv(XUX)\n",
    "vmap_xHx = jax.vmap(lambda x: x @ H @ x)\n",
    "w = vmap_xHx(X.T)\n",
    "assert w.shape == (m,)\n",
    "\n",
    "# <<<< 初期化 <<<<\n",
    "\n",
    "# @jax.jit\n",
    "def one_step_WA_pos(u, L, w):\n",
    "    # １ステップWolfe-Atwoodアルゴリズム。最大化方向の更新\n",
    "    # 1. iの計算\n",
    "    i = w.argmax()\n",
    "    cur_epsilon = ((w - n) / n)[i]\n",
    "    w_i = w[i]\n",
    "    x_i = X[:, i]\n",
    "    u_i = u[i]\n",
    "\n",
    "    # 3. 更新幅の計算とuの更新\n",
    "    lam = (w_i - n) / ((n - 1) * w_i)\n",
    "    lam = jnp.maximum(lam, -u_i)\n",
    "    u = u.at[i].set(u_i + lam)\n",
    "    u = u / (1 + lam)\n",
    "\n",
    "    # 5. wの更新\n",
    "    z_i = jnp.linalg.solve(L, x_i)\n",
    "    H_x_i = jnp.linalg.solve(L.T, z_i)\n",
    "    xx2 = (H_x_i @ X) ** 2\n",
    "    w_updated = (1 + lam) * (w - (lam / (1 + lam * w_i) * xx2))\n",
    "\n",
    "    # 4. 下三角行列の更新\n",
    "    L = cholupdate(L, x_i, lam)\n",
    "\n",
    "    return cur_epsilon, u, L, w_updated\n",
    "\n",
    "\n",
    "# @jax.jit\n",
    "def one_step_WA_neg(u, L, w):\n",
    "    # １ステップWolfe-Atwoodアルゴリズム。最悪値の除去方向の更新\n",
    "    # 1. iの計算\n",
    "    w_ignore_u0 = w + (u == 0) * jnp.inf\n",
    "    i = w_ignore_u0.argmin()\n",
    "    cur_epsilon = ((n - w) / n)[i]\n",
    "    w_i = w[i]\n",
    "    x_i = X[:, i]\n",
    "    u_i = u[i]\n",
    "\n",
    "    # 3. 更新幅の計算とuの更新\n",
    "    lam = (w_i - n) / ((n - 1) * w_i)\n",
    "    lam = jnp.maximum(lam, -u_i)\n",
    "    u = u.at[i].set(u_i + lam)\n",
    "    u = u / (1 + lam)\n",
    "\n",
    "    # 5. wの更新\n",
    "    z_i = jnp.linalg.solve(L, x_i)\n",
    "    H_x_i = jnp.linalg.solve(L.T, z_i)\n",
    "    xx2 = (H_x_i @ X) ** 2\n",
    "    w_updated = (1 + lam) * (w - (lam / (1 + lam * w_i) * xx2))\n",
    "\n",
    "    # 4. 下三角行列の更新\n",
    "    L = cholupdate(L, x_i, lam)\n",
    "\n",
    "    return cur_epsilon, u, L, w_updated\n",
    "\n",
    "\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import HTML\n",
    "from celluloid import Camera\n",
    "\n",
    "\n",
    "def plot_ellipsoid(cur_epsilon, w, u):\n",
    "    # データの描画\n",
    "    x, y = Y\n",
    "    p1 = plt.scatter(x, y, s=3, color=\"lightblue\")\n",
    "\n",
    "    delta = cur_epsilon\n",
    "    thresh = ((delta ** 2) * (n ** 2)) / 4\n",
    "    thresh = 1 + delta * n / 2 - jnp.sqrt(delta * n - delta + thresh)\n",
    "    thresh = n * thresh\n",
    "    essential_data = w >= thresh\n",
    "\n",
    "    # Core-setの描画\n",
    "    x, y = Y[:, essential_data]\n",
    "    num_coreset = essential_data.sum()\n",
    "    p2 = plt.scatter(x, y, s=30, color=\"red\")\n",
    "\n",
    "    # 楕円の描画\n",
    "    y_ = Y @ u\n",
    "    theta = np.linspace(0, 2 * np.pi, 100)\n",
    "    x = np.cos(theta)\n",
    "    y = np.sin(theta)\n",
    "    XUX = X @ jnp.diag(u) @ X.T\n",
    "    H = jnp.linalg.inv(XUX)\n",
    "    H = H[:-1, :-1]\n",
    "    L_plot = np.linalg.cholesky(H)\n",
    "    L_plot_T = np.linalg.inv(L_plot).T\n",
    "    x2, y2 = [], []\n",
    "    for _x, _y in zip(x, y):\n",
    "        v = y_ + np.sqrt(n) * np.dot(L_plot_T, np.array([_x, _y]))\n",
    "        x2.append(v[0])\n",
    "        y2.append(v[1])\n",
    "    p3 = plt.scatter(x2, y2, s=30, color=\"green\")\n",
    "\n",
    "    plt.legend([p1, p2, p3], [f\"Dataset: size={m}\", f\"Core set: size={num_coreset}\", \"Minimum Volume Ellipsoid\"])\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "camera = Camera(fig)\n",
    "\n",
    "i = w.argmax()\n",
    "j = w.argmin()\n",
    "cur_epsilon_pos = ((w - n) / n)[i]\n",
    "cur_epsilon_neg = ((n - w) / n)[i]\n",
    "count = 0\n",
    "while True:\n",
    "    if count < 30:\n",
    "        plot_ellipsoid(cur_epsilon_pos, w, u)\n",
    "        camera.snap()\n",
    "\n",
    "    cur_epsilon_pos, u, L, w = one_step_WA_pos(u, L, w)\n",
    "    cur_epsilon_neg, u, L, w = one_step_WA_neg(u, L, w)\n",
    "    if max(cur_epsilon_pos, cur_epsilon_neg) <= epsilon:\n",
    "        break\n",
    "    count += 1\n",
    "\n",
    "animation = camera.animate()\n",
    "animation.save(\"./figs/MVEE-coreset_KY_init.gif\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![fig](figs/MVEE-coreset.gif)\n",
    "![fig](figs/MVEE-coreset_KY_init.gif)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "8eaa1d0eb28930fc7ef49205b5d1723439e8a6970d242658c88c7a4f508517fb"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}