Unsymmetric Multifrontal LU Decomposition (`sksparse.umfpack`)¶

The sksparse.umfpack module provides an interface to the SuiteSparse UMFPACK package, which computes basic linear algebra operations for sparse matrices.

The main function of this module is to compute the LU factorization of a sparse matrix \(A\) with a fill-reducing permutation \(P\) and \(Q\), and row-scaling \(R\), such that:

\[LU = PRAQ.\]

These factors can then be used to solve linear systems of the form \(Ax = b\), or \(A^{\top} x = b\).

The umfpack module exposes most of the capabilities of the UMFPACK package including:

Computation of the LU factors 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\).
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.

Quickstart¶

If \(A\) is a sparse matrix then the following code computes the LU decomposition of \(A\):

from sksparse.umfpack import umf_factor
L, U, p, q, r = umf_factor(A)
# L @ U == (r[:, np.newaxis] * A).tocsc()[p[:, np.newaxis], q]

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 square linear system:

from sksparse.umfpack import umf_factor
A = ...  # a sparse matrix
b = ...  # right-hand side
x = umf_solve(A, b)

Example¶

To see how to use the LU 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: umfpack_example.py
#  Created: 2025-11-17 09:35
# =============================================================================

"""An example of using UMFPACK to compute the umf_factor factorization of a
sparse matrix.
"""

from pathlib import Path

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

from sksparse.umfpack import umf_factor

# 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

# compute the LU factorization
L, U, p, q, rscale = umf_factor(A, ordering_method="none")
Lp, Up, pp, qp, rscalep = umf_factor(A, ordering_method="amd")

# Make sure the factorization is correct
PRAQ = (rscale[:, np.newaxis] * A).tocsc()[p[:, np.newaxis], q]
PRAQp = (rscalep[:, np.newaxis] * A).tocsc()[pp[:, np.newaxis], qp]
assert_allclose((L @ U).toarray(), PRAQ.toarray(), atol=1e-9)
assert_allclose((Lp @ Up).toarray(), PRAQp.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("UMFPACK Example: LU Factors of west0479")

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

ax = axs[0, 1]
ax.spy(PRAQ, markersize=1)
ax.set_title(r"Permuted Matrix $PAQ$")

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

ax = axs[1, 1]
ax.spy(Lp + Up, markersize=1)
ax.set_title(r"Permuted LU Factors $L_p + U_p$")
ax.set_xlabel(f"{Lp.nnz + Up.nnz:,} non-zeros")

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

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

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

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

Function Interface¶

For users who want to directly compute the factorization without needing to manipulate the UMFFactor object, the umfpack module provides the umf_factor() function that performs both the symbolic analysis and the numeric factorization in one step. It returns the factorization object as a single output, or the individual factors and permutations can be unpacked as multiple outputs.

To directly solve a linear system without explicitly computing the factorization, the umf_solve() function can be used. This function performs the factorization internally and returns the solution to the linear system.

Object Interface¶

For more advanced usage, users can instantiate the UMFFactor class. This class can be instantiated directly using its constructor, or more conveniently using the umf_factor() function.

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

The umf_factor() function performs both the symbolic analysis and the numeric factorization in one step, and returns an instance of the UMFFactor class. The resulting UMFFactor object can then be used to solve linear systems using its UMFFactor.solve() method. The UMFFactor.factorize() method can be called again to factor a new matrix with the same sparsity pattern.

Customization and Statistics¶

The umf_factor() accepts several optional parameters that control the behavior of the factorization, such as the fill-reducing ordering method to use. These options are stored in a UMFControl object that can be accessed via the UMFFactor.control attribute. The UMFControl object can also be instantiated directly to customize the factorization options before passing it to the UMFFactor constructor or the umf_factor() function. Any options not explicitly set by the user will be given the UMFPACK default values.

After a factorization is computed, statistics about the factorization can be accessed via the UMFFactor.info attribute, which is an instance of the UMFInfo class. This object contains various performance metrics and statistics about the factorization process.

Exceptions and Warnings¶

Warnings issued by UMFPACK are converted into Python warnings of type UMFPACKWarning. 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 UMFPACK or by our wrapper code are converted into exceptions of type UMFPACKError or an appropriate subclass. See the Warnings and Exceptions for details.

Unsymmetric Multifrontal LU Decomposition (sksparse.umfpack)¶