docs/source/example_notebooks/flattening_and_folding_psd_matrix.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Flattening and Folding With Covariance Matrices. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import paragami"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we will consider flattening and folding a simple symmetric positive semi-definite matrix:\n",
    "\n",
    "$$\n",
    "A = \\left[\n",
    "\\begin{matrix}\n",
    "a_{11} & a_{12} & a_{13}  \\\\\n",
    "a_{21} & a_{22} & a_{23}  \\\\\n",
    "a_{31} & a_{32} & a_{33}  \\\\\n",
    "\\end{matrix}\n",
    "\\right].\n",
    "$$\n",
    "\n",
    "Of course, symmetry and positive semi-definiteness impose constraints on the entries $a_{ij}$ of $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flattening and Folding."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In the Original Space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us first consider how to represent $A$ as a vector, which we call simply *flattening*, and then as an unconstrained vector, which we call *free flattening*.\n",
    "\n",
    "When a parameter is flattened, it is simply re-shaped as a vector.  Every number that was in the original parameter will occur exactly once in the flattened shape.  (In the present case of a matrix, this is exactly the same as ``np.flatten``.)\n",
    "\n",
    "$$\n",
    "A = \\left[\n",
    "\\begin{matrix}\n",
    "a_{11} & a_{12} & a_{13}  \\\\\n",
    "a_{21} & a_{22} & a_{23}  \\\\\n",
    "a_{31} & a_{32} & a_{33}  \\\\\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\xrightarrow{flatten}\n",
    "A_{flat} = \n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "a_{flat,1} \\\\\n",
    "a_{flat,2} \\\\\n",
    "a_{flat,3} \\\\\n",
    "a_{flat,4} \\\\\n",
    "a_{flat,5} \\\\\n",
    "a_{flat,6} \\\\\n",
    "a_{flat,7} \\\\\n",
    "a_{flat,8} \\\\\n",
    "a_{flat,9} \\\\\n",
    "\\end{matrix}\\right]\n",
    "=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "a_{11} \\\\\n",
    "a_{12} \\\\\n",
    "a_{13} \\\\\n",
    "a_{21} \\\\\n",
    "a_{22} \\\\\n",
    "a_{23} \\\\\n",
    "a_{31} \\\\\n",
    "a_{32} \\\\\n",
    "a_{33} \\\\\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Converting to and from $A$ and $A_{flat}$ can be done with the `flatten` method of a `paragami.PSDSymmetricMatrixPattern` pattern.  \n",
    "\n",
    "For the moment, because we are flattening, not free flattening, we use the option `free=False`.  We will discuss the `free=True` option shortly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Now, a_flat contains the elements of a exactly as shown in the formula above.\n",
      "\n",
      "a:\n",
      "[[1.2712968  0.74536048 0.33203184]\n",
      " [0.74536048 1.91869072 0.4602062 ]\n",
      " [0.33203184 0.4602062  1.58338   ]]\n",
      "\n",
      "a_flat:\n",
      "[1.2712968  0.74536048 0.33203184 0.74536048 1.91869072 0.4602062\n",
      " 0.33203184 0.4602062  1.58338   ]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# A sample positive semi-definite matrix.\n",
    "a = np.eye(3) + np.random.random((3, 3))\n",
    "a = 0.5 * (a + a.T)\n",
    "\n",
    "# Define a pattern and fold.\n",
    "a_pattern = paragami.PSDSymmetricMatrixPattern(size=3)\n",
    "a_flat = a_pattern.flatten(a, free=False)\n",
    "\n",
    "print('Now, a_flat contains the elements of a exactly as shown in the formula above.\\n')\n",
    "print('a:\\n{}\\n'.format(a))\n",
    "print('a_flat:\\n{}\\n'.format(a_flat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also convert from $A_{flat}$ back to $A$ by 'folding'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Folding the flattened value recovers the original matrix.\n",
      "\n",
      "a:\n",
      "[[1.2712968  0.74536048 0.33203184]\n",
      " [0.74536048 1.91869072 0.4602062 ]\n",
      " [0.33203184 0.4602062  1.58338   ]]\n",
      "\n",
      "a_fold:\n",
      "[[1.2712968  0.74536048 0.33203184]\n",
      " [0.74536048 1.91869072 0.4602062 ]\n",
      " [0.33203184 0.4602062  1.58338   ]]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print('Folding the flattened value recovers the original matrix.\\n')\n",
    "a_fold = a_pattern.fold(a_flat, free=False)\n",
    "print('a:\\n{}\\n'.format(a))\n",
    "print('a_fold:\\n{}\\n'.format(a_fold))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, flattening and folding perform checks to make sure the result is a valid instance of the parameter type -- in this case, a symmetric positive definite matrix. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The diagonal of a positive semi-definite matrix must not be less than 0, and folding checks this when ``validate=True``, which it is by default.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A bad folded value: [-1  0  0  0  0  0  0  0  0]\n",
      "Folding with a_pattern raised the following ValueError:\n",
      "Diagonal is less than the lower bound 0.0.\n"
     ]
    }
   ],
   "source": [
    "a_flat_bad = np.array([-1, 0, 0,  0, 0, 0,  0, 0, 0])\n",
    "print('A bad folded value: {}'.format(a_flat_bad))\n",
    "try:\n",
    "    a_fold_bad = a_pattern.fold(a_flat_bad, free=False)\n",
    "except ValueError as err:\n",
    "    print('Folding with a_pattern raised the following ValueError:\\n{}'.format(err))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If ``validate_value`` is ``False``, folding will produce an invalid matrix without an error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Folding a non-pd matrix with validate=False:\n",
      "[[-1  0  0]\n",
      " [ 0  0  0]\n",
      " [ 0  0  0]]\n"
     ]
    }
   ],
   "source": [
    "a_fold_bad = a_pattern.fold(a_flat_bad, free=False, validate_value=False)\n",
    "print('Folding a non-pd matrix with validate=False:\\n{}'.format(a_fold_bad))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, it will not produce a matrix of the wrong shape even when validate is ``False``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A very bad folded value: [1 0 0].\n",
      "Folding with a_pattern raised the following ValueError:\n",
      "Wrong length for PSDSymmetricMatrix flat value.\n"
     ]
    }
   ],
   "source": [
    "a_flat_very_bad = np.array([1, 0, 0])\n",
    "print('A very bad folded value: {}.'.format(a_flat_very_bad))\n",
    "try:\n",
    "    a_fold_very_bad = a_pattern.fold(a_flat_very_bad, free=False, validate_value=False)\n",
    "except ValueError as err:\n",
    "    print('Folding with a_pattern raised the following ValueError:\\n{}'.format(err))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can always check validity of a folded value with the ``validate_folded`` method of a pattern, which returns a boolean and an error message."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Valid: True.\tMessage: \n",
      "Valid: False.\tMessage: Diagonal is less than the lower bound 0.0.\n"
     ]
    }
   ],
   "source": [
    "valid, msg = a_pattern.validate_folded(a_fold)\n",
    "print('Valid: {}.\\tMessage: {}'.format(valid, msg))\n",
    "\n",
    "valid, msg = a_pattern.validate_folded(a_fold - 10 * np.eye(3))\n",
    "print('Valid: {}.\\tMessage: {}'.format(valid, msg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In an Unconstrained Space: \"Free\" Flattening and Folding."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ordinary flattening converts a 3x3 symmetric PSD matrix into a 9-d vector.  However, as seen above, not every 9-d vector is a valid 3x3 symmetric positive definite matrix.  It is useful to have an \"free\" flattened representation of a parameter, where every finite value of the free flattened vector corresponds is guaranteed valid.\n",
    "\n",
    "To accomplish this for a symmetric positive definite matrix, we consider the Cholesky decomposition $A_{chol}$. This is an lower-triangular matrix with positive diagonal entries such that $A = A_{chol} A_{chol}^T$.  By taking the log of the diagonal of $A_{chol}$ and stacking the non-zero entries, we can construct a 6-d vector, every value of which corresponds to a symmetric PSD matrix.\n",
    "\n",
    "$$\n",
    "% A \\xrightarrow{\\textrm{free flatten}} A_{freeflat} \\quad\\quad \\textrm{where} \\\\\n",
    "A \\xrightarrow{}\n",
    "A_{chol} = \n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "\\alpha_{11} & 0 & 0  \\\\\n",
    "\\alpha_{21} & \\alpha_{22} & 0  \\\\\n",
    "\\alpha_{31} & \\alpha_{32} & \\alpha_{33}  \\\\\n",
    "\\end{matrix}\n",
    "\\right] \\xrightarrow{}\n",
    "A_{freeflat} =\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "\\log(\\alpha_{11}) \\\\\n",
    "\\alpha_{21} \\\\\n",
    "\\alpha_{31} \\\\\n",
    "\\log(\\alpha_{22})\\\\\n",
    "\\alpha_{32} \\\\\n",
    "\\log(\\alpha_{33})\n",
    "\\end{matrix}\n",
    "\\right].\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The details of the freeing transform aren't important to the end user, as `paragami` takes care of the transformation behind the scenes with the option `free=True`.  We denote the flattened $A$ in the free parameterization as $A_{freeflat}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The free flat value ``a_freeflat`` is not immediately recognizable as ``a``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a:\n",
      "[[1.2712968  0.74536048 0.33203184]\n",
      " [0.74536048 1.91869072 0.4602062 ]\n",
      " [0.33203184 0.4602062  1.58338   ]]\n",
      "\n",
      "a_freeflat:\n",
      "[0.12001874 0.66106306 0.19659043 0.2944803  0.21814513 0.18546238]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "a_freeflat = a_pattern.flatten(a, free=True)\n",
    "print('a:\\n{}\\n'.format(a))\n",
    "print('a_freeflat:\\n{}\\n'.format(a_freeflat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, it transforms correctly back to a when folded.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a:\n",
      "[[1.2712968  0.74536048 0.33203184]\n",
      " [0.74536048 1.91869072 0.4602062 ]\n",
      " [0.33203184 0.4602062  1.58338   ]]\n",
      "\n",
      "a_fold:\n",
      "[[1.2712968  0.74536048 0.33203184]\n",
      " [0.74536048 1.91869072 0.4602062 ]\n",
      " [0.33203184 0.4602062  1.58338   ]]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "a_freefold = a_pattern.fold(a_freeflat, free=True)\n",
    "print('a:\\n{}\\n'.format(a))\n",
    "print('a_fold:\\n{}\\n'.format(a_freefold))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Any length-six vector will free fold back to a valid PSD matrix up to floating point error.  Let's draw 100 random vectors, fold them, and check that this is true.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Draw random free vectors and confirm that they are positive semi definite.\n",
    "def assert_is_pd(mat):\n",
    "    eigvals = np.linalg.eigvals(mat)\n",
    "    assert np.min(eigvals) >= -1e-8\n",
    "\n",
    "for draw in range(100):\n",
    "    a_rand_freeflat = np.random.normal(scale=2, size=(6, ))\n",
    "    a_rand_fold = a_pattern.fold(a_rand_freeflat, free=True)\n",
    "    assert_is_pd(a_rand_fold)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}