Orthonormalization.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Matplotlib is building the font cache; this may take a moment.\n"
     ]
    }
   ],
   "source": [
    "#math and linear algebra stuff\n",
    "import numpy as np\n",
    "\n",
    "#plots\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['figure.figsize'] = (15.0, 15.0)\n",
    "#mpl.rc('text', usetex = True)\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Orthonormalization /whitening process\n",
    "This notebooks exposes how to perform orthogonalization on a set of vector, such that the resulting method is fast and accurate. This topic is mentionned in the notebook \"StatisticalTest\" where we investigate the use of whitening and decorrelation transforms\n",
    "\n",
    "## Orthonormalization through Gram-Schmidt\n",
    "\n",
    "The gram-schmidt process is the most reknown method to orthogonormalize a set of vector, see [this page](https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process).\n",
    "\n",
    "The idea, is to loop over the set of vector $\\{\\vec{v_0}, \\vec{v_1}, \\dots, \\vec{v_{n-1}} \\}$, project (by the mean of a simple dot product) each new vector $\\vec{v_k}$ onto the set of previously orthonormalized vector $\\{\\vec{e_0}, \\vec{e_1}, \\dots, \\vec{e_{k-1}} \\}$ that one can substract the linearly dependent components:\n",
    "\n",
    "\\begin{align*}\n",
    "  \\vec{e_k'} = \\vec{v_k} - \\sum_{i=0}^{k-1} \\left( \\vec{v_k} \\cdot \\vec{e_i} \\right)  \\vec{e_i}\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "Once $\\vec{e_k'}$ is built, one can normalize it : $\\vec{e_k}=\\frac{\\vec{e_k'}}{\\|\\vec{e_k'}\\|}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Gram-Schmidt method for basis orthonormalization\n",
    "\"\"\"\n",
    "size = 25\n",
    "nbVec = 15\n",
    "V = np.random.rand(size,nbVec)\n",
    "\n",
    "#Naive orthogonalization loop\n",
    "for i in np.arange(nbVec):\n",
    "  #get current vector\n",
    "  ek=V[:,i]\n",
    "  if i>0:\n",
    "    #Compute projection onto other vector\n",
    "    proj = np.dot(ek,V[:,:i])\n",
    "    #orthogonalize\n",
    "    ek = ek - np.dot(V[:,:i], proj)\n",
    "  #Normalize vector and update V\n",
    "  V[:,i] = ek / np.sqrt(np.dot(ek,ek))\n",
    "\n",
    "#check that V is orthonormal\n",
    "assert( np.allclose(np.identity(nbVec),np.dot(V.T,V)) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Orthonormalization using Gram matrix and Choleski\n",
    "\n",
    "Another possible solution is based on the Gram matrix $M_G$ of the set of vectors $\\{\\vec{v_0}, \\vec{v_1}, \\dots, \\vec{v_{n-1}} \\}$, computed as $M_G(i,j)= \\langle \\vec{v_i}, \\vec{v_j} \\rangle$\n",
    "\n",
    "The gram matrix is a symmetric, positive semidefinite matrix, whose entries are actually all computed during the Gram-Schmidt process, however, at a high, level, one can compute $M_G=V^t V$ with a single matrix-matrix multiplication.\n",
    "\n",
    "The interesting part is what happens next. Indeed, one can see, that orthogonalizing a set of vector, amounts to make its gram matrix a diagonal matrix, which can be obtained by using either the eigendecomposition:\n",
    "\n",
    "\\begin{align*}\n",
    "  M_G &= V^t V\\\\\n",
    "  Q D Q^t &= V^t V \\\\\n",
    "  D &= Q^t V^t V Q \\\\\n",
    "  D &= (V Q)^t (V Q)\n",
    "\\end{align*}\n",
    "\n",
    "Orthonormalizing a set of vector, amounts to make its gram matrix equal to the identity matrix, which can be obtained by using  its factorization as a product of square root of $M_G$:\n",
    "\n",
    "\\begin{align*}\n",
    "  M_G &= V^t V\\\\\n",
    "  M_G^{\\frac{1}{2}t} M_G^{\\frac{1}{2}} &= V^t V\\\\\n",
    "  I &= M_G^{-\\frac{1}{2}t} V^t V M_G^{-\\frac{1}{2}}\\\\\n",
    "  I &= (V M_G^{-\\frac{1}{2}})^t (V M_G^{-\\frac{1}{2}})\n",
    "\\end{align*}\n",
    "\n",
    "Such that $(M_G^{\\frac{1}{2}} V)$ is the new orthonormalized set of vector. Now, we are interested in finding, square root of matrix $M_G$, which hopefully is exactly what the Choleski decomposition does. IE from [choleski factorization](https://en.wikipedia.org/wiki/Cholesky_decomposition) we can see that a hermitian positive-definite matrix $A$, admits a decomposition of the form $A = LL^{*}$  where L is a lower triangular matrix with real and positive diagonal entries, and $L^{*}$ denotes the conjugate transpose of $L$. Every Hermitian positive-definite matrix (and thus also every real-valued symmetric positive-definite matrix) has a unique Cholesky decomposition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fa576edf358>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADMlJREFUeJzt3V2sZWV9x/Hvj2HAzjiIiKLCtGBDTKhpCpkQa401pbVI\nCeNFLyC1hWpCvLCVxsQMNam3Whv7kpqaidjSlOiFr8RAFammaVKo4zi8IyClyJS32hQQUmHw34u9\naA7Hc2bOrLfZw/P9JCdn7b2eddZ/1t6/efZeez37SVUhqT3HHOkCJB0Zhl9qlOGXGmX4pUYZfqlR\nhl9qlOGXGmX4pUYZfqlRx865s5NP2lSnb9982Nvdc+uWCaqRXnr+l6d5tn6cjbSdNfynb9/Mv31t\n+2Fv95uv/6UJqpFeem6uGzfc1pf9UqMMv9SoQeFPcn6S7yW5L8musYqSNL3e4U+yCfgk8E7gLOCS\nJGeNVZikaQ3p+c8F7quq+6vqWeBzwM5xypI0tSHhPxX4wYrbD3X3SToKTH7CL8nlSfYk2fP4D5+f\neneSNmhI+PcDKz+0P62770WqandV7aiqHa9+1aYBu5M0piHh/zZwZpIzkhwHXAxcO05ZkqbW+wq/\nqjqQ5P3A14BNwGeq6o7RKpM0qUGX91bVdcB1I9UiaUZe4Sc1yvBLjZp1VN89t27pNUJv04mv6L3P\n5//nid7bSi9l9vxSowy/1CjDLzXK8EuNMvxSowy/1CjDLzXK8EuNMvxSowy/1CjDLzXK8EuNMvxS\no2Yd1dfXkJF5m044of9+n3yy97bSsrPnlxpl+KVGGX6pUUPm6tue5JtJ7kxyR5IPjFmYpGkNOeF3\nAPhgVe1Nsg34TpIbqurOkWqTNKHePX9VPVxVe7vlp4C7cK4+6agxynv+JKcDZwM3j/H3JE1v8Of8\nSV4OfAG4oqp+6oPxJJcDlwO8jC1DdydpJIN6/iSbWQT/mqr64lptVk7UuZnjh+xO0oiGnO0PcBVw\nV1V9YrySJM1hSM//K8DvAr+WZF/3c8FIdUma2JBZev8FyIi1SJqRV/hJjTL8UqOOiiG9QwwZlnvM\ntm29tvvJU0/13qc0F3t+qVGGX2qU4ZcaZfilRhl+qVGGX2qU4ZcaZfilRhl+qVGGX2qU4ZcaZfil\nRhl+qVEv+VF9Q/QdnXfM1q399/n00723lQ6HPb/UKMMvNcrwS40aHP4km5J8N8lXxyhI0jzG6Pk/\nwGKePklHkaEz9pwG/Bbw6XHKkTSXoT3/XwAfAn4yQi2SZjRkuq4Lgceq6juHaHd5kj1J9jzHj/vu\nTtLIhk7XdVGSB4DPsZi26x9WN3KiTmk59Q5/VV1ZVadV1enAxcA/VdW7R6tM0qT8nF9q1CjX9lfV\nt4BvjfG3JM3Dnl9qlOGXGuWQ3gkMGZZ7zJYt/fb5zDO996k22fNLjTL8UqMMv9Qowy81yvBLjTL8\nUqMMv9Qowy81yvBLjTL8UqMMv9Qowy81yvBLjXJU35LpOzovx/f/fsT6sV+s2iJ7fqlRhl9qlOGX\nGjV0uq4Tk3w+yd1J7kryy2MVJmlaQ0/4/SXwj1X120mOA/p9B5Wk2fUOf5JXAG8DLgOoqmeBZ8cp\nS9LUhrzsPwN4HPjbJN9N8ukkW0eqS9LEhoT/WOAc4G+q6mzgaWDX6kZO1CktpyHhfwh4qKpu7m5/\nnsV/Bi/iRJ3SchoyUecjwA+SvLG76zzgzlGqkjS5oWf7/wC4pjvTfz/w+8NLkjSHQeGvqn3AjpFq\nkTQjr/CTGmX4pUY5pPclYsiw3Gw+rt8+n/OarqOZPb/UKMMvNcrwS40y/FKjDL/UKMMvNcrwS40y\n/FKjDL/UKMMvNcrwS40y/FKjDL/UKEf1qffovBzb/+lTBw703lbjsOeXGmX4pUYZfqlRQyfq/KMk\ndyS5Pclnk7xsrMIkTat3+JOcCvwhsKOq3gRsAi4eqzBJ0xr6sv9Y4GeSHMtiht7/HF6SpDkMmbFn\nP/BnwIPAw8ATVfX1sQqTNK0hL/tfCexkMVvv64GtSd69Rjsn6pSW0JCX/b8O/HtVPV5VzwFfBN6y\nupETdUrLaUj4HwTenGRLkrCYqPOuccqSNLUh7/lvZjEt917gtu5v7R6pLkkTGzpR50eAj4xUi6QZ\neYWf1CjDLzXKIb3qzWG5Rzd7fqlRhl9qlOGXGmX4pUYZfqlRhl9qlOGXGmX4pUYZfqlRhl9qlOGX\nGmX4pUYZfqlRhl9qlOGXGmX4pUYZfqlRhwx/ks8keSzJ7SvuOynJDUnu7X6/ctoyJY1tIz3/3wHn\nr7pvF3BjVZ0J3NjdlnQUOWT4q+qfgf9edfdO4Opu+WrgXSPXJWlifd/zn1JVD3fLjwCnjFSPpJkM\nPuFXVQXUeuudqFNaTn3D/2iS1wF0vx9br6ETdUrLqW/4rwUu7ZYvBb4yTjmS5rKRj/o+C/wr8MYk\nDyV5L/BR4DeS3Mtiqu6PTlumpLEdcsaeqrpknVXnjVyLpBl5hZ/UKMMvNcqJOtWOpN92te4n2Uc1\ne36pUYZfapThlxpl+KVGGX6pUYZfapThlxpl+KVGGX6pUYZfapThlxpl+KVGGX6pUY7qUzv6js7r\nOxpwyD5nYM8vNcrwS40y/FKj+k7U+fEkdye5NcmXkpw4bZmSxtZ3os4bgDdV1S8C9wBXjlyXpIn1\nmqizqr5eVQe6mzcBp01Qm6QJjfGe/z3A9SP8HUkzGvQ5f5IPAweAaw7S5nLgcoCXsWXI7iSNqHf4\nk1wGXAic183Uu6aq2g3sBjghJy3vFQ9SY3qFP8n5wIeAX62qZ8YtSdIc+k7U+dfANuCGJPuSfGri\nOiWNrO9EnVdNUIukGXmFn9Qowy81yiG90qEMGZa7xJOD2vNLjTL8UqMMv9Qowy81yvBLjTL8UqMM\nv9Qowy81yvBLjTL8UqMMv9Qowy81yvBLjXJUnzSlvqPzjtnUb7vnD2MX/fYg6Whn+KVGGX6pUb0m\n6lyx7oNJKsnJ05QnaSp9J+okyXbgHcCDI9ckaQa9Jurs/DmLiTuchUc6CvV6z59kJ7C/qm4ZuR5J\nMznsz/mTbAH+mMVL/o20d6JOaQn16fl/HjgDuCXJA8BpwN4kr12rcVXtrqodVbVjM8f3r1TSqA67\n56+q24DXvHC7+w9gR1X914h1SZpY34k6JR3l+k7UuXL96aNVI2k2XuEnNcrwS41KzTAh4P/vLHkc\n+I91Vp8MLNNJw2WrB5avJus5uCNRz89V1as30nDW8B9Mkj1VteNI1/GCZasHlq8m6zm4ZatnNV/2\nS40y/FKjlin8u490AassWz2wfDVZz8EtWz0vsjTv+SXNa5l6fkkzmj38Sc5P8r0k9yXZtcb6JPmr\nbv2tSc6ZsJbtSb6Z5M4kdyT5wBpt3p7kiST7up8/maqebn8PJLmt29eeNdbPdny6/b1xxb99X5In\nk1yxqs2kx2itb5NKclKSG5Lc2/1+5TrbHvT5NmI9H09yd/eYfCnJietse9DHd1ZVNdsPsAn4PvAG\n4DjgFuCsVW0uAK4HArwZuHnCel4HnNMtbwPuWaOetwNfnfEYPQCcfJD1sx2fdR6/R1h8ljzbMQLe\nBpwD3L7ivj8FdnXLu4CP9Xm+jVjPO4Bju+WPrVXPRh7fOX/m7vnPBe6rqvur6lngc8DOVW12An9f\nCzcBJyZ53RTFVNXDVbW3W34KuAs4dYp9jWi247OG84DvV9V6F2pNotb+NqmdwNXd8tXAu9bYdCPP\nt1HqqaqvV9WB7uZNLIa6L7W5w38q8IMVtx/ip8O2kTajS3I6cDZw8xqr39K9nLs+yS9MXEoB30jy\nne6LUFY7IsenczHw2XXWzXmMAE6pqoe75UeAU9Zoc6SO1XtYvDpby6Ee39k4Yw+Q5OXAF4ArqurJ\nVav3Aj9bVT9KcgHwZeDMCct5a1XtT/Ia4IYkd3c9zRGV5DjgIuDKNVbPfYxepKoqyVJ8bJXkw8AB\n4Jp1mizN4zt3z78f2L7i9mndfYfbZjRJNrMI/jVV9cXV66vqyar6Ubd8HbB5yq8qr6r93e/HgC+x\neOm60qzHZ4V3Anur6tHVK+Y+Rp1HX3i70/1+bI02cz+XLgMuBH6nujf4q23g8Z3N3OH/NnBmkjO6\nnuRi4NpVba4Ffq87q/1m4IkVL+9GlSTAVcBdVfWJddq8tmtHknNZHLMfTlTP1iTbXlhmcRJp9XwJ\nsx2fVS5hnZf8cx6jFa4FLu2WLwW+skabjTzfRpHkfBbfZn1RVT2zTpuNPL7zmfsMI4uz1fewOAv7\n4e6+9wHv65YDfLJbfxuLrwibqpa3sngPdiuwr/u5YFU97wfuYHGm+CbgLRPW84ZuP7d0+zyix2dF\nXVtZhPkVK+6b7Rix+E/nYeA5Fu/b3wu8CrgRuBf4BnBS1/b1wHUHe75NVM99LM4vvPA8+tTqetZ7\nfI/Uj1f4SY3yCj+pUYZfapThlxpl+KVGGX6pUYZfapThlxpl+KVG/R/AmPt3GvOQGQAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa57aec3b00>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"\n",
    "Eigen method for basis orthogonalization\n",
    "\"\"\"\n",
    "size = 25\n",
    "nbVec = 15\n",
    "V = np.random.rand(size,nbVec)\n",
    "\n",
    "#Compute Gram matrix\n",
    "Mg = np.dot(V.T,V)\n",
    "\n",
    "#Perform eigendecomposition\n",
    "D,Q=np.linalg.eig(Mg)\n",
    "\n",
    "#Compute Mg2 with the modified vectors V2\n",
    "V2 = np.dot(V,Q)\n",
    "Mg2 = np.dot(V2.T,V2)\n",
    "\n",
    "#check that Mg2 is orthonormal\n",
    "assert( np.allclose(Mg2-Mg2*np.identity(nbVec),np.zeros((nbVec,nbVec))) )\n",
    "\n",
    "#Show matrix Mg2\n",
    "plt.imshow(Mg2,interpolation='none')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fa5769c0630>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADG1JREFUeJzt3V+sZWV5x/Hvr8MgHaQCogjMtIMNIaGmKWRCKDXWlFaR\nEsaLXkBqC9Vk4oUtNCZmKEm91drYP6mpmYgtTQlcKCoxUBippmlSpg7j8HcQkFJgHP7UJqA2FaZ9\nerEXyeF0n5kze6+1zj68309yctbe6z1nPbP2/s2799prnSdVhaT2/NRaFyBpbRh+qVGGX2qU4Zca\nZfilRhl+qVGGX2qU4ZcaZfilRh035sZOO3VDbd2y8Zh/7rEHNg1QjfTG89/8mFfqJ1nN2FHDv3XL\nRv71ri3H/HPvP/OXBqhGeuPZU/eseqwv+6VGGX6pUXOFP8mlSb6b5IkkO/sqStLwZg5/kg3A54AP\nAOcBVyU5r6/CJA1rnpn/QuCJqnqyql4BbgW291OWpKHNE/6zgGeW3H62u0/SOjD4Ab8kO5LsTbL3\nxR/8z9Cbk7RK84T/ILD0Q/vN3X2vU1W7qmpbVW1721s3zLE5SX2aJ/zfBs5JcnaS44Ergdv7KUvS\n0GY+w6+qDif5GHAXsAH4YlU93FtlkgY11+m9VXUHcEdPtUgakWf4SY0y/FKjRr2q77EHNs10hd5d\n398/8za9IlCazplfapThlxpl+KVGGX6pUYZfapThlxpl+KVGGX6pUYZfapThlxpl+KVGGX6pUYZf\natSoV/XNap4r87wiUJrOmV9qlOGXGmX4pUbN06tvS5JvJnkkycNJru2zMEnDmueA32Hg41W1L8lJ\nwH1JdlfVIz3VJmlAM8/8VXWoqvZ1yz8EDmCvPmnd6OU9f5KtwPnAnj5+n6Thzf05f5I3A18Grquq\nl6es3wHsADiBTfNuTlJP5pr5k2xkEvybq+q2aWOWNurcyJvm2ZykHs1ztD/AjcCBqvpsfyVJGsM8\nM/+vAL8D/FqS/d3XZT3VJWlg83Tp/WcgPdYiaUSe4Sc1yvBLjVoXl/TOYy0uB/ZSYK0HzvxSowy/\n1CjDLzXK8EuNMvxSowy/1CjDLzXK8EuNMvxSowy/1CjDLzXK8EuNMvxSo97wV/XNY9ar82wOqvXA\nmV9qlOGXGmX4pUbNHf4kG5J8J8nX+yhI0jj6mPmvZdKnT9I6Mm/Hns3AbwJf6KccSWOZd+b/c+AT\nwP/2UIukEc3Truty4IWquu8o43Yk2Ztk76v8ZNbNSerZvO26rkjyFHArk7Zdf798kI06pcU0c/ir\n6vqq2lxVW4ErgX+sqg/1VpmkQfk5v9SoXs7tr6pvAd/q43dJGoczv9Qowy81ykt6B2BzUK0HzvxS\nowy/1CjDLzXK8EuNMvxSowy/1CjDLzXK8EuNMvxSowy/1CjDLzXK8EuNMvxSo7yqb8HYHFRjceaX\nGmX4pUYZfqlR87brOjnJl5I8muRAkl/uqzBJw5r3gN9fAP9QVb+V5HhgUw81SRrBzOFP8hbgPcA1\nAFX1CvBKP2VJGto8L/vPBl4E/ibJd5J8IcmJPdUlaWDzhP844ALgr6vqfODHwM7lg2zUKS2mecL/\nLPBsVe3pbn+JyX8Gr2OjTmkxzdOo8zngmSTndnddAjzSS1WSBjfv0f7fB27ujvQ/Cfze/CVJGsNc\n4a+q/cC2nmqRNCLP8JMaZfilRnlJ7xuEzUF1rJz5pUYZfqlRhl9qlOGXGmX4pUYZfqlRhl9qlOGX\nGmX4pUYZfqlRhl9qlOGXGmX4pUZ5VZ9sDtooZ36pUYZfapThlxo1b6POP0zycJKHktyS5IS+CpM0\nrJnDn+Qs4A+AbVX1LmADcGVfhUka1rwv+48DfjrJcUw69H5//pIkjWGejj0HgT8FngYOAS9V1d19\nFSZpWPO87D8F2M6kW++ZwIlJPjRlnI06pQU0z8v+Xwf+raperKpXgduAi5cPslGntJjmCf/TwEVJ\nNiUJk0adB/opS9LQ5nnPv4dJW+59wIPd79rVU12SBjZvo85PAp/sqRZJI/IMP6lRhl9qlJf0amY2\nB13fnPmlRhl+qVGGX2qU4ZcaZfilRhl+qVGGX2qU4ZcaZfilRhl+qVGGX2qU4ZcaZfilRnlVn9aE\nzUHXnjO/1CjDLzXK8EuNOmr4k3wxyQtJHlpy36lJdid5vPt+yrBlSurbamb+vwUuXXbfTuCeqjoH\nuKe7LWkdOWr4q+qfgP9cdvd24KZu+Sbggz3XJWlgs77nP72qDnXLzwGn91SPpJHMfcCvqgqoldbb\nqFNaTLOG//kkZwB0319YaaCNOqXFNGv4bweu7pavBr7WTzmSxrKaj/puAf4FODfJs0k+AnwK+I0k\njzNp1f2pYcuU1LejnttfVVetsOqSnmuRNCLP8JMaZfilRnlJr9YVm4P2x5lfapThlxpl+KVGGX6p\nUYZfapThlxpl+KVGGX6pUYZfapThlxpl+KVGGX6pUYZfapRX9akZNgd9PWd+qVGGX2qU4ZcaNWuj\nzs8keTTJA0m+kuTkYcuU1LdZG3XuBt5VVb8IPAZc33NdkgY2U6POqrq7qg53N+8FNg9Qm6QB9fGe\n/8PAnT38Hkkjmutz/iQ3AIeBm48wZgewA+AENs2zOUk9mjn8Sa4BLgcu6Tr1TlVVu4BdAD+TU1cc\nJ2lcM4U/yaXAJ4Bfrar/6rckSWOYtVHnXwEnAbuT7E/y+YHrlNSzWRt13jhALZJG5Bl+UqMMv9Qo\nL+mVjuKN2hzUmV9qlOGXGmX4pUYZfqlRhl9qlOGXGmX4pUYZfqlRhl9qlOGXGmX4pUYZfqlRhl9q\nlFf1SQMauznohe9f/V/Vc+aXGmX4pUYZfqlRMzXqXLLu40kqyWnDlCdpKLM26iTJFuB9wNM91yRp\nBDM16uz8GZPGHXbhkdahmd7zJ9kOHKyq+3uuR9JIjvlz/iSbgD9i8pJ/NeNt1CktoFlm/p8Hzgbu\nT/IUsBnYl+Qd0wZX1a6q2lZV2zbyptkrldSrY575q+pB4O2v3e7+A9hWVf/RY12SBjZro05J69ys\njTqXrt/aWzWSRuMZflKjDL/UqFSNd45OkheBf19h9WnAIh00XLR6YPFqsp4jW4t6fq6q3raagaOG\n/0iS7K2qbWtdx2sWrR5YvJqs58gWrZ7lfNkvNcrwS41apPDvWusCllm0emDxarKeI1u0el5nYd7z\nSxrXIs38kkY0eviTXJrku0meSLJzyvok+ctu/QNJLhiwli1JvpnkkSQPJ7l2ypj3Jnkpyf7u64+H\nqqfb3lNJHuy2tXfK+tH2T7e9c5f82/cneTnJdcvGDLqPpv01qSSnJtmd5PHu+ykr/OwRn2891vOZ\nJI92j8lXkpy8ws8e8fEdVVWN9gVsAL4HvBM4HrgfOG/ZmMuAO4EAFwF7BqznDOCCbvkk4LEp9bwX\n+PqI++gp4LQjrB9t/6zw+D3H5LPk0fYR8B7gAuChJff9CbCzW94JfHqW51uP9bwPOK5b/vS0elbz\n+I75NfbMfyHwRFU9WVWvALcC25eN2Q78XU3cC5yc5IwhiqmqQ1W1r1v+IXAAOGuIbfVotP0zxSXA\n96pqpRO1BlHT/5rUduCmbvkm4INTfnQ1z7de6qmqu6vqcHfzXiaXui+0scN/FvDMktvP8v/Dtpox\nvUuyFTgf2DNl9cXdy7k7k/zCwKUU8I0k93V/CGW5Ndk/nSuBW1ZYN+Y+Aji9qg51y88Bp08Zs1b7\n6sNMXp1Nc7THdzR27AGSvBn4MnBdVb28bPU+4Ger6kdJLgO+CpwzYDnvrqqDSd4O7E7yaDfTrKkk\nxwNXANdPWT32PnqdqqokC/GxVZIbgMPAzSsMWZjHd+yZ/yCwZcntzd19xzqmN0k2Mgn+zVV12/L1\nVfVyVf2oW74D2DjknyqvqoPd9xeArzB56brUqPtniQ8A+6rq+eUrxt5Hnedfe7vTfX9hypixn0vX\nAJcDv13dG/zlVvH4jmbs8H8bOCfJ2d1MciVw+7IxtwO/2x3Vvgh4acnLu14lCXAjcKCqPrvCmHd0\n40hyIZN99oOB6jkxyUmvLTM5iLS8X8Jo+2eZq1jhJf+Y+2iJ24Gru+Wrga9NGbOa51svklzK5K9Z\nX1FVUxvmrfLxHc/YRxiZHK1+jMlR2Bu6+z4KfLRbDvC5bv2DTP5E2FC1vJvJe7AHgP3d12XL6vkY\n8DCTI8X3AhcPWM87u+3c321zTffPkrpOZBLmtyy5b7R9xOQ/nUPAq0zet38EeCtwD/A48A3g1G7s\nmcAdR3q+DVTPE0yOL7z2PPr88npWenzX6ssz/KRGeYaf1CjDLzXK8EuNMvxSowy/1CjDLzXK8EuN\nMvxSo/4PZNLxa6+SdLoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa577dc6748>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"\n",
    "Choleski method for basis orthonormalization\n",
    "\"\"\"\n",
    "size = 25\n",
    "nbVec = 15\n",
    "V = np.random.rand(size,nbVec)\n",
    "\n",
    "#Compute Gram matrix\n",
    "Mg = np.dot(V.T,V)\n",
    "\n",
    "#Perform Choleski factorization\n",
    "L=np.linalg.cholesky(Mg)\n",
    "\n",
    "#Compute Mg2 with the modified vectors V2\n",
    "V2 = np.dot(V,np.linalg.inv(L.T))\n",
    "Mg2 = np.dot(V2.T,V2)\n",
    "\n",
    "#check that Mg2 is orthonormal\n",
    "assert( np.allclose(Mg2-np.identity(nbVec),np.zeros((nbVec,nbVec))) )\n",
    "\n",
    "#Show matrix Mg2\n",
    "plt.imshow(Mg2,interpolation='none')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remark\n",
    "It should be noticed that a Gram matrix is always symmetric positive semi definite, and as such can be diagonalized in an orthonormal basis, which makes $Q^t$ the inverse of $Q$.\n",
    "\n",
    "In the Choleski case, it should be noticed that computing the inverse of a lower triangular matrix is a low complexity operation. Demonstration sketch is as follow:\n",
    "Let $L$ be a lower triangular matrix of size n, we have $L=D(I+N)$, where D is a diagonal matrix, I is the identity matrix, and N is a strictly lower triangular matrix.\n",
    "As such, N is n-nilpotent, ie $N^n=0$, and, using the taylor development:\n",
    "\n",
    "\\begin{align*}\n",
    "  (I-N)^{-1} &= \\sum_{i=0}^{n-1} N^i \\\\\n",
    "  (I+N)^{-1} &= \\sum_{i=0}^{n-1} (-1)^i N^i\n",
    "\\end{align*}\n",
    "\n",
    "We can compute $L^{-1}=(I+N)^{-1} D^{-1}$ quite easily with matrix exponentiation, or just by noticing that inverse of a upper/lower triangular matrix $L$ can be computed by forward/backward subsitution.\n",
    "In the framework of high performance computing, one may be interested in using dtrtri, the level 3 BLAS Triangular Matrix Inversion (TMI)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### General notes on whitening\n",
    "It might be noticed that if you find a $W_W$ such that the following holds:\n",
    "\\begin{align*}\n",
    "  W_W X (W_W X)^T &= W_W X X^T W_W^T \\\\\n",
    "  &= W_W C W_W^T \\\\\n",
    "  &= I\n",
    "\\end{align*}\n",
    "Then, the property will hold for any $W_W'$ such that $W_W' = A W_W$ with $A$ an arbitrary matrix with unitary property, ie $A^{*}A=I$\n",
    "There might then be multiple choice for the final whitening matrix, and one have some additional freedom to choose the whitenning transform such that it exhibit different sort of optimality, such as zca like transform for instance, that optimize for average cross-covariance or average cross-correlation between original and transformed data.\n",
    "\n",
    "More on this in [Optimal whitening and decorrelation](https://arxiv.org/abs/1512.00809)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import IFrame\n",
    "#IFrame(\"doc/Orthonormalization/optimal_whitening_decorrelation.pdf\", width=1200, height=800)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Geometric interpretation of orthonormalization transforms\n",
    "It is interesting to notice that there is a strong link between orthonormalization transforms, and ellipsoid, as, most on the time, one is interested in factorizing either gram or covariance matrix, which are both hermitian matrices (positive semidefinite )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generalizes of orthogonal decomposition from linear spaces to convex cones: Moreau's decomposition\n",
    "\n",
    "Main source for this analysis can be found there:\n",
    "https://www.convexoptimization.com/wikimization/index.php/Convex_cones\n",
    "https://www.convexoptimization.com/wikimization/index.php/Moreau's_decomposition_theorem"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}