CHOLMOD_Test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import division\n",
    "\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.sparse.linalg\n",
    "from scipy.sparse import SparseEfficiencyWarning\n",
    "import warnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dense cholesky\n",
    "cholesky_dense = lambda x: sp.linalg.cholesky(x, lower=True)\n",
    "\n",
    "# sparse cholesky\n",
    "try:\n",
    "    from sksparse.cholmod import cholesky as cholesky_sparse\n",
    "except ImportError: # older versions\n",
    "    from scikits.sparse.cholmod import cholesky as cholesky_sparse\n",
    "    print('old scikt-sparse')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# configs to try: seed=12, density=0.01, mul=0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "density:  0.049594\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<1000x1000 sparse matrix of type '<class 'numpy.float64'>'\n",
       "\twith 49594 stored elements in Compressed Sparse Column format>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seed = 12\n",
    "mul1 = 0.1 # lower values ensure positive definitness\n",
    "mul2 = 0.5\n",
    "\n",
    "F = sp.sparse.rand(1000, 5, density=0.1, format='csc', random_state=seed)\n",
    "S = F.dot(F.T)\n",
    "S *= mul1\n",
    "S = S+S.T\n",
    "\n",
    "with warnings.catch_warnings():\n",
    "    warnings.simplefilter('ignore', category=SparseEfficiencyWarning)\n",
    "    S.setdiag(1 + mul2*np.random.RandomState(seed).rand(S.shape[0])) # keep 1's on diagonal\n",
    "print('density: ', S.data.size / np.prod(S.shape))\n",
    "S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.74522417])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.sparse.linalg.eigsh(S, k=1, return_eigenvectors=False, which='SM')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "L_dense = cholesky_dense(S.A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "L_sparse = cholesky_sparse(S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# check that $LL^{T} = S$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "numpy/scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.9320068704658742e-16"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(S.A - L_dense.dot(L_dense.T)) / np.linalg.norm(S.A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "cholmod"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "cholmod's representation is $LL^T = PSP^T$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.324570115030523e-16"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PtL = L_sparse.apply_Pt(L_sparse.L())\n",
    "np.linalg.norm(S.A - PtL.dot(PtL.T)) / np.linalg.norm(S.A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# check that $LL^Tx = y$ is equiv. to $Sx = y$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "using the fact:\n",
    "$LL^Tx = y$ $\\rightarrow$\n",
    "$\n",
    "\\begin{cases}\n",
    "  L\\omega &= y \\\\\n",
    "  L^Tx    &= \\omega \\\\\n",
    "\\end{cases}\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.RandomState(seed).rand(S.shape[0],)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = S.dot(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## scipy dense"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_dense = sp.linalg.solve_triangular(L_dense, y, lower=True, trans=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "solx_dense = sp.linalg.solve_triangular(L_dense, w_dense, lower=True, trans=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.9004056755538052e-15"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(x - solx_dense) / np.linalg.norm(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## cholmod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_sparse = L_sparse.solve_L(L_sparse.apply_P(y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "solx_sparse = L_sparse.apply_Pt(L_sparse.solve_Lt(w_sparse))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.43370385938462436"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(x - solx_sparse) / np.linalg.norm(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**NOT OK**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "reason: CHOLMOD solves in the form:\n",
    "$LDL^{T} = S$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CHOLMOD solves in the form:\n",
    "$LL^Tx = y$ $\\rightarrow$\n",
    "$\n",
    "\\begin{cases}\n",
    "  L\\omega &= y \\\\\n",
    "  L^Tx    &= \\omega \\\\\n",
    "\\end{cases}\n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## chomod v.0.4.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "new param `use_LDLt_decomposition` (see https://github.com/scikit-sparse/scikit-sparse/issues/43 and https://github.com/scikit-sparse/scikit-sparse/issues/9#issuecomment-326213771)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_sparse = L_sparse.solve_L(L_sparse.apply_P(y), use_LDLt_decomposition=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "solx_sparse = L_sparse.apply_Pt(L_sparse.solve_Lt(w_sparse, use_LDLt_decomposition=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.8111110836750037e-15"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(x - solx_sparse) / np.linalg.norm(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ok"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "also can use natural ordering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "L_sparse_nat = cholesky_sparse(S, ordering_method='natural')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_sparse_nat = L_sparse_nat.solve_L(y, use_LDLt_decomposition=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "solx_sparse_nat = L_sparse_nat.solve_Lt(w_sparse_nat, use_LDLt_decomposition=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0160427115323124e-15"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(x - solx_sparse_nat) / np.linalg.norm(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ok"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## chomod older v.0.4.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "use in the form    \n",
    "$\n",
    "\\begin{cases}\n",
    "  LD\\omega &= y \\\\\n",
    "  L^Tx    &= \\omega \\\\\n",
    "\\end{cases}\n",
    "$\n",
    " or $\\,$\n",
    "$\n",
    "\\begin{cases}\n",
    "  L\\omega &= y \\\\\n",
    "  DL^Tx    &= \\omega \\\\\n",
    "\\end{cases}\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_sparse_old = L_sparse.solve_LD(L_sparse.apply_P(y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "solx_sparse_old = L_sparse.apply_Pt(L_sparse.solve_Lt(w_sparse_old))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.8512584883722028e-15"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(x - solx_sparse_old) / np.linalg.norm(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ok"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_sparse_old = L_sparse.solve_L(L_sparse.apply_P(y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "solx_sparse_old = L_sparse.apply_Pt(L_sparse.solve_DLt(w_sparse_old))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.8512584883722028e-15"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(x - solx_sparse_old) / np.linalg.norm(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ok"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## scipy sparse (*replace cholmod's solve_L with scipy's spsolve_triangular*)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "with warnings.catch_warnings():\n",
    "    warnings.simplefilter('ignore', category=SparseEfficiencyWarning)\n",
    "    w_scipy = sp.sparse.linalg.spsolve_triangular(L_sparse.L(), L_sparse.apply_P(y), lower=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "solx_scipy = L_sparse.apply_Pt(sp.sparse.linalg.spsolve_triangular(L_sparse.L().T, w_scipy, lower=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6355303407984224e-15"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(x - solx_scipy) / np.linalg.norm(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## scipy vs cholmod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.6450588665930786e-16"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(w_sparse-w_scipy) / np.linalg.norm(w_scipy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.389403919881755e-16"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(solx_sparse-solx_scipy) / np.linalg.norm(solx_scipy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1459431697378629"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(w_sparse_old-w_scipy) / np.linalg.norm(w_scipy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.904698267522099e-16"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(solx_sparse_old-solx_scipy) / np.linalg.norm(solx_scipy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9330963700494369"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(w_sparse_nat-w_scipy) / np.linalg.norm(w_scipy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.3533247848572229e-15"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(solx_sparse_nat-solx_scipy) / np.linalg.norm(solx_scipy)"
   ]
  }
 ],
 "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.6.7"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "135px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}