gdalle / sparsematrixcolorings.jl Goto Github PK

The following test is skipped because it fails at the moment:

I think that we should also export these routines.

The revamp in #38 means we no longer support plain color vectors for decompression. This makes DI incompatible with any coloring algorithm that doesn't come from here, and thus with the ADTypes coloring interface.

In interface.jl we have:

function coloring(
    S::AbstractMatrix,
    ::ColoringProblem{:symmetric,:column,:direct},
    algo::GreedyColoringAlgorithm,
)
    ag = adjacency_graph(S)
    color, star_set = star_coloring(ag, algo.order)
    # TODO: handle star_set
    return DefaultColoringResult{:symmetric,:column,:direct}(S, color)
end

function coloring(
    S::AbstractMatrix,
    ::ColoringProblem{:symmetric,:column,:substitution},
    algo::GreedyColoringAlgorithm,
)
    ag = adjacency_graph(S)
    color, tree_set = acyclic_coloring(ag, algo.order)
    # TODO: handle tree_set
    return DefaultColoringResult{:symmetric,:column,:substitution}(S, color)
end

Should we also add:

function coloring(
    S::AbstractMatrix,
    ::ColoringProblem{:symmetric,:row,:direct},
    algo::GreedyColoringAlgorithm,
)
    ag = adjacency_graph(S)
    color, star_set = star_coloring(ag, algo.order)
    # TODO: handle star_set
    return DefaultColoringResult{:symmetric,:row,:direct}(S, color)
end

function coloring(
    S::AbstractMatrix,
    ::ColoringProblem{:symmetric,:row,:substitution},
    algo::GreedyColoringAlgorithm,
)
    ag = adjacency_graph(S)
    color, tree_set = acyclic_coloring(ag, algo.order)
    # TODO: handle tree_set
    return DefaultColoringResult{:symmetric,:row,:substitution}(S, color)
end

Alternative solution:
Should we add a :row_or_column entry to the partition, and only allow it if structure == :symmetric?

I don't know where to find the matrices used in Table 7.1 of "New Acyclic and Star Coloring Algorithms"

To compare with ColPack (https://dl.acm.org/doi/10.1145/2513109.2513110) we need the following orders:

• Natural order (N)
• Random order (R)
• Largest first (LF)
• Smallest last (SL)
• Incidence degree (ID)
• Dynamic largest first (DLF)
• Saturation degree (SD)

#38 is a great PR (congrats to the guy who wrote it) but it leaves aside the question of decompressing when all we have is the StarSet and not a SparseMatrixCSC.
Opening this issue to keep track once I merge

In DifferentiationInterface, the structure and partition can be deduced from the context, but the decompression method (either :direct or :substitution) should be configurable by the user.

In AutoSparse backend you can give a ColoringAlgorithm.

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

There are some matrix types for which we can directly compute the optimal coloring, like (block-)banded matrices.

ArrayInterface.jl does this already but only for column coloring. We could reuse some of that code (if the license is right) and generalize it to other colorings (row, symmetric, bidirectional).

Related:

• #64

This would allow processing B one column or row at a time.

First implementation in #26

Will never be used in practice with GreedyColoringAlgorithm, but one could imagine optimal colorings with JuMP, which would only output a vector of colors.

#71 contains a LinearSystemColoringResult, which is part of the answer in the symmetric case.

At the moment, ADTypes specifies that the result of coloring is just a vector of integer colors.
However, in the case of symmetric coloring, we can compute other things (a star set) that speed up decompression.
More generally, a coloring algorithm could return an AbstractColoringResult, which is then used in a decompression function whose API we need to hash out. Perhaps #35 can help?

At the moment I implement DirectRecover1 from https://www.cs.purdue.edu/homes/apothen/Papers/hessian-coloring-joc.pdf, but DirectRecover2 is more efficient

Once we have solved #44

It is a linear operation so yes. However, it is a linear operation on vectors of values, not matrices.

For column decompression, let $A \in \mathbb{R}^{m \times n}$ be the original matrix and $B \in \mathbb{R}^{m \times c}$ be the compressed one.
Then there exists a matrix $C \in \mathbb{R}^{mn \times mc}$ such that $\mathrm{vec}(A) = C \mathrm{vec}(B)$.
We can further replace $\mathrm{vec}(A)$ with $\mathrm{nonzeros(A)}$ and then $C$ is block-diagonal with exactly one $1$ per row and at most one $1$ per column.

Example: if

$$ A = \begin{pmatrix} a_{11} & 0 & 0 & a_{14} \\ 0 & a_{22} & 0 & a_{24} \\ a_{31} & 0 & a_{33} & 0 \\ a_{41} & 0 & 0 & 0 \end{pmatrix} \quad \text{and} \quad B = \begin{pmatrix} a_{11} & a_{14} \\ a_{22} & a_{24} \\ a_{31} & a_{33} \\ a_{41} & 0 \end{pmatrix} $$

then we have

$$ \begin{pmatrix} a_{11} \\ a_{31} \\ a_{41} \\ a_{22} \\ a_{14} \\ a_{24} \\ a_{33} \end{pmatrix} = \begin{pmatrix} 1 & . & . & . & . & . & . & . \\ . & . & 1 & . & . & . & . & . \\ . & . & . & 1 & . & . & . & . \\ . & 1 & . & . & . & . & . & . \\ . & . & . & . & . & . & 1 & . \\ . & . & . & . & 1 & . & . & . \\ . & . & . & . & . & 1 & . & . \\ \end{pmatrix} \begin{pmatrix} b_{11} \\ b_{12} \\ b_{13} \\ b_{14} \\ b_{21} \\ b_{22} \\ b_{23} \\ b_{24} \end{pmatrix} $$

Is this useful in practice? How does it generalize

• to symmetric matrices?
• to substitution methods where a linear system must be solved?
• to matrices stored as non-sparse formats (e.g. tridiagonal)?

In all cases we have linear functions.

• In practice we only store one triangle of the Hessian of the Lagrangian.
  We could an option as one of the following symbols :full, :lower or :upper.

I compute it anyway when creating the graph

At the moment, our only optimized decompression is for SparseMatrixCSC in :direct mode: we store a vector of compressed_indices such that nonzeros(A) = vec(B)[compressed_indices].
We can probably find a similar optimization for :substitution mode.

What do we want to do for other matrix types, like:

• the ones from LinearAlgebra: Bidiagonal, Tridiagonal, Symmetric, etc.
• BandedMatrices.jl, BlockBandedMatrices.jl

It would be rather tiring to find optimal decompression methods for each of these. My proposal (as a first step) would be to always have a SparseMatrixCSC buffer into which we decompress, and then copy the A_buffer::SparseMatrixCSC into A::SomeWeirdMatrix.
Essentially, it's easier to implement fast copy from SparseMatrixCSC than fast decompression.

Related:

• #65
• #44

This could help taylor decompression preparation. Right now we make the hypothesis that A = similar(S, Float64) or something to that effect.

The coloring algorithms in this package only color rows or columns, but it would be nice to have bicoloring for mixed mode (forward+reverse) Jacobian evaluation.

Related:

• JuliaDiff/SparseDiffTools.jl#25
• JuliaDiff/SparseDiffTools.jl#112

For the following matrix, how can I know if I should use the "compressed columns" in the AD context associated to colors[i] or colors[j] to determine the coefficient H[i,j]?

10×10 SparseMatrixCSC{Float64, Int64} with 58 stored entries:
 2.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  2.0  2.0
 1.0  2.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0  1.0
 1.0   ⋅   2.0   ⋅    ⋅    ⋅    ⋅    ⋅   1.0  1.0
 1.0   ⋅    ⋅   2.0   ⋅    ⋅    ⋅    ⋅   1.0  1.0
 1.0   ⋅    ⋅    ⋅   2.0   ⋅    ⋅    ⋅   1.0  1.0
 1.0   ⋅    ⋅    ⋅    ⋅   2.0   ⋅    ⋅   1.0  1.0
 1.0   ⋅    ⋅    ⋅    ⋅    ⋅   2.0   ⋅   1.0  1.0
 1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   2.0  1.0  1.0
 2.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  2.0  1.0
 2.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  2.0

julia> colors = symmetric_coloring(H, GreedyColoringAlgorithm())
10-element Vector{Int64}:
 1
 2
 2
 2
 2
 2
 2
 2
 3
 4

Right now the problem and algorithm constructors are not type-stable. This must be solved before v0.4

julia> using SparseMatrixColorings, Test

julia> @inferred ColoringProblem(structure=:symmetric, partition=:column)
ERROR: return type ColoringProblem{:symmetric, :column} does not match inferred return type ColoringProblem
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35

julia> @inferred GreedyColoringAlgorithm(decompression=:direct)
ERROR: return type GreedyColoringAlgorithm{:direct, SparseMatrixColorings.NaturalOrder} does not match inferred return type GreedyColoringAlgorithm{_A, SparseMatrixColorings.NaturalOrder} where _A
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35

I didn't implement it because I'm not even sure it is possible. @amontoison what do you think?

What to do when the patterns don't match? Or the matrix types don't match?