Column Approximate Minimum Degree (COLAMD) Ordering (sksparse.colamd)

The sksparse.colamd module provides efficient an implementation of the Column Approximate Minimum Degree (COLAMD) ordering algorithm for sparse matrices.

It exposes the main functions of the COLAMD package, which computes a column ordering \(Q\) of a sparse matrix that minimizes the fill-in of the Cholesky decomposition of \((AQ)^{\top}(AQ)\). The colamd() function is appropriate for use with non-symmetric and non-square matrices, for LU factorization, QR factorization, and other decompositions that require a column ordering.

This module also provides a symmetric variant, symamd(), which computes a permutation P of a symmetric matrix A such that the Cholesky factorization of \(PAP^{\top}\) has less fill-in and requires fewer floating point operations than A. This function assumes that its input is symmetric.

The colamd() and symamd() functions accept both real and complex matrices, in any format supported by scipy.sparse (CSC format is most efficient).

Quickstart

If \(A\) is a sparse matrix, then the following code computes the COLAMD ordering of \(A\):

from sksparse.colamd import colamd
A = ...  # some sparse matrix
q = colamd(A)
AQ = A[:, q]  # permute the columns of A

to give the permuted matrix \(AQ\), where \(Q\) is the permutation matrix corresponding to the ordering \(q\).

We can then continue from above to compute the LU decompositions of the original and permuted matrix, and compare the number of non-zeros in each:

from scipy.sparse.linalg import splu
lu = splu(A)
luq = splu(A[:, q])
L, U = lu.L, lu.U
Lq, Uq = luq.L, luq.U
print("Number of non-zeros in L + U:  ", (L + U).nnz)
print("Number of non-zeros in Lq + Uq:", (Lq + Uq).nnz)

The number of non-zeros in the LU factorization of the permuted matrix should be less than or equal to the number of non-zeros in the LU factorization of the original matrix, but this is not guaranteed.

Example

To see the effects of COLAMD ordering, 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: colamd_example.py
#  Created: 2025-07-31 14:18
# =============================================================================

"""An example of using the COLAMD (Column Approximate Minimum Degree) algorithm
to find a fill-reducing ordering of a sparse matrix.
"""

from pathlib import Path

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

from sksparse.colamd import colamd

# from numpy.testing import assert_allclose, assert_array_equal


# 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

q = colamd(A)  # compute the COLAMD ordering
AQ = A[:, q]   # apply the column ordering to the matrix

# Compute the LU factorization of the original and permuted matrices
lu = splu(A, permc_spec="NATURAL")
L_, U_ = lu.L, lu.U

luq = splu(AQ, permc_spec="NATURAL")
Lq, Uq = luq.L, luq.U

# 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("COLAMD Example: Fill-Reducing Ordering of west0479")

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

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

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

ax = axs[1, 1]
ax.spy(Lq + Uq, markersize=1)
ax.set_title(r"Permuted LU Factors $L_q + U_q$")
ax.set_xlabel(f"{(Lq + Uq).nnz:,} total non-zeros")

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

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

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

COLAMD Example

The number of non-zeros in the LU factorization of the original matrix (left) and the permuted matrix (right) using COLAMD ordering.

COLAMDStats Objects

An COLAMDStats object is a dataclass returned by the colamd() function when the return_info parameter is set to True. It contains information about the ordering, including the return status.

Typically, the COLAMDStats object is unnecessary, and you can just use the permutation vector returned by colamd().

Convenience Methods

The COLAMD package also provides a convenience function, colamd_get_defaults() to get the default control parameters from the COLAMD package. Most users will not need to use this function, as the default control parameters are used automatically by colamd().

Error Handling

Errors raised by the COLAMD package are converted into Python exceptions. See the Exceptions and Warnings for details.