ldl

sksparse.cholmod.ldl(A, beta=0.0, *, lower=True, order='default', sym_kind=None, supernodal_mode=None)[source]

Compute the LDL factorization of a sparse matrix.

This function computes the LDL factorization of a symmetric matrix A:

\[L D L^{\top} = P A P^{\top},\]

where L is a lower triangular matrix with unit diagonal, and D is a diagonal matrix. Only the lower triangular part of A is used. If lower is False, the upper triangular factor R is returned instead, such that:

\[R^{\top} D R = P A P^{\top}.\]

In this case, only the upper 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. 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.

    • natural or None: No permutation is applied (identity permutation).

    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.

  • D (dia_array) – The diagonal matrix D of the factorization, in sparse DIA format. 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

ldl_factor, cholesky, cho_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.5.0.

References

Examples

>>> import numpy as np
>>> from scipy.sparse import coo_array
>>> from sksparse.cholmod import ldl, ldl_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, D, p = ldl(A, order='amd')
>>> L
<Compressed Sparse Column sparse array of dtype 'float64'
        with 30 stored elements and shape (11, 11)>
>>> D
<DIAgonal sparse array of dtype 'float64'
        with 11 stored elements (1 diagonals) and shape (11, 11)>
>>> p
array([ 4,  8,  6,  0,  3,  5,  1,  2,  9, 10,  7])
>>> f = ldl_factor(A, order='amd')
>>> f
CholeskyFactor(N=11, nnz=30, is_ll=False, is_super=False, itype=np.int64,
    dtype=np.float64, order=amd)
>>> Lf, Df = f.get_factor()
>>> np.allclose(L.toarray(), Lf.toarray(), atol=1e-15)
True
>>> np.allclose(D.toarray(), Df.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