[Python, Scipy] Sparse Cholesky decomposition

SparseCholesky.md

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

Thanks for posting this. It's an interesting idea. However, it seems to me that LU.perm_r == np.arange(n) indicates we're missing the mark a little. A major motivation for sparse Cholesky factorizations is to find a permutation P such that PAP^T = GG^T, where G is sparse. This way we don't lose our sparsity when solving systems with G et cetera. Here you seek to get the Cholesky factorization through the LU (LDL) facotrization: PAP^T = LU = LDL^T. However if we're assuming P = I (which is implied by LU.perm_r == np.arange(n)) then there's a good chance L is dense. Furthermoe, if the column permutation matrix is not the identity matrix, then I don't think this will return the actual Cholesky factorization. For starters,

To recap: this will (probably) get you the Cholesky factor of a sparse SPD matrix; however, I think the Cholesky factor will be dense.

Thanks for sharing the code, but I find that when I want to use the lower triangular matrix ,L produced by Cholesky decomposition, I can't get the correct matrix L* because Permutation matrix P. This is my solution to get the correct matrix L*. Just add the parameter permc_spec = "NATURAL" in the sparse.linalg.splu function.

def sparse_cholesky(A):
  sparse_matrix = A.T @ A
  sparse_matrix += 1e-6 * sparse.identity(sparse_matrix.shape[0]) # force the sparse matrix is positive definite
  n = sparse_matrix.shape[0]
  LU = sparse.linalg.splu(sparse_matrix, diag_pivot_thresh = 0.0, permc_spec = "NATURAL") # sparse LU decomposition

  L = LU.L @ sparse.diags(LU.U.diagonal()**0.5)
  
  return L # return L (lower triangular metrix)