Approximate Minimum Degree (AMD) Ordering (sksparse.amd)¶
The sksparse.amd module provides an interface to the Approximate
Minimum Degree (AMD) ordering algorithm for sparse, square
matrices.
It exposes the main function of the AMD package, which
computes a symmetric ordering of a sparse matrix that minimizes the fill-in of
the Cholesky decomposition. The AMD function accepts both real and complex
matrices, in any format supported by scipy.sparse (CSC format is most
efficient).
Quickstart¶
If \(A\) is a sparse, square matrix, then the following code computes the AMD ordering of \(A\):
from sksparse.amd import amd
A = ... # some sparse matrix
p = amd(A)
PAPT = A[p][:, p]
to give the permuted matrix \(PAP^T\), where \(P\) is the permutation matrix corresponding to the ordering \(p\). If \(A\) is not symmetric, then this is the same as AMD computes the ordering of the symbolically symmetric matrix \(A + A^T\).
We can then compute the Cholesky decompositions of the original and permuted
matrix using cholesky(), and compare the number of
non-zeros in each:
from sksparse.cholmod import cholesky
A_factor = cholesky(A)
PAPT_factor = cholesky(PAPT)
L = A_factor.L()
Lp = PAPT_factor.L()
print("Number of non-zeros in L: ", L.nnz)
print("Number of non-zeros in Lp:", Lp.nnz)
The number of non-zeros in the Cholesky factorization of the permuted matrix should be less than or equal to the number of non-zeros in the Cholesky factorization of the original matrix, but this is not guaranteed.
Example¶
To see the effects of AMD ordering, we can load a sparse matrix from the SuiteSparse Matrix Collection and compute its AMD ordering.
#!/usr/bin/env python3
# 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: amd_example.py
# Created: 2025-07-29 12:33
# =============================================================================
"""An example of using the AMD (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 sksparse.amd import amd
from sksparse.cholmod import cholesky
# Load the bcsstk06 matrix (downloaded from the SuiteSparse Matrix Collection:
# <https://sparse.tamu.edu/HB/bcsstk06>)
filepath = Path("data") / "bcsstk06.mtx"
A = mmread(filepath, spmatrix=False).tocsc() # read the matrix
p = amd(A) # compute the AMD ordering
PAPT = A[p][:, p] # apply the ordering to the matrix
# Compute the Cholesky factorization of each matrix
L = cholesky(A, lower=True)
Lp = cholesky(PAPT, lower=True)
# 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("AMD Example: Fill-Reducing Ordering of bcsstk06")
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^T$")
ax = axs[1, 0]
ax.spy(L, markersize=1)
ax.set_title(r"Original Cholesky Factor $L$")
ax.set_xlabel(f"{L.nnz:,} non-zeros")
ax = axs[1, 1]
ax.spy(Lp, markersize=1)
ax.set_title(r"Permuted Cholesky Factor $L_p$")
ax.set_xlabel(f"{Lp.nnz:,} non-zeros")
for ax in axs.flat:
ax.set_aspect("equal")
ax.set_xticks([])
ax.set_yticks([])
plt.show()
fig.savefig("amd_example.svg")
The figure shows the effect of AMD ordering that reduces the fill-in of the Cholesky factorization of a sparse matrix.
The number of non-zeros in the Cholesky factorization of the original matrix (left) and the permuted matrix (right) using AMD ordering.¶
AMDInfo Objects¶
An AMDInfo object is a dataclass returned by the amd() function
when the return_info parameter is set to True. It contains information
about the AMD ordering, including the return status, the number of non-zeros in
the Cholesky factorization, and others.
We can use AMDInfo objects to compare the number of non-zeros in the
Cholesky factorization of the original matrix, without computing it directly:
from sksparse.amd import amd
from sksparse.cholmod import cholesky
A = ... # some sparse matrix
N = A.shape[0]
p, info = amd(A, return_info=True)
PAPT_factor = cholesky(A[p][:, p])
print("nnz in L of A: ", info.Lnz + N)
print("nnz in L of PAPT:", PAPT_factor.L().nnz)
Typically, the AMDInfo object is unnecessary, and you can just use the
permutation vector returned by amd().
Convenience Methods¶
The AMD package also provides a convenience function,
amd_default_control() to get the default control parameters from the AMD
package. Most users will not need to use this function, as the default control
parameters are used automatically by amd().
Error Handling¶
Errors raised by the AMD package are converted into Python exceptions. See the Exceptions and Warnings for details.
References¶
Amestoy, P. R., Davis, T. A., & Duff, I. S. (1996). An approximate minimum degree ordering algorithm. SIAM Journal on Matrix Analysis and Applications, 17(4), 886-905. <https://epubs.siam.org/doi/abs/10.1137/S0895479894278952>.