Cholesky Decomposition (`sksparse.cholmod`)¶

The sksparse.cholmod module provides an interface to the SuiteSparse CHOLMOD package, which computes basic linear algebra operations for sparse, symmetric, positive-definite matrices.

The main function of this module is to compute the Cholesky factor \(L\) of a sparse, symmetric (Hermitian if complex), positive-definite matrix \(A\) with a fill-reducing permutation \(P\), such that:

\[LL^{\top} = PAP^{\top}.\]

For matrices that are symmetric but may be numerically close to semi-definite, the module can compute the LDL factorization:

\[LDL^{\top} = PAP^{\top}.\]

Either of these factors can then be used to solve linear systems of the form \(Ax = b\).

The cholmod module exposes most of the capabilities of the CHOLMOD package including:

Computation of the Cholesky factor with fill-reducing permutation for both real and complex sparse matrices \(A\), in any format supported by scipy.sparse.
An interface for using this decomposition to solve problems of the form \(Ax = b\).
Functions to compute a rank-\(k\) “update” or “downdate” of the Cholesky factor.
The ability to perform the fill-reduction analysis once, and then re-use it to efficiently decompose many matrices with the same pattern of non-zero entries.

This wrapper handles both 32-bit and 64-bit integer types, depending on the input matrix format.

Quickstart¶

If \(A\) is a sparse, symmetric, positive-definite matrix, then the following code computes the Cholesky decomposition of \(A\):

from sksparse.cholmod import cholesky
L = cholesky(A)

or, with a fill-reducing permutation:

L, p = cholesky(A, order="default")

See the example below for a demonstration of the effect of the fill-reducing permutation.

Once the factorization has been computed, it can be used to solve a linear system:

from sksparse.cholmod import ldl_factor
A = ...            # a symmetric, positive-definite sparse matrix
b = ...            # right-hand side
f = ldl_factor(A)  # compute LDL^T factorization
x = f.solve(b)     # solve Ax = b

Examples¶

Cholesky Example¶

To see how to use the Cholesky factorization, we can load a sparse matrix from the SuiteSparse Matrix Collection and compute its ordering.

# Part of the scikit-sparse project.
# Copyright (C) 2025 Bernard Roesler. All rights reserved.
# See pyproject.toml for full author list and LICENSE.txt for license details.
# SPDX-License-Identifier: BSD-2-Clause
#
# =============================================================================
#     File: cholmod_example.py
#  Created: 2025-08-12 20:12
# =============================================================================

"""An example of using CHOLMOD to compute the Cholesky factorization of a
sparse matrix.
"""

from pathlib import Path

import matplotlib.pyplot as plt
from numpy.testing import assert_allclose
from scipy.io import mmread

from sksparse.cholmod import cholesky

# Load the west0479 matrix (downloaded from the SuiteSparse Matrix Collection:
# <https://sparse.tamu.edu/HB/west0479>)
filepath = Path("data") / "west0479.mtx"
A = mmread(filepath, spmatrix=False)  # read the matrix
A = (A.T @ A).tocsc()  # make it symmetric positive definite

# compute the Cholesky factorization
R = cholesky(A)
Rp, p = cholesky(A, order="amd")

PAPT = A[p][:, p]  # apply the ordering to the matrix

# Make sure the factorization is correct
assert_allclose((R.T.conj() @ R).toarray(), A.toarray(), atol=1e-9)
assert_allclose((Rp.T.conj() @ Rp).toarray(), PAPT.toarray(), atol=1e-9)

# Plot the original and permuted matrices
plt.rcParams.update({"font.size": 10})

fig, axs = plt.subplots(num=1, nrows=2, ncols=2, clear=True)
fig.set_size_inches((6, 6), forward=True)
fig.set_constrained_layout(True)
fig.suptitle("CHOLMOD Example: Cholesky Factor of west0479")

ax = axs[0, 0]
ax.spy(A, markersize=1)
ax.set_title(r"Original Matrix $A$")

ax = axs[0, 1]
ax.spy(PAPT, markersize=1)
ax.set_title(r"Permuted Matrix $PAP^{\top}$")

ax = axs[1, 0]
ax.spy(R, markersize=1)
ax.set_title(r"Original Cholesky Factor $R$")
ax.set_xlabel(f"{R.nnz:,} non-zeros")

ax = axs[1, 1]
ax.spy(Rp, markersize=1)
ax.set_title(r"Permuted R Factor $R_p$")
ax.set_xlabel(f"{Rp.nnz:,} non-zeros")

for ax in axs.flat:
    ax.set_aspect("equal")
    ax.set_xticks([])
    ax.set_yticks([])

plt.show()
fig.savefig("cholesky_example.svg")

This figure shows the effect of AMD ordering that reduces the fill-in of the Cholesky factorization of a sparse matrix.

Cholesky Example with AMD Ordering — The number of non-zeros in the Cholesky factorization of the original matrix (left) and the permuted matrix (right) using AMD ordering.¶

Nested Dissection Example¶

To see the effects of nested dissection ordering, we can load a sparse matrix and compute its ordering.

# Part of the scikit-sparse project.
# Copyright (C) 2025 Bernard Roesler. All rights reserved.
# See pyproject.toml for full author list and LICENSE.txt for license details.
# SPDX-License-Identifier: BSD-2-Clause
#
# =============================================================================
#     File: nesdis_example.py
#  Created: 2025-08-22 11:04
# =============================================================================

"""An example of using CHOLMOD to compute the Cholesky factorization of a
sparse matrix with various orderings.
"""

from pathlib import Path

import matplotlib.pyplot as plt
from scipy.io import mmread
from scipy.sparse.linalg import splu

from sksparse.amd import amd
from sksparse.cholmod import nesdis

# Load the west0479 matrix (downloaded from the SuiteSparse Matrix Collection:
# <https://sparse.tamu.edu/HB/west0479>)
filepath = Path("data") / "west0479.mtx"
A = mmread(filepath, spmatrix=False).tocsc()  # read the matrix

# Compute the AMD permutation
p = amd(A)
PAPT = A[p][:, p]

# Compute the nested dissection A.T @ A permutation
q = nesdis(A)
QAQT = A[q][:, q]

# Compute the LU decompositions
lu = splu(A)
lup = splu(PAPT)
luq = splu(QAQT)

LU = lu.L + lu.U
LUp = lup.L + lup.U
LUq = luq.L + luq.U

# Plot the original and permuted matrices
plt.rcParams.update({"font.size": 10})

fig, axs = plt.subplots(num=1, nrows=3, ncols=2, clear=True)
fig.set_size_inches((6, 10), forward=True)
fig.set_constrained_layout(True)
fig.suptitle("CHOLMOD Example: Cholesky Factor of west0479")

ax = axs[0, 0]
ax.spy(A, markersize=1)
ax.set_title(r"Original Matrix $A$")

ax = axs[1, 0]
ax.spy(PAPT, markersize=1)
ax.set_title(r"AMD-Ordered Matrix $PAP^{\top}$")

ax = axs[0, 1]
ax.spy(LU, markersize=1)
ax.set_title(r"LU Factors $LU = A$")
ax.set_xlabel(f"{LU.nnz:,} non-zeros")

ax = axs[1, 1]
ax.spy(LUp, markersize=1)
ax.set_title(r"AMD Factors $LU = PAP^{\top}$")
ax.set_xlabel(f"{LUp.nnz:,} non-zeros")

ax = axs[2, 0]
ax.spy(QAQT, markersize=1)
ax.set_title(r"Nested Dissection Matrix $QAQ^{\top}$")

ax = axs[2, 1]
ax.spy(LUq, markersize=1)
ax.set_title(r"Nesdis Factors $LU = QAQ^{\top}$")
ax.set_xlabel(f"{LUq.nnz:,} non-zeros")

for ax in axs.flat:
    ax.set_aspect("equal")
    ax.set_xticks([])
    ax.set_yticks([])

plt.show()
fig.savefig("nesdis_example.svg")

This figure shows the effect of nested dissection ordering that reduces the fill-in of the LU factorization of a sparse matrix in a case where the AMD order does not help.

LU Example with Nesdis Ordering — The number of non-zeros in the LU factorization of the original matrix and the permuted matrix using AMD and nested dissection ordering.¶

Function Interface¶

For users who want to directly compute the factorization without needing to manipulate the CholeskyFactor object, the cholmod module provides the cholesky() and ldl() functions that perform both the symbolic analysis and the numerical factorization in one step, and return the matrices directly.

Object Interface¶

For more advanced usage, users can instantiate the CholeskyFactor class. This class can be instantiated directly using its constructor, or more conveniently using the cho_factor() or ldl_factor() functions.

When instantiated directly, the constructor performs a symbolic analysis of the matrix, but does not compute the numerical factorization. The numerical factorization is then performed by calling the CholeskyFactor.factorize() method.

The cho_factor() and ldl_factor() functions perform both the symbolic analysis and the numerical factorization in one step, and return an instance of the CholeskyFactor class.

The resulting CholeskyFactor object can then be used to solve linear systems using its CholeskyFactor.solve() method, or to update the factorization in-place using the update(), rowadd(), rowdel(), and resymbol() methods.

The CholeskyFactor.factorize() method can be called again to factor a new matrix with the same sparsity pattern.

Symbolic Analysis¶

In addition to numerical factorization, cholmod provides symbolic operations symbfact(), and etree() that can be used to analyze the structure of the Cholesky factor and to compute fill-reducing permutations.

Graph Partitioning¶

The cholmod module also includes functions for graph partitioning and node reordering, which can be used like the amd and colamd (and their constrained counterparts) modules to reduce fill-in during factorization.

These functions provide a direct interface to the corresponding CHOLMOD functions that are used internally by cholesky() and ldl() when the order argument is specified. The functions bisect(), metis(), and nesdis() can be used to compute fill-reducing orderings, and the SeparatorTree class represents the resulting separator tree.

Exceptions and Warnings¶

Warnings issued by CHOLMOD are converted into Python warnings of type CholmodWarning. The module will also issue a SparseEfficiencyWarning if the input matrix is not a csc_array (note that the csc_matrix class is not supported, as it will be deprecated and is not recommended for use in new code).

Errors detected by CHOLMOD or by our wrapper code are converted into exceptions of type CholmodError or an appropriate subclass. See the Exceptions and Warnings for details.

Cholesky Decomposition (sksparse.cholmod)¶