ellipsoid-basin-math.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = 1000\n",
    "A_diag = np.exp(np.random.normal(size=(d)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(np.float64(1.6949800788650446), np.float64(0.9973603200322416))"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_diag.mean(), np.exp(np.log(A_diag).mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.float64(1.6487212707001282)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.exp(1/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.diag(A_diag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "f(x) = x^T A x \\\\\n",
    "\\text{want: volume where} f(x) < 1 \\\\\n",
    "\\text{axis radii:} A_i r_i^2 = 1 \\\\\n",
    "r_i = A_i^{-1/2}\\\\\n",
    "vol \\sim \\prod_i r_i = \\exp(-1/2 \\sum_i \\log A_i)\\\\\n",
    "logvol = -1/2 \\sum_i \\log A_i + C\\\\\n",
    "\\approx -\\tfrac d2 E[\\log A_i] + C\n",
    "$$\n",
    "\n",
    "Naive estimator:\n",
    "$$\n",
    "u_i \\sim N(0, 1/d)\\\\\n",
    "x_i = r u_i\\\\\n",
    "f(x) = x_i A_ij x_j = \\sum_i A_i x_i^2\\\\\n",
    "\\approx  \\sum_i A_i r^2 / d = 1\\\\\n",
    "r \\approx \\sqrt{d / \\sum_i A_i}\\\\\n",
    "\\approx \\sqrt{1 / E[A_i]}\\\\\n",
    "\\log r \\approx -1/2 \\log E[A_i]\\\\\n",
    "logvol = d \\log r + C \\approx -\\tfrac d2 \\log E[A_i] + C\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAGdCAYAAAA1/PiZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA5hElEQVR4nO3deXyU5b338e9kX0gmG1lJCDuymCBLRFzApiL1uGBrW2st0h6sLbZ6aF04Pa2nzzkWn/rUY63T2vY51tr2qPVRaF1bjChqo0AgQNgjW8hKCJnJvsxczx8xoxGUBCZzz0w+79crr5dz3zczv7mU3F+v+1psxhgjAACAIBBmdQEAAACDRXABAABBg+ACAACCBsEFAAAEDYILAAAIGgQXAAAQNAguAAAgaBBcAABA0IiwugBf83g8qqmpUUJCgmw2m9XlAACAQTDGqKWlRdnZ2QoL++R+lZALLjU1NcrNzbW6DAAAcBaqqqo0ZsyYTzwfcsElISFBUt8XT0xMtLgaAAAwGC6XS7m5ud77+CcJueDS/3goMTGR4AIAQJA50zAPBucCAICgQXABAABBg+ACAACCBsEFAAAEDYILAAAIGgQXAAAQNEImuDgcDk2bNk1z5861uhQAADBMbMYYY3URvuRyuWS32+V0OlnHBQCAIDHY+3fI9LgAAIDQR3ABAABBg+ACAACCBsEFAAAEjZDbZBEAAAyPpzcd1b76Fl01M0tz8lMsqYEeFwAAMCglexv0u3cOa399q2U1EFwAAMCgdPa4JUkxkdbFB4ILAAAYlK4ejyQpJjLcshoILgAAYFA66HEBAADBwvuoKIIeFwAAEOA6e/uCSzSPigAAQKDr/GCMSyzBBQAABDpmFQEAgKDxYXChxwUAAAQwt8eox20kEVwAAECA6+9tkXhU5BMOh0PTpk3T3LlzrS4FAICQMyC4MB363K1cuVK7d+/W5s2brS4FAICQ09nbN6MoKiJMYWE2y+oImeACAACGT1tXryQpPsq63haJ4AIAAAbB2dEjSbLHRlpaB8EFAACckbOd4AIAAIJEf49LIsEFAAAEOh4VAQCAoEGPCwAACBp761yS6HEBAAABrq2rV+t310uSLpmUZmktBBcAAPCpap0d8hgpITpCF00guAAAgABW09wpScpOirW4EoILAAA4gzpnX3DJSoqxuBKCCwAAOIM6V19wyUwkuAAAgAB3orVLkpQ2KtriSgguAADgDE60dUuSUuKjLK6E4AIAAM7gZDvBBQAABIkTrQQXAAAQBA7Ut2hvXYskggsAAAhwG/Y1eP95/Oh4CyvpQ3ABAACfqH/xuW9cPE5xUREWV0NwAQAAn6K6uUOSlJ9mfW+LRHABAACfwBijA/V941tyAmDVXEmyvs8HAAAEnHpXpx54Za8On2hXXFS4CnOTrS5JUoD2uCxdulTJycn6whe+YHUpAACMSPe/tEdrt1VLkr5/xZSAmFEkBWhwueOOO/Tkk09aXQYAACOSMUb/eL9RknTvkqn6+sXjLK7oQwEZXBYuXKiEhASrywAAYET6t3UVamztVlxUuJYvyLe6nAF8Hlw2btyoq6++WtnZ2bLZbFq3bt0p1zgcDuXn5ysmJkZFRUXatGmTr8sAAABn4WRbt/703lFJ0q2Xjld0RLjFFQ3k8+DS1tamgoICORyO055/5plntGrVKt13333aunWrCgoKtHjxYjU0NJz2egAA4D9v7j8uSYoKD9OdxZMtruZUPp9VtGTJEi1ZsuQTzz/00ENasWKFli9fLkl67LHH9NJLL+nxxx/XvffeO+TP6+rqUldXl/e1y+UaetEAAEDt3b36yct7JElXF2RbXM3p+XWMS3d3t8rKylRcXPxhAWFhKi4uVmlp6Vm955o1a2S3270/ubm5vioXAIARxbGhUg0tXcpMjNE9V06xupzT8mtwaWxslNvtVkZGxoDjGRkZqqur874uLi7WDTfcoJdfflljxoz51FCzevVqOZ1O709VVdWw1Q8AQKiqbu7QbzcekiT9+NrpSk8MjAXnPi4gF6B77bXXBn1tdHS0oqOjh7EaAABC3ys7a9Xt9mjO2GRdMS3jzH/AIn7tcUlLS1N4eLjq6+sHHK+vr1dmZqY/SwEAAB/xSkXfk48rZ2TKZrNZXM0n82twiYqK0uzZs1VSUuI95vF4VFJSovnz5/uzFAAA8IGnNh1V2ZGTCg+zBeyg3H4+f1TU2tqqyspK7+tDhw6pvLxcKSkpysvL06pVq7Rs2TLNmTNH8+bN08MPP6y2tjbvLKOz5XA45HA45Ha7z/UrAAAwYhhj9N9v941tueWifGUE6NiWfjZjjPHlG77xxhtatGjRKceXLVumJ554QpL06KOP6sEHH1RdXZ0KCwv1yCOPqKioyCef73K5ZLfb5XQ6lZiY6JP3BAAgVDW0dGre/SUKs0nl912hxJhIS+oY7P3b5z0uCxcu1Jmy0O23367bb7/d1x8NAACGqKqpXZKUZY+1LLQMRUDuVQQAAPzj6AfBJS8lzuJKBofgAgDACHaosS+45KbEWlzJ4IRMcHE4HJo2bZrmzp1rdSkAAASFtq5e/eL1A5Kk88ckWVvMIIVMcFm5cqV2796tzZs3W10KAABB4alNR9U/LPWyyaOtLWaQQia4AACAwTPG6MUdtZKkK6ZlKJcxLgAAIFA9v7Va5VXNstmk710RmBsqnk5A7lUEAACGh9tj9FzZMf3klT2SpNsXTdSUzASLqxo8ggsAACPIL14/oIdf6xuQOyl9lFYummhxRUMTMo+KmFUEAMCn6+716A+lRyRJ1xVm66lbL1RMZLjFVQ2Nz5f8txpL/gMAcHovbK/Rd57apozEaL1zz+WKCA+c/ovB3r8Dp2IAADBsPrqZ4pfm5AZUaBmK4KwaAAAMyTd+v0XlVc2KDLfpS/PyrC7nrDE4FwCAENbW1avP/+of2lvXIkn6wefOU05ScCzvfzr0uAAAEML+vrvOG1ruuXKqblkwzuKKzg3BBQCAEPVqRZ3+5ZntkqSvzR+rby2cYHFF5y5kggvToQEAGOi3bx2UJEWFh+m2y4I/tEhMhwYAIOT0uj36z5f26Il/HJYkvbv6M8q0x1hb1BkM9v7N4FwAAEJIRbVT//Hibr13qEmS9I2LxwV8aBkKggsAACGiurlDX/p1qdq63ZKkn37+fN0wZ4zFVfkWwQUAgBDwl/Jq3fF0uff1M7deqKLxqdYVNExCZnAuAAAj1dET7QNCyyM3zgrJ0CLR4wIAQFB7eWet7nluhyQpJjJMG+9epPSE0BnT8nEEFwAAgtSuGqfuena72rrdyk+N06NfuSCkQ4sUQsHF4XDI4XDI7XZbXQoAAMOqsqFFv3vnsJ7dckzdbo+mZibohe9crMgg3ThxKFjHBQCAILG7xqX/fvuQntt6zHvs8qnp+tkNBUqOj7KwsnPHOi4AAAQ5Y4z+tqtOz2+tVmVDqw42tnnPTc1M0LcWTtA1Bdmy2WwWVulfBBcAAAJQZ49bK57corcONHqPRYbbdMW0TH31wrGaPyE0Zw2dCcEFAIAAY4zRnU+X660DjQoPs+mbl47X/Ampmp5tV0qQPxI6VwQXAAACzJOlR/TqrjpJ0t2Lp+ibIbJBoi+E/vBjAACCyPvHW3XfX3dJkq4pyNbyBeMsriiwEFwAAAgQO441a8nP3/K+/sFV5ykqglv1R/GoCACAAFDn7NQ1j74jSRqXFq/bLhuvjMTQXkzubIRMcGEBOgBAMDLGaMO+Bv1w3S7vscdvmatxafEWVhW4WIAOAAA/q27u0LGmdlU3d+hP7x1V2ZGT3nP/9aUCLZ01xsLqrMECdAAABKCfv3ZA//Xa/gHHbLa+gbifv2CMLpmUZlFlwYHgAgCAHzjbe/TjF3fp+a3VkqSE6AhNz0nUpPQEfX72GBXmJllbYJAguAAAMMy6et36lz+X6/W9DZKkq2Zm6dGvzBpRS/X7CsEFAIBhtO3oSd373E7tq29RZLhNv7l5jhZNTbe6rKBFcAEAwMec7T3aW+fSs2XHtG5btXo9RvbYSDm+coEuZgzLOSG4AADgI1VN7frzlir96o331ev5cNLupZNH65EvFyopbmTvM+QLBBcAAM6Sx2O0s9qpHcea9fy2am072uw9lxgTocunputrF+VrVm4S41l8hOACAMBZuue5HXq27NiAY/PGpeimojxdU5BNWBkGBBcAAIaoqqldD792QM9t7QsthblJ+uy0DF03K0c5SbEWVxfaCC4AAAxBZUOLrn30HbV1920xs2z+WP342hkWVzVyEFwAABikk23d+vJv3lVbt1tjkmP1f24o0IXjU60ua0QJmb2yHQ6Hpk2bprlz51pdCgAgBFU1teubfyxTY2u3JOnXN88mtFiATRYBADiDhpZOXfFfG9Xc3qPIcJt+9sVCXVOQbXVZIYVNFgEA8IG/lFfrv9bvV3N7j0ZFR+i5b12kKZkJVpc1YhFcAAD4mJNt3Xps4/t692CTtlc1e48/+IXzCS0WI7gAAPCB0vdP6IUdNfrLtmrvrCFJmjM2Wf+5dIamZjIEwWoEFwDAiOb2GG08cFzPlR3TSztr1T/yc2pmgpYvyNcFecmalEEvS6AguAAARqS2rl69tKNWv3rzfR1qbPMeXzw9Q9cW5ujK6ZkKC2Pl20BDcAEAjChuj9Fjb76vR1+vVEdP3+Og6IgwfW5mlr40N5cpzgGO4AIAGDGOnWzXH0qP6NcbD0qSMhKjdcmk0br7yilKT4ixuDoMBsEFABDSjp5oV0WNU28dOK6nNlV5j99z5VTddtl4NkIMMgQXAEDIOXi8Vc9tPabXdjdoX33LgHPj0+L1xbm5uu2yCRZVh3NBcAEABL3jLV366/Ya7TzWrN21Lu2vb/WeCw+zaUaOXTlJMfrC7DG6fGqGhZXiXBFcAABBq97Vqftf2qO/bq855dzFE9N0w5wxWjg5Xfa4SAuqw3AguAAAgobHY7SrxqWNB45r06EmvV3ZKLenb+GVqZkJumpmlqbnJGpall2ZdgbbhiKCCwAg4BhjtGFfg/bUtuhkW7dqnZ063tKlyuOtamrrHnDtxPRRunvxFF0xPdOiauFPBBcAgOU8HqP9DS3acvikdtU4Vfr+CR0+0X7aa0dFR2j+hFQtmJCqwrxkFeYm+bdYWIrgAgDwu+5ej3bXurS9qlm7apx6taJOrs7eAdfYbNIlk0ZramaCMhJjlJkYo6ykGM3MsSsyPMyiymG1kAkuDodDDodDbrf7zBcDAPxuT61Lb+4/rrIjJ7XtaLMaW7sGnI+LCldhbpJm5SVpckaC5uanKDsp1qJqEahsxvRvJxUaXC6X7Ha7nE6nEhPZxRMArNLZ49Y7lY1au61alQ2t2ls3cD2V6IgwnT/GrgvGJqtoXIounTRaEfSkjFiDvX+HTI8LACAwNLZ26Scv7dHa8mp9/H+NF0xM1cLJ6ZqQHq8Lx6cqLorbEIaG/2IAAGfF4zGqqHFqd03fgm9HTrSpxtmpyoYW9bj7EktcVLj+6fwsFZ+XoUkZCRqXFm9x1Qh2BBcAwKA1tnapotqpkj0NWr+7XnWuztNed/4Yu75/xRRdMimNvYDgUwQXAMBpGWNU1dShLUeatPlwk7YdbT5lnMqo6AjNykvSpPQETUiPV3ZSrMalxiufnhUME4ILAGCAls4e7Tjm1C/fqNQ7lScGnLPZpLyUOF00IVVXTMvU/AmpiokMt6hSjEQEFwCA15Olh/Xg3/ap5SNrqhTmJqloXIpm5SVpTn6K0kZFW1ghRjqCCwCMcMdOtuuNfcf1328f0qHGNklSekK0CnOTtHzBOM2fkGpxhcCHCC4AMMJ09ri1YW+D1pVXa3etS1VNHd5zNpt06yXj9f3FU1idFgGJ4AIAIeztA416bU+9DjS0qKWzV62dvap1dqqj58NVxsPDbCrMTdKiKaP15Xl5PApCQCO4AEAI8XiMTrR1q7yqWWu3HdPLO+tOe11OUqwWTR2txdMzVZibpISYSD9XCpwdggsABLF6V6fePtCof7x/QrtrXTrc2DagN0WSrjo/S5dOSlPaqGiNio6QPS5Sk9MTFBbG+ioIPgQXAAhCL+2o1V+3V+uNfcfV1es55fz4tHhdOCFVS2flaG5+igUVAsOD4AIAAcwYo8Mn2nWosW+Twi2HT2pXjVP1rg93Vj4vK1Fz85NVNC5VUzJHKS8lXlERDKxFaCK4AEAAqnV2aNvRZj2/9Zhe29NwyvnwMJu+Mi9P1xZma/bYZJbVx4hBcAEAi3k8RtXNHXr34Am9UlGnAw0tA6YoR4TZNPmDDQovGJus88fYNTk9QfY4BtRi5CG4AICfeTxGx0526O3KRv2/sipVVLvU7R44TiU8zKapmQmalZekr144VlMzEy2qFggsBBcA8ANXZ49e2Vmrdw82acO+BjW39ww4Hxlu09TMRBXmJumK6RmalZesUdH8igY+jr8VAOBjxvQ9+nm1ok776lq0r75FFdVOeczA62bkJOriiaN1dUGWJmcksFItMAgEFwDwgU2HmvTm/gZtPdKs7cea1d7tPuWabHuM/qkgW/MnpOrCcamKjWJXZWCoCC4AMARNbd062tSuelenjp3sUFVTuzYfbtKuGtcp187Msesz56VrUnrfWJXspFgLKgZCC8EFAD5Fg6tTf91eo3cqG7WvrkU1zs5PvPbSyaP12fPSdeH4VOWmxCkmkh4VwNcILgAgqaWzR2VHTqr04AnVOzvV0NKlPbUunfzYIFpJyrLHKD0hWmOS45SdFKOxqfGaPyFVE0aPsqByYGQhuAAYsWqaO7R2W7Ve2F6jyoZW9X589OwHCsbYdcX0TM0bl6IpmQlKZENCwDIBGVxefPFFfe9735PH49E999yjf/7nf7a6JABBzOMx2nLkpPbXt6i8qll761w6drLjlCnJWfYYXTpptCamj1JcdLimZSUqJzlWo0dFszItECACLrj09vZq1apV2rBhg+x2u2bPnq2lS5cqNTXV6tIABJHuXo+e2VKll3fUqvJ4q463dJ32uqJxKfrC7DFaMDFNWfYYAgoQ4AIuuGzatEnTp09XTk6OJGnJkiX6+9//rhtvvNHiygAEOo/H6GR7t96ubNRPX92n6uYPl81PiI5Q0fgUZdljtWBimsamxik3JY5F3oAg4/O/sRs3btSDDz6osrIy1dbWau3atbruuusGXONwOPTggw+qrq5OBQUF+sUvfqF58+ZJkmpqaryhRZJycnJUXV3t6zIBBDmPx2h/Q4sO1LfqtT31Kn3/hE60dcv9kXEqiTERurYwR9cWZmtGjp1ZPkAI8HlwaWtrU0FBgb7+9a/r+uuvP+X8M888o1WrVumxxx5TUVGRHn74YS1evFj79u1Tenq6r8sBEEJqnR16bXe9Xt/boC1HTqqls/e012XbY3Tp5NH60dXTFBdFjwoQSnz+N3rJkiVasmTJJ55/6KGHtGLFCi1fvlyS9Nhjj+mll17S448/rnvvvVfZ2dkDeliqq6u9vTGn09XVpa6uD59du1ynLgIFIDh1dLu15UiTGlu7VOfs0i83VKql68OwEhlu03lZiZqXn6IrpmcqPzVOyfFRLJ0PhDC//q9Id3e3ysrKtHr1au+xsLAwFRcXq7S0VJI0b948VVRUqLq6Wna7Xa+88op++MMffuJ7rlmzRj/+8Y+HvXYAw8cYo8qGVjW0dGnb0ZPaWe3U/vpWHTnRdsr+PmNT4/TFObmaNy5FM3n8A4w4fg0ujY2NcrvdysjIGHA8IyNDe/fu7SsoIkI/+9nPtGjRInk8Ht19992fOqNo9erVWrVqlfe1y+VSbm7u8HwBAOekratXe+tatL++RQc/mOmzt65Fja3damw9/ayfLHuMxo+OV2p8tCZnjNI/XzKesAKMYAH58Peaa67RNddcM6hro6OjFR0dPcwVATgbHo/RuvJqbT58UrtqnNpb16LuXs9pr42JDFNeSpzGpsbrwvGpmpqZoMkZCUobFcUUZQBefg0uaWlpCg8PV319/YDj9fX1yszM9GcpAIbJ8ZYuvbanXu8ePKHS90+o4WPrp6SNitZ5WQmalJ6gtIQojU+LV0ZijCZnJCieqckAzsCvvyWioqI0e/ZslZSUeKdIezwelZSU6Pbbbz+n93Y4HHI4HHK7T91KHsDwaW7v1vvH27S9qlnryqu145hzwPn4qHB9pShPUzITNTUzQdOzE+lBAXDWfB5cWltbVVlZ6X196NAhlZeXKyUlRXl5eVq1apWWLVumOXPmaN68eXr44YfV1tbmnWV0tlauXKmVK1fK5XLJbref69cA8Ama2rr11oHj2na0WfvqWvTeoROnDKA9f4xdCyePVtH4VF2Ql6zYKMakAPANnweXLVu2aNGiRd7X/QNnly1bpieeeEJf+tKXdPz4cf3oRz9SXV2dCgsL9eqrr54yYBeA9foGz7p0vKVLda5OvbC9VntqT11yIMseo7GpcVowIU3Xzx6jnKRYC6oFMBLYjDGn3w41SPX3uDidTiUmJlpdDhA0unrdqj7ZoYoal16tqNWb+46rrfv0j17Py0pUYW6SCsbYNS07UeePSfJvsQBCzmDv34yEA0ag7l6Pthxu0uET7TrR2qV99S36++760874GZ8W790heUL6KN0wZ4zSE2IsqBoAQii4MDgXOD1jzAcLuzVrw94GbTnSpKqmDnW7Tw0pMZFhmpyRoAvyknXFtAydl5Wo5PgoC6oGgNPjUREQgtweo7+UV+v5rdXadKjptCElMSZCs/KSlZkYozHJsZqTn6ILx6cw4weAJXhUBIwwbV29qqh2qs7VqSdLj6jsyMkB5yeMjtclk0brssmjNTF9lLKTYhUeRkgBEFwILkCQqmxo1Z+3VKnsyEk1tHSqztmpHveHHaijoiN0y0X5um5WtrLssSzuBiAk8JsMCAL9M352Vjv1fkOrSg+e0ObDJ0+5LjU+Svlp8ZqZY9fyBfkamxpvQbUAMHwILkAAOt7SpZI99Xp6c5Wqmtp1oq37tNcVjUvR9RfkaGJ6354+uclxCuPxD4AQFjLBhVlFCEYej9G++ha9e/CE6l1damjpVPXJDm0+3HTKarTREWE6LyvRu8/PJZPSNCkjwZrCAcAizCoCLODs6NF3n9qm0vdPnHbGjySNTY3TF+fkauGU0cq2xyopLpIZPwBCFrOKgADT0tmjw43tWr+7Tr8vPSJnR48kKcwmFY1L1dSsBKUnxCg9IVrjR8drVl6yxRUDQOAhuADDxNnRo02HmrTt6Entr29Ryd4GfbR/MzclVg9cf77mjUtRZHiYdYUCQBAhuAA+0tnj1o5jTp1o7dLLFXV6taJ2wPRkSUqOi9SsvGRdW5itq2ZmKYLAAgBDQnABzoGzo0e7a1zafqxZv//HYdU6Owecz0+N04XjU5WbEqdLJqWxGSEAnCOCCzBIHo9RnatTVU3ten1vg8qrmrXlyEm5PzL9JyU+SmNT45SdFKulhTn6zHnpDKgFAB8KmeDCdGgMl8qGFv15yzG9tKNW1c0dp5zPSYrV9OxETc1M0DcuHi97XKQFVQLAyMB0aOA0jOlbX+XpTVX6felh76DaiDCbspNiNX50vJbMyNSsvGRNZi0VADhnTIcGzkKds1M7q536/T8O6+3KRu/xiyem6foLcvS5mVmKiQy3sEIAGNkILhjxOnvcenP/cT309/3aV9/iPW6z9QWWFZeM1yWT0hirAgABgOCCEcvZ3qM1r+zRn7dUeZfXt9mkqZmJykmK1TcuHqf5E1KtLRIAMADBBSNOd69HT5Ye1n+t36+27r7B3HFR4Vo6K0ffu2KKUuKjLK4QAPBJCC4Iec72HlUeb9GJ1m6daOvWs1uqtPVos6S+jQt/cNV5uvnCsTwKAoAgQHBByHF29Oj9462qae7Q/ys7pjf3H9fp5s599cI83XXFVKYvA0AQCZngwjouI1tLZ4/W767XH9894u1N+agse4zSE2OUGh+l3ORYfX72GFaxBYAgxDouCFq9bo8ONrbpbxV1emrTUdV8ZLn99IRo5aXEKcMeoy/OydVlk0dbWCkA4ExYxwUhp8ft0d7aFr3zfqOOnGjTX8trvINrJWl0QrSWzsrRjfPyNC4t3sJKAQDDheCCgLan1qVtR5tVUePU2q3V6ugZ+CgwPipcY5LjdHVBlr52Ub4SYxivAgChjOCCgNPa1as/b67S63sbBqxeK/UFldn5KZqcPkoLJqXpskmjFRbGbCAAGCkILggIG/Y26B/vN2pvXYv+8f6JATsuF4yxa1Zesi6dnKZFU9htGQBGMoILLONs79ErFbX62646bdh3fMC58WnxuunCsSoal6IZOXaLKgQABBqCCyzxyzcq9ejrlWr/yODa62flaEpmgubkp2hWbhKPgAAApyC4wG/cHqPyqma9WlGr3751SJI0NjVO1xZka/6ENPYFAgCcEcEFw8oYo21VzdpV49Jjb7yv6uYO77lvXjpedy2eoojwMAsrBAAEk5AJLqycG3icHT1a9Uy5SvY2eI8lxkSoIDdJCyam6dZLxvM4CAAwJKycC5/p6nXrQH2rdte69Oa+43p1V53cHqPwMJtm5yXroomp+ualExQbFW51qQCAAMPKufAbY4ze2H9c/7a2YsCjIKlv9+Vf3nSBPnNehkXVAQBCCcEFZ+1kW7c2Hjiu32w8qF01Lkl9j4LGpcXrgrHJuiAvWYumpmtUNP+ZAQB8gzsKzsrG/ce14skt6ur1eI999cI8rV5ynuIJKgCAYcIdBoP29oFGPbrhgE60dutAQ6skKTMxRp+dlqGbLszT1EzGFAEAhhfBBWfU1evWbzce1M9LDqjH/ZGl+HOT9NSKIsVF8Z8RAMA/uOPgU63fXa8frN2phpYuSdKSGZm6+cKxSoyN1HlZiQpnOjMAwI8ILhigu9ej/fUtqnV2qqqpXf/rxd2SpIgwm75z+SStXDSBBeMAAJYhuEDOjh69tKNWh0+06eWdtTp2cuCU5ix7jNavuozZQQAAy3EnGsGMMXpqU5Xuf2m32j6y2WF8VLgmpI9SZmKMpmfbdWNRLqEFABAQuBuNUMYY/a8Xd+t37xyWJKXER+magmyNHx2vz18whinNAICAFDJ3J/YqGrz3j7fqtxsP6unNVZL61l/50T9NV1QEY1cAAIGNvYpGkC2Hm/T8tmo9temo+v+tf3vhBN195VRrCwMAjHjsVQSvimqnnt58VH9896j32PzxqfrczEzdVDTWwsoAABgagkuIO9zYpqW/fMe7cNwFeUn68tw83TBnjGw21mABAAQXgksIe/fgCd3x9Db1uI1sNunnX56lq8/PIrAAAIIWwSXEtHT2aHuVU+t31+nJd4/IGCltVJSe+eZ8TRg9yuryAAA4JwSXEPK3XXW6/X+2DthP6KqZWbp/6QwlxUVZWBkAAL5BcAkRvW6PHvlgE8Qse4yKxqXo0smjdV1hjsLYTwgAECIILiHikdcrtavGpaiIMK1buUAZiTFWlwQAgM8RXIKYx2O0dlu1Xqmo1Wt7GiRJP1k6k9ACAAhZBJcg5PEYPb+tWn9494i2VzV7j3/z0vH6/AU51hUGAMAwI7gEmfbuXv3bugo9v7VakhQbGa6vX5yvSyeNVtH4VIurAwBgeBFcgoDHY/TCjhr96d2jqqhxqr3bLZtN+u7lk/SVojweDQEARgyCS4B7ruyY/vere9XQ0uU9lpMUq7uvnKJrC3ksBAAYWQguAazO2anVa3equ9ejmMgwrVw4UZ+dnqGJo0cpIpydnAEAIw/BJUAZY3T/y3vU3evR1MwErf32AsVGhVtdFgAAliK4BCC3x+j//H2fXtheo/Awm/7zuhmEFgAARHAJKMYYvbanQT/7+z7trWuRJN1/3QzNyU+xuDIAAAJDyAQXh8Mhh8Mht9ttdSlnbf3uet36hzJJ0qjoCN2zZKq+PC/P4qoAAAgcNmOMOfNlwcPlcslut8vpdCoxMdHqcgatu9ejG35dqu1VzZqcMUp//EaR0pnmDAAYIQZ7/2ZqSgDweIzueHqbtlc1KyE6Qo/fMpfQAgDAaYTMo6JgZYzRst9t0lsHGhURZtMvvjJLY5LjrC4LAICARHCxkMdj5NhQqbcONEqSvr94ihZOSbe4KgAAAhfBxUL/+9W9+vXGg5Kkz07L0DcvHW9xRQAABDaCiwXcHqP/+9ZBb2i5a/EU3XrpeNlsNosrAwAgsBFc/OztA4267Y9lau3qlSRdPjVdKxdNtLgqAACCA8HFj8qOnNQ3/7BFbd19a80sX5Cve66canFVAAAED4KLn/S4Pbrtj2Vq63YrKS5SL9x+sXJTmD0EAMBQEFz8wBijX5Qc0PGWLiXEROi1VZcpbVS01WUBABB0WIDOD/66vUaPvF4pSfru5ZMILQAAnCWCyzBrbO3ST1/dJ0lauWiCVjDlGQCAs8ajomH0h9LD+o8X96jb7VFiTIS+vZDZQwAAnAuCyzAwxuj7z+7Qc1uPSZJm5tj179dMU3w0zQ0AwLngTupjmw836fvPbteRE+2SpFWfnazvXD6RxeUAAPABgosPHTnRphVPblFze4+iwsP079dM11eK8qwuCwCAkEFw8aHH3nxfze09Gp0QrVfvuESpzB4CAMCnmFXkQ3vrWiRJ91w5ldACAMAwILj4yF/Kq7XtaLOkvsG4AADA9wguPtDj9mj18zslSTlJsRo/Ot7iigAACE0El3Pk8RiteXmv2rvdstmk9asuVWQ4zQoAwHBgcO456HF7dOuTW7Rh33FJ0rUF2YqLokkBABgu3GXPwZ5alze03LV4ipYvyLe2IAAAQlxAPtNYunSpkpOT9YUvfMHqUj5VvatLkjQ1M0ErF02ktwUAgGEWkMHljjvu0JNPPml1GWdU5+qUJOWlxFlcCQAAI0NABpeFCxcqISHB6jLOqOGD4JKRGGNxJQAAjAxDDi4bN27U1VdfrezsbNlsNq1bt+6UaxwOh/Lz8xUTE6OioiJt2rTJF7UGnIYPHhWlJ7DYHAAA/jDk4NLW1qaCggI5HI7Tnn/mmWe0atUq3Xfffdq6dasKCgq0ePFiNTQ0eK8pLCzUjBkzTvmpqak5+29igZPt3ZKklFFRFlcCAMDIMOTRpEuWLNGSJUs+8fxDDz2kFStWaPny5ZKkxx57TC+99JIef/xx3XvvvZKk8vLys6v2NLq6utTV1eV97XK5fPbeZ9Lc0SNJSooluAAA4A8+HePS3d2tsrIyFRcXf/gBYWEqLi5WaWmpLz/Ka82aNbLb7d6f3NzcYfmc03G2fxBc4iL99pkAAIxkPg0ujY2NcrvdysjIGHA8IyNDdXV1g36f4uJi3XDDDXr55Zc1ZsyYTw09q1evltPp9P5UVVWddf1D1dzR96jIHktwAQDAHwJy4ZHXXntt0NdGR0crOtqawbHNH/S4EFwAAPAPn/a4pKWlKTw8XPX19QOO19fXKzMz05cfZbnG1i519Xok8agIAAB/8WlwiYqK0uzZs1VSUuI95vF4VFJSovnz5/vyoyz3ZOkRSX27QY+KDsiOKwAAQs6Q77itra2qrKz0vj506JDKy8uVkpKivLw8rVq1SsuWLdOcOXM0b948Pfzww2pra/POMhouDodDDodDbrd7WD+n35bDTZKkFZeMk81m88tnAgAw0tmMMWYof+CNN97QokWLTjm+bNkyPfHEE5KkRx99VA8++KDq6upUWFioRx55REVFRT4p+ExcLpfsdrucTqcSExOH5TM8HqOCH/9dLV29evm7l2ha9vB8DgAAI8Vg799DDi6Bzh/BpbKhVcUPvamYyDBV/PtiRYQH5M4JAAAEjcHev7njnoXyqmZJ0swcO6EFAAA/4q57Ft6pbJQkzR6bYnElAACMLCETXBwOh6ZNm6a5c+cO6+d093q0YV/fvksLp4we1s8CAAADhUxwWblypXbv3q3NmzcP6+e8e/CEmtt7lJ4Qrbn59LgAAOBPIRNc/KV/R+hJGaMUHsY0aAAA/IngMkS97r5JWOFhNB0AAP7G3XeIej19y/xH0tsCAIDfEVyGqOeDHpeIcIILAAD+FjLBxV+zitye/uASMk0HAEDQCJm7r79mFfW4+x4VRfCoCAAAvwuZ4OIvvf09LgzOBQDA77j7DlH/o6JIxrgAAOB3BJch6n9UxBouAAD4H8FliPrXcYlkcC4AAH4XMndff80q+nCMCz0uAAD4W8gEF3/NKurtf1TEGBcAAPwuZIKLv/T3uEQyqwgAAL/j7jtE/Uv+s3IuAAD+R3AZov7BuYxxAQDA/wguQ/ThXkU0HQAA/sbdd4i8j4rocQEAwO8ILkPEdGgAAKwTMsHFb+u49G+yyKMiAAD8LmTuvv5bx4W9igAAsErIBBd/6X9UFM46LgAA+B133yHqH5xLjwsAAP5HcBki73RoelwAAPA77r5D5O6fVUSPCwAAfkdwGSLvrCKmQwMA4HcElyFi5VwAAKzD3XeI3N7doelxAQDA3wguQ9TzwayicIILAAB+FzLBxX8r5/KoCAAAq4TM3ddfK+d6HxUxqwgAAL8LmeDiLz1uHhUBAGAVgssQ9Xp7XGg6AAD8jbvvEPWwjgsAAJYhuAyRd+VclvwHAMDvuPsO0YeziuhxAQDA3wguQ9S/jgvBBQAA/yO4DIHHY2T6Olx4VAQAgAW4+w5Bf2+LRI8LAABWILgMQf/4FkmKpMcFAAC/4+47BP1ruEj0uAAAYIWQCS7+2Kuo1/2RR0Ws4wIAgN+FTHDxx15F/T0u4WE22WwEFwAA/C1kgos/9HoXnyO0AABgBYLLEPSy3D8AAJYiuAxBj3fVXJoNAAArcAceArd3Z2h6XAAAsALBZQj6d4YO51ERAACWILgMQS87QwMAYCnuwENw7GS7JB4VAQBgFYLLELyys06SND3bbnElAACMTASXITjU2CZJWjorx+JKAAAYmQguQ1Dr7JAk5STHWlwJAAAjE8FlkDq63TrZ3iNJyk4iuAAAYAWCyyAdb+mSJMVEhikxJsLiagAAGJkILoPU7XZLkmIiw9lgEQAAixBcBumDtefYpwgAAAsRXAap18OquQAAWI3gMkhuVs0FAMByIXMXdjgcmjZtmubOnTss79+/3D89LgAAWCdkgsvKlSu1e/dubd68eVje/8MeF4ILAABWCZngMtx63fS4AABgNYLLILl5VAQAgOUILoPUP6sogp2hAQCwDMFlkD7scaHJAACwCnfhQeplcC4AAJYjuAwSY1wAALAewWWQ6HEBAMB6BJdBcrPkPwAAliO4DFL/Oi70uAAAYB2CyyAxqwgAAOtxFx4kxrgAAGA9gssgeXtcWIAOAADLEFwGiR4XAACsR3AZJGYVAQBgPYLLINHjAgCA9Qgug+R2M6sIAACrcRceJHpcAACwHsFlkNirCAAA6xFcBqmX4AIAgOUILoPUP6uIR0UAAFiH4DJI9LgAAGA9gssguRmcCwCA5QIuuFRVVWnhwoWaNm2azj//fD377LNWlyTpoz0uAddkAACMGBFWF/BxERERevjhh1VYWKi6ujrNnj1bn/vc5xQfH29pXQsnj1ZyXKTm5idbWgcAACNZwAWXrKwsZWVlSZIyMzOVlpampqYmy4PLFdMzdcX0TEtrAABgpBvyc4+NGzfq6quvVnZ2tmw2m9atW3fKNQ6HQ/n5+YqJiVFRUZE2bdp0VsWVlZXJ7XYrNzf3rP48AAAILUMOLm1tbSooKJDD4Tjt+WeeeUarVq3Sfffdp61bt6qgoECLFy9WQ0OD95rCwkLNmDHjlJ+amhrvNU1NTfra176m3/zmN2fxtQAAQCiyGWPMWf9hm01r167Vdddd5z1WVFSkuXPn6tFHH5UkeTwe5ebm6jvf+Y7uvffeQb1vV1eXPvvZz2rFihW6+eabz3htV1eX97XL5VJubq6cTqcSExOH/qUAAIDfuVwu2e32M96/fTpFpru7W2VlZSouLv7wA8LCVFxcrNLS0kG9hzFGt9xyiy6//PIzhhZJWrNmjex2u/eHx0oAAIQunwaXxsZGud1uZWRkDDiekZGhurq6Qb3HO++8o2eeeUbr1q1TYWGhCgsLtXPnzk+8fvXq1XI6nd6fqqqqc/oOAAAgcAXcrKKLL75Yng+W1x+M6OhoRUdHD2NFAAAgUPi0xyUtLU3h4eGqr68fcLy+vl6ZmUwlBgAA58anwSUqKkqzZ89WSUmJ95jH41FJSYnmz5/vy48CAAAj0JAfFbW2tqqystL7+tChQyovL1dKSory8vK0atUqLVu2THPmzNG8efP08MMPq62tTcuXL/dp4R/ncDjkcDjkdruH9XMAAIB1hjwd+o033tCiRYtOOb5s2TI98cQTkqRHH31UDz74oOrq6lRYWKhHHnlERUVFPin4TAY7nQoAAASOwd6/z2kdl0BEcAEAIPhYso4LAADAcCK4AACAoBFw67icrf7Bub29vZL6upwAAEBw6L9vn2kES8iNcTl27BjL/gMAEKSqqqo0ZsyYTzwfcsHF4/GopqZGCQkJstlsPnvf/s0bq6qqGPQ7zGhr/6Cd/YN29h/a2j+Gq52NMWppaVF2drbCwj55JEvIPCrqFxYW9qlJ7VwlJibyF8JPaGv/oJ39g3b2H9raP4ajne12+xmvYXAuAAAIGgQXAAAQNAgugxQdHa377ruPnaj9gLb2D9rZP2hn/6Gt/cPqdg65wbkAACB00eMCAACCBsEFAAAEDYILAAAIGgQXAAAQNAgug+RwOJSfn6+YmBgVFRVp06ZNVpcUVNasWaO5c+cqISFB6enpuu6667Rv374B13R2dmrlypVKTU3VqFGj9PnPf1719fUDrjl69KiuuuoqxcXFKT09XXfddZd3fyqc6oEHHpDNZtOdd97pPUY7+0Z1dbW++tWvKjU1VbGxsZo5c6a2bNniPW+M0Y9+9CNlZWUpNjZWxcXFOnDgwID3aGpq0k033aTExEQlJSXpG9/4hlpbW/39VQKW2+3WD3/4Q40bN06xsbGaMGGC/uM//mPAXja089nZuHGjrr76amVnZ8tms2ndunUDzvuqXXfs2KFLLrlEMTExys3N1U9/+tNzL97gjJ5++mkTFRVlHn/8cbNr1y6zYsUKk5SUZOrr660uLWgsXrzY/O53vzMVFRWmvLzcfO5znzN5eXmmtbXVe81tt91mcnNzTUlJidmyZYu58MILzUUXXeQ939vba2bMmGGKi4vNtm3bzMsvv2zS0tLM6tWrrfhKAW/Tpk0mPz/fnH/++eaOO+7wHqedz11TU5MZO3asueWWW8x7771nDh48aP72t7+ZyspK7zUPPPCAsdvtZt26dWb79u3mmmuuMePGjTMdHR3ea6688kpTUFBg3n33XfPWW2+ZiRMnmhtvvNGKrxSQ7r//fpOammpefPFFc+jQIfPss8+aUaNGmZ///Ofea2jns/Pyyy+bH/zgB+b55583kszatWsHnPdFuzqdTpORkWFuuukmU1FRYZ566ikTGxtrfv3rX59T7QSXQZg3b55ZuXKl97Xb7TbZ2dlmzZo1FlYV3BoaGowk8+abbxpjjGlubjaRkZHm2Wef9V6zZ88eI8mUlpYaY/r+ooWFhZm6ujrvNb/61a9MYmKi6erq8u8XCHAtLS1m0qRJZv369eayyy7zBhfa2Tfuuecec/HFF3/ieY/HYzIzM82DDz7oPdbc3Gyio6PNU089ZYwxZvfu3UaS2bx5s/eaV155xdhsNlNdXT18xQeRq666ynz9618fcOz66683N910kzGGdvaVjwcXX7XrL3/5S5OcnDzg98Y999xjpkyZck718qjoDLq7u1VWVqbi4mLvsbCwMBUXF6u0tNTCyoKb0+mUJKWkpEiSysrK1NPTM6Cdp06dqry8PG87l5aWaubMmcrIyPBes3jxYrlcLu3atcuP1Qe+lStX6qqrrhrQnhLt7Ct//etfNWfOHN1www1KT0/XrFmz9Nvf/tZ7/tChQ6qrqxvQzna7XUVFRQPaOSkpSXPmzPFeU1xcrLCwML333nv++zIB7KKLLlJJSYn2798vSdq+fbvefvttLVmyRBLtPFx81a6lpaW69NJLFRUV5b1m8eLF2rdvn06ePHnW9YXcJou+1tjYKLfbPeCXuCRlZGRo7969FlUV3Dwej+68804tWLBAM2bMkCTV1dUpKipKSUlJA67NyMhQXV2d95rT/XvoP4c+Tz/9tLZu3arNmzefco529o2DBw/qV7/6lVatWqV//dd/1ebNm/Xd735XUVFRWrZsmbedTteOH23n9PT0AecjIiKUkpJCO3/g3nvvlcvl0tSpUxUeHi632637779fN910kyTRzsPEV+1aV1encePGnfIe/eeSk5PPqj6CC/xu5cqVqqio0Ntvv211KSGnqqpKd9xxh9avX6+YmBirywlZHo9Hc+bM0U9+8hNJ0qxZs1RRUaHHHntMy5Yts7i60PHnP/9Zf/rTn/Q///M/mj59usrLy3XnnXcqOzubdh7BeFR0BmlpaQoPDz9l1kV9fb0yMzMtqip43X777XrxxRe1YcMGjRkzxns8MzNT3d3dam5uHnD9R9s5MzPztP8e+s+h71FQQ0ODLrjgAkVERCgiIkJvvvmmHnnkEUVERCgjI4N29oGsrCxNmzZtwLHzzjtPR48elfRhO33a743MzEw1NDQMON/b26umpiba+QN33XWX7r33Xn35y1/WzJkzdfPNN+tf/uVftGbNGkm083DxVbsO1+8SgssZREVFafbs2SopKfEe83g8Kikp0fz58y2sLLgYY3T77bdr7dq1ev3110/pPpw9e7YiIyMHtPO+fft09OhRbzvPnz9fO3fuHPCXZf369UpMTDzlJjJSfeYzn9HOnTtVXl7u/ZkzZ45uuukm7z/TzuduwYIFp0zn379/v8aOHStJGjdunDIzMwe0s8vl0nvvvTegnZubm1VWVua95vXXX5fH41FRUZEfvkXga29vV1jYwNtUeHi4PB6PJNp5uPiqXefPn6+NGzeqp6fHe8369es1ZcqUs35MJInp0IPx9NNPm+joaPPEE0+Y3bt3m1tvvdUkJSUNmHWBT/etb33L2O1288Ybb5ja2lrvT3t7u/ea2267zeTl5ZnXX3/dbNmyxcyfP9/Mnz/fe75/mu4VV1xhysvLzauvvmpGjx7NNN0z+OisImNoZ1/YtGmTiYiIMPfff785cOCA+dOf/mTi4uLMH//4R+81DzzwgElKSjJ/+ctfzI4dO8y111572umks2bNMu+99555++23zaRJk0b8NN2PWrZsmcnJyfFOh37++edNWlqaufvuu73X0M5np6WlxWzbts1s27bNSDIPPfSQ2bZtmzly5Igxxjft2tzcbDIyMszNN99sKioqzNNPP23i4uKYDu0vv/jFL0xeXp6Jiooy8+bNM++++67VJQUVSaf9+d3vfue9pqOjw3z72982ycnJJi4uzixdutTU1tYOeJ/Dhw+bJUuWmNjYWJOWlma+973vmZ6eHj9/m+Dy8eBCO/vGCy+8YGbMmGGio6PN1KlTzW9+85sB5z0ej/nhD39oMjIyTHR0tPnMZz5j9u3bN+CaEydOmBtvvNGMGjXKJCYmmuXLl5uWlhZ/fo2A5nK5zB133GHy8vJMTEyMGT9+vPnBD34wYHot7Xx2NmzYcNrfycuWLTPG+K5dt2/fbi6++GITHR1tcnJyzAMPPHDOtduM+cgShAAAAAGMMS4AACBoEFwAAEDQILgAAICgQXABAABBg+ACAACCBsEFAAAEDYILAAAIGgQXAAAQNAguAAAgaBBcAABA0CC4AACAoEFwAQAAQeP/AzIw8h1Z5htYAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.sort(A_diag))\n",
    "plt.yscale('log')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Spherical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.6938282186236218 0.09976510158544534 1.6949800788650446\n",
      "0.769348754034222 0.022423761619716504 -0.26263728867901676 0.029233769165452724\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(np.float64(-262.6372886790168),\n",
       " np.float64(29.233769165452724),\n",
       " np.float64(0.2923376916545272))"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_samples = 10_000\n",
    "vecs = np.random.normal(size=(n_samples, d))\n",
    "vecs /= np.linalg.norm(vecs, axis=1, keepdims=True)\n",
    "unscaled_f = np.sum(A_diag * vecs**2, axis=1)\n",
    "print(unscaled_f.mean(), unscaled_f.std(), A_diag.mean())\n",
    "r_squared = 1 / unscaled_f\n",
    "r = np.sqrt(r_squared)\n",
    "print(r.mean(), r.std(), np.log(r).mean(), np.log(r).std())\n",
    "logvol = d * np.log(r)\n",
    "logvol.mean(), logvol.std(), logvol.std() / np.sqrt(n_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Axis-aligned"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(np.float64(1.3215850330529575), np.float64(517.390097144195))"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "logvol = -d/2 * np.log(A_diag)\n",
    "logvol.mean(), logvol.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.float64(-263.83549394068416)"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "-d/2 * np.log(np.mean(A_diag))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.float64(1.3215850330529575)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "-d/2 * np.log(A_diag).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-250.0"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# large d limit of naive spherical estimator\n",
    "-d/2 * (1/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.0"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# large d limit of axis-aligned estimator\n",
    "-d/2 * 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "jax311",
   "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}