SDiagonalizability.jl – Private API

AbstractSDiagResult

Abstract base type for all $S$-diagonalizability and $S$-bandwidth problem results.

Interface

Concrete subtypes of AbstractSDiagResult must implement parametric types

• A<:Union{AbstractGraph,AbstractMatrix{<:Integer}};
• B<:Tuple{Vararg{Integer}}; and
• C<:Union{Nothing,Eigen},

alongside the following fields:

• network::A: the network whose S-bandwidth is investigated;
• S::B: the set from whose entries we are allowed to construct eigenvectors;
• diagonalization::C: an S-diagonalization of the matrix representation of the network, if it satisfies the specified S-bandwidth constraints; otherwise, nothing.

ArbitraryGraphLaplacian <: ClassifiedLaplacian

[TODO: Write here]

Supertype Hierarchy

ArbitraryGraphLaplacian <: ClassifiedLaplacian

ClassifiedLaplacian

An abstract type representing a classified Laplacian matrix of an undirected graph.

Interface

Concrete subtypes of _TypedLaplacian must implement the following fields:

• matrix::Matrix{Int}: the Laplacian matrix of the graph.

CompleteGraphLaplacian <: ClassifiedLaplacian

[TODO: Write here]

Supertype Hierarchy

CompleteGraphLaplacian <: ClassifiedLaplacian

EfficiencyWarning(msg)

[TODO: Write here]

EmptyGraphLaplacian <: ClassifiedLaplacian

[TODO: Write here]

Supertype Hierarchy

EmptyGraphLaplacian <: ClassifiedLaplacian

KOrthogonality

An abstract type representing the property of k-orthogonality of a collection of vectors.

Recall that an (ordered) collection of vectors $v₁, v₂, ..., vₙ$ is said to be $k$-orthogonal if we have the inner product $⟨`vᵢ, vⱼ⟩ = 0$ whenever $|i - j| ≥ k$ (i.e., if every pair of vectors at least $k$ indices apart is orthogonal). This is equivalent to the vectors' Gram matrix having bandwidth at most k, where we define the bandwidth of a matrix $A$ to be the minimum integer $k ∈ \{1, 2, …, n\}$ such that $Aᵢⱼ = 0$ whenever $|i - j| ≥ k$ [JP25, p. 313]. (Note that many texts instead define matrix bandwidth using zero-based indexing—that is, with the condition $|i - j| > k$ [Maf14, p. 186].)

This type is used as a template for concretely defined properties corresponding to specific values of $k$. In the context of the overarching $S$-bandwidth algorithm, we perform a different depth-first search for each family of values of $k$ on our "tree" of $S$-eigenvectors to determine whether there exists a $k$-orthogonal collection of them.

Interface

Concrete subtypes of KOrthogonality must implement the following fields:

• k::Int: the $k$-orthogonality parameter. Must be a positive integer.

References

• [JP25]: N. Johnston and S. Plosker. Laplacian {−1,0,1}- and {−1,1}-diagonalizable graphs. Linear Algebra and its Applications 704, 309–39 (2025). https://doi.org/10.1016/j.laa.2024.10.016.
• [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.

NotImplementedError{Nothing}(f, subtype, abstracttype)
NotImplementedError{Symbol}(f, arg, subtype, abstracttype)

An exception indicating that a function lacks dispatch to handle a specific argument type.

Semantically, this differs from MethodError in that it connotes a developer-side failure to implement a method rather than erroneous user input. Throughout this package, it is often used to warn when an existing function with multiple dispatch on some abstract type is called on a newly created subtype for which no method has been defined.

Fields

• f::Function: the function called.
• arg::Symbol: the name of the argument with the unsupported type, if the function has multiple arguments. If the function has only one argument, this field should be set to nothing.
• subtype::Type: the type of the argument. May be the actual concrete type or some intermediate supertype. (For instance, if the relevant input has concrete type A with hierarchy A <: B <: C and the abstracttype field is C, then both A and B are perfectly valid choices for subtype.)
• abstracttype::Type: the abstract type under which the argument is meant to fall.

Supertype Hierarchy

NotImplementedError <: Exception

Constructors

• NotImplementedError(::Function, ::Type, ::Type): constructs a new NotImplementedError instance for a single-argument function. Throws an error if the second type is not abstract or the first type is not a subtype of the second.
• NotImplementedError(::Function, ::Symbol, ::Type, ::Type): constructs a new NotImplementedError instance for a multi-argument function. Throws an error if the second type is not abstract or the first type is not a subtype of the second.

NullGraphLaplacian <: ClassifiedLaplacian

[TODO: Write here]

Supertype Hierarchy

NullGraphLaplacian <: ClassifiedLaplacian

Orthogonality <: KOrthogonality

The property of pairwise orthogonality for a collection of vectors.

Recall that a collection of vectors $v₁, v₂, ..., vₙ$ is said to be pairwise orthogonal if we have the inner product $⟨vᵢ, vⱼ⟩ = 0$ whenever $|i - j| ≥ 1$ (i.e., if every pair of vectors is orthogonal). This is equivalent to the vectors' Gram matrix being diagonal.

Fields

• k::Int: the $k$-orthogonality parameter; always necessarily $1$.

Supertype Hierarchy

Orthogonality <: KOrthogonality

Constructors

• Orthogonality(): constructs a new Orthogonality object with k = 1.

QuasiOrthogonality <: KOrthogonality

The property of quasi-orthogonality for a collection of vectors.

Recall that an (ordered) collection of vectors $v₁, v₂, ..., vₙ$ is said to be quasi-orthogonal if we have the inner product $⟨vᵢ, vⱼ⟩ = 0$ whenever $|i - j| ≥ 2$ (i.e., if every pair of vectors at least $2$ indices apart is orthogonal). This is equivalent to the vectors' Gram matrix being tridiagonal [JP25, p. 313].

Fields

• k::Int: the $k$-orthogonality parameter; always necessarily $2$.

Supertype Hierarchy

QuasiOrthogonality <: KOrthogonality

Constructors

• QuasiOrthogonality(): constructs a new QuasiOrthogonality object with k = 2.

References

• [JP25]: N. Johnston and S. Plosker. Laplacian {−1,0,1}- and {−1,1}-diagonalizable graphs. Linear Algebra and its Applications 704, 309–39 (2025). https://doi.org/10.1016/j.laa.2024.10.016.

SBandMinimizationResult{A,B,C,D} <: AbstractSDiagResult

[TODO: Write here]

Supertype Hierarchy

SBandMinimizationResult <: AbstractSDiagResult

SBandRecognitionResult{A,B,C} <: AbstractSDiagResult

[TODO: Write here]

Supertype Hierarchy

SBandRecognitionResult <: AbstractSDiagResult

SDiagonalizabilityResult{A,B,C} <: AbstractSDiagResult

[TODO: Write here]

Supertype Hierarchy

SDiagonalizabilityResult <: AbstractSDiagResult

SSpectra{A,B,C,D,E}

[TODO: Write here]

struct SpectrumIntegralResult

Data on whether a matrix is spectrum integral (i.e., whether its eigenvalues are integers).

This struct also contains a map from each eigenvalue to its multiplicity, provided that the eigenvalues are indeed all integers. (Otherwise, the associated field is simply nothing.)

Fields

• matrix::Matrix{Int}: the matrix whose eigenvalues and their integrality are of interest.
• spectrum_integral::Bool: whether the eigenvalues of matrix are all integers.
• multiplicities::Union{Nothing,OrderedDict{Int,Int}}: a map from each eigenvalue to its multiplicity, sorted first by ascending multiplicity then by ascending eigenvalue. (This field is nothing if and only if spectrum_integral is false, since we cannot map non-integer eigenvalues to data.)

Notes

If an undirected graph with integer edge weights is $\{-1, 0, 1\}$-diagonalizable (or, more restrictively, $\{-1, 1\}$-diagonalizable), then its Laplacian matrix has integer eigenvalues [JP25, p. 312]. Hence, validating Laplacian integrality serves as a useful screening step in this package's principal $S$-bandwidth minimization algorithm.

References

• [JP25]: N. Johnston and S. Plosker. Laplacian {−1,0,1}- and {−1,1}-diagonalizable graphs. Linear Algebra and its Applications 704, 309–39 (2025). https://doi.org/10.1016/j.laa.2024.10.016.

WeakOrthogonality <: KOrthogonality

The property of "weak orthogonality" for a collection of vectors.

In particular, "weak orthogonality" is an ad hoc term used to refer to the property of $k$-orthogonality for $k > 2$. Recall that an (ordered) collection of vectors $v₁, v₂, ..., vₙ$ is said to be k-orthogonal if we have the inner product $⟨vᵢ, vⱼ⟩ = 0$ whenever $|i - j| ≥ k$ (i.e., if every pair of vectors at least $k$ indices apart is orthogonal). This is equivalent to the vectors' Gram matrix having bandwidth at most $k$.

Fields

• k::Int: the $k$-orthogonality parameter. Must be greater than $2$.

Supertype Hierarchy

WeakOrthogonality <: KOrthogonality

Notes

The term "weak orthogonality" is not standard terminology in the literature, but it is used here to emphasize the weaker nature of this property compared to orthogonality and quasi-orthogonality. It is an ad hoc term coined for this module and is not intended to be formally introduced in the broader literature.

_QOBasisSearchNode

[TODO: Write here]

rank(A::QRPivoted{<:Any, T}; atol::Real=0, rtol::Real=min(n,m)*ϵ) where {T}

Compute the numerical rank of the QR factorization A by counting how many diagonal entries of A.factors are greater than max(atol, rtol*Δ₁) where Δ₁ is the largest calculated such entry. This is similar to the rank(::AbstractMatrix) method insofar as it counts the number of (numerically) nonzero coefficients from a matrix factorization, although the default method uses an SVD instead of a QR factorization. Like rank(::SVD), this method also re-uses an existing matrix factorization.

Computing rank via QR factorization should almost always produce the same results as via SVD, although this method may be more prone to overestimating the rank in pathological cases where the matrix is ill-conditioned. It is also worth noting that it is generally faster to compute a QR factorization than an SVD, so this method may be preferred when performance is a concern.

atol and rtol are the absolute and relative tolerances, respectively. The default relative tolerance is n*ϵ, where n is the size of the smallest dimension of A and ϵ is the eps of the element type of A.

_assert_graph_has_defined_s_bandwidth(g) -> Nothing

Validate that g is an undirected graph with no self-loops.

The property of $S$-bandwidth is only defined, by [JP25]'s conventions, for undirected graphs without self-loops, so this function checks that g satisfies these conditions.

Arguments

• g::AbstractGraph: the graph of interest.

Returns

• nothing: if the check is passed, no output is produced.

Throws

• DomainError: if g is directed or has self-loops.

Examples

The Petersen graph passes the check:

julia> using Graphs

julia> g = PetersenGraph()
{10, 15} undirected simple Int64 graph

julia> isnothing(SDiagonalizability._assert_graph_has_defined_s_bandwidth(g))
true

This random tournament digraph fails the check:

julia> using Graphs

julia> g = random_tournament_digraph(3; seed=13)
{3, 3} directed simple Int64 graph

julia> SDiagonalizability._assert_graph_has_defined_s_bandwidth(g)
ERROR: DomainError with SimpleDiGraph{Int64}(3, [[2, 3], [3], Int64[]], [Int64[], [1], [1, 2]]):
S-bandwidth is not defined for directed graphs
[...]

This multigraph with self-loops fails the check:

julia> using Graphs

julia> g = SimpleGraph(3)
{3, 0} undirected simple Int64 graph

julia> add_edge!(g, 1, 1)
true

julia> SDiagonalizability._assert_graph_has_defined_s_bandwidth(g)
ERROR: DomainError with SimpleGraph{Int64}(1, [[1], Int64[], Int64[]]):
S-bandwidth is not defined for multigraphs; got a graph with self-loops
[...]

Notes

At first blush, it may seem as though the choice of DomainError over something like ArgumentError (or even simply the return of a boolean) constitutes poor design. However, this is informed by the simple ad hoc use of this function to validate inputs for other functions requiring Laplacian matrices. Certainly, this function is never meant to be publicly exposed on its own.

References

• [JP25]: N. Johnston and S. Plosker. Laplacian {−1,0,1}- and {−1,1}-diagonalizable graphs. Linear Algebra and its Applications 704, 309–39 (2025). https://doi.org/10.1016/j.laa.2024.10.016.

_assert_matrix_is_undirected_laplacian(L) -> Nothing

Validate that L is the Laplacian matrix of an undirected, possibly weighted graph.

Arguments

• L::AbstractMatrix{<:Integer}: the purported Laplacian matrix.

Returns

• nothing: if the check is passed, no output is produced.

Throws

• DomainError: if L is not symmetric or has nonzero row sums.

Examples

The Laplacian matrix of the (undirected) star graph on $5$ vertices passes the check:

julia> L = [ 4  -1  -1  -1  -1;
            -1   1   0   0   0;
            -1   0   1   0   0;
            -1   0   0   1   0;
            -1   0   0   0   1]
5×5 Matrix{Int64}:
  4  -1  -1  -1  -1
 -1   1   0   0   0
 -1   0   1   0   0
 -1   0   0   1   0
 -1   0   0   0   1

julia> isnothing(SDiagonalizability._assert_matrix_is_undirected_laplacian(L))
true

The adjacency matrix of the (undirected) cycle graph on $4$ vertices is symmetric but has nonzero row sums, so it fails the check:

julia> A = [0  1  0  1;
            1  0  1  0;
            0  1  0  1;
            1  0  1  0]
4×4 Matrix{Int64}:
 0  1  0  1
 1  0  1  0
 0  1  0  1
 1  0  1  0

julia> SDiagonalizability._assert_matrix_is_undirected_laplacian(A)
ERROR: DomainError with [0 1 0 1; 1 0 1 0; 0 1 0 1; 1 0 1 0]:
Matrix has nonzero row sums; cannot be an (undirected) Laplacian
[...]

Both the in-degree and out-degree Laplacian matrices of this random tournament digraph have zero row sums but are not symmetric, so they fail the check. (These are the two standard ways of extending the concept of the Laplacian to directed graphs [VL20, p. 196].)

julia> using Graphs

julia> G = random_tournament_digraph(3; seed=87)
{3, 3} directed simple Int64 graph

julia> L_in = laplacian_matrix(G; dir=:in)
3×3 SparseArrays.SparseMatrixCSC{Int64, Int64} with 5 stored entries:
  ⋅  ⋅   ⋅
 -1  2  -1
 -1  ⋅   1

julia> SDiagonalizability._assert_matrix_is_undirected_laplacian(L_in)
ERROR: DomainError with sparse([2, 3, 2, 2, 3], [1, 1, 2, 3, 3], [-1, -1, 2, -1, 1], 3, 3):
Matrix is not symmetric; cannot be an (undirected) Laplacian
[...]

julia> L_out = laplacian_matrix(G; dir=:out)
3×3 SparseArrays.SparseMatrixCSC{Int64, Int64} with 5 stored entries:
 2  -1  -1
 ⋅   ⋅   ⋅
 ⋅  -1   1

julia> SDiagonalizability._assert_matrix_is_undirected_laplacian(L_out)
ERROR: DomainError with sparse([1, 1, 3, 1, 3], [1, 2, 2, 3, 3], [2, -1, -1, -1, 1], 3, 3):
Matrix is not symmetric; cannot be an (undirected) Laplacian
[...]

Notes

If edges are to be bidirectional, then L must be symmetric. L must also have zero row sums, since the $(i, i)$-th entry is the weighted degree of node $i$ (the sum of all incident edges' weights) and the $(i, j)$-th entry for $i ≠ j$ is the negation of the weight of edge $(i, j)$ (or simply $0$, if no such edge exists).

Given the highly optimized, lazy, zero-allocation implementation of LinearAlgebra.issymmetric, the symmetry check is performed first. (Both steps are $O(n²)$ in the worst case, but testing for symmetry is far more performant in practice.) This also allows us to (also lazily) check for nonzero column sums rather than nonzero row sums (since these are equivalent for symmetric matrices) in the second step, taking advantage of Julia's column-major storage model.

At first blush, it may seem as though the choice of DomainError over something like ArgumentError (or even simply the return of a boolean) constitutes poor design. However, this is informed by the simple ad hoc use of this function to validate inputs for other functions requiring Laplacian matrices. Certainly, this function is never meant to be publicly exposed on its own.

References

• [VL20]: J. J. Veerman and R. Lyons. A Primer on Laplacian Dynamics in Directed Graphs. Nonlinear Phenomena in Complex Systems 23, 196–206 (2020). https://doi.org/10.33581/1561-4085-2020-23-2-196-206.

_classified_laplacian_01neg_spectra(CL) -> SSpectra

[TODO: Write here. Also, comment inline]

_classified_laplacian_1neg_spectra(CL) -> SSpectra

[TODO: Write here. Also, comment inline]

_classified_laplacian_s_spectra(CL, S) -> SSpectra

[TODO: Write here]

_extract_independent_cols(A) -> Matrix{Int}

Return a (not necessarily unique) independent spanning subset of the columns of A.

Computing a rank-revealing (pivoted) QR decomposition of A, the scaling coefficients from the orthogonalization process are used to determine the rank (rather than recompute it with an SVD), while the pivots are used to extract a spanning set of independent columns.

The rank-revealing Businger–Golub QR algorithm is used for the pivoting strategy, appending the "most independent" column with respect to the current set of pivots at each step via Householder transformations [BG65, pp. 269–70].

Arguments

• A::AbstractMatrix{T<:Integer}: the matrix whose independent columns to extract.

Returns

• ::Matrix{Int}: a spanning set of independent columns of A.

Examples

Observe how columns with greater Euclidean norms are given priority in the pivot ordering:

julia> A = [3  0  0  0  2  1   5   0
            0  3  0  0  2  1  -5   0
            0  0  3  0  2  1   5   4
            0  0  0  3  2  1   0  -4
            0  0  0  0  0  0   0   0]
5×8 Matrix{Int64}:
 3  0  0  0  2  1   5   0
 0  3  0  0  2  1  -5   0
 0  0  3  0  2  1   5   4
 0  0  0  3  2  1   0  -4
 0  0  0  0  0  0   0   0

julia> SDiagonalizability._extract_independent_cols(A)
5×4 Matrix{Int64}:
  5   0  2  3
 -5   0  2  0
  5   4  2  0
  0  -4  2  0
  0   0  0  0

Notes

Since we already need a pivoted QR decomposition to identify independent columns of A (or, rather, to order the columns in such a way that the first rank(A) ones are guaranteed to be independent), it makes sense to use data from the resulting factorization object to compute the rank of A rather than compute a separate SVD. We thus count the nonzero scaling coefficients—that is, the diagonal entries of the R matrix in A = QR—to determine the rank, similarly to how we count the nonzero singular values in an SVD.

It is worth noting that we manually specify a higher relative tolerance for this rank computation. Further discussion can be found in the _rank_rtol documentation, but in short, a critical part of the formula for LinearAlgebra.rank's default rtol uses the minimum dimension of the input matrix. This may result in rank overestimation for tall-and-skinny and short-and-fat matrices (precisely the type we expect to encounter when dealing with all $\{-1, 0, 1\}$-eigenvectors of a Laplacian matrix, which is the intended use case of this helper function in this package). Our replacement tolerance, on the other hand, is a widely accepted standard in numerical analysis which uses the maximum dimension instead [PTVF07, p. 795].

References

• [BG65]: P. Businger and G. H. Golub. Linear Least Squares Solutions by Householder Transformations. Numerische Mathematik 7, 269–76 (1965). https://doi.org/10.1007/BF01436084.

• [PTVF07]: W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery. Numerical Recipes: The Art of Scientific Computing. 3rd Edition (Cambridge University Press, Cambridge, UK, 2007). ISBN: 978-0-521-88068-8. https://dl.acm.org/doi/10.5555/1403886.

_find_basis_idxs_with_prop(
    curr_idxs, prop::Orthogonality, col_space, col_rank, depth
) -> Union{Nothing,AbstractVector{Int}}
_find_basis_idxs_with_prop(
    curr_idxs, prop::QuasiOrthogonality, col_space, col_rank, depth, nodes, union_find
) -> Union{Nothing,AbstractVector{Int}}
_find_basis_idxs_with_prop(
    curr_idxs, prop::WeakOrthogonality, col_space, col_rank, depth, gram_matrix
) -> Union{Nothing,AbstractVector{Int}}

[TODO: Write here]

_find_basis_with_property(col_space::T, col_rank, prop) -> Union{Nothing,T}

[TODO: Write here]

_laplacian_01neg_spectra(L) -> SSpectra

[TODO: Write here]

_laplacian_1neg_spectra(L) -> SSpectra

[TODO: Write here]

_pot_kernel_01neg_eigvecs(n) -> Iterators.Flatten{<:Base.Generator}

Lazily compute all potential kernel $\{-1, 0 ,1\}$-eigenvectors of an n×n Laplacian.

Each vector is normalized so that its first nonzero entry is $1$, enforcing pairwise linear independence between all generated vectors.

Arguments

• n::Integer: the order of the Laplacian matrix of some undirected graph for which to find potential kernel $\{-1, 0, 1\}$-eigenvectors.

Returns

• ::Iterators.Flatten{<:Base.Generator}: a lazily evaluated iterator over all $\{-1, 0, 1\}$-vectors in $ℝⁿ$, unique up to span. Eltype is Vector{Int}.

Examples

Generate all potential kernel $\{-1, 0, 1\}$-eigenvectors of a $3×3$ Laplacian matrix:

julia> stack(SDiagonalizability._pot_kernel_01neg_eigvecs(3))
3×13 Matrix{Int64}:
  1   1   1   1  1  1   1  1  1   0  0  0  0
 -1  -1  -1   0  0  0   1  1  1   1  1  1  0
 -1   0   1  -1  0  1  -1  0  1  -1  0  1  1

Notes

The number of potential kernel $\{-1, 0, 1\}$-eigenvectors (unique up to span) of an $n×n$ Laplacian matrix is equal to $(3ⁿ - 1) / 2$. See also the relevant OEIS sequence [Slo10].

Regrettably, the implementation here is rather clunky and unidiomatic, but it is worth noting that eigenvector generation is one of two major bottlenecks in the overall $S$-bandwidth minimization algorithm. Given how much potential there is for optimization in this piece of code, we thus prioritize performance over readability in this particular case, making every effort to include inline comments wherever clarification may be needed.

References

• [Slo10]: N. J. Sloane, a(n) = (3^n - 1)/2. Entry A003462 (2010). Accessed: 2025-05-22. https://oeis.org/A003462.

_pot_kernel_1neg_eigvecs(n) -> Base.Generator

Lazily compute all potential kernel $\{-1, 1\}$-eigenvectors of an n×n Laplacian.

Each vector is normalized so that its first entry is $1$, enforcing pairwise linear independence between all generated vectors.

Arguments

• n::Integer: the order of the Laplacian matrix of some undirected graph for which to find potential kernel $\{-1, 1\}$-eigenvectors.

Returns

• ::Base.Generator: a lazily evaluated iterator over all $\{-1, 1\}$-vectors in $ℝⁿ$, unique up to span. Eltype is Vector{Int}.

Examples

Generate all potential kernel $\{-1, 1\}$-eigenvectors of a $5×5$ Laplacian matrix:

julia> stack(SDiagonalizability._pot_kernel_1neg_eigvecs(5))
5×16 Matrix{Int64}:
  1   1   1   1   1   1   1   1   1   1   1   1   1   1   1  1
 -1  -1  -1  -1  -1  -1  -1  -1   1   1   1   1   1   1   1  1
 -1  -1  -1  -1   1   1   1   1  -1  -1  -1  -1   1   1   1  1
 -1  -1   1   1  -1  -1   1   1  -1  -1   1   1  -1  -1   1  1
 -1   1  -1   1  -1   1  -1   1  -1   1  -1   1  -1   1  -1  1

Notes

The number of potential kernel $\{-1, 1\}$-eigenvectors (unique up to span) of an $n×n$ Laplacian matrix is equal to $0$ for $n = 0$ and $2ⁿ$ for $n > 0$. See also the relevant OEIS sequence [Slo14] for the $n > 0$ case.

References

• [Slo14]: N. J. Sloane, Powers of 2: a(n) = 2^n. Entry A000079 (2014). Accessed: 2025-08-13. https://oeis.org/A000079.

_pot_nonkernel_01neg_eigvecs(n) -> Iterators.Flatten{<:Base.Generator}

Lazily compute all potential non-kernel $\{-1, 0, 1\}$-eigenvectors of an n×n Laplacian.

Each vector is normalized so that its first nonzero entry is $1$, enforcing pairwise linear independence between all generated vectors. Since all Laplacian matrices have pairwise orthogonal eigenspaces and the all-ones vector is always in the kernel, every non-kernel $\{-1, 0, 1\}$-eigenvector must also have an equal number of $-1$'s and $1$'s.

Arguments

• n::Integer: the order of the Laplacian matrix of some undirected graph for which to find potential non-kernel $\{-1, 0, 1\}$-eigenvectors.

Returns

• ::Iterators.Flatten{<:Base.Generator}: a lazily evaluated iterator over all $\{-1, 0, 1\}$-vectors in $ℝⁿ$ orthogonal to the all-ones kernel vector, unique up to span. Eltype is Vector{Int}.

Examples

Generate all potential non-kernel $\{-1, 0, 1\}$-eigenvectors of a $4×4$ Laplacian matrix:

julia> stack(SDiagonalizability._pot_nonkernel_01neg_eigvecs(4))
4×9 Matrix{Int64}:
  1   1   1   1   1   1   0   0   0
 -1   0   0  -1  -1   1   1   1   0
  0  -1   0  -1   1  -1  -1   0   1
  0   0  -1   1  -1  -1   0  -1  -1

Notes

The number of potential non-kernel $\{-1, 0, 1\}$-eigenvectors (unique up to span) of an $n×n$ Laplacian matrix is, by non-trivial combinatorial arguments, equal to the number of humps in all Motzkin paths of length $n$. See also the relevant OEIS sequence [Deu21].

Regrettably, the implementation here is rather clunky and unidiomatic, but it is worth noting that eigenvector generation is one of two major bottlenecks in the overall $S$-bandwidth minimization algorithm. Given how much potential there is for optimization in this piece of code, we thus prioritize performance over readability in this particular case, making every effort to include inline comments wherever clarification may be needed.

References

• [Deu21]: E. Deutsch. Number of humps in all Motzkin paths of length n. Entry A097861 (2021). Accessed: 2025-05-22. https://oeis.org/A097861.

_pot_nonkernel_1neg_eigvecs(n) -> Base.Generator

Lazily compute all potential non-kernel $\{-1, 1\}$-eigenvectors of an n×n Laplacian.

Each vector is normalized so that its first entry is $1$, enforcing pairwise linear independence between all generated vectors. Since all Laplacian matrices have pairwise orthogonal eigenspaces and the all-ones vector is always in the kernel, every non-kernel $\{-1, 1\}$-eigenvector must half exactly n / 2 $-1$'s and n / 2 $1$'s. (As a direct corollary, if n is odd, an empty iterator is returned.)

Arguments

• n::Integer: the order of the Laplacian matrix of some undirected graph for which to find potential non-kernel $\{-1, 1\}$-eigenvectors.

Returns

• ::Base.Generator: a lazily evaluated iterator over all\{-1, 1\}-vectors inℝⁿorthogonal to the all-ones kernel vector, unique up to span. Eltype isVector{Int}`.

Examples

Generate all potential non-kernel $\{-1, 1\}$-eigenvectors of a $6×6$ Laplacian matrix:

julia> stack(SDiagonalizability._pot_nonkernel_1neg_eigvecs(6))
6×10 Matrix{Int64}:
  1   1   1   1   1   1   1   1   1   1
 -1  -1  -1  -1  -1  -1   1   1   1   1
 -1  -1  -1   1   1   1  -1  -1  -1   1
 -1   1   1  -1  -1   1  -1  -1   1  -1
  1  -1   1  -1   1  -1  -1   1  -1  -1
  1   1  -1   1  -1  -1   1  -1  -1  -1

Notes

The number of potential non-kernel $\{-1, 1\}$-eigenvectors (unique up to span) of an $n×n$ Laplacian matrix is equal to $n$ for $n ≤ 1$ and, by non-trivial combinatorial arguments, the number of Motzkin $(n + 1)$-paths with exactly one flat step for $n > 1$. See also the relevant OEIS sequence [Sut14] for the $n > 1$ case.

References

• [Sut14]: A. V. Sutherland. The number of Motzkin n-paths with exactly one flat step. Entry A138364 (2014). Accessed: 2025-08-13.

function _rank_rtol(A::AbstractMatrix{<:Real}) -> Float64

Return a reasonable relative tolerance for computing matrix rank via SVD or QRD.

The output is intended to be passed as a keyword argument to LinearAlgebra.rank. LinearAlgebra.eigtype is used to determine the appropriate machine epsilon for the element type of A when eltype(A) is not an AbstractFloat.

Arguments

• A::AbstractMatrix{<:Real}: the matrix for which to compute a tolerance.

Returns

• tol::Float64: a reasonable relative tolerance for computing matrix rank via SVD or QRD. This scales proportionally to the maximum dimension of A.

Notes

LinearAlgebra.rank's default rtol of min(m,n) * ϵ for computing the rank of an $m×n$ matrix may result in overestimating rank when $|m - n| ≫ 0$, since condition number (which determines how numerically stable SVD and QRD are) grows with both dimensions [CD05, p. 603]. Given that we often deal with short-and-fat matrices in this package (particularly when processing all $\{-1, 0, 1\}$-eigenvectors of a Laplacian matrix), we turn instead to the same relative tolerance used by MATLAB's and NumPy's rank functions—max(m,n) * ϵ [MAT25; Num25]. (Indeed, this is a widely adopted standard across the field of numerical analysis [PTVF07, p. 795].)

References

• [CD05]: Z. Chen and J. Dongarra. Condition Numbers of Gaussian Random Matrices. SIAM Journal on Matrix Analysis and Applications 27, 603–20 (2005). https://doi.org/10.1137/040616413.
• [MAT25]: MATLAB Developers, rank. MATLAB reference documentation – R2025a (2025). Accessed: 2025-05-29. https://www.mathworks.com/help/matlab/ref/rank.html.
• [Num25]: NumPy Developers, numpy.linalg.matrix_rank. NumPy reference documentation – v2.2 (2025). Accessed: 2025-05-22. https://numpy.org/doc/stable/reference/generated/numpy.linalg.matrix_rank.html.
• [PTVF07]: W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery. Numerical Recipes: The Art of Scientific Computing. 3rd Edition (Cambridge University Press, Cambridge, UK, 2007). ISBN: 978-0-521-88068-8. https://dl.acm.org/doi/10.5555/1403886.

_sort_tuple(S::T) where {T<:Tuple{Vararg{Integer}}} -> T

Return a sorted version of the input tuple.

Arguments

• S::T<:Tuple{Vararg{Integer}}: the tuple to sort.

Returns

• ::T: a sorted version of S.

Examples

julia> S = (-1, 1, 0)
(-1, 1, 0)

julia> SDiagonalizability._sort_tuple(S)
(-1, 0, 1)

check_spectrum_integrality(A) -> SpectrumIntegralResult

Check whether the eigenvalues of A are integers (up to floating-point error).

If the eigenvalues are indeed all integers, then an eigenvalue-multiplicity map is constructed as well.

Arguments

• A::AbstractMatrix{<:Integer}: the matrix whose eigenvalues to check for integrality.

Returns

• ::SpectrumIntegralResult: a struct containing the following fields:
  • matrix: the A matrix, copied to avoid shared mutability.
  • spectrum_integral: whether the eigenvalues of A are integers.
  • multiplicities: a map from each eigenvalue to its multiplicity, sorted first by ascending multiplicity then by ascending eigenvalue. (This field is nothing if the eigenvalues are not all integers.)

Examples

Confirm that the rotation matrix by $π/2$ radians counterclockwise is not spectrum integral (rather, it has eigenvalues $±i$ [Joy15, p. 1]):

julia> R = Int8.([0 -1; 1 0])
2×2 Matrix{Int8}:
 0  -1
 1   0

julia> res = SDiagonalizability.check_spectrum_integrality(R);

julia> res.matrix
2×2 Matrix{Int64}:
 0  -1
 1   0

julia> res.spectrum_integral
false

julia> isnothing(res.multiplicities)
true

Confirm that the adjacency matrix of the Petersen graph is spectrum integral, with correct eigenvalues and multiplicities of $\{3: 1, -2: 4, 1: 5\}$ [Fox09, p. 2]:

julia> using Graphs

julia> G = smallgraph(:petersen)
{10, 15} undirected simple Int64 graph

julia> A = adjacency_matrix(G)
10×10 SparseArrays.SparseMatrixCSC{Int64, Int64} with 30 stored entries:
 ⋅  1  ⋅  ⋅  1  1  ⋅  ⋅  ⋅  ⋅
 1  ⋅  1  ⋅  ⋅  ⋅  1  ⋅  ⋅  ⋅
 ⋅  1  ⋅  1  ⋅  ⋅  ⋅  1  ⋅  ⋅
 ⋅  ⋅  1  ⋅  1  ⋅  ⋅  ⋅  1  ⋅
 1  ⋅  ⋅  1  ⋅  ⋅  ⋅  ⋅  ⋅  1
 1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1  1  ⋅
 ⋅  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1  1
 ⋅  ⋅  1  ⋅  ⋅  1  ⋅  ⋅  ⋅  1
 ⋅  ⋅  ⋅  1  ⋅  1  1  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  1  ⋅  1  1  ⋅  ⋅

julia> res = SDiagonalizability.check_spectrum_integrality(A);

julia> res.matrix
10×10 Matrix{Int64}:
 0  1  0  0  1  1  0  0  0  0
 1  0  1  0  0  0  1  0  0  0
 0  1  0  1  0  0  0  1  0  0
 0  0  1  0  1  0  0  0  1  0
 1  0  0  1  0  0  0  0  0  1
 1  0  0  0  0  0  0  1  1  0
 0  1  0  0  0  0  0  0  1  1
 0  0  1  0  0  1  0  0  0  1
 0  0  0  1  0  1  1  0  0  0
 0  0  0  0  1  0  1  1  0  0

julia> res.spectrum_integral
true

julia> res.multiplicities
OrderedCollections.OrderedDict{Int64, Int64} with 3 entries:
  3  => 1
  -2 => 4
  1  => 5

Notes

If an undirected graph with integer edge weights is $\{-1, 0, 1\}$-diagonalizable (or, more restrictively, $\{-1, 1\}$-diagonalizable), then its Laplacian matrix has integer eigenvalues [JP25, p. 312]. Hence, validating Laplacian integrality serves as a useful screening step in this package's principal $S$*-bandwidth minimization algorithm.

References

• [Fox09]: J. Fox. Lecture 19: The Petersen graph and Moore graphs. Lecture notes, MAT 307: Combinatorics (2009). Accessed: 2025-07-25. https://math.mit.edu/~fox/MAT307.html.
• [Joy15]: D. Joyce. Rotations and complex eigenvalues. Lecture notes, Math 130: Linear Algebra (2015). http://aleph0.clarku.edu/~ma130/complexeigen.pdf.
• [JP25]: N. Johnston and S. Plosker. Laplacian {−1,0,1}- and {−1,1}-diagonalizable graphs. Linear Algebra and its Applications 704, 309–39 (2025). https://doi.org/10.1016/j.laa.2024.10.016.

classify_k_orthogonality(k)

Classifies the k-orthogonality property based on the given k parameter.

When searching for a $k$-orthogonal $S$-basis of a given Laplacian eigenspace, the family of values to which our $k$ parameter belongs informs our choice of algorithm.

Arguments

• k::Integer: the k-orthogonality parameter to classify. Must be a positive integer.

Returns

• ::KOrthogonality: An instance of a concrete KOrthogonality subtype corresponding to k.

Throws

• DomainError: if $k$ is not a positive integer.

Examples

julia> SDiagonalizability.classify_k_orthogonality(1)
SDiagonalizability.Orthogonality(1)

julia> SDiagonalizability.classify_k_orthogonality(2)
SDiagonalizability.QuasiOrthogonality(2)

julia> SDiagonalizability.classify_k_orthogonality(3)
SDiagonalizability.WeakOrthogonality(3)

julia> SDiagonalizability.classify_k_orthogonality(13)
SDiagonalizability.WeakOrthogonality(13)

classify_laplacian(L)

Classify the Laplacian matrix L and wrap it in a ClassifiedLaplacian object.

It is first verified that L is indeed a Laplacian matrix by _assert_matrix_is_undirected_laplacian, which throws a DomainError otherwise. It is then classified based on any properties which may be exploited in computing data on its \{-1, 0, 1\}-spectrum.

Arguments

• L::AbstractMatrix{<:Integer}: the Laplacian matrix to classify.

Returns

• CL::ClassifiedLaplacian: the Laplacian wrapped in a concrete ClassifiedLaplacian subtype associated with the category of the graph represented by L. In the case of complete graphs, CL also contains data on the (necessarily uniform) weight of each edge.

Examples

Correctly recognizes the Laplacian matrix of the null graph:

julia> using Graphs

julia> G = Graph(0)
{0, 0} undirected simple Int64 graph

julia> L = laplacian_matrix(G)
0×0 SparseArrays.SparseMatrixCSC{Int64, Int64} with 0 stored entries

julia> SDiagonalizability.classify_laplacian(L)
SDiagonalizability.NullGraphLaplacian(Matrix{Int64}(undef, 0, 0))

Correctly recognizes the Laplacian matrix of the empty graph on $3$ nodes:

julia> using Graphs

julia> G = Graph(3)
{3, 0} undirected simple Int64 graph

julia> L = laplacian_matrix(G)
3×3 SparseArrays.SparseMatrixCSC{Int64, Int64} with 0 stored entries:
 ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅

julia> SDiagonalizability.classify_laplacian(L)
SDiagonalizability.EmptyGraphLaplacian([0 0 0; 0 0 0; 0 0 0])

Correctly recognizes the Laplacian matrix of the complete graph on $4$ nodes:

julia> using Graphs

julia> G = complete_graph(4)
{4, 6} undirected simple Int64 graph

julia> L = laplacian_matrix(G)
4×4 SparseArrays.SparseMatrixCSC{Int64, Int64} with 16 stored entries:
  3  -1  -1  -1
 -1   3  -1  -1
 -1  -1   3  -1
 -1  -1  -1   3

julia> SDiagonalizability.classify_laplacian(L)
SDiagonalizability.CompleteGraphLaplacian([3 -1 -1 -1; -1 3 -1 -1; -1 -1 3 -1; -1 -1 -1 3], 1)

Correctly recognizes the Laplacian matrix of this random generic graph:

julia> using Graphs

julia> G = erdos_renyi(5, 8; seed=87)
{5, 8} undirected simple Int64 graph

julia> L = laplacian_matrix(G)
5×5 SparseArrays.SparseMatrixCSC{Int64, Int64} with 21 stored entries:
  3   ⋅  -1  -1  -1
  ⋅   2  -1  -1   ⋅
 -1  -1   4  -1  -1
 -1  -1  -1   4  -1
 -1   ⋅  -1  -1   3

julia> SDiagonalizability.classify_laplacian(L)
SDiagonalizability.ArbitraryGraphLaplacian([3 0 … -1 -1; 0 2 … -1 0; … ; -1 -1 … 4 -1; -1 0 … -1 3])

Notes

Right now, the possible types of Laplacian matrices (or, to be more precise, the graphs they represent) are:

• NullGraphLaplacian: the (unique) graph with no nodes.
• EmptyGraphLaplacian: any graph with no edges on $n ≥ 1$ nodes.
• CompleteGraphLaplacian: any graph on $n ≥ 2$ nodes where every pair of nodes is connected by an edge and all edges possess the same weight.
• ArbitraryGraphLaplacian: any graph on $n ≥ 3$ nodes with no particular strcuture of relevance in the context of investigating $\{-1, 0, 1\}$-spectra.

find_k_orthogonal_basis(
    col_space, col_rank, k, pre_basis=nothing
) -> Union{Nothing,AbstractMatrix{<:Integer}}

Find a k-orthogonal spanning subset of the column set of a matrix.

This function does not allow for any k-orthogonal basis of col_space (in which case it would be trivial to construct a pairwise orthogonal basis via a QR decomposition). Instead, it restricts its search to spanning subsets comprised exclusively of columns already in col_space, returning nothing if no such basis exists.

Arguments

• col_space::T<:AbstractMatrix{<:Integer}: the matrix whose column space is searched for a k-orthogonal basis.
• col_rank::Integer: the rank of the column space of col_space, pre-computed for efficiency.
• k::Integer: the minimum $k$-orthogonality parameter of the desired basis. Must be a positive integer.
• pre_basis::Union{Nothing,AbstractMatrix{<:Integer}}=nothing: an optional precomputed submatrix of col_space whose columns form a spanning subset of the column space of col_space. In case col_space has at most k columns, existing_basis is guaranteed to already be a k-orthogonal basis and thus is used to skip unnecessary computations.

Returns

• ::Union{Nothing,AbstractMatrix{<:Integer}}: a k-orthogonal spanning subset of the columns of col_space, if one exists; otherwise, nothing.

laplacian_s_spectra(L, S) -> SSpectra

[TODO: Write here. Also, comment inline]

pot_kernel_s_eigvecs(n, S) -> Union{Base.Generator,Base.Iterators.Flatten{<:Base.Generator}}

Lazily compute all potential kernel S-eigenvectors of an n×n Laplacian.

The only accepted values of S are (-1, 0, 1), (-1, 1), and permutations thereof (e.g., (1, -1, 0), which will be sorted internally to (-1, 0, 1) anyway), in line with the sets studied by [JP25].

Each vector is normalized so that its first nonzero entry is $1$, enforcing pairwise linear independence between all generated vectors.

Arguments

• n::Integer: the order of the Laplacian matrix of some undirected graph for which to find potential kernel S-eigenvectors.
• S::Tuple{Vararg{Integer}}: the set of possible entries of the potential eigenvectors, provided in the form of a tuple of unique Integer's. Must be either (-1, 0, 1), (-1, 1), or a permutation thereof.

Returns

• gen::Union{Base.Generator,Base.Iterators.Flatten{<:Base.Generator}}: a lazily evaluated iterator over all S-vectors in $ℝⁿ$, unique up to span. Eltype is Vector{Int}.

Throws

• DomainError: if n is negative.
• ArgumentError: if S is not (-1, 0, 1), (-1, 1), or a permutation thereof.

Examples

Generate all potential kernel $\{-1, 0, 1\}$-eigenvectors of a $4×4$ Laplacian matrix:

julia> stack(SDiagonalizability.pot_kernel_s_eigvecs(4, (-1, 0, 1)))
4×40 Matrix{Int64}:
  1   1   1   1   1   1   1   1   1  …   0   0  0  0   0  0  0   0  0  0  0
 -1  -1  -1  -1  -1  -1  -1  -1  -1      1   1  1  1   1  1  1   0  0  0  0
 -1  -1  -1   0   0   0   1   1   1     -1   0  0  0   1  1  1   1  1  1  0
 -1   0   1  -1   0   1  -1   0   1      1  -1  0  1  -1  0  1  -1  0  1  1

Generate all potential kernel $\{-1, 1\}$-eigenvectors of a $5×5$ Laplacian matrix, with the S set provided in a different order:

julia> stack(SDiagonalizability.pot_kernel_s_eigvecs(5, (1, -1)))
5×16 Matrix{Int64}:
  1   1   1   1   1   1   1   1   1   1   1   1   1   1   1  1
 -1  -1  -1  -1  -1  -1  -1  -1   1   1   1   1   1   1   1  1
 -1  -1  -1  -1   1   1   1   1  -1  -1  -1  -1   1   1   1  1
 -1  -1   1   1  -1  -1   1   1  -1  -1   1   1  -1  -1   1  1
 -1   1  -1   1  -1   1  -1   1  -1   1  -1   1  -1   1  -1  1

Only the first two entry sets are supported, so the following command throws an error:

julia> stack(SDiagonalizability.pot_kernel_s_eigvecs(6, (1, 2)))
ERROR: ArgumentError: Unsupported entry set S: (1, 2)
[...]

References

• [JP25]: N. Johnston and S. Plosker. Laplacian {−1,0,1}- and {−1,1}-diagonalizable graphs. Linear Algebra and its Applications 704, 309–39 (2025). https://doi.org/10.1016/j.laa.2024.10.016.

pot_nonkernel_s_eigvecs(n, S) -> Union{Base.Generator,Base.Iterators.Flatten{<:Base.Generator}}

Lazily compute all potential non-kernel S-eigenvectors of an n×n Laplacian.

The only accepted values of S are (-1, 0, 1), (-1, 1), and permutations thereof (e.g., (1, -1, 0), which will be sorted internally to (-1, 0, 1) anyway), in line with the sets studied by [JP25].

Arguments

• n::Integer: the order of the Laplacian matrix of some undirected graph for which to find potential non-kernel S-eigenvectors.
• S::Tuple{Vararg{Integer}}: the set of possible entries of the potential eigenvectors, provided in the form of a tuple of unique Integer's. Must be either (-1, 0, 1), (-1, 1), or a permutation thereof.

Returns

• gen::Union{Base.Generator,Base.Iterators.Flatten{<:Base.Generator}}: a lazily evaluated iterator over all S-vectors in $ℝⁿ$ orthogonal to the all-ones kernel vector, unique up to span. Eltype is Vector{Int}.

Throws

• DomainError: if n is negative.
• ArgumentError: if S is not (-1, 0, 1), (-1, 1), or a permutation thereof.

Examples

Generate all potential non-kernel $\{-1, 0, 1\}$-eigenvectors of a $5×5$ Laplacian matrix, with the S set provided in a different order:

julia> stack(SDiagonalizability.pot_nonkernel_s_eigvecs(5, (1, -1, 0)))
5×25 Matrix{Int64}:
  1   1   1   1   1   1   1   1   1  …   0   0   0   0   0   0   0   0   0
 -1   0   0   0  -1  -1  -1  -1  -1      1   1   1   1   1   1   0   0   0
  0  -1   0   0  -1  -1   0   0   1     -1   0   0  -1  -1   1   1   1   0
  0   0  -1   0   0   1  -1   1  -1      0  -1   0  -1   1  -1  -1   0   1
  0   0   0  -1   1   0   1  -1   0      0   0  -1   1  -1  -1   0  -1  -1

Generate all potential non-kernel $\{-1, 1\}$-eigenvectors of a $6×6$ Laplacian matrix:

julia> stack(SDiagonalizability.pot_nonkernel_s_eigvecs(6, (-1, 1)))
6×10 Matrix{Int64}:
  1   1   1   1   1   1   1   1   1   1
 -1  -1  -1  -1  -1  -1   1   1   1   1
 -1  -1  -1   1   1   1  -1  -1  -1   1
 -1   1   1  -1  -1   1  -1  -1   1  -1
  1  -1   1  -1   1  -1  -1   1  -1  -1
  1   1  -1   1  -1  -1   1  -1  -1  -1

Only the first two entry sets are supported, so the following command throws an error:

julia> stack(SDiagonalizability.pot_nonkernel_s_eigvecs(7, (0, 3)))
ERROR: ArgumentError: Unsupported entry set S: (0, 3)
[...]

References

• [JP25]: N. Johnston and S. Plosker. Laplacian {−1,0,1}- and {−1,1}-diagonalizable graphs. Linear Algebra and its Applications 704, 309–39 (2025). https://doi.org/10.1016/j.laa.2024.10.016.