Sparse QR Decomposition (`sksparse.spqr`)¶

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

The main function of this module is to compute the QR factorization of a sparse matrix \(A\) with a fill-reducing permutation \(E\) such that:

\[QR = AE\]

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

The spqr module exposes most of the capabilities of the SPQR package including:

Computation of the QR 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 QR decomposition of \(A\):

from sksparse.spqr import spqr_factor
Q, R, p = spqr(A)
# Q @ R == A[:, p]

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.spqr import spqr_factor
A = ...  # a sparse matrix
b = ...  # right-hand side
x = spqr_solve(A, b)

Example¶

To see how to use the QR 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: spqr_example.py
#  Created: 2025-11-17 11:58
# =============================================================================

"""An example of using SPQR to compute the QR 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.spqr import spqr

# 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.tocsc()

# compute the LU factorization
Q, R, _ = spqr(A, order="fixed")
Qp, Rp, p = spqr(A, order="colamd")

# Permute A for visualization
AE = A[:, p]

# Make sure the factorization is correct
assert_allclose((Q @ R).toarray(), A.toarray(), atol=1e-9)
assert_allclose((Qp @ Rp).toarray(), AE.toarray(), atol=1e-9)

# 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("SPQR Example: QR 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(AE, markersize=1)
ax.set_title(r"Permuted Matrix $AE$")

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

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

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

ax = axs[2, 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("spqr_example.svg")

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

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

Function Interface¶

For users who just need to compute the QR factorization and return the factors as sparse matrices, the spqr() function can be used. This function takes a sparse matrix as input and returns the \(Q\) and \(R\) factors as sparse matrices, along with the fill-reducing column permutation vector. There is also a mode argument that is similar to the one in scipy.linalg.qr(), which allows users to specify whether they want the full or reduced factors, just the \(R\) factor, or the Householder form of the \(Q\) factor.

Typically, the Householder form is substantially less memory-intensive than the explicit \(Q\) factor, especially for large sparse matrices. The function spqr_qmult() can then be used to efficiently apply the \(Q\) factor to a dense matrix or vector without explicitly forming \(Q\).

To directly solve a linear system without explicitly computing the factorization, the spqr_solve() function can be used. This function performs the factorization internally and returns the solution to the linear system. It can be used on non-square matrices as well. If \(A\) is a matrix of shape \(M\) by \(N\), solving an overdetermined system (\(M > N\)) will return the least-squares solution, while solving an underdetermined system (\(M < N\)) will return the minimum-norm solution.

Object Interface¶

For more advanced usage, users can instantiate the SPQRFactor class. This class can be instantiated directly using its constructor, or more conveniently using the spqr_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 SPQRFactor.factorize() method.

The spqr_factor() function performs both the symbolic analysis and the numeric factorization in one step, and returns an instance of the SPQRFactor class. The resulting SPQRFactor object can then be used to solve linear systems using its SPQRFactor.solve() method. The SPQRFactor.factorize() method can be called again to factor a new matrix with the same sparsity pattern. The SPQRFactor class also provides a method to apply the \(Q\) factor to a dense matrix or vector, using the SPQRFactor.qmult() method.

Currently, it is not possible to extract the \(Q\) and \(R\) factors as sparse matrices from the SPQRFactor object. Users who need the explicit factors should use the spqr() function instead.

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

Exceptions and Warnings¶

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

Sparse QR Decomposition (sksparse.spqr)¶