Scipy does not currently provide a routine for cholesky decomposition of a sparse matrix, and one have to rely on another external package such as scikit.sparse for the purpose. Here I implement cholesky decomposition of a sparse matrix only using scipy functions. Our implementation relies on sparse LU deconposition.
The following function receives a sparse symmetric positive-definite matrix A and returns a spase lower triangular matrix L such that A = LL^T.
from scipy.sparse import linalg as splinalg
import scipy.sparse as sparse
import sys
def sparse_cholesky(A): # The input matrix A must be a sparse symmetric positive-definite.
n = A.shape[0]
LU = splinalg.splu(A,diag_pivot_thresh=0) # sparse LU decomposition
if ( LU.perm_r == np.arange(n) ).all() and ( LU.U.diagonal() > 0 ).all(): # check the matrix A is positive definite.
return LU.L.dot( sparse.diags(LU.U.diagonal()**0.5) )
else:
sys.exit('The matrix is not positive definite')
Helpful documentation for the SuperLU object returned from splu:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.SuperLU.html