aseyboldt/covadapt_subspace.ipynb

## covadapt_subspace.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c9a75c9f-f473-4746-b868-098da478e91f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import linalg, special, stats\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "072271ca-a114-44b3-9fb8-c9a758618426",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "k = 50\n",
    "\n",
    "def make_cov(N, eigdist, diagdist):\n",
    "    eigs = eigdist.rvs(N)\n",
    "    diag = diagdist.rvs(N)\n",
    "    vecs = stats.ortho_group(N).rvs()\n",
    "    return np.diag(diag) @ (vecs @ np.diag(eigs) @ vecs.T) @ np.diag(diag)\n",
    "\n",
    "points = np.random.randn(k, N)\n",
    "grads = np.random.randn(k, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "09842717-b41d-40ae-90a1-3161be93c677",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(43)\n",
    "cov = make_cov(N, stats.lognorm(s=2, scale=1), stats.lognorm(s=1, scale=1))\n",
    "#sns.clustermap(cov, center=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2c4364d3-c82a-44cb-ab21-36ae72f56a31",
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_value_grads(cov, k):\n",
    "    draws = np.random.multivariate_normal(np.zeros(len(cov)), cov, size=k)\n",
    "    grads = - linalg.solve(cov, draws.T, assume_a=\"pos\")\n",
    "    return draws, grads.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "41f9aeeb-b148-4ab2-a8d7-b0a6f56baecb",
   "metadata": {},
   "outputs": [],
   "source": [
    "points, grads = draw_value_grads(cov, k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "af6ccd3c-e9d6-4013-b96f-58c5596eb2f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plt.plot(np.log(linalg.eigvalsh(np.cov(grads.T))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3dba24ef-eb98-4afe-ac2a-043790f4dccc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f7d6069d7c0>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAiDElEQVR4nO3deXxV9Z3/8deHhAQIYc1CEgj7LrJF0GJdUFEolWrVqtXWZX6ov7YzrTPT0XEe/bU+fjPtdGynnTpqqbXWpa5VsYoKtsqiIiTKvshO9oUACWTP/fz+yHV+kSYs3pvc5N738/HII+eec3K+n28I75x87/ecY+6OiIhEvx6RLkBERDqHAl9EJEYo8EVEYoQCX0QkRijwRURiRHykCziZlJQUHzFiRKTLEBHpNvLy8ircPbWtbV068EeMGEFubm6kyxAR6TbM7EB72zSkIyISIxT4IiIxQoEvIhIjFPgiIjFCgS8iEiMU+CIiMUKBLyISIxT4IiJdyIptpTyyck+HHFuBLyLShazYVsLv3tvXIcdW4IuIdCHH65tJSuyYmyAo8EVEupBj9U30VeCLiES/4wp8EZHYcKy+SUM6IiKxQEM6IiIxIBBwKo7VMygpoUOOH3Lgm9l4M9vQ6qPKzL57wj4XmdnRVvv8INR2RUSiTcHhWuoaA4xN69shxw/57wZ33wlMAzCzOKAQeLmNXVe7+8JQ2xMRiVYFh2sAyB7cp0OOH+4hnUuAPe7e7hNXRESkbWXV9QCkJffqkOOHO/CvB55pZ9t5ZrbRzN4ws8ntHcDMFptZrpnllpeXh7k8EZGuq/zTwO+X2CHHD1vgm1kCcCXwQhubPwKGu/tU4FfAK+0dx92XuHuOu+ekprb5HF4RkaiUe6CS1OREkrvBLJ35wEfuXnriBnevcvdjweVlQE8zSwlj2yIi3d6Wwiq+MHowZtYhxw9n4N9AO8M5ZjbEgj0ws1nBdg+FsW0RkW7N3SmrrmNI/44Zv4cwzNIBMLM+wGXAHa3W3Qng7o8A1wB3mVkTUAtc7+4ejrZFRKJBxbEGGpudIf26eOC7ew0w+IR1j7RafhB4MBxtiYhEoyfXtkxunJ49sMPa0JW2IiIR5u48t/4g549JYerQ/h3WjgJfRCTCdpcdo7SqnoVnZ3TYG7agwBcRibhVuyoAOH9sx05eVOCLiERQdV0jv39/P+PTkxk6sGNuqfApBb6ISAS9s7Ocg5U13DN/Qoe3pcAXEYmgNbvKSYzv0eHDOaDAFxGJmOKjtbz0USFfO2cYPeM6Po4V+CIiEfL8+gKaAs5tc0Z2SnsKfBGRCKg83sBvVu/lsknpjEhJ6pQ2FfgiIhHwyMo91DQ08f3Lx3damwp8EZFOdqy+iT/mFTBv0hDGpid3WrsKfBGRTvb8+nwOHW/gtvM7Z+z+Uwp8EZFOVHC4hv94ayfnjBjIrJGDOrVtBb6ISCdasmovtY3N3L/orE5vW4EvItJJ8g5U8uTaA1w+OZ0JQzpv7P5TCnwRkU6QX1nDbY/nMqRfLx64dmqH3hWzPQp8EZEOVtPQxI2PrqU54Dxx2yySe/WMSB0d82h0EREBoL6pmcVP5JFfWctDX5/RqdMwT6QzfBGRDvR/lm5lze4Kfnz1FBZMyYhoLQp8EZEO8sGeQ7yQV8AtXxjBDbOyI12OAl9EpCPkHajktsfXMyolie9dOi7S5QBhCnwz229mm81sg5nltrHdzOy/zGy3mW0ysxnhaFdEpCt6Z0cZ33xsPYOSEnj6f82mf5/IvEl7onC+aXuxu1e0s20+MDb4MRt4OPhZRCRq1Dc1808vbuKVDUWMTevLY7ecQ1pyr0iX9T86a5bOIuAJd3dgrZkNMLMMdy/upPZFRDpUQ1OAu5/byOubi/nbuWP41twxJMbHRbqszwjXGL4Dy80sz8wWt7E9C8hv9boguO6vmNliM8s1s9zy8vIwlSci0nFKq+q466k8Xt9czD8vmMDd88Z3ubCH8J3hz3H3IjNLA1aY2Q53X9Vqe1uXlHlbB3L3JcASgJycnDb3ERHpKgoO13DTox9SeKSWe+ZPYPEFoyNdUrvCEvjuXhT8XGZmLwOzgNaBXwAMa/V6KFAUjrZFRCLl4KEavvzgGhqbAzy7+FxmDu/cu1+eqZCHdMwsycySP10G5gFbTtjtVeAbwdk65wJHNX4vIt3Z1qKj3PL4OgIB5+m/md3lwx7Cc4afDrwcvBFQPPAHd3/TzO4EcPdHgGXAAmA3UAPcGoZ2RUQ6nbvz2zX7+Ndl2+mbEM+vbpzO9OyBkS7rtIQc+O6+F5jaxvpHWi078K1Q2xIRiaSDh2q4/7VtvL29lMsmpfOTq6cwuG9ipMs6bbp5mojIafj44GG++dg6Gpud+xZM5G++ODIitzgOhQJfROQk3J23tpZy9/Mb6N+7J7+79RwmDOkX6bI+FwW+iEg7jtQ0sPjJPNbtq2Rcel+eun02af26zpWzZ0qBLyJyAnfnmXX5PLxyN6VH67lvwUS+8YXhXfJiqjOhwBcRaaWsqo6H3t3D4+/vZ9qwAfzHNVM5d9TgSJcVFgp8EZGg5VtL+M4zH1PfFOCGWcP4169MoUeP7vXG7Mko8EUk5uVX1vDA8p0s3VDEqNQkltw8kzFpkXsUYUdR4ItIzAoEnGfX5/PjN7bT2BzgjgtG8XeXjqVPQnRGY3T2SkTkJNyd59bn859vf0JpVT2zRg7i/kWTu+10y9OlwBeRmLKjpIofLN3Kun2VTM8ewD8vmMiVUzO73UVUn4cCX0RiQsWxeh5/bz+PrNxDcq94fnz1FK7LGUZcFL0peyoKfBGJakdrGnlo5W5+t2Y/Dc0B5k5I42fXTmVgUkKkS+t0CnwRiUrVdY38fMUnvJhbwLGGJq6YPIQ7LxzN1GEDIl1axCjwRSSqBALOb1bv5Ter93HoeD1XTs3kzgtHMzEjut+QPR0KfBGJCoGA8+bWEn69ai8b848wc/hAHrxxetRcJRsOCnwR6daq6hpZ+nEhv//gALvLjpHRvxc/vnoK184cSnxcyA/1iyoKfBHplvaUH+OptQd4Ma+A6romJgxJ5t+umsLXzomtmTdnQoEvIt1KfmUN//3Obl7IKyDOjHmT01l8wSimZPWPibn0oVDgi0i3UHC4hp+8sYM3t5TQw4ybzx3Ot+eOIaUbPWIw0hT4ItJlHa9vYvm2Ep5bn8/avZXE9zCuzRnG314yhoz+vSNdXrcTcuCb2TDgCWAIEACWuPsvT9jnImApsC+46iV3vz/UtkUk+jQHnDW7K3j5owLe2lpKbWMzQwf25nuXjuPqGVkMG9Qn0iV2W+E4w28C/t7dPzKzZCDPzFa4+7YT9lvt7gvD0J6IRKGm5gBvbi3hoXf2sK24in694rlqRhZXTc9iZvbAqLovfaSEHPjuXgwUB5erzWw7kAWcGPgiIp/R0BTgvT0VvLG5mBXbSjlc08iolCR++tWzWTQ9s9s/UrCrCesYvpmNAKYDH7ax+Twz2wgUAf/g7lvbOcZiYDFAdnZ2OMsTkS6irrGZF/IKeGzNPvZVHCc5MZ65E9NYMCWDSyema1plBwlb4JtZX+CPwHfdveqEzR8Bw939mJktAF4BxrZ1HHdfAiwByMnJ8XDVJyKR9enY/PPr83l3ZxnHG5oZn57Mr26YzmWT0unVU2fzHS0sgW9mPWkJ+6fd/aUTt7f+BeDuy8zsITNLcfeKcLQvIl2Tu7OjpJrlW0t58aN88itr6ZsYz6LpWcw/awjnj0nR3PlOFI5ZOgb8Ftju7j9vZ58hQKm7u5nNAnoAh0JtW0S6prKqOlbtquD37+9nc+FRAGaNHMQ9V0zkkolpOpuPkHCc4c8BbgY2m9mG4Lp/BrIB3P0R4BrgLjNrAmqB691dwzUiUaShKcAHew/x5Af7WflJOY3NTkrfRH745UksnJqpC6S6gHDM0lkDnPRvMnd/EHgw1LZEpGtpDjird5Xz7s5y/rSxiEPHG0hNTuSmc4dzXc4wxqcnazplF6IrbUXkjJVW1fFCbj7PrMun8EgtifE9uGh8Kl+ZlsWF41Ppk6Bo6Yr0ryIip6W8up7nc/N5dv1B8itrAZgzZjD/8qWJzJ2Ypjnz3YACX0TaVVZdxx/zCnlzSzGbCo/i3hLyN587nEsmpjM6tW+kS5QzoMAXkc8oPlrLM+vyWbmz7H9Cfnx6Mt+9ZBxfOjuDMWkK+e5KgS8ilFXX8d7uCpZtLuHP20tx4Jzhg/jO3LFcOTVTIR8lFPgiMcjdKTxSy8pPynnl40LW7z8MQErfBBZfMJqvz87WXSmjkAJfJIZU1TWyYmspz+Xms25fJQCjUpP4h3njuHBcGpMz+2kaZRRT4ItEMXfnk9JjLN9awjs7y9hYcJTmgDN0YG/+Yd445k0ewti0vrq9QYxQ4ItEoV2l1bz0cSF/2lhEweGWKZRThw3grgtHc/GENGZkD1DIxyAFvkiUKK+u59WNRbz8cQFbCquI62F8cWwK35k7hvPHppI1QI8EjHUKfJFurPBILcs2tTw8JO/gYZoDzpSs/vxg4SSunKb718hnKfBFuhF3Z2dpNcs2FbN8Wyk7S6txhwlDkrnjglFcPSOLMWnJkS5TuigFvkgX1xxwPtx7iBXbS3lnRxn7D9VgBueOHMz3Lh3H/LOGMDZdIS+npsAX6YJqG5rZVHCEFdtK+dOmIkqr6unVswdnDx3AHReOZu6ENNL79Yp0mdLNKPBFuojahmaWbS7m7e2lrN5VwbH6JnrGGReNT2PRtEwunajHAEpoFPgiEVRV18g7O8pYvq2UlTvLOVbfRGb/Xlw2KZ2FZ2eQM3wQ/fv0jHSZEiUU+CKdrORoHSu2l7J8awlr9x4KPhkqgS9PzWDh2ZmcN2qwrnaVDqHAF+kEBYdreOmjQt7cUsK24ioARqYkcduckcybnM60YQOJU8hLB1Pgi3SQyuMNvL29lKUbCnl/zyHcYUb2AP7x8vFcPrnlXvK62lU6kwJfJIz2VRxnxbaSlguhDhwm4DB8cB/+7pKxLJiSwThNn5QIUuCLhKChKUDugUre3VnOqk/K2VFSDcDEjH58e+5Y5k1KZ3JmP53JS5cQlsA3syuAXwJxwKPu/pMTtltw+wKgBrjF3T8KR9sine1obSNrdlWwYlsJf95RRnVdE3E9jJzhA/nBwknMm5zO0IG6l7x0PSEHvpnFAf8NXAYUAOvN7FV339Zqt/nA2ODHbODh4GeRLs/dOVhZw6pdFby7o4zVuytoaAqQlBDHvMlDuHxyOuePTaVvov5glq4tHD+hs4Dd7r4XwMyeBRYBrQN/EfCEuzuw1swGmFmGuxeHoX2RsAsEnG3FVby3u4KnPjxAfmXLLYaHD+7DjbOy+fLUDKYOHUB8XI8IVypy+sIR+FlAfqvXBfz12Xtb+2QBfxX4ZrYYWAyQnZ0dhvJETk8g4OQdPMzz6/NZvauCkqo6ACZl9OP+RaOYPXIw49I1s0a6r3AEfls//f459mlZ6b4EWAKQk5PT5j4i4dIccD7cd4g3Npfw1tYSyqrr6ZMQxwVjU7lsUjo5IwYyfHBSpMsUCYtwBH4BMKzV66FA0efYR6RTNDUHWLe/khdzC1i9u4Ly6pYbk82dkMb8szK4eEKaxuMlKoXjp3o9MNbMRgKFwPXAjSfs8yrw7eD4/mzgqMbvpTM1NAV4f08F7+4s5+WPCzla20hyr3guGp/GvEnpXDIxjT4JCnmJbiH/hLt7k5l9G3iLlmmZj7n7VjO7M7j9EWAZLVMyd9MyLfPWUNsVOZWGpgDv7ixjxbZS3tpaQlVdE/E9jIsnpHHV9CwuHp9G7wTdfVJiR1hOadx9GS2h3nrdI62WHfhWONoSOZny6nqWbyvh3Z3lvL+7guMNzfRNjGfepHQWTs1g9sjBJGm4RmKUfvKl26s83sCyzcUs3VDI+v2HAUjvl8hVM7K4cFwaF4xLITFeZ/IiCnzplg4dq+ft7aX8ZUcZqz6poLaxmZEpSdx92TjmTU5nXFqybjEscgIFvnQbFcfqeeXjQt7dWc7avYdoCjiZ/Xvx1ZlZfC0nm8mZ/RTyIiehwJcu7cCh46zYVsqft5eRd+AwDc0Bxqcnc9v5I1k0LZNJGboxmcjpUuBLl3OkpoHfrtnH65uK2VtxHIAJQ5K5Zc4IrssZxpi0vhGuUKR7UuBLlxAIOGt2V/D4+/tZ9Uk5TQHni2NT+MZ5w7l4QpqudhUJAwW+RIy7s6Wwipc/LuSNLcUUH60jvV8it58/kkXTspiU2S/SJYpEFQW+dLptRVU8uXY/qz6poPBILQlxPbhgXCr/ePl45p+VoYuhRDqIAl86RU1DE6t3VfBCbgF/2VFKfFwPLpmQxh0XjmLh2ZkMSkqIdIkiUU+BLx2mvqmZ9fsO89h7+1izq4KG5gCDkxK448LR3DZnJKnJiZEuUSSmKPAlrNydrUVVPLX2AK9vLqa6ronePeO46dzhzJ2QxuxRg+iph4aIRIQCX8Kiqq6RP20s4qm1B9leXEWfhDiuOGsI8yYNYc6YwST36hnpEkVingJfPre6xmaWbyvlrS0lvLm1hOaAMzGjHz/88iS+dHamhmxEuhgFvpyxTQVHeD43n9c3FXO4ppGUvgncNDubRdOzmD5sgK58FemiFPhyWo7WNvLG5mJ+/8EBthdX0atnDy4en8bXZw/nC6MH6x42It2AAl/aFQg4a/cd4rn1+by2qZjmgDM6NYkfLJzEVdOzGKiplCLdigJf/kppVR0v5hXw3Pp8DlbW0Cchjq/PzmbRtExmZA/UkI1IN6XAFwCaA87b20t5fn0+7+wsI+Bw7qhB3H3ZOC6fPERXv4pEAQV+jDtS08DSDUU8tfYAu8qOkZqcyB0Xjua6nGGMTNENy0SiiQI/Brk7H+cf4em1B3ltUxH1TQHOyurHf90wnQVnDSFeF0aJRCUFfgw5WtPIUx8e4LVNxWwvriIpIY5rZg7lhlnZnJXVP9LliUgHCynwzew/gC8DDcAe4FZ3P9LGfvuBaqAZaHL3nFDalTOTX1nD87n5PP7+fqrrmpic2Y//+5Wz+Mr0LPom6ne+SKwI9X/7CuBed28ys38H7gX+qZ19L3b3ihDbk9NU39TMm1tKeG59Pu/vOYQZXDoxnbsvG8fEDN1nXiQWhRT47r681cu1wDWhlSOhyt1fydINRby2qYjDNY0MG9Sbuy8bx1dnDiVrQO9IlyciERTOv+dvA55rZ5sDy83MgV+7+5L2DmJmi4HFANnZ2WEsL3oFAs7avYd4dM0+/rKjjN4947h4Qio3zMpmzugUXQUrIsBpBL6ZvQ0MaWPTfe6+NLjPfUAT8HQ7h5nj7kVmlgasMLMd7r6qrR2DvwyWAOTk5Php9CFm1TY08/LHhTz23j52lx0jOTGeuy4azXfmjqFPgsbmReSzTpkK7n7pybab2TeBhcAl7t5mQLt7UfBzmZm9DMwC2gx8ObX8yhqeWnuAZ9fnc7S2kcmZ/fj5dVNZMCWDXj11gZSItC3UWTpX0PIm7YXuXtPOPklAD3evDi7PA+4Ppd1YVV5dz7+/uYOXPirAzLhi8hBumTOCnOG63YGInFqof/c/CCTSMkwDsNbd7zSzTOBRd18ApAMvB7fHA39w9zdDbDem5FfW8Ku/7OKVj4sAuHXOSG4/fySZehNWRM5AqLN0xrSzvghYEFzeC0wNpZ1YlV9Zw0Pv7uaF3AJ69DC+OnMot58/kjFpfSNdmoh0Q3pnrwv6pLSax9bs48W8AnqYcePsbP73RWMY0r9XpEsTkW5Mgd+F7Cqt5r5XtrBuXyUJ8T34+uxs7rxoNBn9NXQjIqFT4HcBq3eV88jKPby3+xCDkxK4d/4Ers0ZxiA9YEREwkiBH0F7y4/x65V7eS43n7TkRL518WiuPyebYYP6RLo0EYlCCvwIOF7fxE/f3MEf1h3EzLh1zgjunT+RhHjdllhEOo4CvxMFAs7DK/fw6Oq9HK5p5IZZw7j7svGkJidGujQRiQEK/E7QHHCWbijkmXUHWb//MHMnpHHXRaM5Z8SgSJcmIjFEgd/BDh6q4ZbH17G3/DiZ/Xvxb1dN4YZZw3RlrIh0OgV+B2lsDvDU2gM88NZOHPjF16axaFqmgl5EIkaBH2buzupdFfz0rR1sKaxiRvYAfnrN2YxJS450aSIS4xT4YXSkpoEfL9vBc7n5DEpK4KGvz2D+WUN0Vi8iXYICP0zW7avk7uc3UHC4lpvOzebe+RNJ0vNiRaQLUSKFqLqukZ8t/4Sn1h6gb694/njXecwcrtk3ItL1KPBDsGZXBf/0x00UH63lupxh/P08zakXka5Lgf85HK9v4kd/2srzuQWMSknihTu/wMzhAyNdlojISSnwz9DqXeV8/8VNFB+t46szhvKvV52lxwqKSLegwD8D+ZU13PFkHhn9e/HoN3K4ZGKaZuCISLehwD9NeQcO8/0XNxJnxpO3z9bjBUWk21Hgn0Ig4Dyyag+/WLGL1OREHrl5psJeRLolBf5J1DQ08Y8vbuL1TcVcND6VX3xtGgP66KEkItI9hXQDdjP7oZkVmtmG4MeCdva7wsx2mtluM7snlDY7S31TM9c8/AHLNhdz7/wJ/O6WcxT2ItKtheMM/z/d/YH2NppZHPDfwGVAAbDezF51921haLtDVNc1cvvjuWwrruLBG6ez8OzMSJckIhKyznjE0ixgt7vvdfcG4FlgUSe0+7kcq2/i+iVrWbe/kgeunaqwF5GoEY7A/7aZbTKzx8ysrauPsoD8Vq8LguvaZGaLzSzXzHLLy8vDUN7pO17fxJ1P5rGjpJpfXj+Na2YO7dT2RUQ60ikD38zeNrMtbXwsAh4GRgPTgGLgZ20doo113l577r7E3XPcPSc1NfX0ehEGNQ1N3PFkHmt2V/AvX5rIomnt/k4SEemWTjmG7+6Xns6BzOw3wGttbCoAhrV6PRQoOq3qOtGPXt3Gmt0V3DN/ArfOGRnpckREwi7UWToZrV5eBWxpY7f1wFgzG2lmCcD1wKuhtBtuz6w7yHO5+dx10WjuvHB0pMsREekQoc7S+amZTaNliGY/cAeAmWUCj7r7AndvMrNvA28BccBj7r41xHbDwt35yZs7+PXKvXxxbArfu3RcpEsSEekwIQW+u9/czvoiYEGr18uAZaG01RFWbCvl1yv3cuPsbO6/cjLxcZ0xaUlEJDJiNuEqjzfw4zd2MDIlSWEvIjEhJm+tUFZVx7W//oCSo3X87pZzFPYiEhNiMvB/+eddFB2p5YnbZnPe6MGRLkdEpFPE3KntU2sP8PSHB7kuZ5jCXkRiSkwFfsWxen71l13MyB7Aj66cHOlyREQ6VcwEvrvzjd+u40hNI397yViN24tIzImZMfzXNxezrbiKn1w9hYvGp0W6HBGRThcTp7nuzgNv7WRyZj/dEE1EYlZMBP7Woir2H6rh5nOHayhHRGJWTKTf65uLiethXD55SKRLERGJmKgPfHfn9U3FzBmTwsAkPaJQRGJX1Af+lsIqDlbW8KUpOrsXkdgW9YH/2uYi4nsY8yYp8EUktkV14Dc0BfjThiK+oOEcEZHoDvyNBUcoOlrHtZqKKSIS3YFfWlUHwLj05AhXIiISeVEd+GVV9QCkJidGuBIRkciL6sAvPFJLQnwPBvTuGelSREQiLqoDf2P+Ec7K7EePHhbpUkREIi5qA7+usZlNBUeZOXxgpEsREekSojbwNxUcpaE5wDkjBkW6FBGRLiGk2yOb2XPA+ODLAcARd5/Wxn77gWqgGWhy95xQ2j0dSzcUEt/DFPgiIkEhBb67f+3TZTP7GXD0JLtf7O4VobR3JjbkH9EFVyIirYRlSMfMDLgOeCYcxwuH0qp6sgb0inQZIiJdRrjG8L8IlLr7rna2O7DczPLMbPHJDmRmi80s18xyy8vLP1cxjc0BDh2vJzVZgS8i8qlTDumY2dtAW3ceu8/dlwaXb+DkZ/dz3L3IzNKAFWa2w91XtbWjuy8BlgDk5OT4qeprS+HhWtxh6MDen+fLRUSi0ikD390vPdl2M4sHrgZmnuQYRcHPZWb2MjALaDPww2FvxTEARqcmdVQTIiLdTjiGdC4Fdrh7QVsbzSzJzJI/XQbmAVvC0G679pYfB2BkSt+ObEZEpFsJR+BfzwnDOWaWaWbLgi/TgTVmthFYB7zu7m+God127a04zoA+PRmkGToiIv8jpGmZAO5+SxvrioAFweW9wNRQ2zkT5dX1DOmnN2xFRFqLyittq2ob6acbpomIfEZ0Bn5dE/16KfBFRFqLysCvrmukX++QR6tERKJKVAZ+VW2jzvBFRE4QlYE/d0IaZw/tH+kyRES6lKgc9/jF9dMjXYKISJcTlWf4IiLy1xT4IiIxQoEvIhIjFPgiIjFCgS8iEiMU+CIiMUKBLyISIxT4IiIxwtw/11MEO4WZlQMHPueXpwAVYSynO1CfY4P6HP1C6e9wd09ta0OXDvxQmFmuu+dEuo7OpD7HBvU5+nVUfzWkIyISIxT4IiIxIpoDf0mkC4gA9Tk2qM/Rr0P6G7Vj+CIi8lnRfIYvIiKtKPBFRGJE1AW+mV1hZjvNbLeZ3RPpesLFzIaZ2Ttmtt3MtprZ3wXXDzKzFWa2K/h5YKuvuTf4fdhpZpdHrvrQmFmcmX1sZq8FX0d1n81sgJm9aGY7gv/e58VAn78X/LneYmbPmFmvaOuzmT1mZmVmtqXVujPuo5nNNLPNwW3/ZWZ22kW4e9R8AHHAHmAUkABsBCZFuq4w9S0DmBFcTgY+ASYBPwXuCa6/B/j34PKkYP8TgZHB70tcpPvxOft+N/AH4LXg66juM/B74G+CywnAgGjuM5AF7AN6B18/D9wSbX0GLgBmAFtarTvjPgLrgPMAA94A5p9uDdF2hj8L2O3ue929AXgWWBThmsLC3Yvd/aPgcjWwnZb/KItoCQiCn78SXF4EPOvu9e6+D9hNy/enWzGzocCXgEdbrY7aPptZP1qC4bcA7t7g7keI4j4HxQO9zSwe6AMUEWV9dvdVQOUJq8+oj2aWAfRz9w+8Jf2faPU1pxRtgZ8F5Ld6XRBcF1XMbAQwHfgQSHf3Ymj5pQCkBXeLlu/FL4DvA4FW66K5z6OAcuB3wWGsR80siSjus7sXAg8AB4Fi4Ki7LyeK+9zKmfYxK7h84vrTEm2B39ZYVlTNOzWzvsAfge+6e9XJdm1jXbf6XpjZQqDM3fNO90vaWNet+kzLme4M4GF3nw4cp+VP/fZ0+z4Hx60X0TJ0kQkkmdlNJ/uSNtZ1qz6fhvb6GFLfoy3wC4BhrV4PpeVPw6hgZj1pCfun3f2l4OrS4J95BD+XBddHw/diDnClme2nZXhurpk9RXT3uQAocPcPg69fpOUXQDT3+VJgn7uXu3sj8BLwBaK7z5860z4WBJdPXH9aoi3w1wNjzWykmSUA1wOvRrimsAi+E/9bYLu7/7zVpleBbwaXvwksbbX+ejNLNLORwFha3uzpNtz9Xncf6u4jaPm3/Iu730R097kEyDez8cFVlwDbiOI+0zKUc66Z9Qn+nF9Cy3tU0dznT51RH4PDPtVmdm7we/WNVl9zapF+57oD3glfQMsMlj3AfZGuJ4z9Op+WP902ARuCHwuAwcCfgV3Bz4Nafc19we/DTs7gnfyu+AFcxP+fpRPVfQamAbnBf+tXgIEx0OcfATuALcCTtMxOiao+A8/Q8h5FIy1n6rd/nj4COcHv0x7gQYJ3TDidD91aQUQkRkTbkI6IiLRDgS8iEiMU+CIiMUKBLyISIxT4IiIxQoEvIhIjFPgiIjHi/wGDPuL+glf25QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.log(linalg.eigvalsh(cov)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "061618c7-9dd5-404c-a8da-b8f419be4a0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "points -= points.mean(0, keepdims=True)\n",
    "grads -= grads.mean(0, keepdims=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2417b107-ec82-4149-b758-ab72763e7bac",
   "metadata": {},
   "outputs": [],
   "source": [
    "span = np.concatenate([grads, points]).T\n",
    "subspace, svdvals, _ = linalg.svd(span, full_matrices=False)\n",
    "subspace = subspace.T\n",
    "\n",
    "# Remove zero values from subtracting mean\n",
    "subspace = subspace[:-2, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f73f83aa-1dcc-49bc-9c50-80235a0062fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "points_sub = points @ subspace.T\n",
    "grads_sub = grads @ subspace.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "5d556afa-3035-41bb-8fbf-ccdfec43c399",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_140260/1413332209.py:1: RuntimeWarning: invalid value encountered in log\n",
      "  plt.plot(np.log(linalg.eigvalsh(np.cov(points_sub.T))[11:]))\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f7d605ac130>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD6CAYAAABEUDf/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZR0lEQVR4nO3de3Bc93ne8e+7VywWAAEQ4P1OkbIoxVIkiLLjNpYlOZZbW2qTuqFd13LqjMatM4477jRWNW2mF814kk6mmWnThnXcqqltjepGkcaxY0tKXCeWZJmS7EjUlTeJICkSIHEhFsBe3/6xCxCmQJHUOcv9cfF8Znawe85izysKz8vDd3+LY+6OiIi0p0SrCxARkeZRkxcRaWNq8iIibUxNXkSkjanJi4i0MTV5EZE2FkuTN7OvmtkJM3thwbZ+M3vUzF5rfO2L41giInLhLI518mb2i8AU8L/c/ZrGtt8BTrn7l83sS0Cfu//W273OwMCAb9q0KXI9IiJLyTPPPDPq7oOL7UvFcQB3/4GZbTpr853AzY379wPfB962yW/atIk9e/bEUZKIyJJhZq+fa18zZ/Ir3f0YQOPriiYeS0REFtHyN17N7G4z22Nme0ZGRlpdjohIW2lmkz9uZqsBGl9PLPYkd9/t7kPuPjQ4uOhISURE3qFmNvlHgLsa9+8CHm7isUREZBFxLaH8BvAkcKWZDZvZZ4AvAx80s9eADzYei4jIJRTX6pqPn2PXrXG8voiIvDMtf+NVRESaJ5YzeRERuXilSo29Ryd49o1xVvV08HffvTr2Y6jJi4hcAuVqjQMjBV48NsHeI5M8d3ic549MUKrUALjj2jVq8iIioZstV3nt+BQHRqc4NDrNwdEp9o1M8erxqfmGnkkluGZND596z0au39jH9Rv6WLWsoyn1qMmLiER0ZHyGv3j5BH/58gme2D/KbLnezM1gzbIcWwbzfPoXNnHV6m52rF7GlsE86eSleUtUTV5E5AIUK1XenJhleGyGg6MFDo0WOHSywP6RAgdHCwBs6O9k140buGlzP1sGu9i4vJOOdLKldavJi8iS5e6MT5c5fnqW45NFTkzOcrJQYmy6xFihxKlCmdGpIkfGZxg5XfyZ782mEmxanmf7yi4+sXMDH3jXCrYO5jGzFv3XLE5NXkSWjGrN2Xt0gif2n+SJ/SfZc+gU06XqW56XSSXo78zQ25lmoCvLLVeuYE1vjjW9HaztzbFpIM+qng4SibAa+mLU5EXksuTuVGrObLnKbLnGbLnKG6emee34afaNTPHa8SlOFkrUak7Nnao744Uyp4sVALat6OJXrl/H5oE8K3s6WNmTZUV3BwPdGXLpZHBn5O+UmryItIy7c7pYYeR0kZlSlWKlRrFS/zo5U+bkVIlThVJ9hFIocaoxRhmbLjE+XaZSW/yiRz0dKbat7Gb7yi4SZiTMSCaMfDbJjZv6ee/W5azobs5qltCoyYvIO1ap1hibLjNdqlAoVpkuVZhuNOtytX4rVmoUihVOFeoNe2y6xOhUiROT9Tn4TPmt45KFEgb9+cz87YoVXfTlM/Tm0nRmkmRTSTrSCbLpJOt6c1yxsovBrmzbnIlHpSYvssSVKjXGpkvzTXjhbWy6RKFYnW/YpUqNQqnCyakSo1NFxqbLF3wcM+jrPNOsf25dL7d2Z1nZk2WwO0tnJkVHOkk2lSCTStDTkWZ5PsOyXPqymH2HSk1e5DJVrTmlSr3xFiv1ufTYdOnMrVDm9GyFqWKZqWKlcb/C1OyZ+5MzZ2bUZzODZbk0+UyKTCpBOmmkkwly6SRbB7u4aUs/y/NZlndl6Mqm6Myk6Mwk6cwk6UgnSSfPfE8+m2JZLk1SzfqSU5MXuURqNWdyttw4Qy4zMVOfK49P15twde4Nwlr9drrRhCdmykzOVigUK0wXK0yXq0wXq5SqtQs6bmcmSVc2RVdHiq5siu6OFMu7OunKpuv38xn68pm3fO3NpUldog/sSPOoyYtcIHdnsnEGXKnWKFedSq3G5EyFI+PTDJ+aYXhshmOTs8yWq40z7PpZ9uRMmbHpMtVzvFE4xwySjTcJuztS9HSk6cmlWZZLs643R65xppzLJMml6/PobCpBNp0gm0rS15mmtzNDX2eavs4M3R0pNeolTk1elpRqY8ndTLnKbLnKxEyZ0akSI6eLjE4VGZ8uU67W6k285hTLNU6cnuXo+AzHJmYXXVO90GB3ltXLOsilk3R3pBhoNOGeXPpnzpSXNZpwby5Nb2eafDZFKmF6s1BipyYvl43ZcrX+MfITBQ6MTM2POKruuMNMqcrJQolThSKnCiXGZ8pUq/X9c6OQcvXtz6QzjTlyqvE1k0ww2J1l+8pu3r99BWt6O+pnx4kEqQXz5nV9Odb25lr+EXaRs6nJS9O5OzPl6vz8eXzmras46iMQb6zicErV+odbio2z7kKxyrGJGRZOOzrSifr6ZzPMoCOdpD+fYXlXhp/r66U3lyadTJBMQCJRf97ccrtcJklHqn62PdidZaCrvsIjn1UkpL3oJ1pi9fKbkzzw9GFOnJ7lzYn6OuiRqeL8r1hdzLJcmq5sfQVHKlE/O06nEnSkEvTlM6xO1efQ6/s72bqii62DeTYP5OnM6MdX5HyUEonVV//6IN98Znj+o+I7N/ezojs7v1qjt7P+RuLcWum+zswl+5WrIkuRmrzEarZcY0N/J49/8eZWlyIi6ELeErNSpUYmpR8rkVAojRKrUlVNXiQkSqPEqlSpkdGMXSQYSqPESuMakbAojRKrYrVGJqUPBImEQk1eYqVxjUhYlEaJValSJatxjUgwlEaJlVbXiIRFaZRYaVwjEhalUWKl1TUiYVEaJVZq8iJhURolVprJi4RFaZTY1Gr1i3JoJi8SDqVRYjN3YWmdyYuEQ2mU2Mw1ea2TFwmH0iixmbv6k87kRcKhNEps5pu8ZvIiwWh6Gs3sdjN7xcz2mdmXmn08aR2dyYuEp6lpNLMk8F+ADwM7gI+b2Y5mHlNaR2+8ioSn2WncCexz9wPuXgIeAO5s8jGlRTSuEQlPs9O4Fji84PFwY5u0oaLGNSLBaXYabZFt/jNPMLvbzPaY2Z6RkZEmlyPNpJm8SHiancZhYP2Cx+uAowuf4O673X3I3YcGBwebXI40k9bJi4Sn2Wn8MbDNzDabWQbYBTzS5GNKi5yZyevyfyKhSDXzxd29Yma/AXwXSAJfdfe9zTymtI7GNSLhaWqTB3D3bwPfbvZxpPVK1SqgJi8SEqVRYqMzeZHwKI0SG62TFwmP0iix0Tp5kfAojRIbLaEUCY/SKLHRuEYkPEqjxKZUqZFKGInEYh90FpFWUJOX2JQquoi3SGiUSIlNqaomLxIaJVJiU6rUNI8XCYwSKbHRuEYkPEqkxEbjGpHwKJESG41rRMKjREpsStWaPgglEhglUmKjmbxIeJRIiY2avEh4lEiJTamqmbxIaJRIiU2pUiOtJi8SFCVSYqNxjUh4lEiJTVFNXiQ4SqTERksoRcKjREps9GEokfAokRIbzeRFwqNESmz0u2tEwqNESiyqNadaczLJZKtLEZEF1OQlFvPXd9WZvEhQlEiJhZq8SJiUSIlFsVoF1ORFQqNESizmzuSzWkIpEhQlUmKhcY1ImJRIiUWpqiYvEiIlUmIxfyavcY1IUJRIiYXGNSJhUiIlFmryImFSIiUWRc3kRYKkREosNJMXCZMSKbGYXyevM3mRoCiREgvN5EXCpERKLLROXiRMkRJpZh8zs71mVjOzobP23WNm+8zsFTP7ULQyJXSayYuEKRXx+18Afhn4w4UbzWwHsAu4GlgDPGZm2929GvF4EiiNa0TCFCmR7v6Su7+yyK47gQfcvejuB4F9wM4ox5KwaVwjEqZmJXItcHjB4+HGNmlTRY1rRIJ03nGNmT0GrFpk173u/vC5vm2RbX6O178buBtgw4YN5ytHAlWq1MgkE5gt9r9eRFrlvE3e3W97B687DKxf8HgdcPQcr78b2A0wNDS06F8EEr5SRRfxFglRs1L5CLDLzLJmthnYBjzdpGNJAErVqpq8SICiLqH8+2Y2DLwX+DMz+y6Au+8FHgReBP4c+JxW1rS3uXGNiIQl0hJKd38IeOgc++4D7ovy+nL50LhGJExKpcSiVFWTFwmRUimx0LhGJExKpcSiqHGNSJCUSomFZvIiYVIqJRalak2/S14kQEqlxEIzeZEwKZUSC41rRMKkVEostIRSJExKpcRC4xqRMCmVEguNa0TCpFRKLNTkRcKkVEosiprJiwRJqZTI3J1SpUZWM3mR4CiVElm5Wr/Wi87kRcKjVEpkuoi3SLiUSomspIt4iwRLqZTI5pt8KtniSkTkbGryEtmZJq8fJ5HQKJUSWalav3yvmrxIeJRKiayombxIsJRKiWxuXKPfJy8SHqVSItNMXiRcSqVEpnXyIuFSKiUyrZMXCZdSKZFpXCMSLqVSItO4RiRcSqVEpiWUIuFSKiUyLaEUCZdSKZFpJi8SLqVSItNMXiRcSqVEpiWUIuFSKiWyUqVGwiClJi8SHKVSIivpIt4iwVIyJbJSpaZRjUiglEyJrH4mr6tCiYRITV4iq5/JW6vLEJFFqMlLZKWKZvIioVIyJTI1eZFwRUqmmf2umb1sZn9jZg+ZWe+CffeY2T4ze8XMPhS5UgmWVteIhCtqMh8FrnH3dwOvAvcAmNkOYBdwNXA78Admpnfm2pRW14iEK1Iy3f177l5pPHwKWNe4fyfwgLsX3f0gsA/YGeVYEi6Na0TCFWcy/wnwncb9tcDhBfuGG9vewszuNrM9ZrZnZGQkxnLkUilqCaVIsFLne4KZPQasWmTXve7+cOM59wIV4Gtz37bI832x13f33cBugKGhoUWfI2HTuEYkXOdt8u5+29vtN7O7gI8At7r7XJMeBtYveNo64Og7LVLCVqpU9bvkRQIVdXXN7cBvAXe4+/SCXY8Au8wsa2abgW3A01GOJeHS6hqRcJ33TP48/jOQBR41M4Cn3P2z7r7XzB4EXqQ+xvmcu1cjHksCpXGNSLgiNXl3v+Jt9t0H3Bfl9eXyoNU1IuFSMiUyNXmRcCmZEplm8iLhUjIlklrNKVddM3mRQCmZEoku4i0SNiVTIplr8lonLxImJVMiKVV0Ji8SMiVTIplv8prJiwRJyZRIdCYvEjYlUyLRG68iYVMyJRKNa0TCpmRKJEWNa0SCpmRKJJrJi4RNyZRItE5eJGxKpkRyZiavy/+JhEhNXiLRuEYkbEqmRFKq1q8FoyYvEiYlUyLRmbxI2JRMiUTr5EXCpmRKJFonLxI2JVMi0RJKkbApmRKJxjUiYVMyJZJSpUYqYSQS1upSRGQRavISSamii3iLhEzplEhKVTV5kZApnRJJqVLTPF4kYEqnRKJxjUjYlE6JpKhxjUjQlE6JROMakbApnRJJqVLTB6FEAqZ0SiSayYuETemUSLSEUiRsSqdEopm8SNiUTolE4xqRsCmdEkl9XKPru4qESk1e3rGZUpVCsaJxjUjAUq0uQC4vb07M8vjLx3n8pRP8cN8oxUqNlT3ZVpclIuegJi9va3y6xFMHTvLE/vpt34kpANb35/jETRu47aqV3LS5v8VVisi5RGryZvbvgTuBGnAC+LS7H23suwf4DFAFPu/u341Yq1wCEzNlnj54iqcOnOSpAyd58dgk7pBLJ9m5uZ9/cMM6bnnXCrat6MJMv0NeJHRRz+R/193/NYCZfR74N8BnzWwHsAu4GlgDPGZm2929GvF4EoPZcpUf7hvl/706wvHJWSZnKkzOlpmYKXNkfAb3+jVbb9jQxxdu3c77rljOu9f1ahWNyGUoUpN398kFD/OAN+7fCTzg7kXgoJntA3YCT0Y5nly8Ws0ZnSoyPD7Dq2+e5vGXT/BXr40wW66RzyRZ19dJTy7Fqp4OrlzZzfr+Tt67dTnXre+lI61VMyKXu8gzeTO7D/gUMAF8oLF5LfDUgqcNN7ZJE1RrzvHJWQ6NFtg/WuDAyBQHRgq8cWqaI2Mz8xfbBljbm+MfDq2vz9K39JPV8keRtnbeJm9mjwGrFtl1r7s/7O73Avc2ZvC/Afw2sNiw1hfZhpndDdwNsGHDhgute8koV2v85PA4LxyZoFCsMFWsMl2qcHq2wtHxGY6Mz/DmxCyV2pk/3lw6yeaBPDtW9/BLO1ayti/H2t4cG5fn2TqY1yxdZAk5b5N399su8LW+DvwZ9SY/DKxfsG8dcPQcr78b2A0wNDS06F8ES4G7M1WsMFYoMzZd4m+OTPCDV0d4cv9JpoqV+edlkgk6s0m6svURyw0b+1jbm2NNb45Ny/NsGcyzqqdDF9YWESD66ppt7v5a4+EdwMuN+48AXzez36P+xus24Okox7rczZarjE4VOTlVqo9WThY4OFrgwEiB109Oc7JQpFz92b/j1vbm+Oi1a3j/9gFu2NjPslxab36KyEWJOpP/spldSX0J5evAZwHcfa+ZPQi8CFSAzy2llTVHx2d45vUxnn1jjGdfH2PfiSkKpbf+5/fnM2weyPO+KwZY0ZOlvzNDb2d6fvvmAY1WRCSaqKtrfuVt9t0H3Bfl9UMzW67y08PjHBmfmZ+HH5uY5fRshanZClPFCqdny0zO1scruXSSa9cv42ND6xnszjLQlWF5Pstgd5aNyzvp7cy0+L9IRNqdPvF6AQ6OFvjaU6/zzWeHGZ8uz28f6MqwalkHy3Jpluc76epI0ZVNsXkgz9DGft61upu0fq+LiLSQmvwi3J3hsRl+dPAUf/rcEf563yiphPGhq1fxy9evZfNAnjW9Oa0jF5Hgqck3FIoV/vQnR3hy/0l+fOgUxyeLAKxZ1sEXP7idX71xPSt6OlpcpYjIxVnyTX6sUOL+Jw/xP584xPh0mdXLOrhp83Ju3NTHjZv72b6iW8sRReSytWSb/ORsmd9/7DW+8fQbTJeq3HbVSv7pzVu5YWNfq0sTEYnNkmzyY4USn/rq07x4bJI7rl3DZ9+/lStXdbe6LBGR2C25Jj86VeSTX/kRB0YK7P7HN3DrVStbXZKISNMsqSZ/fHKWT/z3pzgyPsMffXqIv71tsNUliYg01ZJp8vtOTPGZ+3/M6Oki9//aTm7asrzVJYmINF1bN/n9I1N85/ljfOeFN9l7dJLujhR//Os3cf0GvbkqIktD2zX5as359vPH+MMf7OeFI/Vrmly/oZd7/85VfPTaNaxaprXuIrJ0tE2TL1drPPTcEf7b9/dzYLTA1sE8v/3RHdx+zSpWL8u1ujwRkZZoiyb/08Pj/LOvPcuR8Rl2rO7hD/7R9dx+9Sp9iElElry2aPJzF8v4D3/vGm6+clC/nldEpKEtmvyyzjR//JmbWl2GiEhw9HtwRUTamJq8iEgbU5MXEWljavIiIm1MTV5EpI2pyYuItDE1eRGRNqYmLyLSxszdW13DPDMbAV5vdR1nGQBGW13EBVKtzXG51Hq51AmqNW4b3X3RC2QE1eRDZGZ73H2o1XVcCNXaHJdLrZdLnaBaLyWNa0RE2piavIhIG1OTP7/drS7gIqjW5rhcar1c6gTVesloJi8i0sZ0Ji8i0sbU5BvMbL2Z/aWZvWRme83sNxvb+83sUTN7rfG15VcBN7MOM3vazH7aqPXfhlrrHDNLmtlzZvatxuMgazWzQ2b2vJn9xMz2NLaFWmuvmX3TzF5u/Ny+N8RazezKxp/n3G3SzL4QaK3/vJGpF8zsG42sBVfnxVCTP6MCfNHdrwLeA3zOzHYAXwIed/dtwOONx61WBG5x92uB64Dbzew9hFnrnN8EXlrwOORaP+Du1y1YNhdqrb8P/Lm7vwu4lvqfb3C1uvsrjT/P64AbgGngIQKr1czWAp8Hhtz9GiAJ7CKwOi+au+u2yA14GPgg8AqwurFtNfBKq2s7q85O4FngplBrBdZRD8ctwLca20Kt9RAwcNa24GoFeoCDNN5XC7nWs+r7JeCHIdYKrAUOA/3Ur5r3rUa9QdV5sTedyS/CzDYBPw/8CFjp7scAGl9XtLC0eY3xx0+AE8Cj7h5srcB/Av4lUFuwLdRaHfiemT1jZnc3toVY6xZgBPgfjTHYV8wsT5i1LrQL+EbjflC1uvsR4D8CbwDHgAl3/x6B1Xmx1OTPYmZdwP8FvuDuk62u51zcver1f/6uA3aa2TUtLmlRZvYR4IS7P9PqWi7Q+9z9euDD1Ed2v9jqgs4hBVwP/Fd3/3mgQOBjBDPLAHcA/6fVtSymMWu/E9gMrAHyZvbJ1lYVnZr8AmaWpt7gv+buf9LYfNzMVjf2r6Z+5hwMdx8Hvg/cTpi1vg+4w8wOAQ8At5jZ/ybMWnH3o42vJ6jPjXcSZq3DwHDjX3AA36Te9EOsdc6HgWfd/XjjcWi13gYcdPcRdy8DfwL8AuHVeVHU5BvMzIA/Al5y999bsOsR4K7G/buoz+pbyswGzay3cT9H/YfzZQKs1d3vcfd17r6J+j/V/8LdP0mAtZpZ3sy65+5Tn8e+QIC1uvubwGEzu7Kx6VbgRQKsdYGPc2ZUA+HV+gbwHjPrbPSDW6m/mR1anRdFH4ZqMLO/BfwV8DxnZsf/ivpc/kFgA/Ufgo+5+6mWFNlgZu8G7qf+7n8CeNDd/52ZLSewWhcys5uBf+HuHwmxVjPbQv3sHerjkK+7+30h1gpgZtcBXwEywAHg12j8PBBerZ3U39Tc4u4TjW3B/bk2liP/KvXVds8Bvw50EVidF0NNXkSkjWlcIyLSxtTkRUTamJq8iEgbU5MXEWljavIiIm1MTV5EpI2pyYuItDE1eRGRNvb/AcOsiLJW/VbxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.log(linalg.eigvalsh(np.cov(points_sub.T))[11:]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "ceadacb7-6493-41fd-b16a-a04dd8fb3243",
   "metadata": {},
   "outputs": [],
   "source": [
    "Up, sp, _ = linalg.svd(points_sub.T, full_matrices=False)\n",
    "Ug, sg, _ = linalg.svd(grads_sub.T, full_matrices=False)\n",
    "\n",
    "# remove one value due to zero eigenvalue from mean subtraction\n",
    "Up = Up[:, :-1]\n",
    "sp = sp[:-1]\n",
    "Ug = Ug[:, :-1]\n",
    "sg = sg[:-1]\n",
    "\n",
    "sp /= np.sqrt(k)\n",
    "sg /= np.sqrt(k)\n",
    "sp = sp ** 2\n",
    "sg = sg ** 2\n",
    "\n",
    "others_p = np.exp(np.log(sp.min()) - np.log(sg.min()))\n",
    "others_g = 1 / others_p\n",
    "\n",
    "#others_p = others_g = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "2b94688a-dea1-4279-9928-ede55c2f552c",
   "metadata": {},
   "outputs": [],
   "source": [
    "epsilon = 1e-6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "d1ca02e2-f19d-45fb-84bb-a81fe116c40a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(98, 49)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Up.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "48bc4a7e-fcbd-4784-be5f-d731d2d93a0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "Cp = (\n",
    "    Up @ (np.diag(sp) - epsilon * others_p * np.eye(k - 1)) @ Up.T\n",
    "    + epsilon * others_p * np.eye(2 * k - 2)\n",
    ")\n",
    "Cg = (\n",
    "    Ug @ (np.diag(sg) - epsilon * others_g * np.eye(k - 1)) @ Ug.T\n",
    "    + epsilon * others_g * np.eye(2 * k - 2)\n",
    ")\n",
    "\n",
    "Cg_sqrt = (\n",
    "    Ug @ (np.diag(np.sqrt(sg)) - np.sqrt(epsilon * others_g) * np.eye(k - 1)) @ Ug.T\n",
    "    + np.sqrt(epsilon * others_g) * np.eye(2 * k - 2)\n",
    ")\n",
    "Cg_invsqrt = (\n",
    "    Ug @ (np.diag(1 / np.sqrt(sg)) - 1 / np.sqrt(epsilon * others_g) * np.eye(k - 1)) @ Ug.T\n",
    "    + 1 / np.sqrt(epsilon * others_g) * np.eye(2 * k - 2)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "9993a64d-a809-4726-a3cf-3d6d70942805",
   "metadata": {},
   "outputs": [],
   "source": [
    "#assert np.allclose(Cg_sqrt, linalg.sqrtm(Cg))\n",
    "#assert np.allclose(Cg_invsqrt, linalg.inv(linalg.sqrtm(Cg)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "c0ed684e-306c-460f-8c32-612b201b2132",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean = Cg_invsqrt @ linalg.sqrtm(Cg_sqrt @ Cp @ Cg_sqrt) @ Cg_invsqrt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "02f90fe0-6428-4cd1-86b1-021b1bfcfb72",
   "metadata": {},
   "outputs": [],
   "source": [
    "#linalg.eigvalsh(mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "3304b761-1053-422a-b397-da9d1854407d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f7d60524550>]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAApK0lEQVR4nO3deXxd5X3n8c+jfd9lrZblfQNvCGOzJSFAwGwNEzIQsk9L1rbpJG3SMJNpk05enWnaTAgJhEAaaLY2Cy1JSUIgCzHYgG0M3lckW9a+S/dKd/3NH+dijCMby7rSle79vl8vvXTvOUf3PI8tPd9znvOc5zgzQ0REUk9aogsgIiKJoQAQEUlRCgARkRSlABARSVEKABGRFJWR6AKcTUVFhTU2Nia6GCIis8b27dt7zKzyXLad0QHQ2NjItm3bEl0MEZFZwznXcq7bqgtIRCRFKQBERFKUAkBEJEUpAEREUpQCQEQkRSkARERSlAJARCRFKQBEzlGXv4v/8/z/IRQNJbooInGhABA5B8+eeJbbfnobPz70Yw72HUx0cUTiQgEgchbhaJh7dtzDh5/8MGU5Zfzghh+wsmJlooslEhczeioIkUTa07uHz2/5PHt793Lr4lv5zPrPkJuRm+hiicSNAkDkNL6Qj3tfvJfv7f8eZTllfOlNX+JtjW9LdLFE4k4BIBIzEhzh+/u/zyN7H2EwMMg7l76TP1v3ZxRlFSW6aCJTQgEgAvzbgX/jnhfvYTAwyJX1V/LR1R9VX78kPQWApLxvvvxN7nnxHjbUbOAT6z6hhl9ShgJAUtp9L93H13d+nRsX3MgXLvsCGWn6k5DUMelhoM65pc65nad8DTnnPnHaNm92zg2ess3nJrtfkcnwhXz8wwv/wNd3fp2bF97M3132d2r8JeVM+jfezA4AawCcc+nACeDRcTb9vZndONn9iUzGYGCQ7+3/Ht/Z+x2GgkPctuQ2/seG/0Ga0y0xknrifcjzVuCImZ3zI8lEpsOe3j386OCPePzo4/jDft4y9y3cteouLqi4INFFE0mYeAfA7cD3z7Buo3PuJaAN+JSZ7RlvI+fcXcBdAA0NDXEunqSK0fAoL3a9yPPtz7P5xGYO9B8gNyOX6xqv487ld7K0bGmiiyiScM7M4vNBzmXhNe4rzazztHVFQNTMRpxzm4CvmNniN/rMpqYm00Ph5VwNB4f57fHf8svmX/Js27OEoiEyXAYXVl7IpvmbuGHBDRRmFSa6mCJTyjm33cyazmXbeJ4BXA/sOL3xBzCzoVNeP+6c+7pzrsLMeuK4f0lRZsa393ybe1+8l2A0SHV+Nbcvu53Lai9j7Zy15GXmJbqIIjNSPAPgDs7Q/eOcqwY6zcycc+vxRh/1xnHfkqLC0TBffO6L/PDgD7lq7lV84IIPsKpylS7qipyDuASAcy4PuAb40CnLPgxgZvcD7wA+4pwLA6PA7RavvidJKZFohL6xPsbCY4xGRvny9i+z+cRm/vjCP+ZP1/6pGn6RCYhLAJiZHyg/bdn9p7y+F7g3HvuS1BS1KL945Rfc8+I9nBg5cXJ5ukvncxs/x21Lbktg6URmJ935IjNS1KL0jfXR4evg+PBxHt7zMHt697C0dCmfWf8ZCjILyM7IZkHxApaULkl0cUVmJQWAzAhmxqGBQzzX/hxb27eyvXM7vpDv5Prq/Gq+ePkXuWHBDermEYkTBYBMu3A0TM9oD13+LlqGWtjavpUtbVvoHu0GYF7RPDbN38Ti0sVU51VTnV/NwpKFZKVnJbjkIslFASDTom+sjydbnuSJ5ifY1rmNiEVOrivOLmZjzUYurb2UDTUbqCmoSWBJRVKHAkCmxFh4jP19+9nSvoWtbVt5qfslIhahsaiR9654L/WF9VTlVVGdX82ikkWkp6UnusgiKUcBIOdtLDzGCx0vcLD/IM1DzbQMtdDp66Q/0M9oeBQAh2NF+Qo+eMEHeVvj21hSugTnXIJLLiKgAJBzELUox4aO0T3aTf9YP13+Lra0b+H59ucZi4wBUJ5TTmNxI+uq1lGaU0pZThkNhQ2sr15PSU5JYisgIuNSAMi4On2dbGnfwpa2LWxt30rfWN/r1tcV1HHr4lu5sv5KVlWu0hw7IrOQAkCIRCO0DLdwoO8AO7t2srV9K0cHjwLekf3G2o1cUn0JNQU1lGaXUp5bTnlOubpyRGY5BUCKCUfD7Ovdx+7e3RzoO8DB/oMc6j90sisnJz2Hi6ou4tbFt7KhZgOLSxdr3L1IklIAJKnR8CjHh4/T6eukZ7SHDn8Hu7p3saNrx8kbrIqzi1laupR3LHkHy8uXs7R0KQuKF5CZnpng0ovIdFAAzGJmRjgapj/Qz+6e3ezq2cXunt00DzXT4ev4g+0bixq5Yf4NXFxzMWsq11CVV6VuHJEZps8XZH/7EJcuqpjyfSkAZonBwCDbO7fzYteL7OjawdGBo4yGR193Q1WGy2Bx6WIurrqYeUXzmFc0j+r8airzKqnIrSA7PTuBNRCRN7KrdZAPf2c7/mCYzZ++ivzsqW2iFQAzTCgaosvfRbe/m57RHnb37GZr+1b29u7FMDLTMrmg4gJuXngz+Zn55GTkUJBZwIryFSwrW0ZORk6iqyAi5+FH21v57KO7qMjP4tsfWD/ljT8oAGaMY0PH+NcD/8qjhx9lODh8cnmGy2BV5So+svojXFJzCSsrVupIXiSJmBl///P9fOPpo1y6sJyv3rGW8oLp+RtXACTAwf6DPNnyJN2j3QwFhuge7ebFrhfJcBlcPe9qNtZupCK3gsrcShqKGsjPzE90kUVkinzz90f5xtNHefeGBv7mppVkpE/fqDsFwBSKWpQTIycYDAwyGBikeaiZnx75KXt695Dm0ijJLqEoq4ji7GI+uvqjvGPJO6jMq0x0sUVkmjz2UhtffHw/N6yq4fM3X0Ba2vQOylAAxJGZ0TfWx+6e3fzm+G/47fHf0jv2+kcfLyldwqcv/jQ3LLiB0pzSxBRURBIqGjV+d7CbT/3bS6yfX8Y/3rZ62ht/UABMmD/kZ1vnNra0baHT38loeJTR8Cj9Y/20jbSdvKEqPzOfy+suZ0PNBipyKyjOLqYip4L6wnoNvRRJQaPBCI9saWbz4R52Hh9geCzMojkFfPM9TeRkJmY2XAXAWezt3ct3932X/rF+gpEgvpCPA/0HCEVD5KTnUFdQR25GLjkZOSwoXsDldZdTW1DLwpKFXDTnIt1QJSKYGY/v6uB//+de2gbHWFZdyE2ra1k7t4RrV1RTnJe4diIuAeCcawaGgQgQNrOm09Y74CvAJsAPvN/MdsRj31PhlcFXuPfFe3mi5QkKMwuZWzSX7PRsCrIKuHP5nWys3chFVRdpNI6InNWBjmH+9qd7ePZIL8trivjKHWu5uLEs0cU6KZ5nAG8xs54zrLseWBz7ugS4L/Y9IYKRIDu7drK/bz8DgQGGgkP0jfXRNtJGu6+dvrE+cjNy+dCqD/G+le/TTJciMiED/iBf/tVBvvPcMQqyM/jCLSt51yXzSE9AP//ZTFcX0C3AI2ZmwFbnXIlzrsbM2qdj56FoiAN9B9jeuZ3nO57nhY4XTj6wJM2lUZRVRGlOKTX5NSwrW0ZDUQO3LLyF8tzy6SieiCSBzqExnjncwzOHe3lyXyfDYyHuvGQe//2aJZTmz8znWccrAAx4wjlnwDfM7IHT1tcBx0953xpb9gcB4Jy7C7gLoKGh4bwLNBIc4aljT/GL5l+wvXP7yQa/odBr3C+tvZQ1c9ZQnF2s2S5FZEKGx0I8sqWF3+zvonskQPdwAH/Qm5alJC+TKxZX8LG3LGJ5TVGCS3p28QqAy8yszTk3B/iVc26/mT19yvrxzntsvA+KhccDAE1NTeNuczb+kJ+7N9/N061PE4wGqSuo4+2L3s66qnWsnbOWOXlzJvqRIiKA17Xz8LMtPLT5KENjYdY2lLC6voTKwmxqinPYsKCcFTVFCRnSeT7iEgBm1hb73uWcexRYD5waAK3A3FPe1wNt8dj36XIzchkODnPb0tu4fv71rKpYpWGXInLehsdC/Hp/Fz99qY3fHewmFDGuWVHFn121mAvrixNdvEmZdAA45/KBNDMbjr2+Fvj8aZs9BnzcOfcDvIu/g1PV/++c48G3PTgVHy0iSc7MONLtY/OhbnYeH2DXiUGO9vgwg+qiHN63sZH/clH9jO/aOVfxOAOoAh6NHWVnAN8zs1845z4MYGb3A4/jDQE9jDcM9ANx2K+IyKQEw1EOdg7zUusAO48N8MzhHtoGvZs5q4qyubCumJtW13Lpwgqa5pXOmq6dczXpADCzo8DqcZbff8prAz422X2JiMTDy60DPLKlhZ+93MZYKAp4F283zC/n41dVcsXiCuaW5SW4lFNPdwKLSNLpGByjuddHKBIlFIkyPBamfXCM9oFRdh4f4KXWQfKy0nn72jouXVjB6voS5pblptz1QgWAiMxaZkafL0hLn5/mHh87jvXz7JFejnb7xt2+MCeDxvJ8/uamFdx6UT1FOak9XYsCQERmvHAkyq4Tgzx7pJcXmvvoHAow4A/S7w+e7MIByMtKZ/38Mu64uIEVtUVkZ6SRmZ5GXlY61cU5FKZ4g386BYCIJFwoEqWl18exPj89w0F6fAG6hgK09o/S2u/nWJ//5I1WS6sKmVuWxwW1RZTkZVJTnMu88jzmlefTUJZHVoZu7DxXCgARmTZmRs9IkP0dQ+w6McjuE4Ps7xjmWK+fcPT1930WZGdQV5JLfWkuGxaU09RYyoYF5VRM0+MSU4ECQESmzIA/yJYjvTxzpIe9bUMc6fYxOBo6uX5uWS7Lq4u4/oJqFlYW0FiRT2VBNhUF2eRmJWaO/FSiABCRCRsJhDnQMcz+jiFe6fYxNBbCF4gwHAgTCEUIhKP4AmEOd49g5h3Nr6gt4oZVNSyqLGBJVSEX1BVRkjczJ0lLFQoAETmr1n4/W470enfFdvs42j1y8mYpgNzMdEryMsnPziA/O4OcjDQKczKoKMhm04U1XLG4gtVzS8icxoedy7lRAIjISYFwhAMdw+w87t0Z+3xzH6393ky6hdkZLKjM55IF5SyoyGdZTRHLqgupL0298fPJQgEgkgIG/EEOdo7Q7w8yOBpiaDTEYOxrwB+iY3CM4/1+OobGsNi12IqCbC6aV8IfXz6fjQsrWDynIOmmQkh1CgCRJGBmdA8H2N8xTEufn96RAH2+IG0DY+xrH+LEwOgf/IxzUJSTSXFuJtXFOVy6sIL60lyWVBWypqGE2uIcHdknOQWAyCwxEgjT2u+nY3CMjsEx2gZGOR4bJ3+k20efL/i67YtyMphTlMO6eaW8d+M8ltUUUZ6fRXFuJkW5mRRmZ+iIPsUpAERmkNFghMNdI3SPjNE7EqR7JMC+9mF2nxjklZ7XT2+Q5qCm2Bsnf+2KKpZWF7K0upCFlQWU5Wfpoqu8IQWAyDTzB8Mc6Yrd9ToSoHckQOvAKHtODHGoa5jT7oeitjiHC+qKuXVtHfMr86kuyqEq9qW7XmUyFAAicWRmdA4FONI94n11jdA1HGAkEMYXCNM5FPiD/vg0511wXVFbxLUrq1hZW0RVUQ7l+dmUFWRRkK0/U5ka+s0SmaSWXh/3/+4oe9oGOdI1gi82Zw14QyerinMoyM6gIDuDpsY8bq+cy6I5sbteC7MpzcsiXX3xkgAKAJHzNOAP8tVfH+aRLc1kpKXR1FjKbU1zWViZz8I5BSyqLKCyMFsjaWTGUgCInIc+X5Br/ul39PuDvLNpLn9xzRKqinISXSyRCVEAiJyHQ53D9PqCfPWOtdy0ujbRxRE5LxpCIHIeXp2bvr40N8ElETl/kw4A59xc59xvnHP7nHN7nHN/Ps42b3bODTrndsa+PjfZ/Yokki8YBiBfI3RkFovHb28Y+KSZ7XDOFQLbnXO/MrO9p233ezO7MQ77E0k4f8A7A8jTnPUyi036DMDM2s1sR+z1MLAPqJvs54rMZCfPALJ0BiCzV1yvATjnGoG1wHPjrN7onHvJOfdz59zKs3zGXc65bc65bd3d3fEsnkjc+AJeAORl6wxAZq+4BYBzrgD4MfAJMxs6bfUOYJ6ZrQa+Cvz7mT7HzB4wsyYza6qsrIxX8UTiyheMkJnuyM5QAMjsFZcAcM5l4jX+3zWzn5y+3syGzGwk9vpxINM5VxGPfYskgj8QJk/dPzLLxWMUkAMeAvaZ2T+dYZvq2HY459bH9ts72X2LJIovGCFfF4BllovHIcxlwHuAXc65nbFlnwUaAMzsfuAdwEecc2FgFLjdzGyczxKZFfzBMHkaAiqz3KR/g81sM3DWyU7M7F7g3snuS2Sm8AV0BiCzn+4EFjkP/qCuAcjspwAQOQ++QIR8DQGVWU4BIHIedAYgyUABIHIefEGdAcjspwAQOQ+6D0CSgQJAZIKiUYudASgAZHZTAIhM0GjImwlUw0BltlMAiEzQqzOB6kYwme0UACIT9OqzAHQGILOdAkBkgk6eAegisMxyCgCRCXr1ecAaBiqznQJAZIJOPgxGZwAyyykARCZIZwCSLBQAIhP06hmAngcss50CQGSCXj0DyNMoIJnlFAAiEzTy6hmA7gOQWU4BIDJB/mCY9DRHdob+fGR202+wyAT5AhHystKJPeZaZNZSAIhMkD8Y1gVgSQpxCQDn3HXOuQPOucPOuc+Ms9455+6JrX/ZObcuHvsVSQRfMEKehoBKEph0ADjn0oGvAdcDK4A7nHMrTtvsemBx7Osu4L7J7lckUfwBnQFIcojHGcB64LCZHTWzIPAD4JbTtrkFeMQ8W4ES51xNHPYtMu18wYiGgEpSiEcA1AHHT3nfGls20W0AcM7d5Zzb5pzb1t3dHYfiicSXPxjWEFBJCvEIgPGGQth5bOMtNHvAzJrMrKmysnLShROJN39AZwCSHOIRAK3A3FPe1wNt57GNyKzg0yggSRLxCIAXgMXOufnOuSzgduCx07Z5DHhvbDTQBmDQzNrjsG+RaecLaBSQJIdJH8aYWdg593Hgl0A68C0z2+Oc+3Bs/f3A48Am4DDgBz4w2f2KJIKZ4QuGKdA1AEkCcfktNrPH8Rr5U5fdf8prAz4Wj32JJNJYKIqZngUgyUF3AotMwKuPg9SzACQZKABEJuDVB8LrDECSgQJAZAJOngFoGKgkAQWAyAT4YwGQp4vAkgQUACIT4It1AekMQJKBAkBkAk6eAegagCQBBYDIBJw8A9AoIEkCCgCRCfDpDECSiAJAZAJ0BiDJRAEgMgH+YBjnIDdTASCznwJAZAJ8gQj5WRl6ILwkBQWAyAT4g2E9C0CShgJAZAJ8wYieBiZJQwEgMgH+gM4AJHkoAEQmQE8Dk2SiABCZAH9QTwOT5KEAEJkAX0BnAJI8FAAiE+ALRHQNQJKGAkBkAnzBsEYBSdJQAIicIzPDH4xoGghJGpM6lHHO/QNwExAEjgAfMLOBcbZrBoaBCBA2s6bJ7FckEQLhKJGoaSI4SRqTPQP4FXCBma0CDgJ/fZZt32Jma9T4y2zlD+phMJJcJhUAZvaEmYVjb7cC9ZMvksjM5AvocZCSXOJ5DeCDwM/PsM6AJ5xz251zd53tQ5xzdznntjnntnV3d8exeCKT89oZgAJAksMb/iY7554EqsdZdbeZ/Udsm7uBMPDdM3zMZWbW5pybA/zKObffzJ4eb0MzewB4AKCpqcnOoQ4i0+Lkw2B0EViSxBsGgJldfbb1zrn3ATcCbzWzcRtsM2uLfe9yzj0KrAfGDQCRmcof0BmAJJdJdQE5564DPg3cbGb+M2yT75wrfPU1cC2wezL7FZlu4UiUPW2DALoRTJLGZA9l7gWy8bp1ALaa2Yedc7XAg2a2CagCHo2tzwC+Z2a/mOR+RaZc70iA51/p47cHunlibwf9/hDFuZnUluQmumgicTGpADCzRWdY3gZsir0+CqyezH5EplJrv5+tR/t4pWeEfn+Ifl+Qw10jHOoaAaAgO4Orl8/h+gtreNOSSnL0OEhJEurMlJTQNTzGc0f72N7Sz9BoiFDUCIWj7Gkf5HjfKADpaY6S3ExK87OoK8nl7evquGR+Oavqi8lM103zknwUADLrhSJRuocD9IwEGBoNMzgaos8fpKXHR3Ovj8NdIzT3epeo8rLSKcvPIjM9jYw0x7LqIj542Xw2LixnyZxC0tL0rF9JHQoAmRVGgxF6RrxG/mi3j10nBtl9YpBXenz0+oLj/kx2RhqN5fksqy7ijvUNbFhQzsraIjJ0NC8CKABkBhkaC7GvbYi97UPsax/ixMAoHYNjdA0FGA6EX7dtbmY6K2qLuHZlFVVFOcwpzKGiIIuSvCyKcjMoyc1iTmG2juhFzkIBINPOzGju9fNy6wC7Wgc52DXCoc5h2gfHTm5TUZBFQ1keS6sLuWJxJZWF2VQWZFNRmEV9aR4LKwtIV+MuMikKAIm7UCTKke4R9rcPs69jiAMdw5zoHyUcNUKRKIOjIYbHvCP67Iw0Fs0pYMOCchZXFbC8poiVNUXMKcpJcC1Ekp8CQCYsGjWae33sax9mX/sQff4go8EI/mCYY32jHO4aJhTxbgrPTHcsmlPIgsp8sjLSyUx35GdlsLK2iFX1JSypKlCfvEiCKADkjKJRo88fpGckQNvAKDuPDbD9WD87jw3gi02Mlp7mKM3LJDcrnbzMDKqLc3jTkkqW1xSyrLqIBZX5GkIpMkMpAIRBf4hnj/Tw+8M97GkbYtAfZHA0xNBYmEj0temd0tMcy6oLuXVdPRfWF7OipojFVQVkZ+jGKJHZSAGQ5ALhCF1DATqHxjjW56e5109zj4+u4TEGR8MMjYZoHxwlat4dr6vnFtNQVkJJbiYleZlUFGRTWZjNnMJsltcU6Xm4IklEf81JIByJsrttiK1He9nR0k/PSIABv3cz1IA/9Lpt0xzUluRSU5xDXUkuK2qKqC/N5fLFFayZW6LuGpEUogCYZcKRKI+91MYTezrp9QXo9QXpGBw7+bCSBRX51JTkUFuSS2meNxa+qiiHyqJs5pbmMbcsV102IgIoAGY0MyMQjhIIRwmGo/z2QBdf+81hmnv91JfmUl+ay/LqIt60pJKL5pVyyfxyKguzE11sEZklFAAzhJnR7w9xtHuEnccHeP6VPra19NN32jQHK2uL+MZ7LuKa5VW6y1VEJkUBkCCjwQjbWvrYcqSX517p41DnMENjr013MK88j6uWzWF+RT45melkZaTRWJ7H5YsqiD1bQURkUhQA02QkEGbzoW62NfezraWfPW2DhCJGeppjVX0xN6+pZX5FAfMr8lhZW0yV7oQVkSmmAJgikajRPjjKS8cH+dnLbfx6fxeBcJSsjDRW1xfzwcvns2FBORc3llGgoZUikgBqeeJgLBRhf8cw21v62dHSz772IVr7RwlGogBUFGRz+8Vz2XRhDWsbSsnK0FBLEUk8BcB5aBsY5aHNr7Dz+ADH+/x0DQdOrqsryeXCumKuWVnFvLJ8FlcVsK6hVDNXisiMowCYgNZ+P1//7RF+uO04ZnDRvFLetKSSuWV5LJrjNfTVxeq7T0pm0P4SHPg5XPmXkK4/HZn9JvVb7Jz7G+BPgO7Yos+a2ePjbHcd8BUgHXjQzP5+MvudbjuPD/Dg74/y890dpDl4Z9NcPvLmhdSX5iW6aDKVhjugYxc0b4a9/w79zeDSYen1ULsmwYUTmbx4HMZ82cy+dKaVzrl04GvANUAr8IJz7jEz2xuHfcedLxDmib0dNPf4OTEwyoGOYXadGKQwO4MPXtbIBy6bT21JbqKLKfE2OgBtO6B1G7S+AG07wdflrUvLgAVvhis+CUtvgPzyBBZUJH6m4zx2PXDYzI4COOd+ANwCzKgAGB4L8ciWFh78/VH6/SGcgzmF3vQJ/+umFdzWNFejdWa7aASG22HgOAy2wuBx6NwDbS9C35HYRg4ql8Hia6B6FVRf6H3lFCW06CJTIR4t2sedc+8FtgGfNLP+09bXAcdPed8KXHKmD3PO3QXcBdDQ0BCH4p1dOBLl4S0t3PPUIQZHQ1y1bA4fffNCLqwv1pw5s100Avt+Ci9+B3oPeY1+9PXPFqao3uvOWfMuqFsHdRdBTnFCiisy3d4wAJxzTwLV46y6G7gP+AJgse//CHzw9I8Y52dtnGXeCrMHgAcAmpqazrhdPLzQ3Mf//Pfd7O8Y5sollfzltUu5sF5//LNaaMxr7Js3w9b7YKAFShqg/mJY+XbvdXEDFNdDcR1kFya6xCIJ84YBYGZXn8sHOee+CfxsnFWtwNxT3tcDbedUuikyForwv/9zH/+ytYW6klzuf/dFvG1llaZYmE0iYejaA+0vQ89B6DkEPQe8C7Xm3X/B3Evg2r+DZTdAms7mRE432VFANWbWHnv7dmD3OJu9ACx2zs0HTgC3A++azH4no7nHx0e/u4O97UP8t8vn88lrl5CXpb79GS8agRM74PCT0PKM9zrk89alZ0H5Iq+v/sJ3QuVSqFrpfReRM5psy/d/nXNr8Lp0moEPATjnavGGe24ys7Bz7uPAL/GGgX7LzPZMcr/n5ee72vnLH71MRrrjW+9v4qplVYkohpxN0Ae+bhjpgr6j3tF99wFoeRZG+wAHNath7Z1Qv97rty9t1BG+yHmYVACY2XvOsLwN2HTK+8eBP7g/YLpEo8ZXnjrEV546xLqGEu591zoN5ZwJQqPQuRdaNkPzM3B8K4wNvn4blw5l871ROYuvhYVXQV5ZYsorkmSSvu/DHwzzqR++xOO7Orjtonr+7u0XaHRPIkTC3p20LZvh+PPQtQ/6X3mtv758May4xTuaz58DBXOgdL73PiMrkSUXSVpJHQBjoQjveeh5dhzr5+5Ny/njK+brQu9UCgeg97DXZdNz0BtvP9LpffUdheCIt13ZQq+P/sLbYM5yaNgAheMNNBORqZS0AWBmfPbRXWxv6efed63lxlW1iS5ScjGDoTboeBmOPwfHtnoXZiOvToznvEa9oAoKa7wROfMuhXmXQaGuvYjMBEkbAA9tfoWf7DjBX1y9RI1/PJh5Qy0P/RKO/NqbKmG0z1uXlgE1a2D9n0DtWqhY4o3KydJcSSIzWVIGwO8OdvPFx/dx/QXV/OlVixJdnNkj6PfmvxkbgsCQNxqnaz907YX2nTBwzNuucjksv+m1aRKqV6mxF5mFki4A+n1B/vR7O1hSVciXblutB6ePZ2zIm/+mbYfXbdNzCIbb/nAEDgAOyhZ4Qy8v+3NvJE7J1E/RISJTL+kCoDQ/iy/80QWsayglP9UnbwuNwdAJbwK0oTavsT/2rDfF8aujb0obYc5KaLwcimq8ETg5xd7kZ7llULEYMjVkViQZJWULecuaukQXIbGG2uCZe2D7tyE8+tryjFyob/IeaDJ3PdSu05h6kRSWlAGQUsy8/vqew97Qy2Nb4KXve1MnrPqvMP9K78i+sMYbV68x9SISowCYTQLDcPgpOPA4nNjuPcRkbBCiode2Sc+C1XfA5X/h3UErInIGCoCZLBr1Zrw88ms48htvErRIEHJLofEKyK+M9dWXesMuK5Z4ffrpmYkuuYjMAgqAmaT7IOz5iXeRtu8Vb6qEkN9bV7kc1t/lPY927gY9lFxEJk2tSCKZeY8kPPwk7P6xd1etS4sdyc+HBW/yxtkveDMU6WY2EYkvBcB0iIS8I/q+ozDUCkPt3tF982ZvnhzwRuRc9/feU6s0L46ITAMFwFTpOQw7vg0Hfu41/hZ5bZ1L947oG6/wpjde8Gbv8YQiItNIARBP0Sgc+E947hvQ/HtvjpyFV3nTHJcv9i7UFtd7Ux3rASYikmAKgHiIhLw+/M1fhu793lQJb/0crLlT3TkiMmMpAM7XUPtrM2Me/R2MDcCcFfBfHoIVf6RROiIy46mVmihfDzz9Jdj2kDcmv7AWlt0IK26GRddAWlqiSygick4UAOdqdAC23gdb7vXG5q+5EzZ+DCqXgZ4yJiKz0KQCwDn3r8DS2NsSYMDM1oyzXTMwDESAsJk1TWa/02p0wLuou/Vr3rQLy2+Cqz4HlUsSXTIRkUmZVACY2X999bVz7h+B8SaUf9VbzKxnMvubVr5eeO4+eO4BCAx63Txv+itvXnwRkSQQly4g5z1p/Z3AVfH4vIQaaodnY1Mph0Zh+Y3e9Mlq+EUkycTrGsAVQKeZHTrDegOecM4Z8A0ze+BMH+Scuwu4C6ChYRqfPDXcCc/8P9j2LW9Y56p3ejNqVi59wx8VEZmN3jAAnHNPAuMNZr/bzP4j9voO4Ptn+ZjLzKzNOTcH+JVzbr+ZPT3ehrFweACgqanJ3qh8kzbc4T08Zdu3vFE9q++AKz+lqZRFJOm9YQCY2dVnW++cywBuBS46y2e0xb53OeceBdYD4wbAtBnphqf/L2x/GKJh7+EpV34KyhcmtFgiItMlHl1AVwP7zax1vJXOuXwgzcyGY6+vBT4fh/2ev6F2ePhG6G+GNe+Cy/+7jvhFJOXEIwBu57TuH+dcLfCgmW0CqoBHvevEZADfM7NfxGG/52e4Ax6+yQuB9/8nNGxIWFFERBJp0gFgZu8fZ1kbsCn2+igwM4bQjHTFGv82ePeP1PiLSEpLnXkLjr8A33wrDLbCnT+EeZcmukQiIgmV/AEQjcIzX4F/vs6bsuH9P4PGyxJdKhGRhEvuuYDM4Cd/Art/BMtvhpu/CrkliS6ViMiMkNwBsPc/vMb/yr+Ct3xWk7aJiJwiebuAxgbh55+G6lXwpk+r8RcROU3yngE89XnwdcEd39fDWURExpGcZwDHn4cXHoL1H4K6dYkujYjIjJR8ARAJwU8/AYU1cNXdiS6NiMiMlXx9I+ExqFsLS66H7MJEl0ZEZMZKvgDILoRbvpboUoiIzHjJ1wUkIiLnRAEgIpKiFAAiIilKASAikqIUACIiKUoBICKSohQAIiIpSgEgIpKinJklugxn5JzrBlrO88crgJ44Fmc2SdW6p2q9QXVX3V8zz8wqz+WHZ3QATIZzbpuZNSW6HImQqnVP1XqD6q66nx91AYmIpCgFgIhIikrmAHgg0QVIoFSte6rWG1T3VDWpuiftNQARETm7ZD4DEBGRs1AAiIikqKQLAOfcdc65A865w865zyS6PFPJOTfXOfcb59w+59we59yfx5aXOed+5Zw7FPtemuiyTgXnXLpz7kXn3M9i71Ol3iXOuR855/bH/u83plDd/yL2u77bOfd951xOstbdOfct51yXc273KcvOWFfn3F/H2r0Dzrm3ncs+kioAnHPpwNeA64EVwB3OuRWJLdWUCgOfNLPlwAbgY7H6fgZ4yswWA0/F3iejPwf2nfI+Ver9FeAXZrYMWI33b5D0dXfO1QF/BjSZ2QVAOnA7yVv3bwPXnbZs3LrG/u5vB1bGfubrsfbwrJIqAID1wGEzO2pmQeAHwC0JLtOUMbN2M9sRez2M1xDU4dX54dhmDwN/lJACTiHnXD1wA/DgKYtTod5FwJXAQwBmFjSzAVKg7jEZQK5zLgPIA9pI0rqb2dNA32mLz1TXW4AfmFnAzF4BDuO1h2eVbAFQBxw/5X1rbFnSc841AmuB54AqM2sHLySAOQks2lT5f8BfAdFTlqVCvRcA3cA/x7q/HnTO5ZMCdTezE8CXgGNAOzBoZk+QAnU/xZnqel5tX7IFgBtnWdKPc3XOFQA/Bj5hZkOJLs9Uc87dCHSZ2fZElyUBMoB1wH1mthbwkTxdHmcV6+++BZgP1AL5zrl3J7ZUM8Z5tX3JFgCtwNxT3tfjnSImLedcJl7j/10z+0lscadzria2vgboSlT5pshlwM3OuWa8br6rnHPfIfnrDd7veKuZPRd7/yO8QEiFul8NvGJm3WYWAn4CXEpq1P1VZ6rrebV9yRYALwCLnXPznXNZeBdFHktwmaaMc87h9QXvM7N/OmXVY8D7Yq/fB/zHdJdtKpnZX5tZvZk14v0f/9rM3k2S1xvAzDqA4865pbFFbwX2kgJ1x+v62eCcy4v97r8V77pXKtT9VWeq62PA7c65bOfcfGAx8PwbfpqZJdUXsAk4CBwB7k50eaa4rpfjnea9DOyMfW0CyvFGCByKfS9LdFmn8N/gzcDPYq9Tot7AGmBb7P/934HSFKr73wL7gd3AvwDZyVp34Pt41zpCeEf4/+1sdQXujrV7B4Drz2UfmgpCRCRFJVsXkIiInCMFgIhIilIAiIikKAWAiEiKUgCIiKQoBYCISIpSAIiIpKj/D61re0TsVzv2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.log(linalg.eigvalsh(mean)))\n",
    "plt.plot(-np.log(sg))\n",
    "plt.plot(np.log(sp)[::-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "9e494c5e-7bcb-4c3c-a689-92a32a85e2fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "S, U = linalg.eigh(mean)\n",
    "U = subspace.T @ U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "d6b33960-cecb-4ec8-b606-b3131ff5de1b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 98)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "53caf398-c22f-4528-ba60-9359ca64fe4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "full = U @ (np.diag(S) - np.eye(len(mean))) @ U.T + np.eye(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "59b69485-3438-4c63-86fb-ea6e8ca22f45",
   "metadata": {},
   "outputs": [],
   "source": [
    "#linalg.eigvalsh(full)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "0bf4f469-0430-48c9-a0b3-13866a45b1b4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f7d605055b0>]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfi0lEQVR4nO3deXxU9b3/8dcnG5AIhH1LQhLZEQQMq4qK+27Vui+gFuvjam293rre29vf7a+1te2tWltFcSkqaAG3qnVDBRUIS1CWsAayEbIQQkIg63zvHxlbiiDLnGRyJu/n45FHZs7MfM/nO5B3vvnOOedrzjlERMS/osJdgIiIhEZBLiLicwpyERGfU5CLiPicglxExOdiwrHT7t27u9TU1HDsWkTEt1asWFHmnOtx4PawBHlqairLly8Px65FRHzLzHIPtl1TKyIiPqcgFxHxOQW5iIjPKchFRHxOQS4i4nMKchERn1OQi4j4nIJcRKQFFFfW8LsPNpBTusfzthXkIiItYFtZNU8s2EzR7hrP21aQi4i0gF176wBIjI/1vG0FuYhICyivrgega0Kc520ryEVEWsA3I/Iu8QpyERFfKti1l8T4WNrHRnvetoJcRKQFrC7czQl9OzdL2wpyEZFmVtcQYMOOKob369Qs7SvIRUSa2cbiKuobHSP6teIRuZklmtlcM1tvZtlmNtGLdkVEIkFWfgVAs02teLVC0GPA351zV5pZHBDvUbsiIr43b0UBx/dIoH+35onGkEfkZtYJmAzMBHDO1TnnKkJtV0QkEqwp3M2q/AquH98fM2uWfXgxtZIOlALPm1mWmT1rZgkHPsnMppvZcjNbXlpa6sFuRURav1cy82gXE8UVY5KabR9eBHkMMAb4s3NuNFAN3H/gk5xzM5xzGc65jB49vrUItIhIxNlT28CbWYVcfGJfOjfDqfnf8CLIC4AC59zS4P25NAW7iEib9uaqQqrrGrlufEqz7ifkIHfO7QDyzWxwcNOZwLpQ2xUR8bOa+kZmfr6VoX06MTo5sVn35dVRK3cBLwePWMkBpnnUroiIL72eVUhOaTUzb85otg85v+FJkDvnVgEZXrQlIuJ3zjlezyqkf7d4pgzp2ez705mdIiIee2LBZjK3ljN1Umqzj8ZBQS4i4qmtZdX84aONXDaqL1MnpbbIPhXkIiIecc7xuw82EBMdxYMXDm2R0TgoyEVEPLN4y07+9nURd54xgJ4d27fYfhXkIiIeqG8M8PsPN5IYH8v0yektum8FuYiIB578ZDPLc3dx5xkDmmUVoO+iIBcRCdG7q4t47ONNnDW0J7ed2rKjcVCQi4iEpKSyhv/39jrSuifw2DWjw1KDglxE5BjtqW3gxpmZVOyr42cXDyehnVcnyx+d8OxVRMTnGgOOH8/JYnPpHv5yyzhOHtA9bLVoRC4icgx+8uoqPsou4WcXDwtriIOCXETkqGXl7eKtr7bzg1PTuGliarjLUZCLiByN3fvqufevX5EYH8vdZw0KdzmA5shFRI5YbUMjVz+9mJyyambdMp7jwvTh5oE0IhcROQIFu/Zy1dNLWL+jil9fPpJTBoZ3Xnx/rePXiYhIK1axt45rZiyhtKqW3191Ipc340LKx0JBLiLyHZxz3Dfva4ora3jmpgxOH9z8C0UcLQW5iMghlFbV8qt3s3l/bTEPXzi0VYY4KMhFRA6qYm8dlz35BYUV+7j7zIHcekpauEs6JAW5iMgB8nbu5d65X1FSVcPcH04kI7VruEv6Tp4FuZlFA8uBQufcRV61KyLSktYU7ubaGUuoawzw2++f2OpDHLwdkd8NZAOdPGxTRKTFfL6pjNtnLScxPo450yeQ3DU+3CUdEU+OIzezJOBC4Fkv2hMRaWmfbCjh1heX0SEuxlchDt6dEPQH4KdA4FBPMLPpZrbczJaXlpZ6tFsRkdDNzsxj2vPL6JvYgTf+bZKvQhw8mFoxs4uAEufcCjM7/VDPc87NAGYAZGRkuFD3KyISqpLKGn40J4slOeWcMqA7z9yUQYe4ll2mzQtezJGfDFxiZhcA7YFOZvaSc+4GD9oWEWkWn20s5cH5qynavY8fTRnAXWcOJDban1ctCblq59wDzrkk51wqcA2wQCEuIq2Vc455KwqY+nwm0VHGrFvHc885g30b4qDjyEWkDanYW8dds7NYtKmMoX068dcfTmw1VzAMhac9cM59CnzqZZsiIl7I3VnNzc9lUlixj/vOG8IPTk0jxsej8P35/1eRiMhhzF1RwH+/tRYzeHHaOCaFeWk2rynIRSRiFe3ex6Pvb2D+ykKG9O7IMzdl+O7QwiOhIBeRiLQqv4LrnllCQ6Pj9snp/OjMgSREwHz4wURmr0SkzXLO8fbXRdz716/onhDH7OkT6N8tIdxlNSsFuYhEjMyt5fz2gw1kbi1nUK/j+MPVoyM+xEFBLiIRYGNxFb9+bz0fry+ha0Icv/zeCK4em0x0lIW7tBahIBcRX3ttWT4PvbGa+LgY/uPcwdxycpovT7MPhYJcRHwpu6iS//9ONp9vLmNkUmeenzqWbse1C3dZYaEgFxFf2VvXwB8+2sSzi3JIaBfDj88ayA9PO572sW1rFL4/BbmI+IJzjndWF/E/f1tHcWUtF47ow88vHU73NjoK35+CXERavay8XTwwfzXrd1SR2i2eF6aNbbUr2oeDglxEWq2NxVXMWpzLy0tz6d2pPb+5YiSXj+kXMddI8YqCXERandyd1Tz1WQ5zluURE2VcPTaFBy4YQqf2seEurVVSkItIq7FueyX/+9FGPsouJtqMa8elcO85g+maEBfu0lo1BbmIhF1O6R7+/OkW5q0sIDY6iltOTuPmiamkdIu8C1w1BwW5iITNl5vLeGlpLu+t2UG0GTdO6M89Zw+mc7ymUI6GglxEWtym4ir+551sFm4sJTE+lpsnpnL7aen06dwh3KX5koJcRFqEc47PNpby5CebWbZtF507xHL/+UO4fnwKHfUhZkgU5CLSrJxzfLCumMc/3sTa7ZUkxsfy72cP4rrxKW32lHqvKchFpFnUNwb4OLuEJxY0BXh69wQeuXwEl47q1+YuatXcFOQi4qnKmno+WFvMMwtz2FBcRa9O7fjl90ZwVUaSTuRpJiEHuZklA38BegMBYIZz7rFQ2xURfymprOHFxdt4dtFWahsCpPdI4LFrRnHu8N5t+oJWLcGLEXkD8O/OuZVm1hFYYWYfOufWedC2iLRyVTX1PLFgM7OX5lFV28D5J/TmB5PTGZWUSFQbWdgh3EIOcudcEVAUvF1lZtlAP0BBLhLB8sv38uj7G/hg3Q7qGgKcNbQXd04ZwIh+nTFTgLckT+fIzSwVGA0sPchj04HpACkpKV7uVkRa0OqC3Ty9cAsfrC2m0TmuGNOP68b3Z1RyYrhLa7M8C3IzOw6YB/zYOVd54OPOuRnADICMjAzn1X5FpPnVNQT4KLuYZxblkJVXQaf2MVw1Nonppx6v0+hbAU+C3MxiaQrxl51z871oU0TCb1tZNbMz85i7ooCd1XX07dyehy8cyvdG99Mx4K2IF0etGDATyHbO/T70kkQk3PZflT46yjhraE+uHZfCqQN7tJmV6f3EixH5ycCNwGozWxXc9qBz7l0P2haRFlJeXceH63Ywb0UhmdvKOS64Hua141Lo1al9uMuT7+DFUSufA/oVLeJTTYs4bOGNrO3sq28ktVs8D5w/hCtOStJ6mD6hMztF2qCGxgCvZObx9zU7WJyzk9joKC4c0YdbT0ljeN9OOnzQZxTkIm3I9op9vLw0l1eXFVC2p5ZBvY5j+uR0pk5K1SVkfUxBLhLhnHN8uWUnsxbn8sG6HQBMGdKLa8Ymc+bQnhp9RwAFuUiE2lVdx7yVBbyyNI+csmq6xMdy+2nHc/34FJK66NjvSKIgF4kwK3J38eyiHD7KLqa+0ZHRvwt3nTmA80/oo4tXRSgFuUgE+Oba3y8vzWXRpjIS42O5aWIq389IYkjvTuEuT5qZglzEx/LL9/LqsnxeW55PSVUtvTu1595zBnHzpFQtn9aGKMhFfGZfXSOzlmzjs42lfLllJwacMbjpzMvTB/fQ4g1tkIJcxCcKdu3lxS+3MW9lIeXVdQzqdRx3TRnI1WOT6ZeoQwfbMgW5SCu2e189n20s5a1VhXy6oRSAyYN6cNupaUw6vnuYq5PWQkEu0gqVVtUyf2UBzyzaStmeWrolxDF1Uiq3nJJGX42+5QAKcpFWojHgWLSplJeW5PHJhhIaA45RyYk8euVIJg/SVQfl0BTkImGWU7qHuSsKmL+ykB2VNXRsF8PUSalcOy6ZAT07hrs88QEFuUgYrN9RyXurd/DF5jKW5+4iyuD0wT352cXDmDK0J+1idOKOHDkFuUgLWpVfwV++3Mb8rEKiDAb16sj95w/h8tH96KlrfssxUpCLNLPq2gbeWV3Ea8vyWZ67i7iYKG47JY07Tj9ey6WJJxTkIs2gpr6R99YUsWhjGe+v3UF1XSPpPRJ48IIhXHlSMl0T4sJdokQQBbmIhypr6pm1OJeXluRStLuGLvGxnD+iD9eMTeak/l10yVhpFgpykRAFAo4vtpTxwdpi3sgqpKq2gXFpXfnt909kYno3onTYoDQzBbnIMarYW8dry/N5aUkeeeV7iYtpWi7thgkpjEnR6FtajoJc5CjsqW3g72t28OaqQjK3llPbEGBcalf+49zBnDW0Fx3idNigtDxPgtzMzgMeA6KBZ51zj3jRrkhr4Jxj8ZadzFmWzwfrdlBTHyC5aweuG5/CVRnJDO2j631LeIUc5GYWDTwJnA0UAMvM7C3n3LpQ2xbxknOOgIOGQIDGgKMh4GhsbPreEAjQ0Oj+uT3gqKlv5J3VRby7uoiCXfvo3CGWK8YkcfmYfpo6kVbFixH5OGCzcy4HwMzmAJcCCnIJi+rg9Mf8rAI2Fu+hpr6RmvpG6hvdUbcVZTBlSE9+NGUgl47uqzMupVXyIsj7Afn73S8Axh/4JDObDkwHSElJ8WC3Iv9U1xBgdeFuXluWz9tfb2dvXSP9u8UzZXBPOsRF0z42mrhoIyY6iugoIybK/vE9Jjrqn/ejjeioKGKD94f26URyVy1ULK2bF0F+sL8vvzX0cc7NAGYAZGRkHP3QSOQgSqpqeGlJHrMz8yitqiU+LpqLRvbh+xnJZOi4bWkjvAjyAiB5v/tJwHYP2hU5KOccm0v28MSCzby3poj6Rsf4tK48dMFQzhjSk84dtFaltC1eBPkyYKCZpQGFwDXAdR60K/Iv9tY18PjHm/n7miK27dxLfFw0N0zozxVjkjihX+dwlycSNiEHuXOuwczuBN6n6fDD55xza0OuTISm0feq/ArmrSzgrVXbqaxpYGJ6N249NZ1zhvWil64YKOLNceTOuXeBd71oSwQgd2c176wuYt6KAraUVtMuJopzh/fmpon9yUjtGu7yRFoVndkprUZtQyMLskuYu6KABRtKcA7GpnbhB6emc8HIPnRqr7lvkYNRkEvYrSnczV+X5/P210WUV9fRq1M7/u30AVw9NlmH/okcAQW5hEVdQ4BFm0p5+rMcMreVExcTxdlDe/H9jCROHaiFhkWOhoJcWtSW0j3MWpzL377eTtmeOvp2bs/DFw7lklF96dlRH1yKHAsFuTS7xoDjk/UlvLh4G4s2lREXHcWUIT25ZFRfzh7Wi9joqHCXKOJrCnJpNrk7q3k9q5B5KwvIL99H707tufecQVw9NoUeHbVWpYhXFOTiqYbGAPNWFvCXxbms3V4JwKTju/HA+UM5Z1gvYjT6FvGcglw80dAYYP7KQp5euIUtpdWk90jgvy8exmmDe5LWPSHc5YlENAW5hKS4sob5KwuZsyyP3J17Se+RwFM3nMS5w3vpglUiLURBLsdk/Y5Knlm4lTdXFdIQcIxM6syT143hvBN669BBkRamIJcjVtcQ4I1VhTyzMIdNJXtoHxvFDRP6c1VGMsP6arkzkXBRkMth7dxTy5xl+bz45TZKqmpJ6tKB/7poGJeN7kfXhLhwlyfS5inI5ZCqaup57vNtvJKZS3FlLScmJ/KfFw3jghF9NH0i0oooyOVbtpVV8+G6YmZ+vpXiqhomHd+NGTdmcGJyYrhLE5GDUJDLPyzaVMqfPtnC4pydAJyY1JmnbjyJUQpwkVZNQS58lV/B7z7cyMKNpfTo2I77zhvCxSf2IamLrjwo4gcK8jaqMeD4fHMZsxZv46PsEjp3iOXhC4dy08RU4mJ09qWInyjI25jGgGP+ygJ+/+FGinbXkBgfy73nDGLqyWkc107/HUT8SD+5bUR5dR2zM/OYnZlHwa59pHVP4E/Xj+HMoT1pFxMd7vJEJAQK8gi3pnA3Ly/N5fWsQmrqA0xM78b95w/h3OG9dflYkQgRUpCb2aPAxUAdsAWY5pyr8KAuCdHXBRU8vTCHd74uIj4umotG9uX2yekM7NUx3KWJiMdCHZF/CDzgnGsws18DDwD3hV6WHKsvN5fx1MIcFm4spUNsNHeeMYDpp6Vr4WKRCBZSkDvnPtjv7hLgytDKkWP1/todPPf5VpZuLadHx3bcc/Ygpp6cqgAXaQO8nCO/BXj1UA+a2XRgOkBKSoqHu227nHO8v7aYmZ/nsGzbLpK7duDhC4dy48T++gBTpA05bJCb2UdA74M89JBz7s3gcx4CGoCXD9WOc24GMAMgIyPDHVO18g/ZRZX84p11fLF5J0ldOvCfFw3jxgn9dQy4SBt02CB3zp31XY+b2c3ARcCZzjkFdDNbv6OSGQtzeCOrkI7tY3nogqFMOzlVS6iJtGGhHrVyHk0fbp7mnNvrTUlyMCWVNTy+YBOvLM2jXUw0005O464pA0iM12VkRdq6UOfI/wi0Az4MLuu1xDn3w5Crkn/YvbeemV9s5Y8LNhFwcPmYfvznhcPoouuAi0hQqEetDPCqEPlXu6rreHzBJmZn5lFTH2DyoB48cP4QhvbRSjwi8q90ZmcrU1JVw1Of5vDS0lzqGwNcOSaJqSenMrxv53CXJiKtlIK8ldhX18gLX27jjws2UdMQ4NJRfbn1lDQFuIgcloI8zJxzzFtZyC/fzaa8uo4Tkzrzi8tGMCJJAS4iR0ZBHkbZRZU8sWAT767eweiURB69ciRnDO5JlNbDFJGjoCAPg80lVfzinWw+3VBKXHQU954ziDtOH6AFjUXkmCjIW1BjwDE7M4+fv72WDrHR/OSsQdwwIYVux7ULd2ki4mMK8hZQ3xjg72t28Mh76yms2Mfwvp14Ydo4enRUgItI6BTkzagx4Hgjq5BnFuWwfkcVKV3j+dP1YzhveG/Ng4uIZxTkzaRo9z7unrOKzK3l9O8Wz2+uGMmVJyUpwEXEcwpyj31zadl7XlsFwC8uO4GrxyZrWTURaTYKcg+VV9dxz2ur+HRDKand4nlu6ljSexwX7rJEJMIpyD1Q3xjg9axCHnlvPVU19fz0vMFMm5RGhzgt7iAizU9BHqIvNpdx37yvKdi1j1HJiTxyxQiG9NaFrUSk5SjIj9HOPbU8/MYa3luzg76d2zPz5gydlSkiYaEgPwYLN5Zyz2urqKpp4EdnDuQHp6bRUYsci0iYKMiP0tOfbeFX760ntVs8L0wbxwn9dHErEQkvBfkRqqlv5BfvrOOlJXmMS+vKMzdm0Dleo3ARCT8F+RFYlV/BT+d+xcbiPVw7LoWfXTyM9rE6IkVEWgcF+XcIBBzvrC7intdWkRgfx+PXjuaSE/uGuywRkX+hID+E3fvqmfZ8JivzKhjSuyOzbh2vi1yJSKukID+I3J3VTHt+GXnle/n1FSO4YkwSMTrFXkRaKU/SyczuNTNnZt29aC+cPs4u5rInv6BsTy0v3zaeq8emKMRFpFULeURuZsnA2UBe6OWEj3OOTzeUcsdLK0nu2oGnb8xgQE9dJ0VEWj8vplb+F/gp8KYHbYVFTX0jt89awWcbS+mX2IG//nASXRPiwl2WiMgRCSnIzewSoNA595XZd5+abmbTgekAKSkpoezWU1U19fzXm2v5bGMp954ziNtOTdehhSLiK4cNcjP7COh9kIceAh4EzjmSHTnnZgAzADIyMtxR1NhsNpfs4caZSynaXcNPzhrEnVMGhrskEZGjdtggd86ddbDtZjYCSAO+GY0nASvNbJxzboenVTaD8uo6pr2QSW1DgBk3nsQ5ww/2u0pEpPU75qkV59xqoOc3981sG5DhnCvzoK5mlbdzL7e8uIziylrmTJ/AmJQu4S5JROSYtbnjyL/cUsatLywnJtqYdcs4hbiI+J5nQe6cS/WqreayNGcnNz+XSYfYaObeMYlBvTqGuyQRkZC1mRH52u27mT5rBcld4pl7hw4vFJHI0SZOWSyurOGu2VlERxlP33iSQlxEIkrEB3kg4Hjo9TUUVdTwp+vHMFDTKSISYSJ6aqWmvpFbXljGl1t2Mu3kVCakdwt3SSIinovoIJ/5+Va+3LKT2yenc++5g8NdjohIs4jYIF++rZxH39/A5EE9uO+8IVrdXkQiVkTOkTcGHD97ay19OrfnqRvGKMRFJKJFXJDXNwb48aurWLu9kgcvGEp8XMT+0SEiAkRgkM/8fCtvf7Wd+84bwsVaX1NE2oCICvKtZdU8/vEmTurfhTtOPz7c5YiItIiICfLKmnpufXEZAI9cPiLM1YiItJyImUB+9+sickqreW5qhk76EZE2JSJG5KVVtfzuw4307xbPGYN7Hv4FIiIRJCJG5L96L5vSqqZrix9uyTkRkUjj+xH50pydzF9ZyG2npOkUfBFpk3wf5K9nFZIQF61T8EWkzfJ1kK/I3cWry/O5YEQfrXwvIm2Wr4P85SW5tIuJ4uELh4W7FBGRsPFtkC/fVs78rEJumphK5/jYcJcjIhI2vg3yV5flEx8Xzd1nDgx3KSIiYeXLIK+ubeCd1UVcPLIvCe0i4ghKEZFjFnKQm9ldZrbBzNaa2W+8KOpwsvIq2FvXyAUj+7TE7kREWrWQhrNmdgZwKTDSOVdrZi1yWuVXBRUAjEpObIndiYi0aqGOyO8AHnHO1QI450pCL+nwluTsJL1HAp076ENOEZFQg3wQcKqZLTWzz8xs7KGeaGbTzWy5mS0vLS095h2WVNawaFMZ5w3vfcxtiIhEksNOrZjZR8DBUvOh4Ou7ABOAscBrZpbunHMHPtk5NwOYAZCRkfGtx4/U2qJKAE4b1ONYmxARiSiHDXLn3FmHeszM7gDmB4M708wCQHfg2Ifch5EdDPIhfTo11y5ERHwl1KmVN4ApAGY2CIgDykJs8zutLaykX2IHzY+LiASFGuTPAelmtgaYA9x8sGkVrzQGHF9uKdNVDkVE9hPS4YfOuTrgBo9qOax12yvZtbeeyYO6t9QuRURaPV+d2ZlTtgeAYZofFxH5B18FeWlVLQA9O7UPcyUiIq2Hr4K8uLKG9rFRdGqv66uIiHzDV0F+fI/juOTEvlqXU0RkP74a2l4zLoVrxqWEuwwRkVbFVyNyERH5NgW5iIjPKchFRHxOQS4i4nMKchERn1OQi4j4nIJcRMTnFOQiIj5nzXjV2UPv1KwUyD3Gl3enma953gqpz22D+tw2hNLn/s65by2PFpYgD4WZLXfOZYS7jpakPrcN6nPb0Bx91tSKiIjPKchFRHzOj0E+I9wFhIH63Daoz22D53323Ry5iIj8Kz+OyEVEZD8KchERn/NVkJvZeWa2wcw2m9n94a7HC2aWbGafmFm2ma01s7uD27ua2Ydmtin4vct+r3kg+B5sMLNzw1d9aMws2syyzOxvwfsR3WczSzSzuWa2PvjvPbEN9Pknwf/Xa8xstpm1j7Q+m9lzZlZiZmv223bUfTSzk8xsdfCxx+1olkJzzvniC4gGtgDpQBzwFTAs3HV50K8+wJjg7Y7ARmAY8Bvg/uD2+4FfB28PC/a9HZAWfE+iw92PY+z7PcArwN+C9yO6z8CLwG3B23FAYiT3GegHbAU6BO+/BkyNtD4Dk4ExwJr9th11H4FMYCJgwHvA+Udag59G5OOAzc65HOdcHTAHuDTMNYXMOVfknFsZvF0FZNP0A3ApTT/4BL9fFrx9KTDHOVfrnNsKbKbpvfEVM0sCLgSe3W9zxPbZzDrR9AM/E8A5V+ecqyCC+xwUA3QwsxggHthOhPXZObcQKD9g81H10cz6AJ2cc4tdU6r/Zb/XHJafgrwfkL/f/YLgtohhZqnAaGAp0Ms5VwRNYQ/0DD4tUt6HPwA/BQL7bYvkPqcDpcDzwemkZ80sgQjus3OuEPgtkAcUAbudcx8QwX3ez9H2sV/w9oHbj4ifgvxg80URc+ykmR0HzAN+7Jyr/K6nHmSbr94HM7sIKHHOrTjSlxxkm6/6TNPIdAzwZ+fcaKCapj+5D8X3fQ7OC19K0xRCXyDBzG74rpccZJuv+nwEDtXHkPrupyAvAJL3u59E059pvmdmsTSF+MvOufnBzcXBP7cIfi8Jbo+E9+Fk4BIz20bTFNkUM3uJyO5zAVDgnFsavD+XpmCP5D6fBWx1zpU65+qB+cAkIrvP3zjaPhYEbx+4/Yj4KciXAQPNLM3M4oBrgLfCXFPIgp9MzwSynXO/3++ht4Cbg7dvBt7cb/s1ZtbOzNKAgTR9SOIbzrkHnHNJzrlUmv4dFzjnbiCy+7wDyDezwcFNZwLriOA+0zSlMsHM4oP/z8+k6TOgSO7zN46qj8HplyozmxB8r27a7zWHF+5PfI/y0+ELaDqqYwvwULjr8ahPp9D0J9TXwKrg1wVAN+BjYFPwe9f9XvNQ8D3YwFF8st0av4DT+edRKxHdZ2AUsDz4b/0G0KUN9PnnwHpgDTCLpqM1IqrPwGyaPgOop2lkfeux9BHICL5PW4A/Ejzz/ki+dIq+iIjP+WlqRUREDkJBLiLicwpyERGfU5CLiPicglxExOcU5CIiPqcgFxHxuf8Df6FPWRJDElgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.log(linalg.eigvalsh(cov, full)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "4050720e-69df-47fe-ab07-7d1a6046168b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7319.017302693897"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Distance to true covaritance matrix\n",
    "(np.log(linalg.eigvalsh(full, cov)) ** 2).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "c02974ec-1f60-4d6c-a062-f53ebb252e78",
   "metadata": {},
   "outputs": [],
   "source": [
    "#(np.log(linalg.eigvalsh(full_max, cov)) ** 2).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "cb5582dd-efb3-471d-9508-4089a8d4461c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7140.027986303995"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Distance of perfect draws diag approx\n",
    "(np.log(linalg.eigvalsh(np.diag(points.var(0)), cov)) ** 2).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "ca3c3cbf-e17d-4787-b983-20ce12c234ac",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3951.497131807437"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Diag approx using draws and grads\n",
    "(np.log(linalg.eigvalsh(np.diag(np.sqrt(points.var(0) / grads.var(0))), cov)) ** 2).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "c547d39c-f315-4383-bff2-b9c0f88b7d95",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/adr/mambaforge/envs/pymc-dev/lib/python3.9/site-packages/tqdm/auto.py:22: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
      "  from .autonotebook import tqdm as notebook_tqdm\n",
      "2022-08-19 16:08:35.041350: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcudart.so.11.0'; dlerror: libcudart.so.11.0: cannot open shared object file: No such file or directory\n",
      "2022-08-19 16:08:35.041392: I tensorflow/stream_executor/cuda/cudart_stub.cc:29] Ignore above cudart dlerror if you do not have a GPU set up on your machine.\n",
      "/home/adr/mambaforge/envs/pymc-dev/lib/python3.9/site-packages/pkg_resources/__init__.py:123: PkgResourcesDeprecationWarning: dev is an invalid version and will not be supported in a future release\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "import covadapt\n",
    "import covadapt.spd_manifold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "c3e44ddf-a48d-4a8d-b747-708008ae865f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#solver = covadapt.spd_manifold.ManifoldSolver()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "59a1eb89-7811-4574-9d21-42f974d5620c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#estimator = covadapt.spd_manifold.CovarianceEstimator(solver.estimator, compute_full_matrix=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "6aa516bb-8299-4382-af6b-638d90ca5a7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "optimizer = covadapt.spd_manifold.PerturbingConjugateGradient(\n",
    "    min_step_size=1e-6,\n",
    "    min_gradient_norm=1e-2,\n",
    "    perturbance_cooldown=20,\n",
    "    perturbance=0.01,\n",
    "    #max_iterations=self._maxiter,\n",
    "    perturbance_stepsize_threshold=1e-4,\n",
    "    # beta_type=pymanopt.solvers.conjugate_gradient.BetaTypes.PolakRibiere,\n",
    "    #verbosity=2,\n",
    ")\n",
    "\n",
    "solver = covadapt.spd_manifold.ManifoldSolver(\n",
    "    #n_eigs=self._n_eigs,\n",
    "    n_eigs=20,\n",
    "    stiefel_retraction=True,\n",
    "    exp_stiefel=True,\n",
    "    optimizer_args={\n",
    "        \"verbosity\": 2,\n",
    "        \"min_step_size\": 1e-8,\n",
    "        \"min_gradient_norm\": 1e-2,\n",
    "        #\"max_iterations\": self._maxiter,\n",
    "    },\n",
    "    solver_method=optimizer,\n",
    "    #alpha=self._alpha,\n",
    "    #delta=self._delta,\n",
    "    #beta=self._beta,\n",
    "    #gamma=self._gamma,\n",
    "    gamma=10,\n",
    ")\n",
    "\n",
    "estimator = covadapt.spd_manifold.CovarianceEstimator(solver, compute_full_matrix=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "735806fd-b5a1-46f8-aab9-cfb12e190e16",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 24s, sys: 1.1 s, total: 1min 25s\n",
      "Wall time: 21.8 s\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<style>#sk-container-id-1 {color: black;background-color: white;}#sk-container-id-1 pre{padding: 0;}#sk-container-id-1 div.sk-toggleable {background-color: white;}#sk-container-id-1 label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.3em;box-sizing: border-box;text-align: center;}#sk-container-id-1 label.sk-toggleable__label-arrow:before {content: \"▸\";float: left;margin-right: 0.25em;color: #696969;}#sk-container-id-1 label.sk-toggleable__label-arrow:hover:before {color: black;}#sk-container-id-1 div.sk-estimator:hover label.sk-toggleable__label-arrow:before {color: black;}#sk-container-id-1 div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}#sk-container-id-1 div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}#sk-container-id-1 input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}#sk-container-id-1 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {content: \"▾\";}#sk-container-id-1 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}#sk-container-id-1 div.sk-estimator {font-family: monospace;background-color: #f0f8ff;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;margin-bottom: 0.5em;}#sk-container-id-1 div.sk-estimator:hover {background-color: #d4ebff;}#sk-container-id-1 div.sk-parallel-item::after {content: \"\";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}#sk-container-id-1 div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-serial::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: 0;}#sk-container-id-1 div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;padding-right: 0.2em;padding-left: 0.2em;position: relative;}#sk-container-id-1 div.sk-item {position: relative;z-index: 1;}#sk-container-id-1 div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;position: relative;}#sk-container-id-1 div.sk-item::before, #sk-container-id-1 div.sk-parallel-item::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: -1;}#sk-container-id-1 div.sk-parallel-item {display: flex;flex-direction: column;z-index: 1;position: relative;background-color: white;}#sk-container-id-1 div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}#sk-container-id-1 div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}#sk-container-id-1 div.sk-parallel-item:only-child::after {width: 0;}#sk-container-id-1 div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0 0.4em 0.5em 0.4em;box-sizing: border-box;padding-bottom: 0.4em;background-color: white;}#sk-container-id-1 div.sk-label label {font-family: monospace;font-weight: bold;display: inline-block;line-height: 1.2em;}#sk-container-id-1 div.sk-label-container {text-align: center;}#sk-container-id-1 div.sk-container {/* jupyter's `normalize.less` sets `[hidden] { display: none; }` but bootstrap.min.css set `[hidden] { display: none !important; }` so we also need the `!important` here to be able to override the default hidden behavior on the sphinx rendered scikit-learn.org. See: https://github.com/scikit-learn/scikit-learn/issues/21755 */display: inline-block !important;position: relative;}#sk-container-id-1 div.sk-text-repr-fallback {display: none;}</style><div id=\"sk-container-id-1\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>CovarianceEstimator(compute_full_matrix=True,\n",
       "                    estimator=&lt;covadapt.spd_manifold.ManifoldSolver object at 0x7f7cf31be370&gt;)</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item sk-dashed-wrapped\"><div class=\"sk-label-container\"><div class=\"sk-label sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-1\" type=\"checkbox\" ><label for=\"sk-estimator-id-1\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">CovarianceEstimator</label><div class=\"sk-toggleable__content\"><pre>CovarianceEstimator(compute_full_matrix=True,\n",
       "                    estimator=&lt;covadapt.spd_manifold.ManifoldSolver object at 0x7f7cf31be370&gt;)</pre></div></div></div><div class=\"sk-parallel\"><div class=\"sk-parallel-item\"><div class=\"sk-item\"><div class=\"sk-label-container\"><div class=\"sk-label sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-2\" type=\"checkbox\" ><label for=\"sk-estimator-id-2\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">estimator: ManifoldSolver</label><div class=\"sk-toggleable__content\"><pre>&lt;covadapt.spd_manifold.ManifoldSolver object at 0x7f7cf31be370&gt;</pre></div></div></div><div class=\"sk-serial\"><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-3\" type=\"checkbox\" ><label for=\"sk-estimator-id-3\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">ManifoldSolver</label><div class=\"sk-toggleable__content\"><pre>&lt;covadapt.spd_manifold.ManifoldSolver object at 0x7f7cf31be370&gt;</pre></div></div></div></div></div></div></div></div></div></div>"
      ],
      "text/plain": [
       "CovarianceEstimator(compute_full_matrix=True,\n",
       "                    estimator=<covadapt.spd_manifold.ManifoldSolver object at 0x7f7cf31be370>)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "estimator.fit(points, grads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "94a31b96-d8a7-421c-8b87-0c19d66f530d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3755.4621520250294"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(np.log(linalg.eigvalsh(estimator.covariance_, cov)) ** 2).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "e44c9d2b-7c43-4eec-b76b-6dde47d8fc28",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 4.16194474e+01, -1.01955801e+00,  9.22741697e+00, ...,\n",
       "         2.24696907e+00,  2.65851672e+00, -4.15119837e+00],\n",
       "       [-1.01955801e+00,  9.44641100e+00,  2.93655165e+00, ...,\n",
       "         3.68409027e-01,  1.73689467e+00, -2.56409377e-02],\n",
       "       [ 9.22741697e+00,  2.93655165e+00,  4.98757181e+01, ...,\n",
       "         2.49232691e+00,  3.76440338e+00,  6.31676795e-02],\n",
       "       ...,\n",
       "       [ 2.24696907e+00,  3.68409027e-01,  2.49232691e+00, ...,\n",
       "         4.03760086e+00,  8.59499682e-01, -1.44358900e-01],\n",
       "       [ 2.65851672e+00,  1.73689467e+00,  3.76440338e+00, ...,\n",
       "         8.59499682e-01,  7.07855407e+00, -5.33890328e-01],\n",
       "       [-4.15119837e+00, -2.56409377e-02,  6.31676795e-02, ...,\n",
       "        -1.44358900e-01, -5.33890328e-01,  1.67265877e+01]])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5604fe8-0c7e-4a27-b509-de4462328194",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pymc",
   "language": "python",
   "name": "pymc"
  },
  "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}