cholesky¶

sksparse.cholmod.cholesky(A, beta=0.0, *, lower=False, order=None, sym_kind=None, supernodal_mode=None)[source]¶

Compute the Cholesky factorization of a sparse matrix.

This function computes the Cholesky factorization of a symmetric positive definite matrix A:

\[R^{\top} R = P A P^{\top},\]

where R is an upper triangular matrix. Only the upper triangular part of A is used. If lower is True, the lower triangular factor L is returned instead, such that:

\[L L^{\top} = P A P^{\top}.\]

In this case, only the lower triangular part of A is used.

If beta is a scalar value, compute the factorization of:

\[P A P^{\top} + \beta I,\]

where I is the identity matrix.

Parameters:

A ((N, N) {array_like, sparse array}) – An array convertible to a sparse matrix in Compressed Sparse Column (CSC) format. The matrix must be square and symmetric positive definite. Only the upper or lower triangular part of the matrix is used, and no check is made for symmetry.
beta (float, optional) – The scalar value to add to the diagonal of the matrix before factorization.
lower (bool, optional) – If True, return the lower triangular factor L, otherwise return the upper triangular factor R.
order (None or str in {“default”, “best”, “natural”, “metis”, “nesdis”, “amd”, “colamd”, “postordered”}, optional) – The permutation algorithm to use for the factorization. By default, the natural ordering of the input matrix is used. The other options are:
- default: Use the default method, which first tries AMD, then METIS.
- best: Automatically select the best ordering based on the input.
- metis: Use the METIS library for graph partitioning.
- nesdis: Use the NESDIS library for nested dissection.
- amd: Use the Approximate Minimum Degree (AMD) algorithm.
- colamd: Use the Approximate Minimum Degree (AMD) algorithm for the
  symmetric case, or the COLAMD algorithm for the unsymmetric case (\(A A^{\top}\) or \(A^{\top} A\)).
- postordered: Use natural ordering followed by postordering.
By default, methods other than natural will also be postordered.

Warning

The ordering method best may be quite slow for large matrices, but if the factorization is reused many times, it can be worth it.
sym_kind (str in {“sym”, “row”, “col”}, optional) – The type of factorization for which to analyze the matrix:
- sym: Symmetric factorization. No check is made for symmetry.
- row: Unsymmetric factorization of \(A A^{\top}\).
- col: Unsymmetric factorization of \(A^{\top} A\).
supernodal_mode (str in {“auto”, “simplicial”, “supernodal”}, optional) – The type of factorization to use:
- auto: Automatically select the factorization type.
- simplicial: Use a simplicial factorization.
- supernodal: Use a supernodal factorization.
Note that the simplicial mode may be slow for large matrices.

Returns:

R (csc_array) – The triangular factor of the Cholesky decomposition. The data type will match that of A.
p (ndarray of int, optional) – The permutation vector used in the factorization. Only returned if the ordering is not None.

Raises:

CholmodNotPositiveDefiniteError – If the input matrix is not positive definite.

See also

cho_factor, ldl, ldl_factor

Notes

This function is an interface to the CHOLMOD library, which is part of the SuiteSparse collection by Timothy A. Davis. For more details, see the documentation in the header file [1].

Added in version 0.1.0.

Changed in version 0.5.0: The function now returns the matrix directly instead of a Factor object, and the permutation vector when an ordering method is specified.

References

Examples

>>> import numpy as np
>>> from scipy.sparse import coo_array
>>> from sksparse.cholmod import cholesky, cho_factor
>>> # Create a symmetric positive definite matrix from (Davis, Eqn 2.1)
>>> N = 11
>>> rows = np.array([5, 6, 2, 7, 9, 10, 5, 9, 7, 10, 8, 9, 10, 9, 10, 10])
>>> cols = np.array([0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 7, 7, 9])
>>> rng = np.random.default_rng(56)
>>> vals = rng.random(len(rows), dtype=np.float64)
>>> L = coo_array((vals, (rows, cols)), shape=(N, N))
>>> A = L + L.T   # make it symmetric
>>> A.setdiag(N)  # make it strongly positive definite
>>> A = A.tocsc()
>>> L, p = cholesky(A, order='amd', lower=True)
>>> L
<Compressed Sparse Column sparse array of dtype 'float64'
        with 30 stored elements and shape (11, 11)>
>>> p
array([ 4,  8,  6,  0,  3,  5,  1,  2,  9, 10,  7])
>>> f = cho_factor(A, order='amd', lower=True)
>>> f
CholeskyFactor(N=11, nnz=30, is_ll=True, is_super=False, itype=np.int64,
    dtype=np.float64, order=natural)
>>> np.allclose(L.toarray(), f.get_factor().toarray(), atol=1e-15)
True
>>> np.array_equal(p, f.get_perm())
True
>>> # Solve a linear system
>>> expect_x = np.arange(N, dtype=np.float64)
>>> b = A @ expect_x
>>> x = f.solve(b)
>>> np.allclose(x, expect_x)
True