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 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 taken from the following articles:

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

Alternatives

Citing

Please cite this software using the Zenodo DOI of the version you used. The link https://zenodo.org/doi/10.5281/zenodo.11314275 resolves to the latest version.

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())
GreedyColoringAlgorithm(order=NaturalOrder(); decompression=:direct)
  • order::AbstractOrder: the order in which the columns or rows are colored, which can impact the number of colors.
  • 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{Int}: 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 colors used to color the matrix.

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

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

Return a vector group such that for every 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 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.

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

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

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.DynamicDegreeBasedOrder — Type
DynamicDegreeBasedOrder{degtype,direction}

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

Danger

This order is still experimental and needs more tests, correctness is not yet guaranteed.

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

References

source