Public API Documentation

MatrixBandwidth

Fast algorithms for matrix bandwidth minimization and recognition, written in Julia.

The bandwidth of an $n×n$ matrix $A$ is the minimum non-negative integer $k ∈ \{0, 1, …, n - 1\}$ such that $A[i, j] = 0$ whenever $|i - j| > k$. Reordering the rows and columns of a matrix to reduce its bandwidth has many practical applications in engineering and scientific computing: it can improve performance when solving linear systems, approximating partial differential equations, optimizing circuit layout, and more [Maf14, p. 184]. There are two variants of this problem: minimization, which involves finding a permutation matrix $P$ such that the bandwidth of $PAPᵀ$ is minimized, and recognition, which entails determining whether there exists a permutation matrix $P$ such that the bandwidth of $PAPᵀ$ is at most some fixed non-negative integer (an optimal permutation that fully minimizes the bandwidth of $A$ is not required).

Many matrix bandwidth reduction algorithms exist in the literature, but implementations in the open-source ecosystem are scarce, with those that do exist primarily tackling older, less efficient algorithms. The Boost libraries in C++, the NetworkX library in Python, and the MATLAB standard library all only implement the reverse Cuthill–McKee algorithm from 1971. This gap in the ecosystem not only makes it difficult for theoretical researchers to benchmark and compare new proposed algorithms but also precludes the application of the most performant modern algorithms in real-life industry settings. MatrixBandwidth.jl aims to bridge this gap by presenting a unified interface for matrix bandwidth reduction algorithms in Julia, designed with extensibility to further methods in mind.

The following matrix bandwidth reduction algorithms are currently available:

• Minimization
  • Exact
    • Del Corso–Manzini (Minimization.DelCorsoManzini)
    • Del Corso–Manzini with perimeter search (Minimization.DelCorsoManziniWithPS)
    • Caprara–Salazar-González (Minimization.CapraraSalazarGonzalez)
    • Saxe–Gurari–Sudborough (Minimization.SaxeGurariSudborough)
    • Brute-force search (Minimization.BruteForceSearch)
  • Heuristic
    • Gibbs–Poole–Stockmeyer (Minimization.GibbsPooleStockmeyer)
    • Cuthill–McKee (Minimization.CuthillMcKee)
    • Reverse Cuthill–McKee (Minimization.ReverseCuthillMcKee)
• Recognition
  • Del Corso–Manzini (Recognition.DelCorsoManzini)
  • Del Corso–Manzini with perimeter search (Recognition.DelCorsoManziniWithPS)
  • Caprara–Salazar-González (Recognition.CapraraSalazarGonzalez)
  • Saxe–Gurari–Sudborough (Recognition.SaxeGurariSudborough)
  • Brute-force search (Recognition.BruteForceSearch)

Recognition algorithms determine whether any row-and-column permutation of a matrix induces bandwidth less than or equal to some fixed integer. Exact minimization algorithms always guarantee optimal orderings to minimize bandwidth, while heuristic minimization algorithms produce near-optimal solutions more quickly. Metaheuristic minimization algorithms employ iterative search frameworks to find better solutions than heuristic methods (albeit more slowly); no such algorithms are already implemented, but several (e.g., simulated annealing) are currently under development.

This package also exports several additional core functions, including (but not limited to) bandwidth and profile to compute the original bandwidth and profile of a matrix prior to any reordering.

The full documentation is available at GitHub Pages.

References

• [Maf14]: L. O. Mafteiu-Scai. The Bandwidths of a Matrix. A Survey of Algorithms. Annals of West University of Timisoara - Mathematics and Computer Science 52, 183–223 (2014). https://doi.org/10.2478/awutm-2014-0019.

AbstractAlgorithm

Abstract base type for all matrix bandwidth reduction algorithms.

Interface

Concrete subtypes of AbstractAlgorithm must implement the following methods:

• Base.summary(::T) where {T<:AbstractAlgorithm}: returns a String indicating the name of the algorithm (e.g., "Gibbs–Poole–Stockmeyer").
• _requires_structural_symmetry(::T) where {T<:AbstractAlgorithm}: returns a Bool indicating whether the algorithm requires the input matrix to be structurally symmetric.

Direct subtypes of AbstractAlgorithm must implement the following method:

• _problem(::T) where {T<:AbstractAlgorithm}: returns a Symbol indicating the matrix bandwidth problem tackled by the algorithm (e.g., :minimization).

AbstractResult

Abstract base type for all matrix bandwidth reduction results.

Interface

Concrete subtypes of AbstractResult must implement parametric types

• A<:AbstractAlgorithm;
• M<:AbstractMatrix{<:Number}; and
• O<:Union{Nothing,Vector{Int}},

alongside the following fields:

• algorithm::A: the algorithm used to investigate the bandwidth.
• matrix::M: the matrix whose bandwidth is investigated.
• ordering::O: the corresponding ordering of the rows and columns, if a relevant one is found; otherwise, nothing.

bandwidth(A) -> Int

Compute the bandwidth of A before any permutation of its rows and columns.

The bandwidth of an $n×n$ matrix $A$ is the minimum non-negative integer $k ∈ \{0, 1, …, n - 1\}$ such that $A[i, j] = 0$ whenever $|i - j| > k$.

In contrast to minimize_bandwidth, this function does not attempt to minimize the bandwidth of A by permuting its rows and columns—it simply computes its bandwidth as is.

Arguments

• A::AbstractMatrix{<:Number}: the (square) matrix whose bandwidth is computed.

Returns

• ::Int: the bandwidth of A.

Performance

Given an $n×n$ input matrix, this relatively simple algorithm runs in $O(n²)$ time.

Examples

bandwidth correctly identifies the bandwidth of a pentadiagonal matrix as $2$ and does not attempt to find a minimizing permutation upon shuffling of its rows and columns:

julia> using Random

julia> Random.seed!(242622);

julia> (n, k) = (8, 2);

julia> perm = randperm(n);

julia> A = (!iszero).(MatrixBandwidth.random_banded_matrix(n, k))
8×8 BitMatrix:
 1  0  0  0  0  0  0  0
 0  1  0  1  0  0  0  0
 0  0  0  1  1  0  0  0
 0  1  1  1  0  1  0  0
 0  0  1  0  0  0  0  0
 0  0  0  1  0  0  0  0
 0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0

julia> bandwidth(A)
2

julia> A_shuffled = A[perm, perm]
8×8 BitMatrix:
 0  0  0  0  0  0  1  0
 0  0  0  0  0  0  0  0
 0  0  0  1  0  0  0  0
 0  0  1  0  0  0  1  0
 0  0  0  0  1  0  0  0
 0  0  0  0  0  1  1  0
 1  0  0  1  0  1  1  0
 0  0  0  0  0  0  0  0

julia> bandwidth(A_shuffled)
6

Notes

Some texts define matrix bandwidth to be the minimum non-negative integer $k$ such that $A[i, j] = 0$ whenever $|i - j| ≥ k$ instead, particularly in more mathematically-minded communities. Effectively, this definition treats diagonal matrices as bandwidth $1$, tridiagonal matrices as bandwidth $2$, and so on. Our definition, on the other hand, is more common in computer science contexts, treating diagonal matrices as bandwidth $0$ and tridiagonal matrices as bandwidth $1$. (Both definitions, however, agree that the bandwidth of an empty matrix is simply $0$.)

bandwidth_lower_bound(A) -> Int

Compute a lower bound on the bandwidth of A using [CS05, pp. 359–60]'s results.

A is assumed to be structurally symmetric, since the bound from [CS05, pp.359–60] was discovered in the context of undirected graphs (whose adjacency matrices are symmetric). Since the original algorithm is defined only for connected graphs, we compute the bound on each connected component of the graph represented by A and return the maximum of these.

The bandwidth of an $n×n$ matrix $A$ is the minimum non-negative integer $k ∈ \{0, 1, …, n - 1\}$ such that $A[i, j] = 0$ whenever $|i - j| > k$.

In contrast to minimize_bandwidth, this function does not attempt to truly minimize the bandwidth of A—it simply returns a lower bound on its bandwidth up to symmetric permutation of its rows and columns. This bound is not generally tight, but it indeed matches the true minimum in many non-trivial cases and is easily computable in $O(n³)$ time (dominated by the Floyd–Warshall algorithm call; the core logic itself runs in $O(n²)$ time).

Arguments

• A::AbstractMatrix{<:Number}: the (square) matrix on whose bandwidth a lower bound is to be computed. A must be structurally symmetric (i.e., A[i, j] must be nonzero if and only if A[j, i] is nonzero for $1 ≤ i, j ≤ n$).

Returns

• ::Int: a lower bound on the bandwidth of A. (This bound is tight in many non-trivial cases but not universally so.)

Examples

julia> using Random, SparseArrays, Combinatorics

julia> Random.seed!(21);

julia> (n, p) = (9, 0.4);

julia> A = sprand(n, n, p);

julia> A = A + A' # Ensure structural symmetry;

julia> minimize_bandwidth(A, Minimization.BruteForceSearch())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Brute-force search
 * Approach: exact
 * Minimum Bandwidth: 5
 * Original Bandwidth: 8
 * Matrix Size: 9×9

julia> bandwidth_lower_bound(A) # Always less than or equal to the true minimum bandwidth
4

Notes

Some texts define matrix bandwidth to be the minimum non-negative integer $k$ such that $A[i, j] = 0$ whenever $|i - j| ≥ k$ instead, particularly in more mathematically-minded communities. Effectively, this definition treats diagonal matrices as bandwidth $1$, tridiagonal matrices as bandwidth $2$, and so on. Our definition, on the other hand, is more common in computer science contexts, treating diagonal matrices as bandwidth $0$ and tridiagonal matrices as bandwidth $1$. (Both definitions, however, agree that the bandwidth of an empty matrix is simply $0$.)

References

• [CS05]: A. Caprara and J.-J. Salazar-González. Laying Out Sparse Graphs with Provably Minimum Bandwidth. INFORMS Journal on Computing 17, 356–73 (2005). https://doi.org/10.1287/ijoc.1040.0083.

profile(A) -> Int

Compute the profile of A before any permutation of its rows and columns.

The profile of a structurally symmetric $n×n$ matrix $A$ is traditionally defined as the sum of the distances from each diagonal entry to the leftmost nonzero entry in that row—in other words, $∑ᵢ₌₁ⁿ (i - fᵢ)$, where each $fᵢ$ is the smallest index such that $A[i, fᵢ] ≠ 0$ [Maf14, pp. 187-88]. Generalizing this property to all square matrices, we define the column profile of a matrix to be the sum of the distances from each diagonal entry to the farthest (not necessarily topmost) nonzero entry in that column and the row profile to be the sum of the distances from each diagonal entry to the farthest (not necessarily leftmost) nonzero entry in that row. (Note that both of these properties are equal to traditional matrix profile for structurally symmetric matrices.)

One of the most common contexts in which matrix profile is relevant is sparse matrix storage, where lower-profile matrices occupy less space in memory [Maf14, p.188]. Since the SparseArrays.jl standard library package defaults to compressed sparse column storage over compressed sparse row, we therefore compute column profile by default unless the dimension is otherwise specified.

Arguments

• A::AbstractMatrix{<:Number}: the (square) matrix whose profile is computed.

Keyword Arguments

• dim::Symbol=:col: the dimension along which the profile is computed; must be either :col (the default) or :row.

Returns

• ::Int: the profile of A along the specified dimension.

Performance

Given an $n×n$ input matrix, this relatively simple algorithm runs in $O(n²)$ time.

Examples

profile computes the column profile of a matrix by default:

julia> using Random, SparseArrays

julia> Random.seed!(2287);

julia> (n, p) = (25, 0.05);

julia> A = sprand(n, n, p)
25×25 SparseMatrixCSC{Float64, Int64} with 29 stored entries:
⎡⠀⠀⠀⠀⠀⠀⠐⠀⠒⠀⡀⠀⠀⎤
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠐⠂⎥
⎢⠀⠀⠀⢀⠌⠀⠀⠀⢀⠈⠀⠀⠀⎥
⎢⠠⢄⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠂⡀⠀⠀⠌⠀⠈⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠢⡀⢄⡈⠀⠀⠀⎥
⎣⠀⠀⠀⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⎦

julia> profile(A)
211

The dimension (:row or :col) can also be explicitly specified:

julia> using Random, SparseArrays

julia> Random.seed!(3647);

julia> (n, p) = (25, 0.05);

julia> A = sprand(n, n, p)
25×25 SparseMatrixCSC{Float64, Int64} with 31 stored entries:
⎡⠄⠀⠀⠀⠀⠀⠀⠘⠀⠀⠀⠁⠀⎤
⎢⠄⢀⠀⠀⠁⠀⠀⠀⠀⢀⠀⠀⠀⎥
⎢⠀⢀⡂⠀⠀⠀⠀⠀⠀⠀⠀⡄⠀⎥
⎢⠀⠀⠀⠀⠀⡀⠂⠀⠀⠀⠀⠀⠀⎥
⎢⠁⠀⠁⠀⠁⠀⠀⠁⠄⢀⠈⠀⠀⎥
⎢⠂⠐⠐⠐⠠⠀⠄⠀⠀⠀⠠⣀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎦

julia> profile(A, dim=:row)
147

julia> profile(A, dim=:col)
175

References

• [Maf14]: L. O. Mafteiu-Scai. The Bandwidths of a Matrix. A Survey of Algorithms. Annals of West University of Timisoara - Mathematics and Computer Science 52, 183–223 (2014). https://doi.org/10.2478/awutm-2014-0019.

MatrixBandwidth.Minimization

Exact, heuristic, and metaheuristic algorithms for matrix bandwidth minimization in Julia.

The bandwidth of an $n×n$ matrix $A$ is the minimum non-negative integer $k ∈ \{0, 1, …, n - 1\}$ such that $A[i, j] = 0$ whenever $|i - j| > k$. The matrix bandwidth minimization problem involves finding a permutation matrix $P$ such that the bandwidth of $PAPᵀ$ is minimized.

The following matrix bandwidth minimization algorithms are currently available:

• Exact
  • Del Corso–Manzini (DelCorsoManzini)
  • Del Corso–Manzini with perimeter search (DelCorsoManziniWithPS)
  • Caprara–Salazar-González (CapraraSalazarGonzalez)
  • Saxe–Gurari–Sudborough (SaxeGurariSudborough)
  • Brute-force search (BruteForceSearch)
• Heuristic
  • Gibbs–Poole–Stockmeyer (GibbsPooleStockmeyer)
  • Cuthill–McKee (CuthillMcKee)
  • Reverse Cuthill–McKee (ReverseCuthillMcKee)

Exact minimization algorithms always guarantee optimal orderings to minimize bandwidth, while heuristic minimization algorithms produce near-optimal solutions more quickly. Metaheuristic minimization algorithms employ iterative search frameworks to find better solutions than heuristic methods (albeit more slowly); no such algorithms are already implemented, but several (e.g., simulated annealing) are currently under development.

This submodule is part of the MatrixBandwidth.jl package.

AbstractSolver <: AbstractAlgorithm

Abstract base type for all matrix bandwidth minimization solvers.

Interface

As per the interface of supertype AbstractAlgorithm, concrete subtypes of AbstractSolver must implement the following methods:

• Base.summary(::T) where {T<:AbstractSolver}: returns a String indicating the name of the solver (e.g., "Gibbs–Poole–Stockmeyer").
• _requires_structural_symmetry(::T) where {T<:AbstractSolver}: returns a Bool indicating whether the solver requires the input matrix to be structurally symmetric.

Direct subtypes of AbstractSolver must implement the following method:

• _approach(::T) where {T<:AbstractSolver}: returns a Symbol indicating the category of solver (e.g., :heuristic).

Supertype Hierarchy

AbstractSolver <: AbstractAlgorithm

Notes

To implement a new matrix bandwidth minimization algorithm, define a new concrete subtype of AbstractSolver (or of one of its abstract subtypes like MetaheuristicSolver) then implement a corresponding _minimize_bandwidth_impl(::AbstractMatrix{Bool}, ::NewSolverType) method. Do not attempt to directly implement a new minimize_bandwidth method, as the function contains common preprocessing logic independent of the specific algorithm used.

MinimizationResult{A,M,O} <: AbstractResult

Output struct for matrix bandwidth minimization results.

Fields

• algorithm::A<:AbstractSolver: the solver used to minimize the bandwidth.
• matrix::M<:AbstractMatrix{<:Number}: the original matrix whose bandwidth is minimized.
• ordering::O<:Vector{Int}: the (near-)optimal ordering of the rows and columns.
• bandwidth::Int: the minimized bandwidth of the matrix.
• approach::Symbol: the approach used by the solver. (Should be one of :exact, :heuristic, and :metaheuristic.)

Constructors

• MinimizationResult(algorithm, matrix, ordering, bandwidth): constructs a new MinimizationResult instance with the given fields. The approach field is automatically determined based on the algorithm type.

Supertype Hierarchy

MinimizationResult <: AbstractResult

minimize_bandwidth(A, solver=GibbsPooleStockmeyer()) -> MinimizationResult

Minimize the bandwidth of A using the algorithm defined by solver.

The bandwidth of an $n×n$ matrix $A$ is the minimum non-negative integer $k ∈ \{0, 1, …, n - 1\}$ such that $A[i, j] = 0$ whenever $|i - j| > k$.

This function computes a (near-)optimal ordering $π$ of the rows and columns of $A$ so that the bandwidth of $PAPᵀ$ is minimized, where $P$ is the permutation matrix corresponding to $π$. This is known to be an NP-complete problem; however, several heuristic algorithms such as Gibbs–Poole–Stockmeyer run in polynomial time while still still producing near-optimal orderings in practice. Exact methods like Caprara–Salazar-González are also available, but they are at least exponential in time complexity and thus only feasible for relatively small matrices.

Arguments

• A::AbstractMatrix{<:Number}: the (square) matrix whose bandwidth is minimized.
• solver::AbstractSolver: the matrix bandwidth minimization algorithm to use; defaults to GibbsPooleStockmeyer. (See the Minimization module documentation for a full list of supported solvers.)

Returns

• ::MinimizationResult: a struct containing the algorithm used, the original matrix A, the (near-)optimal ordering of the rows and columns, and the minimized bandwidth.

Examples

Multiple algorithms to minimize the bandwidth of a given matrix are available. In particular, there are exact solvers (which always guarantee optimal solutions), heuristic solvers (which produce near-optimal solutions more quickly than exact methods), and metaheuristic solvers (which employ iterative search frameworks to find better solutions than heuristic methods, but even more slowly).

Certainly, exact solvers will always produce the same optimal bandwidth (but likely different orderings):

julia> using Random, SparseArrays

julia> Random.seed!(38);

julia> (n, p) = (20, 0.05);

julia> A = sprand(n, n, p);

julia> A = A .+ A' # Ensure structural symmetry
20×20 SparseMatrixCSC{Float64, Int64} with 36 stored entries:
⎡⠀⠀⠀⠀⠀⣎⠁⠈⠀⠂⎤
⎢⠀⠀⠀⠀⡀⠐⠀⠀⠀⠠⎥
⎢⡠⢤⢀⠈⠠⠂⠀⠀⠔⠀⎥
⎢⡁⠀⠀⠀⠀⠀⠀⠀⡀⠐⎥
⎣⠠⠀⠀⡀⠐⠁⢀⠈⢐⠔⎦

julia> res_dcm = minimize_bandwidth(A, Minimization.DelCorsoManzini())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Del Corso–Manzini
 * Approach: exact
 * Minimum Bandwidth: 3
 * Original Bandwidth: 17
 * Matrix Size: 20×20

julia> A[res_dcm.ordering, res_dcm.ordering]
20×20 SparseMatrixCSC{Float64, Int64} with 36 stored entries:
⎡⠠⠂⠤⡀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠣⡁⠈⠀⡀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠠⡤⠋⠤⡀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠣⡀⡨⠢⡀⎥
⎣⠀⠀⠀⠀⠀⠀⠈⠢⠁⠀⎦

julia> res_sgs = minimize_bandwidth(A, Minimization.SaxeGurariSudborough())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Saxe–Gurari–Sudborough
 * Approach: exact
 * Minimum Bandwidth: 3
 * Original Bandwidth: 17
 * Matrix Size: 20×20

julia> A[res_sgs.ordering, res_sgs.ordering]
20×20 SparseMatrixCSC{Float64, Int64} with 36 stored entries:
⎡⠀⡠⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠔⣡⠢⡀⠀⠀⠀⠀⎥
⎢⠀⠀⠈⠢⡀⡨⣀⡀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠸⡀⠈⢂⡀⎥
⎣⠀⠀⠀⠀⠀⠀⠈⠰⠄⡡⎦

However, the answers of (meta)heuristic solvers may differ from each other:

julia> using Random, SparseArrays

julia> Random.seed!(405);

julia> (n, p) = (400, 0.002);

julia> A = sprand(n, n, p);

julia> A = A .+ A' # Ensure structural symmetry;

julia> minimize_bandwidth(A, Minimization.GibbsPooleStockmeyer())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Gibbs–Poole–Stockmeyer
 * Approach: heuristic
 * Minimum Bandwidth: 31
 * Original Bandwidth: 380
 * Matrix Size: 400×400

julia> minimize_bandwidth(A, Minimization.ReverseCuthillMcKee())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Reverse Cuthill–McKee
 * Approach: heuristic
 * Minimum Bandwidth: 37
 * Original Bandwidth: 380
 * Matrix Size: 400×400

If no solver is specified, then the heuristic Gibbs–Poole–Stockmeyer algorithm is used by default:

julia> using Random, SparseArrays

julia> Random.seed!(80);

julia> (n, p) = (500, 0.001);

julia> A = sprand(n, n, p);

julia> A = A .+ A' # Ensure structural symmetry
500×500 SparseMatrixCSC{Float64, Int64} with 496 stored entries:
⎡⠀⠀⠁⠐⠀⠀⠀⠄⠠⢐⠄⠀⡈⠀⠀⡆⠠⠑⠀⠰⠀⠀⠀⠐⠀⠄⠀⠀⡠⠀⠀⠀⠆⠐⠀⠄⢠⡈⠈⠀⎤
⎢⢁⠀⠄⠁⠠⠊⠀⠀⠂⠀⡀⠂⠀⢀⠀⠀⠀⠈⠠⠠⠀⠠⠠⠀⠀⡀⠀⠁⡀⠀⡀⠐⠀⠄⠀⠀⠄⠂⡈⠠⎥
⎢⠀⠀⡠⠂⠅⠁⠀⠀⠀⠀⠐⡁⠀⡀⠈⢈⠀⠐⠀⠀⠀⠀⠃⢀⠀⠀⠀⠀⠀⠀⠐⠀⠤⠁⠀⠁⠀⠀⠀⡀⎥
⎢⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠀⠠⠀⠂⠁⠀⠀⡐⠔⠂⠂⠀⠠⠀⠐⠀⠀⠄⠀⠄⠀⠀⡀⠀⢀⠁⠀⎥
⎢⢀⢂⠈⠀⠀⠀⠀⠀⠎⠁⠄⠂⠀⠁⣀⠁⠈⠠⠀⠀⠂⠄⢀⠄⠀⠠⠀⠀⠉⠈⠀⠀⠀⠀⠠⠀⠀⠀⠀⡀⎥
⎢⠀⠁⠠⠈⠔⠠⠠⠀⠠⠁⠊⠄⠀⢀⠀⠀⠠⠀⠠⡌⠀⠈⠀⠀⠀⠠⠈⠀⠀⠀⠂⠀⠀⠀⠀⠡⠀⠠⠠⢂⎥
⎢⠂⠈⠀⢀⠀⠠⠀⠀⠄⠀⠀⢀⡐⠈⠀⠀⠀⠀⠀⠈⠐⠠⠀⠀⠀⠀⠀⡠⠑⠀⠂⠀⠀⠀⡀⠀⡐⠀⠄⠁⎥
⎢⠠⠤⠀⠀⡂⢀⠀⠂⠄⠘⠀⠀⠀⠀⠀⠀⠀⠑⠀⠂⡄⠀⠄⠀⡁⡐⠁⠠⠂⠀⠀⠐⠀⠀⠁⠀⠀⠁⠈⠀⎥
⎢⢄⠂⡀⠀⢀⠀⠌⠀⠂⡀⠀⠂⠀⠀⢄⠀⠀⠀⠀⢀⠀⠀⢀⠄⡀⠂⠂⠀⠅⠀⠀⠀⠀⠁⠀⢀⠄⠀⠀⠀⎥
⎢⢀⡀⠀⡂⠀⠀⠀⠀⠀⠀⡀⠦⡀⠀⠠⠀⠀⢀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⡐⢀⢀⠁⠔⠆⠀⠀⠀⢀⠅⎥
⎢⠀⠀⠀⡀⠀⠀⢐⠌⠈⠄⡀⠀⠐⡀⠀⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠀⠌⡀⠀⠀⠀⠀⠀⠀⠂⎥
⎢⢀⠀⠀⠂⠉⢀⠨⠀⠀⠔⠀⠀⠀⠀⠀⠁⠀⠔⠐⠀⠀⠀⠀⡠⠠⠀⠰⠀⠠⠘⠀⠐⠀⡀⠁⠠⠀⡑⠀⠀⎥
⎢⠀⠄⠀⠠⠀⠀⠀⡀⠀⡀⠀⡀⠀⠀⢁⠨⠠⠈⠀⠀⠀⠀⠀⠂⠀⠀⠁⠀⠁⠂⠄⠀⠠⠀⠀⠀⠀⠀⡀⠀⎥
⎢⠀⠀⠄⠀⠀⠀⢀⠀⠀⠀⠂⠀⠀⡠⠁⡀⠈⠀⠀⠀⠀⠠⠐⠂⠁⠀⠀⠀⢁⠊⠠⠀⢀⠀⠀⠀⠰⠌⠀⠀⎥
⎢⠀⠊⠀⠈⠀⠀⠀⠀⡃⠀⠀⠀⠑⠀⠈⠀⠁⠁⢀⠠⠀⠀⣀⠂⠡⠀⡡⠐⠀⠀⠀⠐⠂⠰⠀⠠⠀⠀⠂⠀⎥
⎢⠀⠀⢀⠈⠐⠀⠀⠁⠀⠀⠈⠀⠈⠀⢀⠀⠀⠀⠀⢐⡀⠄⢀⠀⠀⠁⠀⠂⢀⠀⠤⠃⠀⠐⠀⠀⠀⠈⠀⠀⎥
⎢⢈⠁⠀⠄⠄⠃⠀⠁⠀⠀⠀⠀⠀⠀⠀⠀⠄⠀⢁⠄⠀⠈⠀⠠⠀⠂⠀⠐⢈⡀⢀⠀⠄⠁⠀⠈⠀⠆⠀⠀⎥
⎢⠀⠄⠀⠀⠄⠀⠀⠠⠀⠂⠄⡀⠀⠈⠁⠀⠀⢀⠈⠁⠀⠀⠁⡀⠀⠀⠀⠀⠀⡀⠀⠀⡀⠀⢠⠒⠈⠀⠀⠀⎥
⎢⡀⠲⠠⠁⠀⠀⠀⢀⠀⠀⠀⡀⠐⠈⠄⠀⠀⠁⠀⠀⠀⠀⢄⠠⠀⠀⡐⠆⠀⠀⡀⠀⠠⠄⠂⠀⠤⠃⠀⠁⎥
⎣⠂⠀⠂⡈⠀⠠⠁⠀⠀⠠⠠⢂⠄⠁⠂⠀⠀⠀⠄⠔⠠⠀⠀⠀⠀⠈⠀⠀⠈⠀⠀⠀⠀⠀⠀⠀⠄⠀⠁⠀⎦

julia> res = minimize_bandwidth(A)
Results of Bandwidth Minimization Algorithm
 * Algorithm: Gibbs–Poole–Stockmeyer
 * Approach: heuristic
 * Minimum Bandwidth: 6
 * Original Bandwidth: 487
 * Matrix Size: 500×500

julia> A[res.ordering, res.ordering]
500×500 SparseMatrixCSC{Float64, Int64} with 496 stored entries:
⎡⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠀⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠑⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠐⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠊⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠛⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠻⣢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠛⢄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠱⣦⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⣦⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⢆⡀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⢆⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⎦

Notes

To implement a new matrix bandwidth minimization algorithm, define a new concrete subtype of AbstractSolver (or of one of its abstract subtypes like MetaheuristicSolver) then implement a corresponding _minimize_bandwidth_impl(::AbstractMatrix{Bool}, ::NewSolverType) method. Do not attempt to directly implement a new minimize_bandwidth method, as the function contains common preprocessing logic independent of the specific algorithm used.

Note also that some texts define matrix bandwidth to be the minimum non-negative integer $k$ such that $A[i, j] = 0$ whenever $|i - j| ≥ k$ instead, particularly in more mathematically-minded communities. Effectively, this definition treats diagonal matrices as bandwidth $1$, tridiagonal matrices as bandwidth $2$, and so on. Our definition, on the other hand, is more common in computer science contexts, treating diagonal matrices as bandwidth $0$ and tridiagonal matrices as bandwidth $1$. (Both definitions, however, agree that the bandwidth of an empty matrix is simply $0$.)

MatrixBandwidth.Minimization.Exact

Exact solvers for matrix bandwidth minimization.

Exact methods are those which guarantee an optimal ordering producing the true minimum bandwidth of a matrix. Since bandwidth minimization is an NP-complete problem, existing exact algorithms are, at best, exponential in time complexity—much worse than many polynomial-time heuristic approaches (e.g., Gibbs–Poole–Stockmeyer). Such methods, therefore, are not feasible for large matrices, but they remain useful when precise solutions are required for small-to-medium-sized inputs (say, up to $100×100$).

The following exact matrix bandwidth minimization algorithms are currently available:

• Del Corso–Manzini (DelCorsoManzini)
• Del Corso–Manzini with perimeter search (DelCorsoManziniWithPS)
• Caprara–Salazar-González (CapraraSalazarGonzalez)
• Saxe–Gurari–Sudborough (SaxeGurariSudborough)
• Brute-force search (BruteForceSearch)

This submodule is part of the MatrixBandwidth.Minimization submodule of the MatrixBandwidth.jl package.

BruteForceSearch <: ExactSolver <: AbstractSolver <: AbstractAlgorithm

The simplest exact method for minimizing the bandwidth of a matrix is to iterate over all possible symmetric permutations and compare the bandwidths they induce.

Since $i₁, i₂, … iₙ$ induces the same bandwidth as $iₙ, iₙ₋₁, … i₁$, we restrict our search to orderings such that $i₁ ≤ iₙ$ (with equality checked just in case $n = 1$).

Supertype Hierarchy

BruteForceSearch <: ExactSolver <: AbstractSolver <: MatrixBandwidth.AbstractAlgorithm

Performance

Given an $n×n$ input matrix $A$, this brute-force algorithm runs in $O(n! ⋅ n²)$ time:

• Precisely $n!/2$ permutations are checked (except when $n = 1$, in which case $1! = 1$ permutation is checked). This is, clearly, $O(n!)$.
• For each permutation, the bandwidth function is called on view(A, perm, perm), which takes $O(n²)$ time.
• Therefore, the overall time complexity is $O(n! ⋅ n²)$.

Indeed, due to the need to exhaustively check all permutations, this is close to a lower bound as well on the the algorithm's time complexity. (The only reason we cannot claim to have a precise value for the big-$Θ$ complexity is that the bandwidth function is not exactly $Θ(n²)$, although it is close.)

Examples

The algorithm always iterates over all possible permutations, so it is infeasible to go above $9×9$ or $10×10$ without incurring multiple-hour runtimes. Nevertheless, we see that it is quite effective for, say, $8×8$:

julia> using Random, SparseArrays

julia> Random.seed!(628318);

julia> (n, p) = (8, 0.2);

julia> A = sprand(n, n, p);

julia> A = A + A' # Ensure structural symmetry;

julia> minimize_bandwidth(A, Minimization.BruteForceSearch())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Brute-force search
 * Approach: exact
 * Minimum Bandwidth: 3
 * Original Bandwidth: 6
 * Matrix Size: 8×8

Notes

Brute force is by far the slowest approach to matrix bandwidth minimization and should only be used in very niche cases (like verifying the correctness of other algorithms in unit tests). For $10×10$ matrices, the algorithm already takes several minutes to run (between $2$ to $5$ on most commercial machines) and allocates over $4$ gigabytes of memory. Given the $O(n! ⋅ n²)$ time complexity, minimizing the bandwidth of any $11×11$ matrix would take over an hour.

CapraraSalazarGonzalez <: ExactSolver <: AbstractSolver <: AbstractAlgorithm

The Caprara–Salazar-González minimization algorithm is an exact method for minimizing the bandwidth of a structurally symmetric matrix $A$. For a fixed $k ∈ ℕ$, the algorithm performs a depth-first search of all partial orderings of the rows and columns of $A$, adding indices one at a time. Partial orderings are pruned not only by ensuring that adjacent pairs of currently placed indices are within $k$ of each other but also by employing a branch-and-bound framework with lower bounds on bandwidth compatibility computed via integer linear programming relaxations. This search is repeated with incrementing values of $k$ until a bandwidth-$k$ ordering is found [CS05], with $k$ initialized to some lower bound on the minimum bandwidth of $A$ up to symmetric permutation.

Specifically, this implementation of the Caprara–Salazar-González algorithm uses the $min(α(A), γ(A))$ lower bound from the original paper [CS05, pp. 359–60] as the initial value of $k$. (Further implementation details can be found in the source code for bandwidth_lower_bound.)

As noted above, the Caprara–Salazar-González algorithm requires structurally symmetric input (that is, $A[i, j]$ must be nonzero if and only if $A[j, i]$ is nonzero for $1 ≤ i, j ≤ n$).

Supertype Hierarchy

CapraraSalazarGonzalez <: ExactSolver <: AbstractSolver <: MatrixBandwidth.AbstractAlgorithm

Performance

Given an $n×n$ input matrix, the Caprara–Salazar-González algorithm runs in $O(n! ⋅ n ⋅ Tᵢₗₚ(n, n²))$ time, where $Tᵢₗₚ(n, m)$ is the time taken to solve an integer linear programming (ILP) problem with $O(n)$ variables and $O(m)$ constraints:

• For each underlying "bandwidth ≤ $k$" check, we perform a depth-first search of $O(n!)$ partial orderings.
• At each search node, we solve ILP relaxations with $n$ variables and $O(n²)$ constraints (given by the number of nonzero entries in the computed distance matrix), taking $Tᵢₗₚ(n, n²)$ time. (This dominates the $O(n²)$ auxiliary computations needed to set up the ILP.) Thus, the overall time complexity for each value of $k$ is $O(n! ⋅ Tᵢₗₚ(n, n²))$.
• The difference between the maximum possible bandwidth ($n - 1$) and our initial lower bound grows linearly in $n$, so we run the underlying $O(n! ⋅ Tᵢₗₚ(n, n²))$ recognition algorithm $O(n)$ times.
• Therefore, the overall time complexity is $O(n! ⋅ n ⋅ Tᵢₗₚ(n, n²))$.

Note that $Tᵢₗₚ(n, n²)$ has worst-case complexity $O(2ⁿ)$, although this ultimately depends on the ILP solver used. (Here, we use the HiGHS solver from the HiGHS.jl package.)

Of course, this is all but an upper bound on the time complexity of Caprara–Salazar-González, achieved only in the most pathological of cases. In practice, efficient pruning techniques and compatibility checks—along with [CS05, pp. 359–60]'s relatively tight initial lower bound on the minimum bandwidth—result in approximately exponential growth in time complexity with respect to $n$.

Examples

We verify the optimality of the ordering found by Caprara–Salazar-González for a random $8×8$ matrix via a brute-force search over all possible permutations up to reversal:

julia> using Random, SparseArrays

julia> Random.seed!(5883);

julia> (n, p) = (8, 0.25);

julia> A = sprand(n, n, p);

julia> A = A + A' # Ensure structural symmetry;

julia> res_bf = minimize_bandwidth(A, Minimization.BruteForceSearch())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Brute-force search
 * Approach: exact
 * Minimum Bandwidth: 4
 * Original Bandwidth: 7
 * Matrix Size: 8×8

julia> res_csg = minimize_bandwidth(A, Minimization.CapraraSalazarGonzalez())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Caprara–Salazar-González
 * Approach: exact
 * Minimum Bandwidth: 4
 * Original Bandwidth: 7
 * Matrix Size: 8×8

Notes

For readers of the original paper, what we call the Caprara–Salazar-González algorithm here is designated the LAYOUT_LEFT_TO_RIGHT algorithm in [CS05]. The paper also describes a LAYOUT_BOTH_WAYS algorithm that performs a bidirectional search by adding indices to both the left and right ends of the current partial ordering. However, this version is considerably more complex to implement, and we ran into problems enforcing ILP constraints on node pairs added to opposite ends of the ordering. In any case, computational results demonstrate that neither LAYOUT_LEFT_TO_RIGHT nor LAYOUT_BOTH_WAYS is consistently faster, and the paper states that there is no known heuristic for determining which version will be more performant for a given input [CS05, pp. 368–69]. Therefore, we opt to implement only LAYOUT_LEFT_TO_RIGHT as a matter of practicality, although future developers may wish to extend the interface with LAYOUT_BOTH_WAYS as well.

A final implementation detail worth noting is that we use HiGHS as our solver; it is one of the fastest open-source solvers available for mixed-integer linear programming.

References

• [CS05]: A. Caprara and J.-J. Salazar-González. Laying Out Sparse Graphs with Provably Minimum Bandwidth. INFORMS Journal on Computing 17, 356–73 (2005). https://doi.org/10.1287/ijoc.1040.0083.

DelCorsoManzini <: ExactSolver <: AbstractSolver <: AbstractAlgorithm

The Del Corso–Manzini minimization algorithm is an exact method for minimizing the bandwidth of a structurally symmetric matrix $A$. For a fixed $k ∈ ℕ$, the algorithm invokes a subroutine that determines whether $A$ has bandwidth at most $k$ up to symmetric permutation. This subroutine performs a depth-first search of all partial orderings of the rows and columns of $A$, adding indices one at a time. Partial orderings are pruned not only by ensuring that adjacent pairs of currently placed indices are within $k$ of each other but also by tracking the latest positions at which the remaining indices can be placed. This search is repeated with incrementing values of $k$ until a bandwidth-$k$ ordering is found [DM99], with $k$ initialized to some lower bound on the minimum bandwidth of $A$ up to symmetric permutation.

Specifically, this implementation of the Del Corso–Manzini algorithm uses the $min(α(A), γ(A))$ lower bound from [CS05, pp. 359–60] as the initial value of $k$. (Further implementation details can be found in the source code for bandwidth_lower_bound.) This improves upon the original algorithm, which used the maximum number of nonzero off-diagonal entries in a single row as a lower bound on the minimum bandwidth of $A$ up to symmetric permutation [DM99, p. 192–93].

As noted above, the Del Corso–Manzini algorithm requires structurally symmetric input (that is, $A[i, j]$ must be nonzero if and only if $A[j, i]$ is nonzero for $1 ≤ i, j ≤ n$).

Supertype Hierarchy

DelCorsoManzini <: ExactSolver <: AbstractSolver <: MatrixBandwidth.AbstractAlgorithm

Performance

Given an $n×n$ input matrix, the Del Corso–Manzini algorithm runs in $O(n! ⋅ n³)$ time:

• For each underlying "bandwidth ≤ $k$" check, we perform a depth-first search of $O(n!)$ partial orderings.
• Checking plausibility of each partial ordering takes $O(nk)$ time, resulting in $O(n! ⋅ nk)$ steps for each value of $k$.
• The difference between the maximum possible bandwidth ($n - 1$) and our initial lower bound grows linearly in $n$, so we run the underlying $O(n! ⋅ nk)$ recognition algorithm $O(n)$ times.
• Finally, $∑ₖ₌₀ⁿ⁻¹ nk = O(n³)$, so the overall time complexity is $O(n! ⋅ n³)$.

Of course, this is but an upper bound on the time complexity of Del Corso–Manzini, achieved only in the most pathological of cases. In practice, efficient pruning techniques and compatibility checks—along with [CS05, pp. 359–60]'s relatively tight initial lower bound on the minimum bandwidth—result in approximately exponential growth in time complexity with respect to $n$.

Based on experimental results, the algorithm is feasible for $n×n$ matrices up to $n ≈ 100$ or so.

Examples

We verify the optimality of the ordering found by Del Corso–Manzini for a random $9×9$ matrix via a brute-force search over all possible permutations up to reversal:

julia> using Random, SparseArrays

julia> Random.seed!(0117);

julia> (n, p) = (9, 0.5);

julia> A = sprand(n, n, p);

julia> A = A + A' # Ensure structural symmetry;

julia> res_bf = minimize_bandwidth(A, Minimization.BruteForceSearch())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Brute-force search
 * Approach: exact
 * Minimum Bandwidth: 5
 * Original Bandwidth: 8
 * Matrix Size: 9×9

julia> res_dcm = minimize_bandwidth(A, Minimization.DelCorsoManzini())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Del Corso–Manzini
 * Approach: exact
 * Minimum Bandwidth: 5
 * Original Bandwidth: 8
 * Matrix Size: 9×9

We now generate (and shuffle) a random $40×40$ matrix with minimum bandwidth $10$ using MatrixBandwidth.random_banded_matrix. Del Corso–Manzini then finds a bandwidth-$10$ ordering, which is (we claim) optimal up to symmetric permutation. (In some cases, random_banded_matrix(n, k) does generate matrices with minimum bandwidth < k. Nevertheless, this example demonstrates that Del Corso–Manzini at the very least finds a quite good ordering, even though exact optimality—which is guaranteed by the original paper [DM99]—is not explicitly verified.)

julia> using Random

julia> Random.seed!(0201);

julia> (n, k) = (40, 10);

julia> A = MatrixBandwidth.random_banded_matrix(n, k);

julia> perm = randperm(n);

julia> A_shuffled = A[perm, perm];

julia> bandwidth(A)
10

julia> bandwidth(A_shuffled) # Much larger after shuffling
36

julia> minimize_bandwidth(A_shuffled, Minimization.DelCorsoManzini())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Del Corso–Manzini
 * Approach: exact
 * Minimum Bandwidth: 10
 * Original Bandwidth: 36
 * Matrix Size: 40×40

Notes

For readers of the original paper, what we call the Del Corso–Manzini minimization algorithm here is designated the "MB-ID algorithm" in [DM99, p. 191]. The so-called "MB-PS algorithm," on the other hand, we implement in DelCorsoManziniWithPS.

References

• [CS05]: A. Caprara and J.-J. Salazar-González. Laying Out Sparse Graphs with Provably Minimum Bandwidth. INFORMS Journal on Computing 17, 356–73 (2005). https://doi.org/10.1287/ijoc.1040.0083.
• [DM99]: G. M. Del Corso and G. Manzini. Finding Exact Solutions to the Bandwidth Minimization Problem. Computing 62, 189–203 (1999). https://doi.org/10.1007/s006070050002.

DelCorsoManziniWithPS{D} <: ExactSolver <: AbstractSolver <: AbstractAlgorithm

The Del Corso–Manzini minimization algorithm with perimeter search is an exact method for minimizing the bandwidth of a structurally symmetric matrix $A$. The base Del Corso–Manzini algorithm performs a depth-first search of all partial orderings of the rows and columns of $A$ for some fixed $k ∈ ℕ$, adding indices one at a time. Partial orderings are pruned not only by ensuring that adjacent pairs of currently placed indices are within $k$ of each other but also by tracking the latest positions at which the remaining indices can be placed. This search is repeated with incrementing values of $k$ until a bandwidth-$k$ ordering is found [DM99], with $k$ initialized to some lower bound on the minimum bandwidth of $A$ up to symmetric permutation.

The incorporation of perimeter search to this approach entails precomputing a "perimeter" of $d$-permutations of row indices of $A$, where $d$ is a positive integer passed as a parameter to the solver. Each permutation represents a way to select the last $d$ entries of the ordering, and as the construction of the partial ordering progresses, potential endings are pruned to exclude those incompatible with already placed indices. In addition to pruning a potential ending if it contains indices already placed, compatibility is also checked via precomputed time stamps indicating, for each potential ending, a loose lower bound on the earliest position at which any given index can be placed should said ending be selected.

Like our implementation of the base Del Corso–Manzini algorithm (see DelCorsoManzini), this implementation uses the $min(α(A), γ(A))$ lower bound from [CS05, pp. 359–60] as the initial value of $k$. (Further implementation details can be found in the source code for bandwidth_lower_bound.) This improves upon the original algorithm, which used the maximum number of nonzero off-diagonal entries in a single row as a lower bound on the minimum bandwidth of $A$ up to symmetric permutation [DM99, p. 194].

As noted above, the Del Corso–Manzini algorithm with perimeter search requires structurally symmetric input (that is, $A[i, j]$ must be nonzero if and only if $A[j, i]$ is nonzero for $1 ≤ i, j ≤ n$).

Fields

• depth::D<:Union{Nothing,Int}: the perimeter search depth. If this field is not set (and) thus automatically initialized to nothing), a default depth is automatically computed by Recognition.dcm_ps_optimal_depth as a function of the input matrix every time the solver is passed to MatrixBandwidth.Minimization.minimize_bandwidth. Otherwise, it must be manually set to a positive integer.

Constructors

• DelCorsoManziniWithPS(): constructs a new DelCorsoManziniWithPS instance with the default perimeter search depth initialized to nothing.
• DelCorsoManziniWithPS(depth::Integer): constructs a new DelCorsoManziniWithPS instance with the specified perimeter search depth. depth must be a positive integer.

Supertype Hierarchy

DelCorsoManziniWithPS <: ExactSolver <: AbstractSolver <: MatrixBandwidth.AbstractAlgorithm

Performance

Given an $n×n$ input matrix and perimeter search depth $d$, the Del Corso–Manzini algorithm with perimeter search runs in $O(n! ⋅ nᴰ⁺¹)$ time, where $Dᴰ = max(d, 2)$:

• For each underlying "bandwidth ≤ $k$" check, we perform a depth-first search of $O(n!)$ partial orderings.
• Checking plausibility of each partial ordering takes $O(nk)$ time, and checking compatibility with all size-$d$ LPOs takes $O(nᵈ)$ time. Thus, the overall time complexity for each value of $k$ is $O(n! ⋅ (nᵈ + nk))$.
• The difference between the maximum possible bandwidth ($n - 1$) and our initial lower bound grows linearly in $n$, so we run the underlying $O(n! ⋅ (nᵈ + nk))$ recognition algorithm $O(n)$ times.
• Finally, $∑ₖ₌₀ⁿ⁻¹ (nᵈ + nk) = O(nᵈ⁺¹ + n³)$, so the overall time complexity is $O(n! ⋅ nᴰ⁺¹)$, where $D = max(d, 2)$.

Of course, this is but an upper bound on the time complexity of Del Corso–Manzini with perimeter search, achieved only in the most pathological of cases. In practice, efficient pruning techniques and compatibility checks—along with [CS05, pp. 359–60]'s relatively tight initial lower bound on the minimum bandwidth—result in approximately exponential growth in time complexity with respect to $n$.

Based on experimental results, the algorithm is feasible for $n×n$ matrices up to $n ≈ 100$ or so.

Examples

We verify the optimality of the ordering found by Del Corso–Manzini with perimeter search for a random $9×9$ matrix via a brute-force search over all possible permutations up to reversal. The depth parameter is not explicitly set; instead, some near-optimal value is automatically computed upon the first MatrixBandwidth.Minimization.minimize_bandwidth function call.

julia> using Random, SparseArrays

julia> Random.seed!(548836);

julia> (n, p) = (9, 0.2);

julia> A = sprand(n, n, p);

julia> A = A + A' # Ensure structural symmetry;

julia> res_bf = minimize_bandwidth(A, Minimization.BruteForceSearch())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Brute-force search
 * Approach: exact
 * Minimum Bandwidth: 3
 * Original Bandwidth: 8
 * Matrix Size: 9×9

julia> res_dcm_ps = minimize_bandwidth(A, Minimization.DelCorsoManziniWithPS())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Del Corso–Manzini with perimeter search
 * Approach: exact
 * Minimum Bandwidth: 3
 * Original Bandwidth: 8
 * Matrix Size: 9×9

We now generate (and shuffle) a random $30×30$ matrix with minimum bandwidth $8$ using MatrixBandwidth.random_banded_matrix. Del Corso–Manzini with perimeter search then finds a bandwidth-$8$ ordering, which is (we claim) optimal up to symmetric permutation. (In some cases, random_banded_matrix(n, k) does generate matrices with minimum bandwidth < k. Nevertheless, this example demonstrates that Del Corso–Manzini at the very least finds a quite good ordering, even though exact optimality—which is guaranteed by the original paper [DM99]—is not explicitly verified.) In this case, we set the depth parameter to $4$ beforehand instead of relying on Recognition.dcm_ps_optimal_depth.

julia> using Random

julia> Random.seed!(78779);

julia> (n, k, depth) = (30, 8, 4);

julia> A = MatrixBandwidth.random_banded_matrix(n, k);

julia> perm = randperm(n);

julia> A_shuffled = A[perm, perm];

julia> bandwidth(A)
8

julia> bandwidth(A_shuffled) # Much larger after shuffling
25

julia> minimize_bandwidth(A_shuffled, Minimization.DelCorsoManziniWithPS(depth))
Results of Bandwidth Minimization Algorithm
 * Algorithm: Del Corso–Manzini with perimeter search
 * Approach: exact
 * Minimum Bandwidth: 8
 * Original Bandwidth: 25
 * Matrix Size: 30×30

Notes

For readers of the original paper, what we call the Del Corso–Manzini minimization algorithm with perimeter search here is designated the "MB-PS algorithm" in [DM99. p. 193]. The so-called "MB-ID algorithm," on the other hand, we implement in DelCorsoManzini.

References

• [CS05]: A. Caprara and J.-J. Salazar-González. Laying Out Sparse Graphs with Provably Minimum Bandwidth. INFORMS Journal on Computing 17, 356–73 (2005). https://doi.org/10.1287/ijoc.1040.0083.
• [DM99]: G. M. Del Corso and G. Manzini. Finding Exact Solutions to the Bandwidth Minimization Problem. Computing 62, 189–203 (1999). https://doi.org/10.1007/s006070050002.

ExactSolver <: AbstractSolver <: AbstractAlgorithm

Abstract type for all exact matrix bandwidth minimization solvers.

Exact methods are those which guarantee an optimal ordering producing the true minimum bandwidth of a matrix. Since bandwidth minimization is an NP-complete problem, existing exact algorithms are, at best, exponential in time complexity—much worse than many polynomial-time heuristic approaches (e.g., Gibbs–Poole–Stockmeyer). Such methods, therefore, are not feasible for large matrices, but they remain useful when precise solutions are required for small-to-medium-sized inputs (say, up to $100×100$).

Supertype Hierarchy

ExactSolver <: AbstractSolver <: MatrixBandwidth.AbstractAlgorithm

SaxeGurariSudborough <: ExactSolver <: AbstractSolver <: AbstractAlgorithm

The Saxe–Gurari–Sudborough minimization algorithm is an exact method for minimizing the bandwidth of a structurally symmetric matrix $A$. For a fixed $k ∈ ℕ$, the algorithm invokes a subroutine that determines whether $A$ has bandwidth at most $k$ up to symmetric permutation. This subroutine employs dynamic programming to search over equivalence classes of partial orderings, where two partial orderings of length $l$ are equivalent if they share the same active region. (The active region of a partial ordering is defined as the sequence of the last $min(k, l)$ vertices in the ordering taken together with all dangling edges—edges with one endpoint in the ordering and the other endpoint not yet in the ordering.) It extends these partial layouts one vertex at a time in a breadth-first manner, pruning implausible classes that violate bandwidth-$k$ constraints such as degree bounds on active vertices and excessive numbers of dangling edges [GS84]. This search is repeated with incrementing values of $k$ until a bandwidth-$k$ ordering is found, with $k$ initialized to some lower bound on the minimum bandwidth of $A$ up to symmetric permutation.

Specifically, this implementation of the Saxe–Gurari–Sudborough algorithm uses the $min(α(A), γ(A))$ lower bound from [CS05, pp. 359–60] as the initial value of $k$. (Further implementation details can be found in the source code for bandwidth_lower_bound.)

As noted above, the Saxe–Gurari–Sudborough algorithm requires structurally symmetric input (that is, $A[i, j]$ must be nonzero if and only if $A[j, i]$ is nonzero for $1 ≤ i, j ≤ n$).

Supertype Hierarchy

SaxeGurariSudborough <: ExactSolver <: AbstractSolver <: MatrixBandwidth.AbstractAlgorithm

Performance

Given an $n×n$ input matrix, the Saxe–Gurari–Sudborough algorithm runs in $O(nⁿ⁻¹)$ time:

• For each underlying "bandwidth ≤ k" check, we call the Saxe–Gurari–Sudborough recognition algorithm, which runs in $O(nᵏ)$ time [GS84, p. 531]. (This is an improvement upon the original $O(nᵏ⁺¹)$ Saxe algorithm [Sax80, p. 363].)
• The difference between the maximum possible bandwidth ($n - 1$) and our initial lower bound grows linearly in $n$, so we run the underlying $O(nᵏ)$ recognition algorithm $O(n)$ times.
• Therefore, the overall time complexity is $∑ₖ₌₀ⁿ⁻¹ nᵏ = O(nⁿ⁻¹)$.

Whereas most exact bandwidth minimization algorithms are technically factorial-time (with respect to $n$) in the worst case but practically always approximate exponential time complexity in real life, the $O(nⁿ⁻¹)$ upper bound on Saxe–Gurari–Sudborough is typically a good representation of actual performance in most cases. Indeed, these other types of algorithms tend to outperform Saxe–Gurari–Sudborough for larger $n$, given that their aggressive pruning strategies keep their effective search space very small in practice and $O(2ⁿ)$ ⊂ $O(nⁿ⁻¹)$.

Examples

We verify the optimality of the ordering found by Saxe–Gurari–Sudborough for a random $9×9$ matrix via a brute-force search over all possible permutations up to reversal:

julia> using Random, SparseArrays

julia> Random.seed!(52452);

julia> (n, p) = (9, 0.5);

julia> A = sprand(n, n, p);

julia> A = A + A' # Ensure structural symmetry;

julia> res_bf = minimize_bandwidth(A, Minimization.BruteForceSearch())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Brute-force search
 * Approach: exact
 * Minimum Bandwidth: 5
 * Original Bandwidth: 8
 * Matrix Size: 9×9

julia> res_sgs = minimize_bandwidth(A, Minimization.SaxeGurariSudborough())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Saxe–Gurari–Sudborough
 * Approach: exact
 * Minimum Bandwidth: 5
 * Original Bandwidth: 8
 * Matrix Size: 9×9

We now generate (and shuffle) a random $25×25$ matrix with minimum bandwidth $5$ using MatrixBandwidth.random_banded_matrix. Saxe–Gurari–Sudborough then finds a bandwidth-$4$ ordering, which is (we claim) optimal up to symmetric permutation. (In some cases, random_banded_matrix(n, k) generates matrices with minimum bandwidth < k—this appears to be one such case. Although we do not explicitly verify exact optimality—which is guaranteed by the original paper [GS84]—here via brute-force search, this example demonstrates that Saxe–Gurari–Sudborough at the very least finds a quite good ordering.)

julia> using Random

julia> Random.seed!(937497);

julia> (n, k, p) = (25, 5, 0.25);

julia> A = MatrixBandwidth.random_banded_matrix(n, k; p=p);

julia> perm = randperm(n);

julia> A_shuffled = A[perm, perm];

julia> bandwidth(A)
5

julia> bandwidth(A_shuffled) # Much larger after shuffling
19

julia> minimize_bandwidth(A_shuffled, Minimization.SaxeGurariSudborough())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Saxe–Gurari–Sudborough
 * Approach: exact
 * Minimum Bandwidth: 4
 * Original Bandwidth: 19
 * Matrix Size: 25×25

Notes

The Saxe–Gurari–Sudborough algorithm was originally a bandwidth recognition algorithm, not a minimization one—as previously mentioned, we repurpose it here by repeatedly invoking the original procedure for incrementing values of $k$ until a valid ordering is found. The general family of recognition algorithms to which Saxe–Gurari–Sudborough belongs was conceived as a response to a question posed by [GGJK78, p. 494]: is the "bandwidth ≤ k?" problem NP-complete for arbitrary $k$? [Sax80] answered this question in the negative by providing a $O(nᵏ⁺¹)$ algorithm, constructively proving that the problem is class P. Later, [GS84] improved upon this algorithm by reducing time complexity to $O(nᵏ)$. Whereas the original Saxe algorithm considers extensions of partial orderings with any remaining unplaced vertex (of which there are $O(n)$ at any point in the breadth-first search), the Gurari–Sudborough refinement only considers extensions with vertices reachable by paths beginning with a dangling edge that never again traverse a dangling edge [GS84, pp. 535–36].

References

• [GGJK78]: M. R. Garey, R. L. Graham, D. S. Johnson and D. E. Knuth. Complexity Results for Bandwidth Minimization. SIAM Journal on Applied Mathematics 34, 477–95 (1978). https://doi.org/10.1137/0134037.
• [GS84]: E. M. Gurari and I. H. Sudborough. Improved dynamic programming algorithms for bandwidth minimization and the MinCut Linear Arrangement problem. Journal of Algorithms 5, 531–46 (1984). https://doi.org/10.1016/0196-6774(84)90006-3.
• [Sax80]: J. B. Saxe. Dynamic-Programming Algorithms for Recognizing Small-Bandwidth Graphs in Polynomial Time. SIAM Journal on Algebraic and Discrete Methods 1, 363–69 (1980). https://doi.org/10.1137/0601042.

MatrixBandwidth.Minimization.Heuristic

Heuristic solvers for matrix bandwidth minimization.

Heuristic methods are those which aim to produce near-optimal solutions in a more performant manner than exact methods. While precise bandwidth minimization is NP-complete, many heuristic algorithms (such as Gibbs–Poole–Stockmeyer) run in polynomial time.

Heuristic algorithms differ from metaheuristic ones in that they do not employ higher-level iterative search frameworks (e.g., stochastic techniques) to survey the global search space and escape local minima; instead, they rely on straightforward deterministic procedures to find good solutions in a single pass.

The following heuristic matrix bandwidth minimization algorithms are currently available:

• Gibbs–Poole–Stockmeyer (GibbsPooleStockmeyer)
• Cuthill–McKee (CuthillMcKee)
• Reverse Cuthill–McKee (ReverseCuthillMcKee)

This submodule is part of the MatrixBandwidth.Minimization submodule of the MatrixBandwidth.jl package.

CuthillMcKee <: HeuristicSolver <: AbstractSolver <: AbstractAlgorithm

The Cuthill–McKee algorithm is a heuristic method for minimizing the bandwidth of a structurally symmetric matrix $A$. It considers the graph $G(A)$ whose adjacency matrix is $A$ (ignoring weights and self-loops) and performs a breadth-first search of each connected component of $G(A)$, starting from a low-degree node then visiting its neighbors in order of increasing degree. Particularly effective when $A$ is sparse, this heuristic typically produces an ordering which induces a matrix bandwidth either equal to or very close to the true minimum [CM69, pp. 157–58].

As noted above, the Cuthill–McKee algorithm requires structurally symmetric input (that is, $A[i, j]$ must be nonzero if and only if $A[j, i]$ is nonzero for $1 ≤ i, j ≤ n$).

Fields

• node_finder::Function: a function that selects a node from some connected component of the input matrix from which to start the breadth-first search. If no custom heuristic is specified, this field defaults to bi_criteria_node_finder, which picks a node "farthest" from the others in the component (not necessarily the lowest-degree node).

Supertype Hierarchy

CuthillMcKee <: HeuristicSolver <: AbstractSolver <: MatrixBandwidth.AbstractAlgorithm

Performance

Given an $n×n$ input matrix $A$, the Cuthill–McKee algorithm runs in $O(n²)$ time.

[CG80] provide a linear-time implementation in the number of nonzero entries of $A$, which is still quadratic when $A$ is dense but often much faster when dealing with sparse matrices. However, this would require that $A$ be stored as a graph or a sparse matrix, which runs counter to our desire to provide a bandwidth minimization API for all AbstractMatrix{<:Number} types, including dense matrices. (In the future, however, we may indeed consider supporting this more performant implementation for sparse matrices.)

It was found in [Geo71, pp. 114–15] that reversing the ordering produced by Cuthill–McKee tends to induce a more optimal matrix profile (a measure of how far, on average, nonzero entries are from the diagonal; see also MatrixBandwidth.profile). This so-called reverse Cuthill–McKee variant is preferred in almost all cases—see ReverseCuthillMcKee and the associated method of _minimize_bandwidth_impl for our implementation.

Examples

In the following examples, MatrixBandwidth.random_banded_matrix is used to generate random matrices with minimum bandwidth close to $k$. In some cases, however, the true minimum bandwidth up to symmetric permutation may be even less than $k$, making it hard to verify whether Cuthill–McKee finds a truly optimal ordering or simply a near-optimal one. Nevertheless, the results are still very good in practice.

Cuthill–McKee finds a good ordering for a $30×30$ matrix:

julia> using Random

julia> Random.seed!(13);

julia> (n, k) = (30, 5);

julia> A = MatrixBandwidth.random_banded_matrix(n, k);

julia> perm = randperm(n);

julia> A_shuffled = A[perm, perm];

julia> bandwidth(A)
5

julia> bandwidth(A_shuffled) # Much larger after shuffling
25

julia> minimize_bandwidth(A_shuffled, Minimization.CuthillMcKee())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Cuthill–McKee
 * Approach: heuristic
 * Minimum Bandwidth: 8
 * Original Bandwidth: 25
 * Matrix Size: 30×30

Cuthill–McKee finds a good ordering for a structurally symmetric $183×183$ matrix with multiple (separate) connected components:

julia> using Random, SparseArrays

julia> Random.seed!(37452);

julia> (max_cc_size, max_band, p, num_ccs) = (60, 9, 0.2, 7);

julia> components = Vector{SparseMatrixCSC{Float64, Int64}}(undef, num_ccs);

julia> for i in 1:num_ccs # Some components may themselves be disconnected
           cc_size = rand(1:max_cc_size);
           cc_band = rand(0:min(max_band, cc_size - 1));
           components[i] = sparse(
               MatrixBandwidth.random_banded_matrix(cc_size, cc_band; p=p)
           );
       end

julia> A = blockdiag(components...); # `A` has least 7 connected components

julia> perm = randperm(sum(map(cc -> size(cc, 1), components)));

julia> A_shuffled = A[perm, perm];

julia> res = minimize_bandwidth(A_shuffled, Minimization.CuthillMcKee());

julia> A # The original matrix
276×276 SparseMatrixCSC{Float64, Int64} with 464 stored entries:
⎡⢾⡷⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠘⢻⣲⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠘⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠿⡧⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠉⢯⡷⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠚⣤⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⢻⣶⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠯⡧⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠛⣤⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠛⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⣤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠱⣢⡀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⡢⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢴⣷⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠿⣧⡀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢿⡷⎦

julia> A_shuffled # A far-from-optimal ordering of `A`
276×276 SparseMatrixCSC{Float64, Int64} with 464 stored entries:
⎡⠁⢄⠀⢀⠀⠀⠀⢀⠠⠀⠀⠐⠀⠀⠀⠐⢀⡐⠀⠀⠀⢀⠀⠀⠀⠀⠐⠀⢠⠀⠀⠀⡄⠀⠀⠐⠀⠀⠂⠄⎤
⎢⠀⢀⠱⠂⠀⠀⠀⠈⠀⠀⠀⠀⠀⠀⢨⠀⠀⠀⠀⠀⡀⠁⠠⠀⠘⠀⠀⠡⢀⠈⠀⠀⠀⠀⠀⠀⠄⠀⠁⠁⎥
⎢⠀⠀⠀⠀⠑⢀⠀⠂⠀⠀⠀⠀⢐⠀⠀⠠⠈⠠⠀⠀⠀⠐⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⢢⢀⢀⠀⎥
⎢⠀⢀⡀⠀⠠⠀⠁⠄⠀⠠⠀⠄⠀⠀⠀⠄⠀⠀⠀⠀⢀⠀⠀⢀⠀⠑⠀⠀⠐⠠⠀⠀⠠⠨⠂⠀⠀⠀⠀⠀⎥
⎢⠀⠂⠀⠀⠀⠀⠀⡀⠱⢆⡀⠂⠀⠀⠀⠀⠀⠀⢀⢊⠀⠐⠐⠈⠀⠈⠀⢀⠄⠀⡀⠀⢁⢀⠠⠀⠃⠀⠊⠀⎥
⎢⢀⠀⠀⠀⠀⠀⠀⠄⠠⠈⠑⠀⢀⠐⠀⠌⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠀⢀⠉⢀⠀⠠⠈⠀⠀⣁⠁⎥
⎢⠀⠀⠀⠀⠐⠐⠀⠀⠀⠀⢀⠐⠁⠄⠈⠀⢌⠀⠆⠠⢀⠀⠄⠐⠰⠀⠀⠀⠁⠰⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⢀⠀⠂⠒⠀⡀⠀⠄⠀⠀⡀⠄⠂⠀⠐⢄⠁⢀⠀⠀⠀⡀⠀⠀⠀⠀⡠⠀⠀⠀⠀⠀⠀⠀⠀⢈⠀⠀⠀⠁⎥
⎢⢀⠰⠀⠀⠂⡀⠀⠀⠀⠀⠀⠈⠂⠑⠁⢀⠐⠄⠄⠂⠂⠜⠄⠀⠀⠀⡄⠀⠀⢀⠀⠠⠀⢀⠄⠀⢀⠀⠂⡂⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⡠⢐⠀⠀⠈⡁⠀⠀⠠⠁⠀⢀⠀⠀⠀⠀⡀⠀⠀⢀⠀⠈⠃⠀⠸⠈⠠⠀⠀⠀⢄⠂⎥
⎢⠀⢀⠄⠈⢀⠀⠀⠐⢀⠀⠀⠀⠀⠐⠀⠠⣈⠄⠀⠀⠐⢀⠀⡀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠊⠀⠠⠀⠐⠀⠀⎥
⎢⠀⠀⠀⠂⢀⠀⠀⢀⡐⠀⠀⠀⢀⠁⠀⠀⠀⠁⠀⠀⠀⠠⠄⣥⠉⠈⠀⠀⠀⠀⡀⠀⠀⠀⠀⠀⡀⠀⠀⠀⎥
⎢⠀⠀⠒⠀⠀⠀⢄⠀⡀⠀⠀⠀⠐⠂⠀⠀⠀⠀⠀⠈⠀⠀⡃⠀⠁⢀⠀⠀⢀⡀⢈⠈⠀⠀⠀⠂⠀⠠⠂⠂⎥
⎢⠐⠀⠄⡀⠀⠀⠀⠀⠀⢀⠀⠈⠀⠀⠀⠊⠀⠉⠀⢀⠀⠀⠀⠀⠀⠀⠑⠀⠀⠀⠀⠀⢀⠀⠈⠀⠛⠃⢄⠀⎥
⎢⠀⠒⡀⠐⠀⠀⠐⡀⠀⠁⠀⠀⢁⡀⠀⠀⠀⢀⡀⠀⠀⠀⠀⠀⠀⠰⠀⠀⠀⢄⠀⠰⠀⠠⠠⢀⠀⠀⢂⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⡄⠐⠂⠀⠀⠀⠀⡀⠉⠀⠐⠀⠀⠈⡂⠐⠀⠀⢀⡀⠀⣠⠀⠄⠠⠀⠀⡀⠀⠀⎥
⎢⠀⠉⠀⠀⠀⠀⡀⡂⠁⢐⠀⠐⠀⠀⠀⠀⠀⢀⡒⠂⡠⠀⠀⠀⠀⠀⠀⠐⠀⡀⠀⠄⠑⠄⠀⠀⠀⠀⠀⠀⎥
⎢⢀⠀⠀⠀⠀⡀⠈⠀⠀⠂⡀⠂⠀⠀⡀⢀⠀⠁⠀⠂⠀⡀⠀⠀⠠⠀⠂⠀⠀⢂⠀⠂⠀⠀⠁⢀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠁⠈⢒⠀⠀⠉⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⢀⠀⠀⠈⠀⡀⠿⠀⠀⠀⠀⠠⠀⠀⠀⠀⠱⠆⠀⠀⎥
⎣⠈⠄⠅⠀⠀⠐⠀⠀⠊⠀⠅⠘⠀⠀⠄⠀⠨⠠⠠⠑⠀⠀⠀⠀⠨⠀⠀⠑⠈⠐⠀⠀⠀⠀⠀⠀⠀⠀⠔⢅⎦

julia> A_shuffled[res.ordering, res.ordering] # A near-optimal reordering of `A_shuffled`
276×276 SparseMatrixCSC{Float64, Int64} with 464 stored entries:
⎡⠱⣦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠉⠻⣦⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠘⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⡦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⡦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠺⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠚⣤⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⣤⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⢄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⢀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⢄⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢄⠀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⎦

julia> bandwidth(A)
7

julia> bandwidth(A_shuffled) # Much larger after shuffling
266

julia> res # Even better than the original bandwidth (which was, clearly, not yet optimal)
Results of Bandwidth Minimization Algorithm
 * Algorithm: Cuthill–McKee
 * Approach: heuristic
 * Minimum Bandwidth: 5
 * Original Bandwidth: 266
 * Matrix Size: 276×276

Notes

Note that the node_finder field must be of the form (A::AbstractMatrix{Bool}) -> Integer (i.e., it must take in an boolean matrix and return an integer). If this is not the case, an ArgumentError is thrown upon construction.

References

• [CG80]: W. M. Chan and A. George. A linear time implementation of the reverse Cuthill–McKee algorithm. BIT Numerical Mathematics 20, 8–14 (1980). https://doi.org/10.1007/BF01933580.
• [CM69]: E. Cuthill and J. McKee. Reducing the bandwidth of sparse symmetric matrices. In: Proceedings of the 24th National Conference of the ACM (Brandon Systems Press, 1969); pp. 157–72. https://doi.org/10.1145/800195.805928.
• [Geo71]: J. A. George. Computer Implementation of the Finite Element Method. Ph.D. Thesis, Department of Computer Science, Stanford University (1971). https://apps.dtic.mil/sti/tr/pdf/AD0726171.pdf.

GibbsPooleStockmeyer <: HeuristicSolver <: AbstractSolver <: AbstractAlgorithm

The Gibbs–Poole–Stockmeyer algorithm is a heuristic method for minimizing the bandwidth of a structurally symmetric matrix $A$. It considers the graph $G(A)$ whose adjacency matrix is $A$ (ignoring weights and self-loops) and builds an ordering by identifying a pair of "endpoints" in the graph far from each other, constructing sets of levels from these endpoints, and merging these level structures in such a way that minimizes the size of the largest level in the final combined structure. Based on the classical reverse Cuthill–McKee algorithm [Geo71], this heuristic typically produces an ordering which induces a matrix bandwidth either equal to or very close to the true minimum, with improvements in bandwidths over reverse Cuthill–McKee more noticeable once input size exceeds $400×400$ or so [GPS76, pp. 246–47].

Whereas the original paper outlined a strategy for conditionally reversing the orderings of individual "connected components" (treating the input matrix $A$ as an undirected graph) [GPS76, p. 241], this implementation instead reverses the entire final ordering in every case, similarly to ReverseCuthillMcKee. Conditional reversals are not only more complex to implement but also slightly more time-consuming, with the only benefit being a marginally smaller matrix profile (a measure of how far, on average, nonzero entries are from the diagonal; see also MatrixBandwidth.profile). Since such reversal strategies do not affect matrix bandwidth (the primary focus of this package), we opt instead for the simpler unconditional reversal.

As noted above, the Gibbs–Poole–Stockmeyer algorithm requires structurally symmetric input (that is, $A[i, j]$ must be nonzero if and only if $A[j, i]$ is nonzero for $1 ≤ i, j ≤ n$).

Fields

• node_finder::Function: a function that selects a node from some connected component of the input matrix from which to start the breadth-first search. If no custom heuristic is specified, this field defaults to bi_criteria_node_finder, which picks a node "farthest" from the others in the component (not necessarily the lowest-degree node).

Supertype Hierarchy

GibbsPooleStockmeyer <: HeuristicSolver <: AbstractSolver <: MatrixBandwidth.AbstractAlgorithm

Performance

Given an $n×n$ input matrix, the Gibbs–Poole–Stockmeyer algorithm runs in $O(n²)$ time.

[Lew82] provides a notably faster and more memory-efficient implementation, relying on sparse storage of the input matrix. However, this would run counter to our desire to provide a bandwidth minimization API for all AbstractMatrix{<:Number} types, including dense matrices. (In the future, however, we may indeed consider supporting this more performant implementation for sparse matrices.)

On that note, Gibbs–Poole–Stockmeyer has been found to take considerably less time than reverse Cuthill–McKee when matrices are stored in sparse format [GPS76, pp. 246–47]. The dense-matrix implementations of both algorithms in this package, however, result in reverse Cuthill–McKee consistently outperforming Gibbs–Poole–Stockmeyer in terms of runtime (although Gibbs–Poole–Stockmeyer still typically produces lower-bandwidth orderings for larger matrices). This further motivates the desire to implement a sparse version of both algorithms in the future.

Examples

In the following examples, MatrixBandwidth.random_banded_matrix is used to generate random matrices with minimum bandwidth close to $k$. In some cases, however, the true minimum bandwidth up to symmetric permutation may be even less than $k$, making it hard to verify whether Gibbs–Poole–Stockmeyer finds a truly optimal ordering or simply a near-optimal one. Nevertheless, the results are still very good in practice.

Gibbs–Poole–Stockmeyer finds a good ordering for a $40×40$ matrix:

julia> using Random

julia> Random.seed!(561);

julia> (n, k) = (40, 7);

julia> A = MatrixBandwidth.random_banded_matrix(n, k);

julia> perm = randperm(n);

julia> A_shuffled = A[perm, perm];

julia> bandwidth(A)
7

julia> bandwidth(A_shuffled)
37

julia> minimize_bandwidth(A_shuffled, Minimization.GibbsPooleStockmeyer())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Gibbs–Poole–Stockmeyer
 * Approach: heuristic
 * Minimum Bandwidth: 7
 * Original Bandwidth: 37
 * Matrix Size: 40×40

Gibbs–Poole–Stockmeyer finds a good ordering for a $748×748$ matrix with multiple (separate) connected components:

julia> using Random, SparseArrays

julia> Random.seed!(271828);

julia> (max_cc_size, max_band, p, num_ccs) = (120, 13, 0.3, 11);

julia> components = Vector{SparseMatrixCSC{Float64, Int64}}(undef, num_ccs);

julia> for i in 1:num_ccs # Some components may themselves be disconnected
           cc_size = rand(0:max_cc_size);
           cc_band = rand(1:min(max_band, cc_size - 1));
           components[i] = sparse(
               MatrixBandwidth.random_banded_matrix(cc_size, cc_band; p=p)
           );
       end

julia> A = blockdiag(components...); # `A` has least 8 connected components

julia> perm = randperm(sum(map(cc -> size(cc, 1), components)));

julia> A_shuffled = A[perm, perm];

julia> res = minimize_bandwidth(A_shuffled, Minimization.GibbsPooleStockmeyer());

julia> A # The original matrix
748×748 SparseMatrixCSC{Float64, Int64} with 2526 stored entries:
⎡⢿⣷⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠘⠿⣧⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠘⢿⣷⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠉⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠛⣤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠛⣤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠛⣤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⣤⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⎦

julia> A_shuffled # A far-from-optimal ordering of `A`
748×748 SparseMatrixCSC{Float64, Int64} with 2526 stored entries:
⎡⠰⣦⢎⢪⢐⠆⣗⣔⠆⠀⠀⠠⢃⡦⠵⠸⡐⢌⠴⠤⣤⢅⢒⠰⡢⡄⡰⣄⢂⢊⠎⠀⢝⡀⣼⠤⠅⢒⠰⠢⎤
⎢⡪⣑⡋⢌⠈⢀⣖⣠⢠⡀⡂⠀⠐⠤⣠⠰⣃⠀⢊⢃⡨⡇⡝⠢⣀⠈⢬⢁⠽⡼⠐⡦⠤⠐⠤⠚⠪⢠⠦⢈⎥
⎢⠰⠔⠂⢀⡑⢌⠐⠀⠐⡒⠤⡅⠂⡱⡕⠈⡑⠹⠱⣨⢓⣑⠀⠂⠀⢬⢂⣅⡤⢑⠢⢁⠀⢒⡑⡬⡆⠚⠀⡀⎥
⎢⢙⢽⠘⣹⠐⠀⢱⣶⠂⠂⠈⡵⢰⡤⠖⠰⡂⣠⢨⢰⣛⠀⠰⡌⡐⠄⠠⡦⠜⠃⠐⣦⢂⠄⠷⢅⣉⠰⠿⢴⎥
⎢⠈⠁⠀⠲⢰⠠⠨⠀⣑⣼⠩⡃⠂⢰⢁⠓⢐⣘⣂⢦⠂⡐⠔⢄⡨⣃⠦⢁⢈⠀⡂⠀⡅⠈⠀⢌⠀⡀⠵⠫⎥
⎢⠀⡀⠈⠈⠄⠧⢆⡤⠧⠢⠕⢅⡍⢔⡴⠀⡀⠆⡨⠈⣄⠲⠌⠳⢁⠐⠠⠴⠀⠄⠔⠺⠁⡉⢂⣭⠠⡠⠀⣉⎥
⎢⠩⡴⠐⡄⢌⡠⠐⡶⢈⣀⢃⢍⠋⢄⡖⢩⢌⣖⠈⢀⣴⣙⡀⢓⠁⠠⢬⢈⠅⡤⠅⢰⣁⠌⣌⠆⠸⢀⠒⣁⎥
⎢⣑⡃⢀⡚⡑⠉⢘⡁⢥⠐⠐⠋⡜⣉⠑⢄⠍⣢⣓⢉⠊⢖⠀⢐⡀⠲⡈⢑⠀⠊⠀⠛⠀⢃⠘⠀⢁⣒⢀⢁⎥
⎢⡐⢌⠉⠘⣕⡈⠈⣨⣐⢰⠠⠌⢢⢵⠣⣡⢕⣱⣂⢭⢎⠜⠀⠉⡂⣃⢐⠅⠒⠔⡂⠨⡂⢳⠍⢤⠰⡠⣩⡏⎥
⎢⠐⡇⠮⢐⡑⣢⢂⣒⠨⣜⡂⠊⠂⢀⡝⢘⡌⣜⣱⢞⠂⠾⠊⢔⠙⣢⢭⣹⠑⡑⠀⡉⡄⠁⠀⠴⣇⣃⠃⡑⎥
⎢⠄⢟⠦⠮⢝⢰⠛⠘⢈⠠⢠⡙⣔⢻⢪⢄⣊⠕⣨⡄⢕⣵⢐⠊⠊⣱⡠⣢⠀⠣⠀⣍⠔⠀⢝⢄⣌⡄⡌⠆⎥
⎢⢘⡐⠳⡉⠠⠀⡐⠦⠐⢅⢦⡁⢤⢈⢀⢀⡄⠀⢊⢄⡰⠐⠑⢄⠴⡅⠁⢆⡠⣄⠤⢐⠈⡀⠺⠂⠀⢀⡴⡀⎥
⎢⠈⠮⡀⠘⡀⣄⠐⠌⠦⢪⢁⠐⠁⡀⢠⡈⠬⢨⠳⣠⢎⣠⠔⠧⠛⢄⡀⡹⠡⠌⠃⢀⠣⠁⠉⢀⠑⠁⢦⡉⎥
⎢⠐⢮⠆⢓⠌⢴⠠⡦⠌⢃⢀⡆⡂⢓⢆⢈⠔⠔⣇⣳⠠⣪⠡⢄⣄⡨⠛⢄⠑⢐⠐⡠⠪⠃⢤⢁⡝⠐⡀⢎⎥
⎢⡨⢐⣓⡧⢄⢋⠶⠁⠂⠐⠀⠄⠁⡥⡠⠀⢘⠄⢕⠠⠤⡀⠀⢮⡁⠆⢑⢀⠕⣥⣑⠼⡀⡅⠝⠢⠠⠄⠀⠔⎥
⎢⠊⠁⠰⡤⠌⢂⠰⣤⠈⠈⣰⡁⢁⣁⣤⠀⡈⡈⡄⠠⡄⢤⢀⢃⠉⢀⠐⡠⣑⡜⠑⢄⠀⡼⠠⠲⢈⠠⢃⠁⎥
⎢⠓⠱⢀⠃⢠⢀⠈⠔⡁⠉⡅⠠⡁⠜⠤⢀⢬⣈⠄⠉⠐⠁⠂⠠⠍⠂⠮⠂⠄⠬⣀⡤⠵⢇⠈⠀⠣⠡⠩⡎⎥
⎢⠒⡟⣠⠃⡑⡬⠝⢇⡀⢄⡌⣴⠢⠝⠒⠀⠃⣅⢀⡄⠓⢕⠺⠂⠃⢀⠄⢓⠳⡁⢠⡂⠂⠀⠛⢄⡀⢒⠒⠂⎥
⎢⢡⢁⠊⣂⣨⠉⢃⡘⠀⠠⠀⡢⠒⢂⢡⢰⠐⡢⠭⢹⠂⠽⠀⢀⠕⠀⢓⠉⠀⠆⠂⡐⠍⡂⢠⢈⠁⣤⢀⡈⎥
⎣⠰⡂⡈⢃⠀⠠⢛⣇⡵⡃⡄⢠⠜⢠⠄⢐⡧⠾⢍⠠⠢⠍⠐⠫⡌⠳⡠⢌⢀⠄⠍⠐⡣⠦⠸⠀⡀⠰⡛⢌⎦

julia> A_shuffled[res.ordering, res.ordering] # A near-optimal reordering of `A_shuffled`
748×748 SparseMatrixCSC{Float64, Int64} with 2526 stored entries:
⎡⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠑⣤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠱⢆⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠱⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⢄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢻⣶⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠛⣤⣀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⢻⣶⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠛⣤⡀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⎦

julia> bandwidth(A)
12

julia> bandwidth(A_shuffled) # Much larger after shuffling
731

julia> res # Gets very close to the original bandwidth
Results of Bandwidth Minimization Algorithm
 * Algorithm: Gibbs–Poole–Stockmeyer
 * Approach: heuristic
 * Minimum Bandwidth: 18
 * Original Bandwidth: 731
 * Matrix Size: 748×748

Notes

Note that the node_finder field must be of the form (A::AbstractMatrix{Bool}) -> Integer (i.e., it must take in an boolean matrix and return an integer). If this is not the case, an ArgumentError is thrown upon construction.

References

• [Geo71]: J. A. George. Computer Implementation of the Finite Element Method. Ph.D. Thesis, Department of Computer Science, Stanford University (1971). https://apps.dtic.mil/sti/tr/pdf/AD0726171.pdf.
• [GPS76]: N. E. Gibbs, W. G. Poole Jr. and P. K. Stockmeyer. An Algorithm for Reducing the Bandwidth and Profile of a Sparse Matrix. SIAM Journal on Numerical Analysis 13, 236–50 (1976). https://doi.org/10.1137/0713023.
• [Lew82]: J. G. Lewis. Implementation of the Gibbs–Poole–Stockmeyer and Gibbs–King Algorithms. ACM Transactions on Mathematical Software 8, 180–89 (1982). https://doi.org/10.1145/355993.355998.

HeuristicSolver <: AbstractSolver <: AbstractAlgorithm

Abstract type for all heuristic matrix bandwidth minimization solvers.

Heuristic methods are those which aim to produce near-optimal solutions in a more performant manner than exact methods. While precise bandwidth minimization is NP-complete, many heuristic algorithms (such as Gibbs–Poole–Stockmeyer) run in polynomial time.

Heuristic algorithms differ from metaheuristic ones in that they do not employ higher-level iterative search frameworks (e.g., stochastic techniques) to survey the global search space and escape local minima; instead, they rely on straightforward deterministic procedures to find good solutions in a single pass.

Supertype Hierarchy

HeuristicSolver <: AbstractSolver <: MatrixBandwidth.AbstractAlgorithm

ReverseCuthillMcKee <: HeuristicSolver <: AbstractSolver <: AbstractAlgorithm

The reverse Cuthill–McKee algorithm is a variant of the Cuthill–McKee algorithm—a heuristic method for minimizing the bandwidth of a structurally symmetric matrix $A$. Cuthill–McKee considers the graph $G(A)$ whose adjacency matrix is $A$ (ignoring weights and self-loops) and performs a breadth-first search of each connected component of $G(A)$, starting from a low-degree node then visiting its neighbors in order of increasing degree. Particularly effective when $A$ is sparse, this heuristic typically produces an ordering which induces a matrix bandwidth either equal to or very close to the true minimum [CM69, pp. 157–58]. The reverse Cuthill–McKee algorithm simply reverses the ordering produced by application of Cuthill–McKee; it was found in [Geo71, pp. 114–15] that although the bandwidth remains the same, this tends to produce a more optimal matrix profile (a measure of how far, on average, nonzero entries are from the diagonal; see also MatrixBandwidth.profile).

As noted above, the reverse Cuthill–McKee algorithm requires structurally symmetric input (that is, $A[i, j]$ must be nonzero if and only if $A[j, i]$ is nonzero for $1 ≤ i, j ≤ n$).

Performance

Given an $n×n$ input matrix $A$, the reverse Cuthill–McKee algorithm runs in $O(n²)$ time.

[CG80] provide a linear-time implementation in the number of nonzero entries of $A$, which is still quadratic when $A$ is dense but often much faster when dealing with sparse matrices. However, this would require that $A$ be stored as a graph or a sparse matrix, which runs counter to our desire to provide a bandwidth minimization API for all AbstractMatrix{<:Number} types, including dense matrices. (In the future, however, we may indeed consider supporting this more performant implementation for sparse matrices.)

Fields

• node_finder::Function: a function that selects a node from some connected component of the input matrix from which to start the breadth-first search. If no custom heuristic is specified, this field defaults to bi_criteria_node_finder, which picks a node "farthest" from the others in the component (not necessarily the lowest-degree node).

Supertype Hierarchy

ReverseCuthillMcKee <: HeuristicSolver <: AbstractSolver <: MatrixBandwidth.AbstractAlgorithm

Examples

In the following examples, MatrixBandwidth.random_banded_matrix is used to generate random matrices with minimum bandwidth close to $k$. In some cases, however, the true minimum bandwidth up to symmetric permutation may be even less than $k$, making it hard to verify whether reverse Cuthill–McKee finds a truly optimal ordering or simply a near-optimal one. Nevertheless, the results are still very good in practice.

Reverse Cuthill–McKee finds a good ordering for a $35×35$ matrix:

julia> using Random

julia> Random.seed!(87);

julia> (n, k) = (35, 3);

julia> A = MatrixBandwidth.random_banded_matrix(n, k);

julia> perm = randperm(n);

julia> A_shuffled = A[perm, perm];

julia> bandwidth(A)
3

julia> bandwidth(A_shuffled) # Much larger after shuffling
30

julia> minimize_bandwidth(A_shuffled, Minimization.ReverseCuthillMcKee())
Results of Bandwidth Minimization Algorithm
 * Algorithm: Reverse Cuthill–McKee
 * Approach: heuristic
 * Minimum Bandwidth: 3
 * Original Bandwidth: 30
 * Matrix Size: 35×35

Reverse Cuthill–McKee finds a good ordering for a $235×235$ matrix with multiple (separate) connected components:

julia> using Random, SparseArrays

julia> Random.seed!(5747);

julia> (max_cc_size, max_band, p, num_ccs) = (60, 9, 0.2, 8);

julia> components = Vector{SparseMatrixCSC{Float64, Int64}}(undef, num_ccs);

julia> for i in 1:num_ccs # Some components may themselves be disconnected
           cc_size = rand(0:max_cc_size);
           cc_band = rand(1:min(max_band, cc_size - 1));
           components[i] = sparse(
               MatrixBandwidth.random_banded_matrix(cc_size, cc_band; p=p)
           );
       end

julia> A = blockdiag(components...); # `A` has least 8 connected components

julia> perm = randperm(sum(map(cc -> size(cc, 1), components)));

julia> A_shuffled = A[perm, perm];

julia> res = minimize_bandwidth(A_shuffled, Minimization.ReverseCuthillMcKee());

julia> A # The original matrix
235×235 SparseMatrixCSC{Float64, Int64} with 445 stored entries:
⎡⢾⣳⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠘⢿⡷⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠏⣥⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠉⢴⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠙⠻⢂⣀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠸⣮⣿⣢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠚⢿⡳⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⢰⣶⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠾⡧⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⣢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⡠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢻⠖⣀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⢏⡱⣄⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢮⣷⣄⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢻⣲⣄⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢻⢖⎦

julia> A_shuffled # A far-from-optimal ordering of `A`
235×235 SparseMatrixCSC{Float64, Int64} with 445 stored entries:
⎡⠑⠄⠀⠀⠀⠀⠀⢀⢀⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠄⠀⡐⠀⠐⠀⠂⠀⠀⠀⠀⠀⠀⡂⠀⢀⠄⠁⠠⠐⠀⎤
⎢⠀⠀⠀⢄⡀⠀⢁⠀⠀⠈⠀⠁⠀⠀⠀⢀⠁⠄⠈⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠈⠀⠂⠠⠀⠀⠀⠀⠀⡀⠐⎥
⎢⠀⠀⠀⠈⠁⢀⠀⠑⠀⢀⠁⢀⠀⠈⠀⠘⠌⠀⢀⠀⠄⠀⠂⡄⠄⠁⠀⠀⠈⠀⠀⠀⠀⠀⠀⠀⣂⠀⠀⠀⎥
⎢⠀⢀⠁⠐⢄⠀⠠⢆⠀⠀⠀⠀⠀⠀⠀⢀⠠⡀⠀⠀⠠⠀⠀⠀⠐⠀⠀⡀⠀⢀⠀⠀⠀⠈⠀⡀⠀⠀⠘⠀⎥
⎢⠀⠐⡀⠀⠀⢀⠀⠀⢀⢔⠈⢀⠀⠀⣐⠀⠀⠀⢀⠀⠀⠀⠀⠀⠐⠀⠄⢠⠀⠀⠀⠀⠀⠀⠀⡀⠀⠈⠣⡀⎥
⎢⠀⠁⠄⠀⠁⢀⠀⠀⠂⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⠀⡀⠀⠀⡐⠐⠀⡀⠀⠀⠂⠀⠀⠀⢀⠀⠀⠄⠀⎥
⎢⠀⠀⠀⠀⡀⠀⠀⠀⠀⠀⠀⠀⠑⠀⠀⠀⠀⠀⠀⠀⡠⠀⡀⠀⠀⠄⠀⠀⠠⠀⠀⠠⠀⠀⠀⠀⡠⠄⠀⠄⎥
⎢⠀⠀⠀⢀⣀⠀⠀⢀⠐⠘⠀⠀⠀⠀⠕⢅⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⢀⠐⡀⠀⠈⠀⠂⠀⠀⢀⠀⠃⠀⠄⎥
⎢⠀⠀⠁⠄⠂⠁⠀⠢⠀⠀⠀⠀⠀⠀⠀⠀⠛⢄⢸⠘⠀⠀⠀⠀⠄⠈⠁⠀⠀⠨⠀⠀⢀⠀⠀⠨⠀⠀⠈⠀⎥
⎢⠀⠄⠂⠀⠀⠐⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⣒⠒⠁⠀⢠⠀⠀⠀⠐⠀⠀⠀⢁⠀⠐⠈⠀⠀⠂⠀⠂⡀⠀⠀⎥
⎢⢀⠠⠀⠀⠀⠁⠀⠂⠀⠀⠀⠂⠀⠊⠀⠀⠀⠀⠀⠒⠑⠀⠀⠀⠀⠂⢀⠁⡂⠀⠀⢀⠀⠀⠀⠀⠁⠂⡊⠂⎥
⎢⢀⠀⢀⠀⠈⠤⠀⠀⠀⠀⠀⠈⠀⠈⠀⠠⠀⠀⠀⠀⠀⠀⠁⢀⠅⠀⢀⠀⠀⠀⢀⠀⡁⠀⠀⠀⠀⠠⠀⠀⎥
⎢⠠⠀⠀⠀⠄⠁⠐⠀⠐⠀⢀⠠⠀⠄⠀⠀⡀⠁⠐⠀⠠⠀⠁⠁⠁⠀⢀⠄⠀⠉⠀⠃⠀⠀⠀⠀⠠⢠⠂⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠠⠀⣁⠐⠀⠀⠀⢀⠐⠁⠀⠀⠀⠄⠐⠀⠐⠀⠔⡑⠌⠀⠀⠀⠀⠢⢉⠀⠀⠄⠀⠀⠄⎥
⎢⠀⠀⡀⠀⠂⠀⠀⢀⠀⠀⠀⠈⠀⠂⠀⠈⡀⡀⠁⠐⠈⠈⠀⠀⡄⠀⠀⠀⠄⠅⠀⡀⠀⠀⠀⠠⠀⠰⠀⠂⎥
⎢⠀⠀⠠⠀⠀⠀⠀⠀⠀⠀⠠⠀⠀⡀⠂⠀⠀⠀⡐⠀⠀⢀⠀⠐⠤⠀⠀⠀⠀⠠⠀⢀⠀⠀⠀⠀⠀⠀⠒⢀⎥
⎢⠈⠈⠀⠂⠀⠀⡀⠀⠀⠀⠀⠀⠀⠀⠈⠀⠀⠐⠀⠀⠀⠀⠁⠈⠀⠀⡌⢂⠀⠀⠀⠀⠀⠄⠈⠀⠀⡀⠀⠀⎥
⎢⠀⠔⠀⠀⠀⠀⠀⠠⠀⠠⠀⢀⠀⠀⠀⢀⡀⡀⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠂⠀⠡⠂⠀⠀⡠⠀⎥
⎢⠁⡀⠀⠀⠈⠘⠀⠀⡀⠀⠀⠀⠀⠎⠤⠀⠀⠀⠈⠠⠡⠀⠀⡀⠀⣂⠀⠁⢀⡀⠀⠀⠀⠠⠀⠀⠰⠆⠌⠀⎥
⎣⠐⠀⢀⠈⠀⠀⠒⠀⠉⠢⠀⠁⠀⠄⠀⠄⠂⠀⠀⠀⠪⠈⠀⠀⠈⠀⠀⠄⠠⠀⠘⢀⠀⠀⠀⠊⠂⠁⠐⢀⎦

julia> A_shuffled[res.ordering, res.ordering] # A near-optimal reordering of `A_shuffled`
235×235 SparseMatrixCSC{Float64, Int64} with 445 stored entries:
⎡⠁⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠁⢄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠁⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠑⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢞⡵⡦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠫⢥⣳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠛⣤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠡⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠫⡦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⢻⣶⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢫⣶⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢪⣶⣄⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢾⣳⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⣄⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠿⣧⡀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠱⣦⎦

julia> bandwidth(A)
9

julia> bandwidth(A_shuffled) # Much larger after shuffling
226

julia> res # Gets very close to the original bandwidth
Results of Bandwidth Minimization Algorithm
 * Algorithm: Reverse Cuthill–McKee
 * Approach: heuristic
 * Minimum Bandwidth: 11
 * Original Bandwidth: 226
 * Matrix Size: 235×235

Notes

Note that the node_finder field must be of the form (A::AbstractMatrix{Bool}) -> Integer (i.e., it must take in an boolean matrix and return an integer). If this is not the case, an ArgumentError is thrown upon construction.

See also the documentation for CuthillMcKee—the original (non-reversed) algorithm. (Indeed, the reverse Cuthill–McKee method of _minimize_bandwidth_impl is merely a wrapper around the Cuthill–McKee method.)

References

• [CG80]: W. M. Chan and A. George. A linear time implementation of the reverse Cuthill–McKee algorithm. BIT Numerical Mathematics 20, 8–14 (1980). https://doi.org/10.1007/BF01933580.
• [CM69]: E. Cuthill and J. McKee. Reducing the bandwidth of sparse symmetric matrices. In: Proceedings of the 24th National Conference of the ACM (Brandon Systems Press, 1969); pp. 157–72. https://doi.org/10.1145/800195.805928.
• [Geo71]: J. A. George. Computer Implementation of the Finite Element Method. Ph.D. Thesis, Department of Computer Science, Stanford University (1971). https://apps.dtic.mil/sti/tr/pdf/AD0726171.pdf.

MatrixBandwidth.Minimization.Metaheuristic

Metaheuristic solvers for matrix bandwidth minimization.

Metaheuristic methods are those which employ higher-level iterative search frameworks such as stochastic techniques or nature-inspired processes to survey the global search space and escape local minima. Unlike heuristic methods—which follow fixed deterministic procedures—metaheuristics adaptively refine candidate solutions over multiple iterations. Although metaheuristic approaches are often slower than heuristic ones (but certainly still faster than exact ones), they shine in complex cases where the latter may get trapped in poor-quality local minima.

No metaheuristic matrix bandwidth minimization algorithms are already implemented, but several (e.g., simulated annealing) are currently under development.

This submodule is part of the MatrixBandwidth.Minimization submodule of the MatrixBandwidth.jl package.

AntColony <: MetaheuristicSolver <: AbstractSolver <: AbstractAlgorithm

[TODO: Write here]

Supertype Hierarchy

AntColony <: MetaheuristicSolver <: AbstractSolver <: MatrixBandwidth.AbstractAlgorithm

[TODO: Write here]

GRASP <: MetaheuristicSolver <: AbstractSolver <: AbstractAlgorithm

[TODO: Write here]

Supertype Hierarchy

GRASP <: MetaheuristicSolver <: AbstractSolver <: MatrixBandwidth.AbstractAlgorithm

[TODO: Write here]

GeneticAlgorithm <: MetaheuristicSolver <: AbstractSolver <: AbstractAlgorithm

[TODO: Write here]

Supertype Hierarchy

GeneticAlgorithm <: MetaheuristicSolver <: AbstractSolver <: MatrixBandwidth.AbstractAlgorithm

[TODO: Write here]

MetaheuristicSolver <: AbstractSolver <: AbstractAlgorithm

Abstract type for all metaheuristic matrix bandwidth minimization solvers.

Metaheuristic methods are those which employ higher-level iterative search frameworks such as stochastic techniques or nature-inspired processes to survey the global search space and escape local minima. Unlike heuristic methods—which follow fixed deterministic procedures—metaheuristics adaptively refine candidate solutions over multiple iterations. Although metaheuristic approaches are often slower than heuristic ones (but certainly still faster than exact ones), they shine in complex cases where the latter may get trapped in poor-quality local minima.

Supertype Hierarchy

MetaheuristicSolver <: AbstractSolver <: MatrixBandwidth.AbstractAlgorithm

PSOHC <: MetaheuristicSolver <: AbstractSolver <: AbstractAlgorithm

[TODO: Write here]

Supertype Hierarchy

PSOHC <: MetaheuristicSolver <: AbstractSolver <: MatrixBandwidth.AbstractAlgorithm

[TODO: Write here]

SimulatedAnnealing <: MetaheuristicSolver <: AbstractSolver <: AbstractAlgorithm

[TODO: Write here]

Supertype Hierarchy

SimulatedAnnealing <: MetaheuristicSolver <: AbstractSolver <: MatrixBandwidth.AbstractAlgorithm

[TODO: Write here]

TabuSearch <: MetaheuristicSolver <: AbstractSolver <: AbstractAlgorithm

[TODO: Write here]

Supertype Hierarchy

TabuSearch <: MetaheuristicSolver <: AbstractSolver <: MatrixBandwidth.AbstractAlgorithm

[TODO: Write here]

MatrixBandwidth.Recognition

Algorithms for matrix bandwidth recognition in Julia.

The bandwidth of an $n×n$ matrix $A$ is the minimum non-negative integer $k ∈ \{0, 1, …, n - 1\}$ such that $A[i, j] = 0$ whenever $|i - j| > k$. The matrix bandwidth recognition problem entails determining whether there exists a permutation matrix $P$ such that the bandwidth of $PAPᵀ$ is at most some fixed non-negative integer (an optimal permutation that fully minimizes the bandwidth of $A$ is not required).

The following matrix bandwidth recognition algorithms are currently available:

• Del Corso–Manzini (DelCorsoManzini)
• Del Corso–Manzini with perimeter search (DelCorsoManziniWithPS)
• Caprara–Salazar-González (CapraraSalazarGonzalez)
• Saxe–Gurari–Sudborough (SaxeGurariSudborough)
• Brute-force search (BruteForceSearch)

This submodule is part of the MatrixBandwidth.jl package.

AbstractDecider <: AbstractAlgorithm

Abstract base type for all matrix bandwidth recognition deciders.

Interface

As per the interface of supertype AbstractAlgorithm, concrete subtypes of AbstractDecider must implement the following methods:

• Base.summary(::T) where {T<:AbstractDecider}: returns a String indicating the name of the decider (e.g., "Caprara–Salazar-González").
• _requires_structural_symmetry(::T) where {T<:AbstractDecider}: returns a Bool indicating whether the decider requires the input matrix to be structurally symmetric.

Supertype Hierarchy

AbstractDecider <: AbstractAlgorithm

Notes

To implement a new matrix bandwidth recognition algorithm, define a new concrete subtype of AbstractDecider then implement a corresponding _has_bandwidth_k_ordering_impl(::AbstractMatrix{Bool}, ::Integer, ::NewDeciderType) method. Do not attempt to directly implement a new has_bandwidth_k_ordering method, as the function contains common preprocessing logic independent of the specific algorithm used.

BruteForceSearch <: AbstractDecider <: AbstractAlgorithm

The simplest method for determining, given some fixed $k ∈ ℕ$, whether a matrix has bandwidth at most $k$ up to symmetric permutation is to iterate over all orderings and compute the bandwidth induced by each.

Since $i₁, i₂, … iₙ$ induces the same bandwidth as $iₙ, iₙ₋₁, … i₁$, we restrict our search to orderings such that $i₁ ≤ iₙ$ (with equality checked just in case $n = 1$).

If a bandwidth-$k$ ordering is found, the algorithm breaks early instead of continuing to check subsequent permutations.

Supertype Hierarchy

BruteForceSearch <: AbstractDecider <: AbstractAlgorithm

Performance

Given an $n×n$ input matrix, this brute-force algorithm runs in $O(n! ⋅ n²)$ time:

• Up to $n!/2$ permutations may be checked (except when $n = 1$, in which case $1! = 1$ permutation is checked). This is, clearly, $O(n!)$.
• For each permutation, the bandwidth function is called on view(A, perm, perm), which takes $O(n²)$ time.
• Therefore, the overall time complexity is $O(n! ⋅ n²)$.

Examples

In many cases, the algorithm iterates over all (if $k$ is smaller than the true minimu bandwidth) or almost all (if $k$ is equally to or only slightly larger than the true minimum) possible permutations—in these cases, it is infeasible to go above $9×9$ or $10×10$ without incurring multiple-hour runtimes. (Even when $k$ is considerably larger than the true minimum, it is unlikely that a bandwidth-$k$ ordering will be found in a reasonable time frame.) Nevertheless, we see that it is quite effective for, say, $8×8$:

julia> using Random, SparseArrays

julia> Random.seed!(314159);

julia> (n, p) = (8, 0.5);

julia> A = sprand(n, n, p);

julia> A = A + A' # Ensure structural symmetry;

julia> (k_false, k_true) = (3, 5);

julia> has_bandwidth_k_ordering(A, k_false, Recognition.BruteForceSearch())
Results of Bandwidth Recognition Algorithm
 * Algorithm: Brute-force search
 * Bandwidth Threshold k: 3
 * Has Bandwidth ≤ k Ordering: false
 * Original Bandwidth: 6
 * Matrix Size: 8×8

julia> has_bandwidth_k_ordering(A, k_true, Recognition.BruteForceSearch())
Results of Bandwidth Recognition Algorithm
 * Algorithm: Brute-force search
 * Bandwidth Threshold k: 5
 * Has Bandwidth ≤ k Ordering: true
 * Original Bandwidth: 6
 * Matrix Size: 8×8

Notes

Brute force is by far the slowest approach to matrix bandwidth minimization and should only be used in very niche cases (like verifying the correctness of other algorithms in unit tests). For $10×10$ matrices, the algorithm already takes several minutes to run for difficult values of $k$ (namely, values below or only slightly above the true minimum) and allocates several gigabytes of memory. Given the $O(n! ⋅ n²)$ time complexity, checking "bandwidth ≤ k" would take over an hour for many $11×11$ matrices.

This holds true even when $k$ is considerably larger than the true minimum bandwidth—as long as it remains below the bandwidth induced by the original ordering, it is unlikely that a bandwidth-$k$ ordering will be found early simply by random chance. Additionally, time complexity will remain on the order of $n! ⋅ n²$ in the average case.

See also MatrixBandwidth.Minimization.Exact.BruteForceSearch for the minimization variant of this algorithm (which simply never breaks early, instead iterating over all permutations up to reversal to ensure that the minimum bandwidth is found).

CapraraSalazarGonzalez <: AbstractDecider <: AbstractAlgorithm

The Caprara–Salazar-González recognition algorithm is a method for determining, given some fixed $k ∈ ℕ$, whether a structurally symmetric matrix $A$ has a bandwidth at most $k$ up to symmetric permutation. The algorithm performs a depth-first search of all partial orderings of the rows and columns of $A$, adding indices one at a time. Partial orderings are pruned not only by ensuring that adjacent pairs of currently placed indices are within $k$ of each other but also by employing a branch-and-bound framework with lower bounds on bandwidth compatibility computed via integer linear programming relaxations. This search is repeated with incrementing values of $k$ until a bandwidth-$k$ ordering is found [CS05], with $k$ initialized to some lower bound on the minimum bandwidth of $A$ up to symmetric permutation.

As noted above, the Caprara–Salazar-González algorithm requires structurally symmetric input (that is, $A[i, j]$ must be nonzero if and only if $A[j, i]$ is nonzero for $1 ≤ i, j ≤ n$).

Supertype Hierarchy

CapraraSalazarGonzalez <: AbstractDecider <: AbstractAlgorithm

Performance

Given an $n×n$ input matrix and threshold bandwidth $k$, the Caprara–Salazar-González algorithm runs in $O(n! ⋅ Tᵢₗₚ(n, n²))$ time, where $Tᵢₗₚ(n, m)$ is the time taken to solve an integer linear programming (ILP) problem with $O(n)$ variables and $O(m)$ constraints:

• We perform a depth-first search of $O(n!)$ partial orderings.
• At each search node, we solve ILP relaxations with $n$ variables and $O(n²)$ constraints (given by the number of nonzero entries in the computed distance matrix), taking $Tᵢₗₚ(n, n²)$ time. (This dominates the $O(n²)$ auxiliary computations needed to set up the ILP.)
• Therefore, the overall time complexity is $O(n! ⋅ Tᵢₗₚ(n, n²))$.

Note that $Tᵢₗₚ(n, n²)$ has worst-case complexity $O(2ⁿ)$, although this ultimately depends on the ILP solver used. (Here, we use the HiGHS solver from the HiGHS.jl package.)

Of course, this is all but an upper bound on the time complexity of Caprara–Salazar-González, achieved only in the most pathological of cases. In practice, efficient pruning techniques and compatibility checks result in approximately exponential growth in time complexity with respect to $n$.

Examples

julia> using Random, SparseArrays

julia> Random.seed!(17);

julia> (n, p) = (10, 0.17);

julia> A = sprand(n, n, p);

julia> A = A + A' # Ensure structural symmetry;

julia> (k_false, k_true) = (3, 6);

julia> has_bandwidth_k_ordering(A, k_false, Recognition.CapraraSalazarGonzalez())
Results of Bandwidth Recognition Algorithm
 * Algorithm: Caprara–Salazar-González
 * Bandwidth Threshold k: 3
 * Has Bandwidth ≤ k Ordering: false
 * Original Bandwidth: 9
 * Matrix Size: 10×10

julia> has_bandwidth_k_ordering(A, k_true, Recognition.CapraraSalazarGonzalez())
Results of Bandwidth Recognition Algorithm
 * Algorithm: Caprara–Salazar-González
 * Bandwidth Threshold k: 6
 * Has Bandwidth ≤ k Ordering: true
 * Original Bandwidth: 9
 * Matrix Size: 10×10

Notes

For readers of the original paper, what we call the Caprara–Salazar-González algorithm here is designated the LAYOUT_LEFT_TO_RIGHT algorithm in [CS05]. The paper also describes a LAYOUT_BOTH_WAYS algorithm that performs a bidirectional search by adding indices to both the left and right ends of the current partial ordering. However, this version is considerably more complex to implement, and we ran into problems enforcing ILP constraints on node pairs added to opposite ends of the ordering. In any case, computational results demonstrate that neither LAYOUT_LEFT_TO_RIGHT nor LAYOUT_BOTH_WAYS is consistently faster, and the paper states that there is no known heuristic for determining which version will be more performant for a given input [CS05, pp. 368–69]. Therefore, we opt to implement only LAYOUT_LEFT_TO_RIGHT as a matter of practicality, although future developers may wish to extend the interface with LAYOUT_BOTH_WAYS as well.

Do also note that this algorithm is not the main LAYOUT_LEFT_TO_RIGHT procedure described in the original paper, which actually never explicitly tackles matrix bandwidth recognition [CS05]. However, the LAYOUT_LEFT_TO_RIGHT algorithm presented therein for bandwidth minimization does repeatedly call a recognition subroutine—this is precisely the logic we implement here. (We do, however, also implement said minimization algorithm in MatrixBandwidth.Minimization.Exact.CapraraSalazarGonzalez.)

A final implementation detail worth noting is that we use HiGHS as our solver; it is one of the fastest open-source solvers available for mixed-integer linear programming.

References

• [CS05]: A. Caprara and J.-J. Salazar-González. Laying Out Sparse Graphs with Provably Minimum Bandwidth. INFORMS Journal on Computing 17, 356–73 (2005). https://doi.org/10.1287/ijoc.1040.0083.

DelCorsoManzini <: AbstractDecider <: AbstractAlgorithm

The Del Corso–Manzini recognition algorithm is a method for determining, given some fixed $k ∈ ℕ$, whether a structurally symmetric matrix $A$ has bandwidth at most $k$ up to symmetric permutation. The algorithm performs a depth-first search of all partial orderings of the rows and columns of $A$, adding indices one at a time. Partial orderings are pruned not only by ensuring that adjacent pairs of currently placed indices are within $k$ of each other but also by tracking the latest positions at which the remaining indices can be placed [DM99].

As noted above, the Del Corso–Manzini algorithm requires structurally symmetric input (that is, $A[i, j]$ must be nonzero if and only if $A[j, i]$ is nonzero for $1 ≤ i, j ≤ n$).

Supertype Hierarchy

DelCorsoManzini <: AbstractDecider <: AbstractAlgorithm

Performance

Given an $n×n$ input matrix and threshold bandwidth $k$, the Del Corso–Manzini algorithm runs in $O(n! ⋅ nk)$ time:

• We perform a depth-first search of $O(n!)$ partial orderings.
• Checking plausibility of each partial ordering takes $O(nk)$ time.
• Therefore, the overall time complexity is $O(n! ⋅ nk)$.

Of course, this is but an upper bound on the time complexity of Del Corso–Manzini, achieved only in the most pathological of cases. In practice, efficient pruning techniques and compatibility checks result in approximately exponential growth in time complexity with respect to $n$.

Based on experimental results, the algorithm is feasible for $n×n$ matrices up to $n ≈ 100$ or so.

Examples

julia> using Random, SparseArrays

julia> Random.seed!(7878);

julia> (n, p) = (40, 0.1);

julia> A = sprand(n, n, p);

julia> A = A + A' # Ensure structural symmetry;

julia> (k_false, k_true) = (13, 26);

julia> has_bandwidth_k_ordering(A, k_false, Recognition.DelCorsoManzini())
Results of Bandwidth Recognition Algorithm
 * Algorithm: Del Corso–Manzini
 * Bandwidth Threshold k: 13
 * Has Bandwidth ≤ k Ordering: false
 * Original Bandwidth: 34
 * Matrix Size: 40×40

julia> has_bandwidth_k_ordering(A, k_true, Recognition.DelCorsoManzini())
Results of Bandwidth Recognition Algorithm
 * Algorithm: Del Corso–Manzini
 * Bandwidth Threshold k: 26
 * Has Bandwidth ≤ k Ordering: true
 * Original Bandwidth: 34
 * Matrix Size: 40×40

Notes

For readers of the original paper, what we call the Del Corso–Manzini recognition algorithm here is essentially a wrapper around the underlying AddNode subroutine in what [DM99, p. 191] term the "MB-ID algorithm" for bandwidth minimization (not mere recognition). MB-ID (which we also implement in MatrixBandwidth.Minimization.Exact.DelCorsoManzini) calls this recognition procedure with incrementing values of $k$ until a bandwidth-$k$ ordering is found, with $k$ initialized to some lower bound on the minimum bandwidth of $A$ up to symmetric permutation.

[DM99, p. 193] also describe an "MB-PS algorithm" for bandwidth minimization, which we implement in MatrixBandwidth.Minimization.Exact.DelCorsoManziniWithPS. Similarly, the underlying recognition subroutine for MB-PS is implemented in DelCorsoManziniWithPS.

References

• [DM99]: G. M. Del Corso and G. Manzini. Finding Exact Solutions to the Bandwidth Minimization Problem. Computing 62, 189–203 (1999). https://doi.org/10.1007/s006070050002.

DelCorsoManziniWithPS{D} <: AbstractDecider <: AbstractAlgorithm

The Del Corso–Manzini recognition algorithm with perimeter search is a method for determining, given some fixed $k ∈ ℕ$, whether a structurally symmetric matrix $A$ has bandwidth at most $k$ up to symmetric permutation. The base Del Corso–Manzini algorithm performs a depth-first search of all partial orderings of the rows and columns of $A$, adding indices one at a time. Partial orderings are pruned not only by ensuring that adjacent pairs of currently placed indices are within $k$ of each other but also by tracking the latest positions at which the remaining indices can be placed [DM99].

The incorporation of perimeter search to this approach entails precomputing a "perimeter" of $d$-permutations of row indices of $A$, where $d$ is a positive integer passed as a parameter to the decider. Each permutation represents a way to select the last $d$ entries of the ordering, and as the construction of the partial ordering progresses, potential endings are pruned to exclude those incompatible with already placed indices. In addition to pruning a potential ending if it contains indices already placed, compatibility is also checked via precomputed time stamps indicating, for each potential ending, a loose lower bound on the earliest position at which any given index can be placed should said ending be selected.

As noted above, the Del Corso–Manzini algorithm with perimeter search requires structurally symmetric input (that is, $A[i, j]$ must be nonzero if and only if $A[j, i]$ is nonzero for $1 ≤ i, j ≤ n$).

Fields

• depth::D<:Union{Nothing,Integer}: the perimeter search depth. If this field is not set (and thus automatically initialized to nothing), a default depth is computed by dcm_ps_optimal_depth as a function of the input matrix every time the decider is passed to has_bandwidth_k_ordering as a function of the input matrix. Otherwise, it must be manually set to a positive integer.

Constructors

• DelCorsoManziniWithPS(): constructs a new DelCorsoManziniWithPS instance with the default perimeter search depth initialized to nothing.
• DelCorsoManziniWithPS(depth::Integer): constructs a new DelCorsoManziniWithPS instance with the specified perimeter search depth. depth must be a positive integer.

Supertype Hierarchy

DelCorsoManziniWithPS <: AbstractDecider <: AbstractAlgorithm

Performance

Given an $n×n$ input matrix, perimeter search depth $d$, and threshold bandwidth $k$, the Del Corso–Manzini algorithm with perimeter search runs in $O(n! ⋅ max(nᵈ, nk))$ time:

• We perform a depth-first search of $O(n!)$ partial orderings.
• Checking plausibility of each partial ordering takes $O(nk)$ time, and checking compatibility with all size-$d$ LPOs takes $O(nᵈ)$ time. Thus, the overall time complexity for each value of $k$ is $O(n! ⋅ (nᵈ + nk))$.
• Therefore, the overall time complexity is $O(n! ⋅ max(nᵈ, nk))$.

Of course, this is but an upper bound on the time complexity of Del Corso–Manzini with perimeter search, achieved only in the most pathological of cases. In practice, efficient pruning techniques and compatibility checks result in approximately exponential growth in time complexity with respect to $n$.

Based on experimental results, the algorithm is feasible for $n×n$ matrices up to $n ≈ 100$ or so.

Examples

Here, Del Corso–Manzini with perimeter search ascertains that A random $30×30$ matrix has a minimum bandwidth greater than $9$. The depth parameter is not explicitly set; instead, some near-optimal value is automatically computed upon the first has_bandwidth_k_ordering function call.

julia> using Random, SparseArrays

julia> Random.seed!(5847);

julia> (n, p) = (30, 0.05);

julia> A = sprand(n, n, p);

julia> A = A + A' # Ensure structural symmetry;

julia> k = 6;

julia> has_bandwidth_k_ordering(A, k, Recognition.DelCorsoManziniWithPS())
Results of Bandwidth Recognition Algorithm
 * Algorithm: Del Corso–Manzini with perimeter search
 * Bandwidth Threshold k: 6
 * Has Bandwidth ≤ k Ordering: false
 * Original Bandwidth: 27
 * Matrix Size: 30×30

Now, Del Corso–Manzini with perimeter search recognizes that a random $35×35$ matrix has a minimum bandwidth at most $8$. In this case, we explicitly set the depth parameter to $4$ beforehand instead of relying on Recognition.dcm_ps_optimal_depth.

julia> using Random, SparseArrays

julia> Random.seed!(23552);

julia> (n, p, depth) = (35, 0.02, 4);

julia> A = sprand(n, n, p);

julia> A = A + A' # Ensure structural symmetry;

julia> k = 8;

julia> has_bandwidth_k_ordering(A, k, Recognition.DelCorsoManziniWithPS(depth))
Results of Bandwidth Recognition Algorithm
 * Algorithm: Del Corso–Manzini with perimeter search
 * Bandwidth Threshold k: 8
 * Has Bandwidth ≤ k Ordering: true
 * Original Bandwidth: 32
 * Matrix Size: 35×35

Notes

For readers of the original paper, what we call the Del Corso–Manzini recognition algorithm with perimeter search here is essentially a wrapper around the underlying AddNode1 and Prune subroutines in what [DM99, p. 193] term the "MB-PS algorithm" for bandwidth minimization (not mere recognition). MB-PS (which we also implement in MatrixBandwidth.Minimization.Exact.DelCorsoManziniWithPS) calls this recognition procedure with incrementing values of $k$ until a bandwidth-$k$ ordering is found, with $k$ initialized to some lower bound on the minimum bandwidth of $A$ up to symmetric permutation.

[DM99, p. 191] also describe an "MB-ID algorithm" for bandwidth minimization, which we implement in MatrixBandwidth.Minimization.Exact.DelCorsoManzini. Similarly, the underlying recognition subroutine for MB-ID is implemented in DelCorsoManzini.

References

• [DM99]: G. M. Del Corso and G. Manzini. Finding Exact Solutions to the Bandwidth Minimization Problem. Computing 62, 189–203 (1999). https://doi.org/10.1007/s006070050002.

RecognitionResult{A,M,O} <: AbstractResult

Output struct for matrix bandwidth recognition results.

Fields

• algorithm::A<:AbstractDecider: the decider used to test the bandwidth.
• matrix::M<:AbstractMatrix{<:Number}: the original matrix whose bandwidth is tested.
• ordering::O<:Union{Nothing,Vector{Int}}: an ordering of the rows and columns of matrix inducing a bandwidth at most k, if such an ordering exists; otherwise, nothing.
• k::Integer: the threshold bandwidth against which to test.
• has_ordering::Bool: whether the matrix has an ordering inducing a bandwidth at most k. (This is true if and only if ordering is not nothing.)

Constructors

• RecognitionResult(decider, matrix, ordering, k): constructs a new RecognitionResult instance with the given fields. The has_ordering field is automatically determined based on whether ordering is nothing or a Vector{Int}.

Supertype Hierarchy

RecognitionResult <: AbstractResult

SaxeGurariSudborough <: AbstractDecider <: AbstractAlgorithm

The Saxe–Gurari–Sudborough recognition algorithm is a method for determining, given some fixed $k ∈ ℕ$, whether a structurally symmetric matrix $A$ has bandwidth at most $k$ up to symmetric permutation. The algorithm employs dynamic programming to search over equivalence classes of partial orderings, where two partial orderings of length $l$ are equivalent if they share the same active region. (The active region of a partial ordering is defined as the sequence of the last $min(k, l)$ vertices in the ordering taken together with all dangling edges—edges with one endpoint in the ordering and the other endpoint not yet in the ordering.) It extends these partial layouts one vertex at a time in a breadth-first manner, pruning implausible classes that violate bandwidth-$k$ constraints such as degree bounds on active vertices and excessive numbers of dangling edges [GS84].

As noted above, the Saxe–Gurari–Sudborough algorithm requires structurally symmetric input (that is, $A[i, j]$ must be nonzero if and only if $A[j, i]$ is nonzero for $1 ≤ i, j ≤ n$).

Supertype Hierarchy

SaxeGurariSudborough <: AbstractDecider <: AbstractAlgorithm

Performance

Given an $n×n$ input matrix and threshold bandwidth $k$, the Saxe–Gurari–Sudborough algorithm runs in $O(nᵏ)$ time [GS84, p. 531]. This is an improvement upon the original $O(nᵏ⁺¹)$ Saxe algorithm [Sax80, p. 363]. (Of course, when $k < 3$, then the initial $O(n³)$ bandwidth lower bound computation performed in all has_bandwidth_k_ordering calls dominates the overall complexity, although the constant scaling factor of that subroutine is generally much smaller than that of the algorithm proper).

Whereas most bandwidth recognition algorithms are technically factorial-time (with respect to $n$) in the worst case but practically always approximate exponential time complexity in real life (see: DelCorsoManzini), the $O(nᵏ)$ upper bound on Saxe–Gurari–Sudborough is typically a good representation of actual performance in most cases. Indeed, these other types of algorithms tend to outperform Saxe–Gurari–Sudborough for larger $k$, given that their aggressive pruning strategies keep their effective search space very small in practice.

Examples

julia> using Random, SparseArrays

julia> Random.seed!(274);

julia> (n, p) = (20, 0.08);

julia> A = sprand(n, n, p);

julia> A = A + A' # Ensure structural symmetry;

julia> (k_false, k_true) = (3, 5);

julia> has_bandwidth_k_ordering(A, k_false, Recognition.SaxeGurariSudborough())
Results of Bandwidth Recognition Algorithm
 * Algorithm: Saxe–Gurari–Sudborough
 * Bandwidth Threshold k: 3
 * Has Bandwidth ≤ k Ordering: false
 * Original Bandwidth: 12
 * Matrix Size: 20×20

julia> has_bandwidth_k_ordering(A, k_true, Recognition.SaxeGurariSudborough())
Results of Bandwidth Recognition Algorithm
 * Algorithm: Saxe–Gurari–Sudborough
 * Bandwidth Threshold k: 5
 * Has Bandwidth ≤ k Ordering: true
 * Original Bandwidth: 12
 * Matrix Size: 20×20

Notes

This general family of bandwidth recognition algorithms was conceived as a response to a question posed by [GGJK78, p. 494]: is the "bandwidth ≤ k?" problem NP-complete for arbitrary $k$? [Sax80] answered this question in the negative by providing a $O(nᵏ⁺¹)$ algorithm, constructively proving that the problem is class P. Later, [GS84] improved upon this algorithm by reducing time complexity to $O(nᵏ)$. Whereas the original Saxe algorithm considers extensions of partial orderings with any remaining unplaced vertex (of which there are $O(n)$ at any point in the breadth-first search), the Gurari–Sudborough refinement only considers extensions with vertices reachable by paths beginning with a dangling edge that never again traverse a dangling edge [GS84, pp. 535–36].

References

• [GGJK78]: M. R. Garey, R. L. Graham, D. S. Johnson and D. E. Knuth. Complexity Results for Bandwidth Minimization. SIAM Journal on Applied Mathematics 34, 477–95 (1978). https://doi.org/10.1137/0134037.
• [GS84]: E. M. Gurari and I. H. Sudborough. Improved dynamic programming algorithms for bandwidth minimization and the MinCut Linear Arrangement problem. Journal of Algorithms 5, 531–46 (1984). https://doi.org/10.1016/0196-6774(84)90006-3.
• [Sax80]: J. B. Saxe. Dynamic-Programming Algorithms for Recognizing Small-Bandwidth Graphs in Polynomial Time. SIAM Journal on Algebraic and Discrete Methods 1, 363–69 (1980). https://doi.org/10.1137/0601042.

has_bandwidth_k_ordering(A, k, decider=CapraraSalazarGonzalez()) -> RecognitionResult

Determine whether A has bandwidth at most k using the algorithm defined by decider.

The bandwidth of an $n×n$ matrix $A$ is the minimum non-negative integer $k ∈ \{0, 1, …, n - 1\}$ such that $A[i, j] = 0$ whenever $|i - j| > k$.

Given some fixed non-negative integer k, this function determines (with 100% certainty) whether there exists some ordering $π$ of the rows and columns of $A$ such that the bandwidth of $PAPᵀ$ is at most k, where $P$ is the permutation matrix corresponding to $π$. This is known to be decidable in $O(nᵏ)$ time, although some deciders (e.g., CapraraSalazarGonzalez) run in exponential time instead to produce even quicker runtimes in practice.

If $k ≥ n - 1$, then this function immediately answers in the affirmative, since the maximum possible bandwidth of an $n×n$ matrix is $n - 1$. After this initial check, a preliminary lower bound on the bandwidth is computed in $O(n³)$ time using results from Caprara and Salazar-González (2005). If this lower bound is greater than $k$`

Arguments

• A::AbstractMatrix{<:Number}: the (square) matrix whose bandwidth is tested.
• k::Integer: the threshold bandwidth against which to test.
• decider::AbstractDecider: the matrix bandwidth recognition algorithm to use; defaults to CapraraSalazarGonzalez. (See the Recognition module documentation for a full list of supported deciders.)

Returns

• ::RecognitionResult: a struct containing the algorithm used, the original matrix A, the identified ordering of the rows and columns (if one exists), the threshold bandwidth k, and a boolean indicating whether the ordering exists.

Examples

Multiple algorithms to decide whether a given matrix has bandwidth at most k are available. Naturally, they will always agree, but the final orderings produced (in the case of an affirmative) may differ:

julia> using Random, SparseArrays

julia> Random.seed!(52);

julia> (n, p) = (8, 0.2);

julia> A = sprand(Bool, n, n, p);

julia> A = A .|| A' # Ensure structural symmetry
8×8 SparseMatrixCSC{Bool, Int64} with 22 stored entries:
 1  ⋅  ⋅  1  ⋅  1  1  ⋅
 ⋅  ⋅  ⋅  ⋅  1  1  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1
 1  ⋅  ⋅  ⋅  ⋅  ⋅  1  ⋅
 ⋅  1  ⋅  ⋅  ⋅  ⋅  1  1
 1  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 1  ⋅  ⋅  1  1  ⋅  ⋅  1
 ⋅  ⋅  1  ⋅  1  ⋅  1  1

julia> k = 3;

julia> res_csg = has_bandwidth_k_ordering(A, k, Recognition.CapraraSalazarGonzalez())
Results of Bandwidth Recognition Algorithm
 * Algorithm: Caprara–Salazar-González
 * Bandwidth Threshold k: 3
 * Has Bandwidth ≤ k Ordering: true
 * Original Bandwidth: 6
 * Matrix Size: 8×8

julia> A[res_csg.ordering, res_csg.ordering]
8×8 SparseMatrixCSC{Bool, Int64} with 22 stored entries:
 ⋅  ⋅  1  1  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  1  ⋅  1  ⋅  ⋅  ⋅
 1  1  1  1  ⋅  ⋅  ⋅  ⋅
 1  ⋅  1  ⋅  ⋅  1  1  ⋅
 ⋅  1  ⋅  ⋅  ⋅  1  ⋅  ⋅
 ⋅  ⋅  ⋅  1  1  ⋅  1  ⋅
 ⋅  ⋅  ⋅  1  ⋅  1  1  1
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1  ⋅

julia> res_sgs = has_bandwidth_k_ordering(A, k, Recognition.SaxeGurariSudborough())
Results of Bandwidth Recognition Algorithm
 * Algorithm: Saxe–Gurari–Sudborough
 * Bandwidth Threshold k: 3
 * Has Bandwidth ≤ k Ordering: true
 * Original Bandwidth: 6
 * Matrix Size: 8×8

julia> A[res_sgs.ordering, res_sgs.ordering]
8×8 SparseMatrixCSC{Bool, Int64} with 22 stored entries:
 ⋅  1  ⋅  1  ⋅  ⋅  ⋅  ⋅
 1  ⋅  1  ⋅  1  ⋅  ⋅  ⋅
 ⋅  1  1  ⋅  1  1  ⋅  ⋅
 1  ⋅  ⋅  ⋅  ⋅  ⋅  1  ⋅
 ⋅  1  1  ⋅  ⋅  ⋅  1  1
 ⋅  ⋅  1  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  1  1  ⋅  1  1
 ⋅  ⋅  ⋅  ⋅  1  ⋅  1  ⋅

If no decider is specified, then the Del Corso–Manzini algorithm is used by default:

julia> using Random, SparseArrays

julia> Random.seed!(174);

julia> (n, p, k) = (20, 0.1, 4);

julia> A = sprand(n, n, p);

julia> A = A .+ A' # Ensure structural symmetry;

julia> has_bandwidth_k_ordering(A, k)
Results of Bandwidth Recognition Algorithm
 * Algorithm: Del Corso–Manzini
 * Bandwidth Threshold k: 4
 * Has Bandwidth ≤ k Ordering: false
 * Original Bandwidth: 15
 * Matrix Size: 20×20

Notes

To implement a new matrix bandwidth recognition algorithm, define a new concrete subtype of AbstractDecider then implement a corresponding _has_bandwidth_k_ordering_impl(::AbstractMatrix{Bool}, ::Integer, ::NewDeciderType) method. Do not attempt to directly implement a new has_bandwidth_k_ordering method, as the function contains common preprocessing logic independent of the specific algorithm used.

Note also that some texts define matrix bandwidth to be the minimum non-negative integer $k$ such that $A[i, j] = 0$ whenever $|i - j| ≥ k$ instead, particularly in more mathematically-minded communities. Effectively, this definition treats diagonal matrices as bandwidth $1$, tridiagonal matrices as bandwidth $2$, and so on. Our definition, on the other hand, is more common in computer science contexts, treating diagonal matrices as bandwidth $0$ and tridiagonal matrices as bandwidth $1$. (Both definitions, however, agree that the bandwidth of an empty matrix is simply $0$.)