MatrixBandwidth.jl – Public API

MatrixBandwidth

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

Reordering the rows and columns of a matrix to reduce its bandwidth has many practical applications in engineering and scientific computing. It is a common preprocessing step used to improve performance when solving linear systems, approximating partial differential equations, optimizing circuit layout, and more [Maf14, p. 184].

Recall that 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$. Equivalently, $A$ has bandwidth at most $k$ if all entries above the $k$ᵗʰ superdiagonal and below the $k$ᵗʰ subdiagonal are zero, and $A$ has bandwidth at least $k$ if there exists any nonzero entry in the $k$ᵗʰ superdiagonal or subdiagonal.

The matrix bandwidth minimization problem involves finding a permutation matrix $P$ such that the bandwidth of $PAPᵀ$ is minimized; this is known to be NP-complete. Several heuristic algorithms (such as Gibbs–Poole–Stockmeyer) run in polynomial time while still producing near-optimal orderings in practice, but exact methods (like Caprara–Salazar-González) are at least exponential in time complexity and thus are only feasible for relatively small matrices.

On the other hand, 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 $k ∈ ℕ$—an optimal permutation that fully minimizes the bandwidth of $A$ is not required. Unlike the NP-hard minimization problem, this is decidable in $O(nᵏ)$ time.

The following algorithms are currently supported:

• Minimization
  • Exact
    • Caprara–Salazar-González algorithm (Minimization.CapraraSalazarGonzalez)
    • Del Corso–Manzini algorithm (Minimization.DelCorsoManzini)
    • Del Corso–Manzini algorithm with perimeter search (Minimization.DelCorsoManziniWithPS)
    • Saxe–Gurari–Sudborough algorithm (Minimization.SaxeGurariSudborough)
    • Brute-force search (Minimization.BruteForceSearch)
  • Heuristic
    • Gibbs–Poole–Stockmeyer algorithm (Minimization.GibbsPooleStockmeyer)
    • Cuthill–McKee algorithm (Minimization.CuthillMcKee)
    • Reverse Cuthill–McKee algorithm (Minimization.ReverseCuthillMcKee)
  • Metaheuristic
    • Greedy randomized adaptive search procedure (GRASP) (Minimization.GRASP)
    • Simulated annealing (Minimization.SimulatedAnnealing)
    • Genetic algorithm (Minimization.GeneticAlgorithm)
• Recognition
  • Caprara–Salazar-González algorithm (Recognition.CapraraSalazarGonzalez)
  • Del Corso–Manzini algorithm (Recognition.DelCorsoManzini)
  • Del Corso–Manzini algorithm with perimeter search (Recognition.DelCorsoManziniWithPS)
  • Saxe–Gurari–Sudborough algorithm (Recognition.SaxeGurariSudborough)
  • Brute-force search (Recognition.BruteForceSearch)

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.

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$. Equivalently, $A$ has bandwidth at most $k$ if all entries above the $k$ᵗʰ superdiagonal and below the $k$ᵗʰ subdiagonal are zero, and $A$ has bandwidth at least $k$ if there exists any nonzero entry in the $k$ᵗʰ superdiagonal or subdiagonal.

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 $A$, 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).(random_banded_matrix(8, 2))
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).

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$. Equivalently, $A$ has bandwidth at most $k$ if all entries above the $k$ᵗʰ superdiagonal and below the $k$ᵗʰ subdiagonal are zero, and $A$ has bandwidth at least $k$ if there exists any nonzero entry in the $k$ᵗʰ superdiagonal or subdiagonal.

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

The function correctly computes a bound less than (or equal to) the true minimum bandwidth of a matrix up to symmetric permutation:

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$.)

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 Julia's SparseArrays 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 $A$, 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

random_banded_matrix(n, k; p=0.5, rng=default_rng()) -> Matrix{Float64}

Generate a random n×n structurally symmetric k-banded matrix with band density ≈ p.

By definition of structural symmetry, the $(i, j)$-th entry of the matrix is nonzero if and only if the $(j, i)$-th entry is nonzero as well. All entries from this matrix are from the interval [0, 1]. Entries up to the kᵗʰ superdiagonal and down to the kᵗʰ subdiagonal are nonzero with probability p.

It is also guaranteed that each of these bands (besides the main diagonal) has at least one nonzero entry (even when p is very small), thus ensuring that the matrix has bandwidth precisely k before any reordering. (There may, however, still exist a symmetric permutation inducing a minimum bandwidth less than k, especially for small values of p.)

Arguments

• n::Integer: the order of the matrix to generate. Must be positive.
• k::Integer: the desired matrix bandwidth. Must satisfy 0 ≤ k < n.

Keyword Arguments

• p::Real=0.5: the band density. Must satisfy 0 < p ≤ 1. Defaults to 0.5.
• rng::AbstractRNG=Random.default_rng(): the random number generator to use. Defaults to Random.default_rng().

Returns

• ::Matrix{Float64}: a random n×n matrix with bandwidth exactly k and sparse bands with density p.

Examples

Generate a $6×6$ matrix with bandwidth $1$ and the maximum number of nonzero entries:

julia> using Random

julia> A = random_banded_matrix(6, 1; p=1, rng=MersenneTwister(1228))
6×6 Matrix{Float64}:
 0.310239  0.346413  0.0       0.0        0.0       0.0
 0.509981  0.917073  0.390771  0.0        0.0       0.0
 0.0       0.760045  0.808396  0.0195686  0.0       0.0
 0.0       0.0       0.222338  0.853164   0.806888  0.0
 0.0       0.0       0.0       0.421603   0.132165  0.805813
 0.0       0.0       0.0       0.0        0.305339  0.0799183

julia> bandwidth(A)
1

Generate a $7×7$ matrix with bandwidth $3$ and band density $0.3$:

julia> using Random

julia> A = random_banded_matrix(7, 3; p=0.3, rng=MersenneTwister(0402))
7×7 Matrix{Float64}:
 0.0       0.132699  0.0       0.0       0.0  0.0       0.0
 0.869352  0.0       0.324319  0.926496  0.0  0.0       0.0
 0.0       0.891878  0.0       0.658102  0.0  0.0       0.0
 0.0       0.88859   0.399559  0.0       0.0  0.284285  0.703377
 0.0       0.0       0.0       0.0       0.0  0.0       0.0
 0.0       0.0       0.0       0.489594  0.0  0.0       0.393573
 0.0       0.0       0.0       0.412412  0.0  0.47063   0.0

julia> bandwidth(A)
3

Generate an $8×8$ diagonal (bandwidth $0$) matrix with default band density ($0.5$):

julia> using Random

julia> A = random_banded_matrix(8, 0; rng=MersenneTwister(0102))
8×8 Matrix{Float64}:
 0.0  0.0        0.0       0.0       0.0  0.0      0.0  0.0
 0.0  0.0762399  0.0       0.0       0.0  0.0      0.0  0.0
 0.0  0.0        0.373113  0.0       0.0  0.0      0.0  0.0
 0.0  0.0        0.0       0.726309  0.0  0.0      0.0  0.0
 0.0  0.0        0.0       0.0       0.0  0.0      0.0  0.0
 0.0  0.0        0.0       0.0       0.0  0.41974  0.0  0.0
 0.0  0.0        0.0       0.0       0.0  0.0      0.0  0.0
 0.0  0.0        0.0       0.0       0.0  0.0      0.0  0.293132

julia> bandwidth(A)
0

Notes

Users of the MatrixBandwidth package may find this function useful when generating random test data for whatever frameworks, algorithms, etc. they are implementing.

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$. Equivalently, $A$ has bandwidth at most $k$ if all entries above the $k$ᵗʰ superdiagonal and below the $k$ᵗʰ subdiagonal are zero, and $A$ has bandwidth at least $k$ if there exists any nonzero entry in the $k$ᵗʰ superdiagonal or subdiagonal.

The matrix bandwidth minimization problem involves finding a permutation matrix $P$ such that the bandwidth of $PAPᵀ$ is minimized; this is known to be NP-complete. Several heuristic algorithms (such as Gibbs–Poole–Stockmeyer) run in polynomial time while still producing near-optimal orderings in practice, but exact methods (like Caprara–Salazar-González) are at least exponential in time complexity and thus are only feasible for relatively small matrices.

The following algorithms are currently supported:

• Exact
  • Caprara–Salazar-González algorithm (CapraraSalazarGonzalez)
  • Del Corso–Manzini algorithm (DelCorsoManzini)
  • Del Corso–Manzini algorithm with perimeter search (DelCorsoManziniWithPS)
  • Saxe–Gurari–Sudborough algorithm (SaxeGurariSudborough)
  • Brute-force search (BruteForceSearch)
• Heuristic
  • Gibbs–Poole–Stockmeyer algorithm (GibbsPooleStockmeyer)
  • Cuthill–McKee algorithm (CuthillMcKee)
  • Reverse Cuthill–McKee algorithm (ReverseCuthillMcKee)
• Metaheuristic
  • Greedy randomized adaptive search procedure (GRASP) (GRASP)
  • Simulated annealing (SimulatedAnnealing)
  • Genetic algorithm (GeneticAlgorithm)

This submodule is part of the MatrixBandwidth.jl package.

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$. Equivalently, $A$ has bandwidth at most $k$ if all entries above the $k$ᵗʰ superdiagonal and below the $k$ᵗʰ subdiagonal are zero, and $A$ has bandwidth at least $k$ if there exists any nonzero entry in the $k$ᵗʰ superdiagonal or subdiagonal.

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

[TODO: Add here once more solvers are implemented. For now, refer to the Examples sections of the GibbsPooleStockmeyer, CuthillMcKee, and ReverseCuthillMcKee docstrings.]

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$.)

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 algorithms are currently supported:

• Caprara–Salazar-González algorithm (CapraraSalazarGonzalez)
• Del Corso–Manzini algorithm (DelCorsoManzini)
• Del Corso–Manzini algorithm with perimeter search (DelCorsoManziniWithPS)
• Saxe–Gurari–Sudborough algorithm (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(Bool, 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 bidirectional depth-first search of all partial orderings of the rows and columns of $A$, adding indices one at a time to both the left and right ends. 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 bandwidtth 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

[TODO: Write here]

Examples

[TODO: Write here]

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 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 $A$, 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 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 = 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.

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 $A$ 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 = 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 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 = 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.

SaxeGurariSudborough <: ExactSolver <: AbstractSolver <: AbstractAlgorithm

[TODO: Write here]

Supertype Hierarchy

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

[TODO: Write here]

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 algorithms are currently supported:

• Gibbs–Poole–Stockmeyer algorithm (GibbsPooleStockmeyer)
• Cuthill–McKee algorithm (CuthillMcKee)
• Reverse Cuthill–McKee algorithm (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_selector::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 pseudo_peripheral_node, 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 _bool_minimal_band_ordering 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 = 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: 5
 * 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(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_selector 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.

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" [GPS76, p. 241] (treating the input matrix $A$ as an undirected graph), 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_selector::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 pseudo_peripheral_node, 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 $A$, 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 = 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(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_selector 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.

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_selector::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 pseudo_peripheral_node, 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 = 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(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: 9
 * Original Bandwidth: 226
 * Matrix Size: 235×235

Notes

Note that the node_selector 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 _bool_minimal_band_ordering is merely a wrapper around the Cuthill–McKee method.)

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.

The following metaheuristic algorithms are currently supported:

• Greedy randomized adaptive search procedure (GRASP) (GRASP)
• Simulated annealing (SimulatedAnnealing)
• Genetic algorithm (GeneticAlgorithm)

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

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]

SimulatedAnnealing <: MetaheuristicSolver <: AbstractSolver <: AbstractAlgorithm

[TODO: Write here]

Supertype Hierarchy

SimulatedAnnealing <: 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$. Equivalently, $A$ has bandwidth at most $k$ if all entries above the $k$ᵗʰ superdiagonal and below the $k$ᵗʰ subdiagonal are zero, and $A$ has bandwidth at least $k$ if there exists any nonzero entry in the $k$ᵗʰ superdiagonal or subdiagonal.

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 $k ∈ ℕ$—an optimal permutation that fully minimizes the bandwidth of $A$ is not required. Unlike the NP-hard minimization problem, this is decidable in $O(nᵏ)$ time.

The following algorithms are currently supported:

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

This submodule is part of the MatrixBandwidth.jl package.

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 $A$, 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(Bool, 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 bidirectional depth-first search of all partial orderings of the rows and columns of $A$, adding indices one at a time to both the left and right ends. 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 bandwidtth 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

[TODO: Write here]

Examples

[TODO: Write here]

Notes

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

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 $A$ 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

We demonstrate both an affirmative and a negative result for the Del Corso–Manzini recognition algorithm on a random $40×40$ matrix:

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] terms 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 describes 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.

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 $A$, 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 explitily 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] terms 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 describes 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.

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

[TODO: Write here]

Supertype Hierarchy

SaxeGurariSudborough <: AbstractDecider <: AbstractAlgorithm

[TODO: Write here]

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$. Equivalently, $A$ has bandwidth at most $k$ if all entries above the $k$ᵗʰ superdiagonal and below the $k$ᵗʰ subdiagonal are zero, and $A$ has bandwidth at least $k$ if there exists any nonzero entry in the $k$ᵗʰ superdiagonal or subdiagonal.

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

[TODO: Add here once more deciders are implemented. For now, refer to the Examples sections of the DelCorsoManzini and DelCorsoManziniWithPS docstrings.]

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$.)