10.3_Introduction_to_Frames.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Finite Frames and RIC-POVM's"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maturing out of our reflexive habit of working with orthonormal bases, we can instead work with an overcomplete basis or \"frame,\" which can have as many elements as we like. In fact, with probability 1, any randomly chosen set of vectors ($n \\ge d$) will form a frame.\n",
    "\n",
    "Let's just think about $d=2$ real vector space. Why not a frame made of 5 vectors?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 449,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.1443,  0.0924,  0.7253,  0.4074, -0.2219],\n",
       "       [ 1.2268,  0.234 ,  1.0096,  1.3012, -0.2936]])"
      ]
     },
     "execution_count": 449,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy as sc\n",
    "import scipy.linalg\n",
    "np.set_printoptions(precision=4, suppress=True)\n",
    "\n",
    "d, n = 2, 5\n",
    "R = np.random.randn(d, n); R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can form the Gram matrix, aka the matrix of inner products:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 450,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.5258,  0.2737,  1.1339,  1.5375, -0.3282],\n",
       "       [ 0.2737,  0.0633,  0.3032,  0.3421, -0.0892],\n",
       "       [ 1.1339,  0.3032,  1.5454,  1.6092, -0.4574],\n",
       "       [ 1.5375,  0.3421,  1.6092,  1.8591, -0.4724],\n",
       "       [-0.3282, -0.0892, -0.4574, -0.4724,  0.1355]])"
      ]
     },
     "execution_count": 450,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G = R.T @ R; G"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One condition for everything being kosher is that the Gram matrix ought to be rank-$d$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 451,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4.7074, 0.4216, 0.    , 0.    , 0.    ])"
      ]
     },
     "execution_count": 451,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.round(np.linalg.svd(G)[1], decimals=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We should also check out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 452,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.7706, 1.1721],\n",
       "       [1.1721, 4.3584]])"
      ]
     },
     "execution_count": 452,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R @ R.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll have a frame if the latter has non-0 determinant aka is full rank. This means that the rows of $R$ are linearly independent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 453,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.9847520214675856"
      ]
     },
     "execution_count": 453,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.det(R @ R.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By the way, the technical definition of a frame is a set of vectors $\\{v\\}$, for which there are constants $a, b > 0$ such that for any vector $u$:\n",
    "\n",
    "$$ a ||u||^2 \\leq \\sum_{i} |\\langle u|v_{i}\\rangle|^2 \\leq b||u||^2 $$\n",
    "\n",
    "The $a$ and $b$ are called frame bounds, and this definition is always satisfied if the aforementioned conditions are true. We won't be using this definition too much. \n",
    "\n",
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis and Synthesis\n",
    "\n",
    "Suppose we have some 2D point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 454,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.3576, 0.5659])"
      ]
     },
     "execution_count": 454,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "point = np.random.randn(2); point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can expand it in terms of our frame vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 455,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.6426,  0.1654,  0.8307,  0.882 , -0.2455])"
      ]
     },
     "execution_count": 455,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = R.T @ point; c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " `R.T` is often called the \"analysis\" operator.\n",
    " \n",
    "How do we recover the original point, though? \n",
    "\n",
    "We use the dual frame, obtained by taking the pseudo-inverse of `R`. \n",
    "\n",
    "The inner products with the frame vectors turn out to be the weights on the dual basis vectors, which via summation recover the original point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 456,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.0413,  0.5615],\n",
       "       [ 0.0647,  0.0363],\n",
       "       [ 0.9964, -0.0363],\n",
       "       [ 0.1262,  0.2646],\n",
       "       [-0.314 ,  0.0171]])"
      ]
     },
     "execution_count": 456,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D = np.linalg.pinv(R); D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 457,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.3576, 0.5659])"
      ]
     },
     "execution_count": 457,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum([c[i]*D[i] for i in range(n)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 458,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.3576, 0.5659])"
      ]
     },
     "execution_count": 458,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D.T @ c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`D.T` is often called the \"synthesis\" operator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It works because:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 459,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.],\n",
       "       [-0.,  1.]])"
      ]
     },
     "execution_count": 459,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D.T @ R.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Conversely, we could have expanded in the dual basis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 460,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.0547,  0.0437,  0.3358,  0.1949, -0.1026])"
      ]
     },
     "execution_count": 460,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = D @ point; c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 461,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.3576, 0.5659])"
      ]
     },
     "execution_count": 461,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R @ c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way of thinking about this \"recovery\" procedure: instead, we could take the pseudo-inverse of the Gram matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 462,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.3997, -0.047 , -1.058 ,  0.0172,  0.3365],\n",
       "       [-0.047 ,  0.0055,  0.0631,  0.0178, -0.0197],\n",
       "       [-1.058 ,  0.0631,  0.9942,  0.1161, -0.3135],\n",
       "       [ 0.0172,  0.0178,  0.1161,  0.0859, -0.0351],\n",
       "       [ 0.3365, -0.0197, -0.3135, -0.0351,  0.0989]])"
      ]
     },
     "execution_count": 462,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ginv = np.linalg.pinv(G); Ginv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we have our frame coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 463,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.6426,  0.1654,  0.8307,  0.882 , -0.2455])"
      ]
     },
     "execution_count": 463,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = R.T @ point; c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Acting with $G^{-1}$ will give us the coordinates in the dual basis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 464,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.0547,  0.0437,  0.3358,  0.1949, -0.1026])"
      ]
     },
     "execution_count": 464,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ginv @ c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 465,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.3576, 0.5659])"
      ]
     },
     "execution_count": 465,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R @ Ginv @ c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In other words, using $G^{-1}$, we can avoid explicitly invoking the dual basis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "## Tight Frames\n",
    "\n",
    "We might consider restrictions on our frames. Different kinds of frames are useful for different purposes. \n",
    "\n",
    "For example, everyone loves tight frames. We can get one via a generalization of a polar decomposition for rectangular matrices. Starting with a random $d \\times n$ matrix $R$, we decompose as $R = UH$, where the rows of $U$ are orthonormal, and $H$ is Hermitian. Then we throw out the $H$ part."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 466,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.1471,  0.2691,  0.0433,  0.1094,  0.9445],\n",
       "       [ 0.4285,  0.6161, -0.5379, -0.3818, -0.04  ]])"
      ]
     },
     "execution_count": 466,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R = np.random.randn(d, n)\n",
    "R = sc.linalg.polar(R)[0]; R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 467,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.2053,  0.2244, -0.2369, -0.1797, -0.1561],\n",
       "       [ 0.2244,  0.452 , -0.3198, -0.2058,  0.2296],\n",
       "       [-0.2369, -0.3198,  0.2913,  0.2101,  0.0624],\n",
       "       [-0.1797, -0.2058,  0.2101,  0.1578,  0.1186],\n",
       "       [-0.1561,  0.2296,  0.0624,  0.1186,  0.8937]])"
      ]
     },
     "execution_count": 467,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R.T @ R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The point is that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 468,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1., -0.],\n",
       "       [-0.,  1.]])"
      ]
     },
     "execution_count": 468,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R @ R.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the frame is the same as its dual! In other words, the analysis operator is `R.T` and the synthesis operator is `R`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 469,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.3576, 0.5659])"
      ]
     },
     "execution_count": 469,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 470,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.1899,  0.4449, -0.2889, -0.177 ,  0.3152])"
      ]
     },
     "execution_count": 470,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = R.T @ point; c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 471,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.3576, 0.5659])"
      ]
     },
     "execution_count": 471,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R @ c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More generally, any frame for which $R R^{\\dagger}$ is proportional to the identity is a \"tight frame,\" so that the dual frame is proportional to the original frame. In terms of the \"frame bounds\" from before, a tight frame means that $a=b$. And when $a=b=1$, this is known as a \"Parseval frame.\"\n",
    "\n",
    "Also, for a tight frame, the Gram matrix is a rank-$d$ projector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 472,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.,  1.,  1., -0., -0.])"
      ]
     },
     "execution_count": 472,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.eig(R.T @ R)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unit Norm Tight Frames\n",
    "\n",
    "We might also be interested not just in tight frames, but in \"unit norm tight frames\". As the name suggests, each vector has unit norm. \n",
    "\n",
    "One easy way to get one is to alternate between \"tightening\" and \"normalizing\" our vectors until we're satisfied."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 473,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_funtf(d, n, rtol=1e-15, atol=1e-15):\n",
    "    R = np.random.randn(d, n)\n",
    "    done = False\n",
    "    while not (np.allclose(R @ R.T, (n/d)*np.eye(d), rtol=rtol, atol=atol) and\\\n",
    "               np.allclose(np.linalg.norm(R, axis=0), np.ones(n), rtol=rtol, atol=atol)):\n",
    "        R = sc.linalg.polar(R)[0]\n",
    "        R = np.array([state/np.linalg.norm(state) for state in R.T]).T\n",
    "    return R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 474,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.0479, -0.9568, -0.0103,  0.9641,  0.8079],\n",
       "       [ 0.9989,  0.2906, -0.9999, -0.2657,  0.5893]])"
      ]
     },
     "execution_count": 474,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R = random_funtf(2, 5); R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 475,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.    ,  0.2444, -0.9993, -0.2192,  0.6274],\n",
       "       [ 0.2444,  1.    , -0.2808, -0.9997, -0.6017],\n",
       "       [-0.9993, -0.2808,  1.    ,  0.2558, -0.5976],\n",
       "       [-0.2192, -0.9997,  0.2558,  1.    ,  0.6222],\n",
       "       [ 0.6274, -0.6017, -0.5976,  0.6222,  1.    ]])"
      ]
     },
     "execution_count": 475,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R.T @ R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 476,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2.5, -0. ],\n",
       "       [-0. ,  2.5]])"
      ]
     },
     "execution_count": 476,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R @ R.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that we've had to demand $RR^{\\dagger}$ to be $\\frac{n}{d} I$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 477,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.3576, 0.5659])"
      ]
     },
     "execution_count": 477,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 478,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.5824, -0.1777, -0.5695,  0.1944,  0.6224])"
      ]
     },
     "execution_count": 478,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = R.T @ point; c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 479,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.3576, 0.5659])"
      ]
     },
     "execution_count": 479,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(d/n)*R @ c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equiangular Frames\n",
    "\n",
    "As I say, there are many kinds of frames that people are interested in, which have various geometric and/or information theoretic properties.\n",
    "\n",
    "One of the classics of the genre are the so-called \"equiangular frames,\" which are made up of vectors which form equal angles with each other (or more properly, the absolute values of the inner products are all the same).\n",
    "\n",
    "These are great for many reasons, not least because they provide very resilient representations of data even if you should lose a few components.\n",
    "\n",
    "There are, however, geometrical restrictions on the maximum number of equiangular lines you can have in a given dimension. For example, in $d=2$, you can have a maximum of $3$ equiangular lines. \n",
    "\n",
    "The $3$ are easy to construct:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 480,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.   , -0.5  , -0.5  ],\n",
       "       [ 0.   ,  0.866, -0.866]])"
      ]
     },
     "execution_count": 480,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R = np.array([[r.real, r.imag] for r in [np.exp(2j*np.pi*i/3) for i in range(3)]]).T; R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 481,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1. , -0.5, -0.5],\n",
       "       [-0.5,  1. , -0.5],\n",
       "       [-0.5, -0.5,  1. ]])"
      ]
     },
     "execution_count": 481,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R.T @ R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 482,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.5, 0. ],\n",
       "       [0. , 1.5]])"
      ]
     },
     "execution_count": 482,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R @ R.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 3 lines form the diagonals of a hexagon. \n",
    "\n",
    "<img src=\"img/hexagon.jpg\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "In three dimensions, an equiangular frame is provided by the icosahedron, which has 20 faces and 12 vertices: its six diagonals provide the frame vectors.\n",
    "\n",
    "<img width=300 src=\"img/icosahedron.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 483,
   "metadata": {},
   "outputs": [],
   "source": [
    "def circular_shifts(v): \n",
    "    \"\"\"\n",
    "    Returns all shifted versions of a vector, e.g.:\n",
    "        [1,2,3]\n",
    "        [3,1,2]\n",
    "        [2,3,1]\n",
    "    \"\"\"\n",
    "    shifts = [v]\n",
    "    for i in range(len(v)-1):\n",
    "        u = shifts[-1][:]\n",
    "        u.insert(0, u.pop()) \n",
    "        shifts.append(u)\n",
    "    return shifts\n",
    "\n",
    "def icosahedron_frame():\n",
    "    \"\"\"\n",
    "    Returns equiangular frame based on the icosahedron in three dimensions.\n",
    "    \"\"\"\n",
    "    phi = (1+np.sqrt(5))/2 # golden ratio\n",
    "    # 12 vertices provided by circular shifts of the four vectors\n",
    "    vertices = [np.array(v) for v in \n",
    "                   circular_shifts([0, 1, phi]) + \\\n",
    "                   circular_shifts([0, -1, -phi]) + \\\n",
    "                   circular_shifts([0, 1, -phi]) + \\\n",
    "                   circular_shifts([0, -1, phi])]\n",
    "    # throw out the antipodal vectors\n",
    "    keep = []\n",
    "    for i, a in enumerate(vertices):\n",
    "        for j, b in enumerate(vertices):\n",
    "            if i != j and np.allclose(a, -b) and j not in keep:\n",
    "                keep.append(i)\n",
    "    # normalize\n",
    "    return np.array([e/np.linalg.norm(e) for i, e in enumerate(vertices) if i in keep]).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 484,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.    ,  0.8507,  0.5257,  0.    , -0.8507,  0.5257],\n",
       "       [ 0.5257,  0.    ,  0.8507,  0.5257,  0.    , -0.8507],\n",
       "       [ 0.8507,  0.5257,  0.    , -0.8507,  0.5257,  0.    ]])"
      ]
     },
     "execution_count": 484,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R = icosahedron_frame(); R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 485,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.    ,  0.4472,  0.4472, -0.4472,  0.4472, -0.4472],\n",
       "       [ 0.4472,  1.    ,  0.4472, -0.4472, -0.4472,  0.4472],\n",
       "       [ 0.4472,  0.4472,  1.    ,  0.4472, -0.4472, -0.4472],\n",
       "       [-0.4472, -0.4472,  0.4472,  1.    , -0.4472, -0.4472],\n",
       "       [ 0.4472, -0.4472, -0.4472, -0.4472,  1.    , -0.4472],\n",
       "       [-0.4472,  0.4472, -0.4472, -0.4472, -0.4472,  1.    ]])"
      ]
     },
     "execution_count": 485,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R.T @ R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 486,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2., 0., 0.],\n",
       "       [0., 2., 0.],\n",
       "       [0., 0., 2.]])"
      ]
     },
     "execution_count": 486,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R @ R.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check out <a href=\"examples/frame_viz.py\">frame_viz.py</a> to play around with different kinds of 2D frames! \n",
    "\n",
    "You'll be able to drag around a point in 2D. Press `c` to show standard Cartesian axes, and the projections of the point onto them. Press `f` to show the frame vectors and the projections onto them. Press `d` to show the dual vectors and the projections onto them. Press `q` to show the reconstruction of the point by a weighted sum of frame vectors (provided by the inner products with the dual vectors). Press `p` to show the reconstruction of the point by a weighted sum of dual vectors (provided by the inner products with the frame vectors).\n",
    "\n",
    "Note that to get the projections right, we need to divide frame (dual) coordinates by the norm squared of the corresponding frame (dual) vectors. \n",
    "\n",
    "Because:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 487,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.6091,  0.6271, -0.4837],\n",
       "       [-0.6152, -1.0213, -1.5933]])"
      ]
     },
     "execution_count": 487,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R = np.random.randn(2, 3); R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 488,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.3576, 0.5659])"
      ]
     },
     "execution_count": 488,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to project the point onto the three frame vectors, we can use the projectors formed from the outer products $rr^{T}$, where $r$ is unit norm. This amounts to dividing by $||r||^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 489,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.5008,  0.1915],\n",
       "       [-0.1544,  0.2515],\n",
       "       [ 0.1875,  0.6175]])"
      ]
     },
     "execution_count": 489,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array([np.outer(r/np.linalg.norm(r), r/np.linalg.norm(r))@point for r in R.T])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Meanwhile, we have our inner products:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 490,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.9236, -0.3536, -1.0746])"
      ]
     },
     "execution_count": 490,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = R.T @ point; c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To recover the projections, we multiply each of the inner products by the vector, divided by the norm squared."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 491,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.5008,  0.1915],\n",
       "       [-0.1544,  0.2515],\n",
       "       [ 0.1875,  0.6175]])"
      ]
     },
     "execution_count": 491,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array([c[i]*R.T[i]/np.linalg.norm(R.T[i])**2 for i in range(3)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, in the case of 2D, we have another way of visualizing things. Press `3` to check it out.\n",
    "\n",
    "We can interpret the three inner products with the frame vectors as a point in 3D. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 492,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.9236, -0.3536, -1.0746])"
      ]
     },
     "execution_count": 492,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Moreover, `D` will be a $3 \\times 2$ matrix: we can interpret its 2 columns as specifying two vectors in 3D."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 493,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.4949, -0.0153],\n",
       "       [ 0.3159, -0.3472],\n",
       "       [-0.0114, -0.3991]])"
      ]
     },
     "execution_count": 493,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D = np.linalg.pinv(R); D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can project our 3D point onto these two vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 494,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([-0.5132,  0.3276, -0.0118]), array([-0.031 , -0.7015, -0.8063])]"
      ]
     },
     "execution_count": 494,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[np.outer(D.T[i], D.T[i])/np.linalg.norm(D.T[i])**2 @ c for i in range(2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But this is just the same as weighting each of the two vectors by the Cartesian coordinates of the original point!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 495,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([-0.5132,  0.3276, -0.0118]), array([-0.031 , -0.7015, -0.8063])]"
      ]
     },
     "execution_count": 495,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[point[i]*D.T[i]/np.linalg.norm(D.T[i])**2 for i in range(2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In other words, the three inner products with the frame vector specify a point in 3D. The dual basis specifies two linearly independent axes in 3D. The projection of the 3D point onto these two axes gives two coordinates: these are proportional to the original 2D point in Cartesian coordinates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RIC-POVM's\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now consider that $3$ is also the dimension of the space of $2 \\times 2$ symmetric matrices. \n",
    "\n",
    "Indeed, we might be interested in this fact if we wanted to do real vector space quantum mechanics. \n",
    "\n",
    "In real quantum mechanics, pure states are given by 2D vectors $\\mid \\psi \\rangle = \\begin{pmatrix}x \\\\ y \\end{pmatrix}$ such that $x^2 + y^2 = 1$, aka points on the circle; or better yet: rank-1 density matrices: $\\mid \\psi \\rangle \\langle \\psi \\mid = \\begin{pmatrix}x^2 & xy \\\\ xy & y^2 \\end{pmatrix}$. \n",
    "\n",
    "More generally, a density matrix might be a higher rank statistical mixture. Any symmetric, positive semi-definite matrix with unit trace will do.\n",
    "\n",
    "Now here's the big reveal: it turns out that tight frames are just rank-1 POVM's in disguise. \n",
    "\n",
    "Recall that a POVM, or positive operator valued measure, is a set of positive semi-definite matrices that sum to the identity: $\\sum_{i} E_{i} = I$. Each \"element\" represents an outcome to a generalized quantum measurement that concretely involves entangling your quantum system to an apparatus (with a dimensionality equal to the number of POVM elements) and doing a standard (Von Neumann) measurement on the latter. In the case of a rank-1 POVM, the outcome states of the original sytem will correspond to the normalized POVM elements.\n",
    "\n",
    "A POVM can have any number of elements (well, I suppose at least two!): but we're interested in \"informationally complete\" POVM's, which have as many elements as the dimensions of the relevant operator space, aka the space of density matrices, which for real vector space is $\\frac{d(d+1)}{2}$. (In complex quantum mechanics, we use not symmetric matrices, but Hermitian matrices, which have $d^2$ parameters.) \n",
    "\n",
    "If we have enough elements to span the operator space, then the vector of probabilities for each outcome suffice to completely determine the quantum state. And so, among other things, we can do quantum tomography, treating pure states and mixed states on an equal footing, by repeating one measurement over and over again on many copies of a system, instead of having to do different, complementary measurements.\n",
    "\n",
    "We can upgrade a tight frame into a rank-1 POVM by taking the outer products of the vectors with themselves. (And by the way, it's not at all necessary that POVM elements be rank-1: those we'll be considering here, however, will be.) The POVM condition, as you can see, is satisfied since $RR^{\\dagger}$ amounts to $\\sum_{i} |r\\rangle\\langle r|$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 496,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0.],\n",
       "       [0., 1.]])"
      ]
     },
     "execution_count": 496,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d, n = 2, 3\n",
    "R = sc.linalg.polar(np.random.randn(d, n))[0]\n",
    "R @ R.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 497,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0.],\n",
       "       [0., 1.]])"
      ]
     },
     "execution_count": 497,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E = [np.outer(r, r) for r in R.T]\n",
    "sum(E)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given a random density matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 498,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.123 , -0.3285],\n",
       "       [-0.3285,  0.877 ]])"
      ]
     },
     "execution_count": 498,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def rand_dm(d, n=2):\n",
    "    # n controls the degree of mixedness\n",
    "    A = np.random.randn(d, n)\n",
    "    rho = A @ A.T\n",
    "    return rho/rho.trace()\n",
    "\n",
    "rho = rand_dm(2); rho"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can expand it in terms of POVM elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 499,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.3963, 0.0164, 0.5872])"
      ]
     },
     "execution_count": 499,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = np.array([(e @ rho).trace() for e in E]); p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 500,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9999999999999998"
      ]
     },
     "execution_count": 500,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The components sum to 1: they are the probabilities for the outcomes to this experiment. \n",
    "\n",
    "<hr>\n",
    "\n",
    "Now how do we get the density matrix back? The dual basis, of course! But how do we find it?\n",
    "\n",
    "Given a matrix $A$, denote by $|A\\rangle \\rangle$ the vector you get by reshaping the matrix into a vector. \n",
    "\n",
    "(If you want to be fancy, you can think of this like $\\sum_{i} (A \\otimes I)\\mid i \\rangle\\mid i \\rangle$, or as $A$ acting on the left on the maximally entangled \"cup\" state on the doubled Hilbert space, $H_{d} \\otimes H_{d}$).\n",
    "\n",
    "Then form the \"(operator) frame operator\":\n",
    "\n",
    "$F = \\sum_{i} \\mid E_{i} \\rangle \\rangle \\langle \\langle E_{i} \\mid = \\sum_{i} E_{i} \\otimes E_{i}$.\n",
    "\n",
    "I should note that the right hand side is only true if we have rank-1 POVM elements. In that case:\n",
    "\n",
    "$$ |E_{i}\\rangle\\rangle = (E_{i} \\otimes I) \\sum_{j} |j\\rangle | j \\rangle = (|r_{i}\\rangle\\langle r_{i}| \\otimes I) \\sum_{j} |j\\rangle | j \\rangle = \\sum_{j} |r_{i}\\rangle\\langle r_{i}|j\\rangle | j \\rangle = |r_{i}\\rangle \\sum_{j} \\langle r_{i} | j \\rangle | j \\rangle = |r_{i}\\rangle |r_{i}\\rangle$$\n",
    "\n",
    "So that:\n",
    "\n",
    "$$ \\sum_{i} |E_{i}\\rangle\\rangle \\langle \\langle E_{i} | = \\sum_{i} |r_{i}\\rangle |r_{i}\\rangle \\langle r_{i} | \\langle r_{i} | = \\sum_{i} |r_{i}\\rangle \\langle r_{i}| \\otimes | r_{i} \\rangle \\langle r_{i} | = \\sum_{i} E_{i} \\otimes E_{i}$$\n",
    "\n",
    "In the second line, we use the fact that $| a \\rangle | b \\rangle \\langle c | \\langle d| = |a\\rangle\\langle c| \\otimes |b\\rangle\\langle d|$.\n",
    "\n",
    "Finally, take the (pseudo) inverse of $F$. Then, $F^{-1}\\mid E_{i}\\rangle \\rangle = \\mid \\tilde{E}_{i} \\rangle \\rangle$, where $\\tilde{E}$ is the dual basis. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 501,
   "metadata": {},
   "outputs": [],
   "source": [
    "F = sum([np.kron(e, e) for e in E])\n",
    "Finv = np.linalg.pinv(F)\n",
    "D = [(Finv @ e.reshape(d**2)).reshape(d,d) for e in E]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we can recover the original density matrix by weighting the elements of the dual basis by the probabilities and summing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 502,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.123 , -0.3285],\n",
       "       [-0.3285,  0.877 ]])"
      ]
     },
     "execution_count": 502,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rho"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 503,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.123 , -0.3285],\n",
       "       [-0.3285,  0.877 ]])"
      ]
     },
     "execution_count": 503,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum([p[i]*D[i] for i in range(n)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Conversely, we could expand in dual elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 504,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.1974, -0.0119,  1.5699])"
      ]
     },
     "execution_count": 504,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q = np.array([(e @ rho).trace() for e in D]); q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And recover the original density matrix by weighting POVM elements and summing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 505,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.123 , -0.3285],\n",
       "       [-0.3285,  0.877 ]])"
      ]
     },
     "execution_count": 505,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum([q[i]*E[i] for i in range(n)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By the way, the dual basis of an informationally complete POVM <a href=\"https://physics.stackexchange.com/questions/588467/dual-of-ic-povm-cannot-be-positive\">doesn't form a POVM itself</a>. Indeed, its elements must be indefinite!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, instead of working with the dual basis directly, we could play a similar trick from before. \n",
    "\n",
    "Instead of taking $G^{-1}$, now we take $(G \\circ G )^{-1}$ where $\\circ$ is the Hadamard (element-wise) square. Instead of expanding in the dual basis, we can act with $(G \\circ G )^{-1}$ on the vector of probabilities to get the proper weights for POVM elements required to recover $\\rho$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 506,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.123 , -0.3285],\n",
       "       [-0.3285,  0.877 ]])"
      ]
     },
     "execution_count": 506,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rho"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 507,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.3963, 0.0164, 0.5872])"
      ]
     },
     "execution_count": 507,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 508,
   "metadata": {},
   "outputs": [],
   "source": [
    "G = R.T @ R\n",
    "Gsq_inv = np.linalg.inv(G**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 509,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.1974, -0.0119,  1.5699])"
      ]
     },
     "execution_count": 509,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q = Gsq_inv @ p; q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now get our density matrix back:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 510,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.123 , -0.3285],\n",
       "       [-0.3285,  0.877 ]])"
      ]
     },
     "execution_count": 510,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum([q[i]*E[i] for i in range(n)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "But we can do a little bit more, which will prove useful.\n",
    "\n",
    "We could form a matrix of conditional probabilities supposing we repeat our POVM twice. In other words, a matrix of conditional probabilities: \n",
    "\n",
    "$$ Pr(E_{i}|E_{j}) = tr(E_{i}\\frac{E_{j}}{tr(E_{j})})$$\n",
    "\n",
    "After the first measurement, we end up in some $\\frac{E_{j}}{tr(E_{j})}$, and then we look at the probabilities that we get some $E_{i}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 511,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.4485, 0.0208, 0.3867],\n",
       "       [0.0447, 0.9637, 0.0255],\n",
       "       [0.5069, 0.0155, 0.5879]])"
      ]
     },
     "execution_count": 511,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E_given_E = np.array([[(E[i]@(E[j]/E[j].trace())).trace() for j in range(n)] for i in range(n)]); E_given_E"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that this related to $(G \\circ G )$, except for the fact that we divide by $tr(E_{j})$. \n",
    "\n",
    "We can take its inverse to obtain the \"Born matrix\" $\\Phi$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 512,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 8.7087, -0.0956, -5.7242],\n",
       "       [-0.2054,  1.0407,  0.09  ],\n",
       "       [-7.5032,  0.0549,  6.6341]])"
      ]
     },
     "execution_count": 512,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phi = np.linalg.inv(E_given_E); phi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What's the point? If we take our density matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 513,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.123 , -0.3285],\n",
       "       [-0.3285,  0.877 ]])"
      ]
     },
     "execution_count": 513,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rho"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And its probabilities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 514,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.3963, 0.0164, 0.5872])"
      ]
     },
     "execution_count": 514,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can form:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 515,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.0885, -0.0114,  0.9229])"
      ]
     },
     "execution_count": 515,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q = phi @ p; q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 516,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 516,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is nice, because by acting on our vector of probabilities with $\\Phi$, we get back a vector of *quasi-probabilities*: while the components might be negative, they nevertheless sum to 1. And this will prove to be useful.\n",
    "\n",
    "We can then get our density matrix back by summing over normalized elements of our POVM, weighted by the quasi-probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 517,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.123 , -0.3285],\n",
       "       [-0.3285,  0.877 ]])"
      ]
     },
     "execution_count": 517,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum([q[i]*E[i]/E[i].trace() for i in range(n)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Moreover, we can also obtain this matrix $\\Phi$ from the dual elements. Actually, it's nicer to renormalize the dual elements first so that they have the same weights as the POVM elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 518,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0.],\n",
       "       [0., 1.]])"
      ]
     },
     "execution_count": 518,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D = [E[i].trace()*D[i] for i in range(n)]\n",
    "sum(D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 519,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 8.7087, -0.0956, -5.7242],\n",
       "       [-0.2054,  1.0407,  0.09  ],\n",
       "       [-7.5032,  0.0549,  6.6341]])"
      ]
     },
     "execution_count": 519,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D_given_D = np.array([[(D[i]@(D[j]/D[j].trace())).trace() for j in range(n)] for i in range(n)]); D_given_D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Expanding in renormalized dual elements gives us our quasiprobabilities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 520,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.0885, -0.0114,  0.9229])"
      ]
     },
     "execution_count": 520,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q = np.array([(d @ rho).trace() for d in D]); q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 521,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.123 , -0.3285],\n",
       "       [-0.3285,  0.877 ]])"
      ]
     },
     "execution_count": 521,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum([q[i]*E[i]/E[i].trace() for i in range(n)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And symmetrically:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 522,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.123 , -0.3285],\n",
       "       [-0.3285,  0.877 ]])"
      ]
     },
     "execution_count": 522,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum([p[i]*D[i]/D[i].trace() for i in range(n)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deforming the Law of Total Probability\n",
    "\n",
    "Okay, but why should we care?\n",
    "\n",
    "We're going to consider two measurements on our quantum system, $A$ and $B$, where $A$ is an IC-POVM, and $B$ could be a standard Von Neumann measurement or a POVM, doesn't matter. Let's say it's an IC-POVM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 523,
   "metadata": {},
   "outputs": [],
   "source": [
    "d, n = 2, 3\n",
    "A = [np.outer(r, r) for r in sc.linalg.polar(np.random.randn(d, n))[0].T]\n",
    "B = [np.outer(r, r) for r in sc.linalg.polar(np.random.randn(d, n))[0].T]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 524,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.123 , -0.3285],\n",
       "       [-0.3285,  0.877 ]])"
      ]
     },
     "execution_count": 524,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rho"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can of course nail down our quantum state by its probabilities on measurement $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 525,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.4322, 0.525 , 0.0428])"
      ]
     },
     "execution_count": 525,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_A = np.array([(a @ rho).trace() for a in A]); p_A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're going to compare two sequences of events.\n",
    "\n",
    "First: we measure $A$, and then measure $B$. Before we look at the outcome for $A$, however, what probabilities should we assign to the outcomes of $B$?\n",
    "\n",
    "To figure this out, we could use a matrix of conditional probabilities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 526,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.111 , 0.0003, 0.1581],\n",
       "       [0.0468, 0.8323, 0.4232],\n",
       "       [0.8423, 0.1674, 0.4187]])"
      ]
     },
     "execution_count": 526,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B_given_A = np.array([[(b@a/a.trace()).trace() for a in A] for b in B]); B_given_A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are the conditional probabilities for each of $B$'s outcomes, given each outcome state of $A$. This gives us conditional probabilities $Pr(B_{i}|A_{j})$. \n",
    "\n",
    "We end up with a column stochastic matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 527,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 1., 1.])"
      ]
     },
     "execution_count": 527,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(B_given_A, axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then use the law of total probability: $Pr(B_{i}) = \\sum_{j} Pr(B_{i}|A_{j}) Pr(A_{j})$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 528,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.0549, 0.4753, 0.4698])"
      ]
     },
     "execution_count": 528,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_AB = B_given_A @ p_A; p_AB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 529,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9999999999999991"
      ]
     },
     "execution_count": 529,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(p_AB)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check this. After the measurement of $A$, we should regard our system's state as a statistical mixture of $A$'s outcomes, weighted by the probability of their occurance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 530,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.2658, -0.0957],\n",
       "       [-0.0957,  0.7342]])"
      ]
     },
     "execution_count": 530,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rho_post_A = sum([p_A[i]*a/a.trace() for i, a in enumerate(A)]); rho_post_A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then get the $B$ probabilities for this density matrix by taking inner products with $B$'s elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 531,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.0549, 0.4753, 0.4698])"
      ]
     },
     "execution_count": 531,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_AB2 = np.array([(b @ rho_post_A).trace() for b in B]); p_AB2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 532,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9999999999999993"
      ]
     },
     "execution_count": 532,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(p_AB2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get the same probabilities. \n",
    "\n",
    "<hr>\n",
    "\n",
    "To recap: we did measurement of $A$, but didn't look at the answer. We just wanted the probabilities for the future measurement $B$. So we used the law of total probability: $p(B_{i}) = \\sum_{j} p(B_{i}|A_{j}) p(A_{j})$.\n",
    "\n",
    "Okay. \n",
    "\n",
    "The second experiment is that we leave $A$ hypothetical, but instead just do the $B$ measurement. We know what the answer should be: we just expand in terms of $B$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 533,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.0465, 0.3126, 0.641 ])"
      ]
     },
     "execution_count": 533,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_B = np.array([(b @ rho).trace() for b in B]); p_B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But since our probabilities on $A$ completely nail down the state, there ought to be a way of getting these latter probabilities in terms of them...\n",
    "\n",
    "It's easy! We simply form the \"Born matrix\" $\\Phi$ from before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 534,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi = np.linalg.inv([[(a1 @ a0/a0.trace()).trace() for a0 in A] for a1 in A])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And instead of doing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 535,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.0549, 0.4753, 0.4698])"
      ]
     },
     "execution_count": 535,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B_given_A @ p_A "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 536,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.0465, 0.3126, 0.641 ])"
      ]
     },
     "execution_count": 536,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B_given_A @ phi @ p_A "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence the name the \"Born matrix\": it allows us to recover the Born rule. Notice how everything is still formulated in terms of a) the conditional probability matrix for B's outcomes given A's outcomes, and b) the probability vector for A's outcomes. \n",
    "\n",
    "One can now see the usefulness of:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 537,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.8023,  0.4678, -0.2701])"
      ]
     },
     "execution_count": 537,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q = phi @ p_A; q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 538,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 538,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fact that $\\Phi$ takes us to quasiprobabilities makes sure we play nicely with the conditional probability matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way of thinking about this is to form a \"conditional quasiprobability matrix\":\n",
    "\n",
    "$$ Qr(B_{i}|\\tilde{A}_{j}) = tr(B_{i}\\frac{\\tilde{A}_{j}}{tr(\\tilde{A}_{j})})$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 539,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.1535, -0.0528,  0.1832],\n",
       "       [-0.8035,  1.2063,  0.6183],\n",
       "       [ 1.65  , -0.1535,  0.1985]])"
      ]
     },
     "execution_count": 539,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "FA = sum([np.outer(e.reshape(d**2), e.reshape(d**2)) for e in A])\n",
    "FAinv = np.linalg.pinv(FA)\n",
    "DA = [e.trace()*(FAinv @ e.reshape(d**2)).reshape(d,d) for e in A]\n",
    "B_given_DA = np.array([[(B[i]@(DA[j]/DA[j].trace())).trace() for j in range(n)] for i in range(n)]); B_given_DA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 540,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.0465, 0.3126, 0.641 ])"
      ]
     },
     "execution_count": 540,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B_given_DA @ p_A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 541,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 541,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(B_given_DA, B_given_A@phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perhaps more concretely:\n",
    "\n",
    "If we perform POVM $A$, then:\n",
    "\n",
    "$$ \\rho \\rightarrow \\rho^{\\prime} = \\sum_{i} Pr(A_{i}|\\rho)\\frac{A_{i}}{tr(A_{i})}$$\n",
    "\n",
    "If we *don't* perform POVM $A$, then:\n",
    "\n",
    "$$ \\rho \\rightarrow \\sum_{i} Pr(A_{i}|\\rho)\\frac{\\tilde{A}_{i}}{tr(\\tilde{A}_{i})}$$\n",
    "\n",
    "Which is just $\\rho$! \n",
    "\n",
    "If we perform POVM $A$, and then POVM $B$:\n",
    "\n",
    "$$ \\rho \\rightarrow \\sum_{i} \\Big{(} \\sum_{j} Pr(B_{i}|A_{j}) Pr(A_{j}|\\rho) \\Big{)} \\frac{B_{i}}{tr(B_{i})} = \\sum_{i} Pr(B_{i}|\\rho^{\\prime}) \\frac{B_{i}}{tr(B_{i})}$$\n",
    "\n",
    "If we *don't* perform POVM $A$, and then do perform POVM $B$:\n",
    "\n",
    "$$ \\rho \\rightarrow  \\sum_{i} \\big{(} \\sum_{j} Qr(B_{i}|\\tilde{A}_{j}) Pr(A_{j}|\\rho) \\big{)} \\frac{B_{i}}{tr(B_{i})} =  \\sum_{i} Pr(B_{i}|\\rho)  \\frac{B_{i}}{tr(B_{i})} $$\n",
    "\n",
    "And if we did neither:\n",
    "\n",
    "$$ \\rho \\rightarrow  \\sum_{i} \\big{(} \\sum_{j} Qr(B_{i}|\\tilde{A}_{j}) Pr(A_{j}|\\rho) \\big{)} \\frac{\\tilde{B}_{i}}{tr(\\tilde{B}_{i})} =  \\sum_{i} Pr(B_{i}|\\rho)  \\frac{\\tilde{B}_{i}}{tr(\\tilde{B}_{i})} $$\n",
    "\n",
    "Which is just $\\rho$!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "## Quantumness \n",
    "\n",
    "So in the first case, we did $A$, and afterwards we described our system as a mixed state of A's outcomes, weighted by the probability of them happening. No \"coherence\" between outcomes: there's just \"classical ignorance.\" And we used the old school law of total probability to get the probabilities for the outcomes of $B$.\n",
    "\n",
    "In the second case, we didn't do $A$: and so \"quantum coherence\" was preserved. But we could still get the probabilities for a $B$ measurement from the probabilities of the (hypothetical) $A$ measurement.\n",
    "\n",
    "One way of thinking of this is that quantum mechanics encourages us to modify the law of total probability. Classically speaking, we could use the law of total probability whether or not we did an intermediate measurement of $A$ before $B$. But quantum mechanically, it matters. If we don't do $A$, we have to stick the Born matrix $\\Phi$ in between the conditional probability matrix and the probability vector. In matrix notation, instead of $\\vec{Pr(B_{i})} = \\hat{Pr(B_{i}|A_{j})}\\vec{Pr(A_{j})}$, we have $\\vec{Pr(B_{i})} = \\hat{Pr(B_{i}|A_{j})}\\hat{\\Phi}\\vec{Pr(A_{j})}$.\n",
    "\n",
    "Now, if $\\Phi$ were the identity, then there would be no difference, as it were, between classical and quantum. No coherence, everything having \"objective properties,\" law of total probability always valid, etc. The lesson of quantum mechanics is that: $\\Phi$ is not the identity.\n",
    "\n",
    "So you might wonder: what measurements are associated with a $\\Phi$ such that $\\Phi$ is *as close to the identity as possible*?\n",
    "\n",
    "With this in mind, we form a quantity we'll call \"quantumness\": $Q = ||\\Phi - I||$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 542,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quantumness(E):\n",
    "    phi = np.linalg.inv([[(e1@e0/e0.trace()).trace() for e0 in E] for e1 in E])\n",
    "    return np.linalg.norm(phi - np.eye(len(E)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 543,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.0756601741411402, 19.00145583246328)"
      ]
     },
     "execution_count": 543,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "quantumness(A), quantumness(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This quantity diagnoses the extent to which the law of total probability is deformed.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now it's known that SIC-POVM's aka \"symmetrically informationally complete POVM's,\" which correspond to equiangular lines, minimize quantumness. And in fact, the proof that they do so doesn't depend on the particulars of which norm you choose to measure the difference between $\\Phi$ and the identity. \n",
    "\n",
    "In other words, and it makes sense, SIC-POVM's, whose elements all make the same angle with each other, have $\\Phi$ matrices as close to the identity as possible. They are the quantum measurements that are \"irreducibly quantum.\"\n",
    "\n",
    "Now much of the study of SIC-POVM's takes place in the context of *real* quantum mechanics, I mean, complex vector space quantum mechanics. In this case, SIC-POVM's would correspond to *complex equiangular lines*, which are a bit more difficult to visualize. And to span the space of density matrices, you'd need $d^2$ of them for dimension $d$.\n",
    "\n",
    "For example, in the case of $d=2$, the four elements of a SIC-POVM make up the vertices of a regular tetrahedron on the Bloch sphere. \n",
    "\n",
    "It's been conjectured, and widely believed, on the basis of a great deal of theoretical and computational effort, that SIC-POVM's exist in all complex dimensions. And conditional on this being true, then we know the irreducibly quantum measurements in any dimension: they are the SIC's.\n",
    "\n",
    "<hr>\n",
    "\n",
    "Here's the thing, though. \n",
    "\n",
    "As I've said: in real quantum mechanics, RIC-POVM's (real informationally complete POVM's) require $n=\\frac{d(d+1)}{2}$ elements. As we've seen, this is fine for $d=2$: in this case, $n=3$, and we have the \"Mercedes-Benz\" frame. So too, all is well in $d=3$. Then $n=6$, and we have our icosahedron.\n",
    "\n",
    "But then the shit hits the fan. \n",
    "\n",
    "Here's a <a href=\"https://oeis.org/A002853\">table</a> listing the known maximum number of equiangular lines in different real dimensions compared to the dimensionality of the symmetric operator space. \n",
    "\n",
    "I've starred the ones where SIC-POVM's can exist, namely: 2, 3, 7, 23 (so far as we know).\n",
    "\n",
    "| d | max| n |\n",
    "|---|-----|----------------|\n",
    "| 2 |  3  |  3* |\n",
    "| 3 |  6  |  6*|\n",
    "| 4 |  6  |  10|\n",
    "| 5 | 10  |  15|\n",
    "| 6 | 16  |  21|\n",
    "| 7 | 28 |  28*|\n",
    "| 8 | 28 | 36|\n",
    "| 9 | 28 | 45|\n",
    "| 10 | 28 | 55|\n",
    "| 11 | 28 | 66|\n",
    "| 12 | 28 | 78|\n",
    "| 13 | 28 | 91|\n",
    "|14 | 28| 105|\n",
    "|15 | 36 | 120|\n",
    "|16 | 40 | 136|\n",
    "|17 | 48-49| 153|\n",
    "|18 | 56-60| 171|\n",
    "|19 | 72-74| 190|\n",
    "|20 | 90-94| 210|\n",
    "|21 | 126 | 231|\n",
    "|22 | 176 | 253|\n",
    "|23 | 276 | 276*|\n",
    "\n",
    "(Incidentally, the famous Gerzon bound says there definitely can't be more that $\\frac{d(d+1)}{2}$ equiangular lines in $\\mathbb{R}^n$).\n",
    "\n",
    "You might wonder: given that (unlike in the complex case), we definitely know there aren't SIC's in all real dimensions, what's the *next best thing*? Specifically, what are the RIC's with lowest quantumness in each real dimension where a SIC doesn't exist? \n",
    "\n",
    "It's an open question!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "## Bibliography"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://www.math.auckland.ac.nz/~waldron/Preprints/Frame-book/draft26-7.17.pdf\">\n",
    "An Introduction to Finite Tight Frames</a>\n",
    "\n",
    "<a href=\"https://faculty.sites.iastate.edu/esweber/files/inline-files/book_final.pdf\">Frames for Undergraduates</a>\n",
    "\n",
    "<a href=\"http://www1.chapman.edu/~mbvajiac/conferences/chapman_lec03_slides.pdf\">Slides on Frame Theory</a>\n",
    "    \n",
    "<a href=\"https://arxiv.org/abs/1311.3219\">New bounds for equiangular lines</a>\n",
    "\n",
    "<a href=\"https://www.lib.uidaho.edu/digital/etd/items/oldroyd_idaho_0089e_10817.html\">Generalizations and Approximations of Equiangular Tight Frames</a>\n",
    "\n",
    "<a href=\"https://arxiv.org/abs/1812.08762\">The Varieties of Minimal Tomographically Complete Measurements</a>\n",
    "\n",
    "<a href=\"https://arxiv.org/abs/1805.08721\">Symmetric Informationally Complete Measurements Identify the Irreducible Difference between Classical and Quantum Systems</a>\n",
    "\n",
    "<a href=\"https://arxiv.org/abs/quant-ph/0310075\">Symmetric Informationally Complete Quantum Measurements</a>\n",
    "\n",
    "<a href=\"https://arxiv.org/abs/1911.05809\">Sporadic SICs and Exceptional Lie Algebras</a>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "VPython",
   "language": "python",
   "name": "vpython"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}