GitHub

API reference

The docstrings on this page define the public API of the package.

SparseMatrixColorings.SparseMatrixColorings — Module
SparseMatrixColorings

SparseMatrixColorings.jl

Build Status Stable Documentation Dev Documentation Coverage Code Style: Blue arXiv DOI

Coloring algorithms for sparse Jacobian and Hessian matrices.

Getting started

To install this package, run the following in a Julia Pkg REPL:

pkg> add SparseMatrixColorings

Background

The algorithms implemented in this package are described in the following preprint:

and inspired by previous works:

Some parts of the articles (like definitions) are thus copied verbatim in the documentation.

Alternatives

Citing

Please cite this software using the provided CITATION.cff file or the .bib entry below:

@unpublished{montoison2025revisitingsparsematrixcoloring,
      title={Revisiting Sparse Matrix Coloring and Bicoloring}, 
      author={Alexis Montoison and Guillaume Dalle and Assefaw Gebremedhin},
      year={2025},
      eprint={2505.07308},
      archivePrefix={arXiv},
      primaryClass={math.NA},
      url={https://arxiv.org/abs/2505.07308}, 
}

The link https://zenodo.org/doi/10.5281/zenodo.11314275 resolves to the latest version on Zenodo.

Exports

source

Main function

SparseMatrixColorings.coloring — Function
coloring(
    S::AbstractMatrix,
    problem::ColoringProblem,
    algo::GreedyColoringAlgorithm;
    [decompression_eltype=Float64, symmetric_pattern=false]
)

Solve a ColoringProblem on the matrix S with a GreedyColoringAlgorithm and return an AbstractColoringResult.

The result can be used to compress and decompress a matrix A with the same sparsity pattern as S. If eltype(A) == decompression_eltype, decompression might be faster.

For a :nonsymmetric problem (and only then), setting symmetric_pattern=true indicates that the pattern of nonzeros is symmetric. This condition is weaker than the symmetry of actual values, so it can happen for some Jacobians. Specifying it allows faster construction of the bipartite graph.

Example

julia> using SparseMatrixColorings, SparseArrays

julia> S = sparse([
           0 0 1 1 0 1
           1 0 0 0 1 0
           0 1 0 0 1 0
           0 1 1 0 0 0
       ]);

julia> problem = ColoringProblem(; structure=:nonsymmetric, partition=:column);

julia> algo = GreedyColoringAlgorithm(; decompression=:direct);

julia> result = coloring(S, problem, algo);

julia> column_colors(result)
6-element Vector{Int64}:
 1
 1
 2
 1
 2
 3

julia> collect.(column_groups(result))
3-element Vector{Vector{Int64}}:
 [1, 2, 4]
 [3, 5]
 [6]

See also

source
SparseMatrixColorings.ColoringProblem — Type
ColoringProblem{structure,partition}

Selector type for the coloring problem to solve, enabling multiple dispatch.

It is passed as an argument to the main function coloring.

Constructors

ColoringProblem{structure,partition}()
ColoringProblem(; structure=:nonsymmetric, partition=:column)
  • structure::Symbol: either :nonsymmetric or :symmetric
  • partition::Symbol: either :column, :row or :bidirectional
Warning

The second constructor (based on keyword arguments) is type-unstable.

Link to automatic differentiation

Matrix coloring is often used in automatic differentiation, and here is the translation guide:

matrixmodestructurepartitionimplemented
Jacobianforward:nonsymmetric:columnyes
Jacobianreverse:nonsymmetric:rowyes
Jacobianmixed:nonsymmetric:bidirectionalyes
Hessian-:symmetric:columnyes
Hessian-:symmetric:rowno
source
SparseMatrixColorings.GreedyColoringAlgorithm — Type
GreedyColoringAlgorithm{decompression} <: ADTypes.AbstractColoringAlgorithm

Greedy coloring algorithm for sparse matrices which colors columns or rows one after the other, following a configurable order.

It is passed as an argument to the main function coloring.

Constructors

GreedyColoringAlgorithm{decompression}(order=NaturalOrder(); postprocessing=false)
GreedyColoringAlgorithm(order=NaturalOrder(); postprocessing=false, decompression=:direct)
  • order::Union{AbstractOrder,Tuple}: the order in which the columns or rows are colored, which can impact the number of colors. Can also be a tuple of different orders to try out, from which the best order (the one with the lowest total number of colors) will be used.
  • postprocessing::Bool: whether or not the coloring will be refined by assigning the neutral color 0 to some vertices.
  • decompression::Symbol: either :direct or :substitution. Usually :substitution leads to fewer colors, at the cost of a more expensive coloring (and decompression). When :substitution is not applicable, it falls back on :direct decompression.
Warning

The second constructor (based on keyword arguments) is type-unstable.

ADTypes coloring interface

GreedyColoringAlgorithm is a subtype of ADTypes.AbstractColoringAlgorithm, which means the following methods are also applicable:

See their respective docstrings for details.

See also

source
SparseMatrixColorings.ConstantColoringAlgorithm — Type
ConstantColoringAlgorithm{partition}  <: ADTypes.AbstractColoringAlgorithm

Coloring algorithm which always returns the same precomputed vector of colors. Useful when the optimal coloring of a matrix can be determined a priori due to its specific structure (e.g. banded).

It is passed as an argument to the main function coloring, but will only work if the associated problem has :nonsymmetric structure. Indeed, for symmetric coloring problems, we need more than just the vector of colors to allow fast decompression.

Constructors

ConstantColoringAlgorithm{partition}(matrix_template, color)
ConstantColoringAlgorithm(matrix_template, color; partition=:column)
  • partition::Symbol: either :row or :column.
  • matrix_template::AbstractMatrix: matrix for which the vector of colors was precomputed (the algorithm will only accept matrices of the exact same size).
  • color::Vector{<:Integer}: vector of integer colors, one for each row or column (depending on partition).
Warning

The second constructor (based on keyword arguments) is type-unstable.

We do not necessarily verify consistency between the matrix template and the vector of colors, this is the responsibility of the user.

Example

julia> using SparseMatrixColorings, LinearAlgebra

julia> matrix_template = Diagonal(ones(Bool, 5))
5×5 Diagonal{Bool, Vector{Bool}}:
 1  ⋅  ⋅  ⋅  ⋅
 ⋅  1  ⋅  ⋅  ⋅
 ⋅  ⋅  1  ⋅  ⋅
 ⋅  ⋅  ⋅  1  ⋅
 ⋅  ⋅  ⋅  ⋅  1

julia> color = ones(Int, 5)  # coloring a Diagonal is trivial
5-element Vector{Int64}:
 1
 1
 1
 1
 1

julia> problem = ColoringProblem(; structure=:nonsymmetric, partition=:column);

julia> algo = ConstantColoringAlgorithm(matrix_template, color; partition=:column);

julia> result = coloring(similar(matrix_template), problem, algo);

julia> column_colors(result)
5-element Vector{Int64}:
 1
 1
 1
 1
 1

ADTypes coloring interface

ConstantColoringAlgorithm is a subtype of ADTypes.AbstractColoringAlgorithm, which means the following methods are also applicable (although they will error if the kind of coloring demanded not consistent):

source

Result analysis

SparseMatrixColorings.AbstractColoringResult — Type
AbstractColoringResult{structure,partition,decompression}

Abstract type for the result of a coloring algorithm.

It is the supertype of the object returned by the main function coloring.

Type parameters

Combination between the type parameters of ColoringProblem and GreedyColoringAlgorithm:

  • structure::Symbol: either :nonsymmetric or :symmetric
  • partition::Symbol: either :column, :row or :bidirectional
  • decompression::Symbol: either :direct or :substitution

Applicable methods

Warning

Unlike the methods above, the concrete subtypes of AbstractColoringResult are not part of the public API and may change without notice.

source
SparseMatrixColorings.ncolors — Function
ncolors(result::AbstractColoringResult)

Return the number of different non-zero colors used to color the matrix.

For bidirectional partitions, this number is the sum of the number of non-zero row colors and the number of non-zero column colors.

source
SparseMatrixColorings.column_groups — Function
column_groups(result::AbstractColoringResult)

Return a vector group such that for every non-zero color c, group[c] contains the indices of all columns that are colored with c.

source
SparseMatrixColorings.row_groups — Function
row_groups(result::AbstractColoringResult)

Return a vector group such that for every non-zero color c, group[c] contains the indices of all rows that are colored with c.

source
SparseMatrixColorings.sparsity_pattern — Function
sparsity_pattern(result::AbstractColoringResult)

Return the matrix that was initially passed to coloring, without any modifications.

Note

This matrix is not necessarily a SparseMatrixCSC, nor does it necessarily have Bool entries.

source

Decompression

SparseMatrixColorings.compress — Function
compress(A, result::AbstractColoringResult)

Compress A given a coloring result of the sparsity pattern of A.

  • If result comes from a :column (resp. :row) partition, the output is a single matrix B compressed by column (resp. by row).
  • If result comes from a :bidirectional partition, the output is a tuple of matrices (Br, Bc), where Br is compressed by row and Bc by column.

Compression means summing either the columns or the rows of A which share the same color. It is undone by calling decompress or decompress!.

Example

julia> using SparseMatrixColorings, SparseArrays

julia> A = sparse([
           0 0 4 6 0 9
           1 0 0 0 7 0
           0 2 0 0 8 0
           0 3 5 0 0 0
       ]);

julia> result = coloring(A, ColoringProblem(), GreedyColoringAlgorithm());

julia> collect.(column_groups(result))
3-element Vector{Vector{Int64}}:
 [1, 2, 4]
 [3, 5]
 [6]

julia> B = compress(A, result)
4×3 Matrix{Int64}:
 6  4  9
 1  7  0
 2  8  0
 3  5  0

See also

source
SparseMatrixColorings.decompress — Function
decompress(B::AbstractMatrix, result::AbstractColoringResult{_,:column/:row})
decompress(Br::AbstractMatrix, Bc::AbstractMatrix, result::AbstractColoringResult{_,:bidirectional})

Decompress B (or the tuple (Br,Bc)) into a new matrix A, given a coloring result of the sparsity pattern of A. The in-place alternative is decompress!.

Compression means summing either the columns or the rows of A which share the same color. It is done by calling compress.

Warning

For some coloring variants, the result object is mutated during decompression.

Example

julia> using SparseMatrixColorings, SparseArrays

julia> A = sparse([
           0 0 4 6 0 9
           1 0 0 0 7 0
           0 2 0 0 8 0
           0 3 5 0 0 0
       ]);

julia> result = coloring(A, ColoringProblem(), GreedyColoringAlgorithm());

julia> collect.(column_groups(result))
3-element Vector{Vector{Int64}}:
 [1, 2, 4]
 [3, 5]
 [6]

julia> B = compress(A, result)
4×3 Matrix{Int64}:
 6  4  9
 1  7  0
 2  8  0
 3  5  0

julia> decompress(B, result)
4×6 SparseMatrixCSC{Int64, Int64} with 9 stored entries:
 ⋅  ⋅  4  6  ⋅  9
 1  ⋅  ⋅  ⋅  7  ⋅
 ⋅  2  ⋅  ⋅  8  ⋅
 ⋅  3  5  ⋅  ⋅  ⋅

julia> decompress(B, result) == A
true

See also

source
SparseMatrixColorings.decompress! — Function
decompress!(
    A::AbstractMatrix, B::AbstractMatrix,
    result::AbstractColoringResult{_,:column/:row}, [uplo=:F]
)

decompress!(
    A::AbstractMatrix, Br::AbstractMatrix, Bc::AbstractMatrix
    result::AbstractColoringResult{_,:bidirectional}
)

Decompress B (or the tuple (Br,Bc)) in-place into A, given a coloring result of the sparsity pattern of A. The out-of-place alternative is decompress.

Note

In-place decompression is faster when A isa SparseMatrixCSC.

Compression means summing either the columns or the rows of A which share the same color. It is done by calling compress.

For :symmetric coloring results (and for those only), an optional positional argument uplo in (:U, :L, :F) can be passed to specify which part of the matrix A should be updated: the Upper triangle, the Lower triangle, or the Full matrix. When A isa SparseMatrixCSC, using the uplo argument requires a target matrix which only stores the relevant triangle(s).

Warning

For some coloring variants, the result object is mutated during decompression.

Example

julia> using SparseMatrixColorings, SparseArrays

julia> A = sparse([
           0 0 4 6 0 9
           1 0 0 0 7 0
           0 2 0 0 8 0
           0 3 5 0 0 0
       ]);

julia> result = coloring(A, ColoringProblem(), GreedyColoringAlgorithm());

julia> collect.(column_groups(result))
3-element Vector{Vector{Int64}}:
 [1, 2, 4]
 [3, 5]
 [6]

julia> B = compress(A, result)
4×3 Matrix{Int64}:
 6  4  9
 1  7  0
 2  8  0
 3  5  0

julia> A2 = similar(A);

julia> decompress!(A2, B, result)
4×6 SparseMatrixCSC{Int64, Int64} with 9 stored entries:
 ⋅  ⋅  4  6  ⋅  9
 1  ⋅  ⋅  ⋅  7  ⋅
 ⋅  2  ⋅  ⋅  8  ⋅
 ⋅  3  5  ⋅  ⋅  ⋅

julia> A2 == A
true

See also

source
SparseMatrixColorings.decompress_single_color! — Function
decompress_single_color!(
    A::AbstractMatrix, b::AbstractVector, c::Integer,
    result::AbstractColoringResult, [uplo=:F]
)

Decompress the vector b corresponding to color c in-place into A, given a :direct coloring result of the sparsity pattern of A (it will not work with a :substitution coloring).

  • If result comes from a :nonsymmetric structure with :column partition, this will update the columns of A that share color c (whose sum makes up b).
  • If result comes from a :nonsymmetric structure with :row partition, this will update the rows of A that share color c (whose sum makes up b).
  • If result comes from a :symmetric structure with :column partition, this will update the coefficients of A whose value is deduced from color c.
Warning

This function will only update some coefficients of A, without resetting the rest to zero.

For :symmetric coloring results (and for those only), an optional positional argument uplo in (:U, :L, :F) can be passed to specify which part of the matrix A should be updated: the Upper triangle, the Lower triangle, or the Full matrix. When A isa SparseMatrixCSC, using the uplo argument requires a target matrix which only stores the relevant triangle(s).

Warning

For some coloring variants, the result object is mutated during decompression.

Example

julia> using SparseMatrixColorings, SparseArrays

julia> A = sparse([
           0 0 4 6 0 9
           1 0 0 0 7 0
           0 2 0 0 8 0
           0 3 5 0 0 0
       ]);

julia> result = coloring(A, ColoringProblem(), GreedyColoringAlgorithm());

julia> collect.(column_groups(result))
3-element Vector{Vector{Int64}}:
 [1, 2, 4]
 [3, 5]
 [6]

julia> B = compress(A, result)
4×3 Matrix{Int64}:
 6  4  9
 1  7  0
 2  8  0
 3  5  0

julia> A2 = similar(A); A2 .= 0;

julia> decompress_single_color!(A2, B[:, 2], 2, result)
4×6 SparseMatrixCSC{Int64, Int64} with 9 stored entries:
 ⋅  ⋅  4  0  ⋅  0
 0  ⋅  ⋅  ⋅  7  ⋅
 ⋅  0  ⋅  ⋅  8  ⋅
 ⋅  0  5  ⋅  ⋅  ⋅

julia> A2[:, [3, 5]] == A[:, [3, 5]]
true

See also

source

Orders

SparseMatrixColorings.AbstractOrder — Type
AbstractOrder

Abstract supertype for the vertex order used inside GreedyColoringAlgorithm.

In this algorithm, the rows and columns of a matrix form a graph, and the vertices are colored one after the other in a greedy fashion. Depending on how the vertices are ordered, the number of colors necessary may vary.

Options

source
SparseMatrixColorings.RandomOrder — Type
RandomOrder(rng=default_rng(), seed=nothing)

Instance of AbstractOrder which sorts vertices using a random permutation, generated from rng with a given seed.

  • If seed = nothing, the rng will never be re-seeded. Therefore, two consecutive calls to vertices(g, order) will give different results.
  • Otherwise, the rng will be re-seeded before each sample. Therefore, two consecutive calls to vertices(g, order) will give the same result.
Warning

Do not use a seed with the default_rng(), otherwise you will affect the global state of your program. If you need reproducibility, create a new rng specifically for your RandomOrder. The package StableRNGs.jl offers random number generators whose behavior is stable across Julia versions.

source
SparseMatrixColorings.DynamicDegreeBasedOrder — Type
DynamicDegreeBasedOrder{degtype,direction}(; reproduce_colpack=false)

Instance of AbstractOrder which sorts vertices using a dynamically computed degree.

This order works by assigning vertices to buckets based on their dynamic degree, and then updating buckets iteratively by transfering vertices between them.

Type parameters

  • degtype::Symbol: can be :forward (for the forward degree) or :back (for the back degree)
  • direction::Symbol: can be :low2high (if the order is defined from lowest to highest, i.e. 1 to n) or :high2low (if the order is defined from highest to lowest, i.e. n to 1)

Concrete variants

Settings

  • reproduce_colpack::Bool: whether to manage the buckets in the exact same way as the original ColPack implementation.
    • When reproduce_colpack=true, we always append and remove vertices at the end of a bucket (unilateral).
    • When reproduce_colpack=false (the default), we can append and remove vertices either at the start or at the end of a bucket (bilateral).

Allowing modifications on both sides of a bucket enables storage optimization, with a single fixed-size vector for all buckets instead of one dynamically-sized vector per bucket. As a result, the default setting reproduce_colpack=false is slightly more memory-efficient.

References

source
SparseMatrixColorings.PerfectEliminationOrder — Type
PerfectEliminationOrder(elimination_algorithm=CliqueTrees.MCS())

Instance of AbstractOrder which computes a perfect elimination ordering when the underlying graph is chordal. For non-chordal graphs, it computes a suboptimal ordering.

The elimination_algorithm must be an instance of CliqueTrees.EliminationAlgorithm.

Warning

This order can only be applied for symmetric or bidirectional coloring problems. Furthermore, its theoretical guarantees only hold for decompression by substitution.

Danger

This order is implemented as a package extension and requires loading CliqueTrees.jl.

References

source