Chordal coloring algorithm.

[samuelsonric]

Chordal graphs can be minimally colored in linear time, and every such coloring is acyclic. This paper discusses chordal graphs in §4 and suggests that they arise in practice. This patch adds an algorithm ChordalColoringAlgorithm that handles this case. Here are benchmarks for acyclic coloring of a band graph (defined in the paper).

using BandedMatrices, BenchmarkTools, Random, SparseArrays, SparseMatrixColorings

n = 10000 # vertices
m = 400   # bandwidth
matrix = sparse(brand(n, n, div(m, 2), 0))
matrix .+= matrix'
perm = randperm(n)
permute!(matrix, perm, perm)
problem = ColoringProblem(; structure=:symmetric, partition=:column)

greedy coloring

julia> algorithm = GreedyColoringAlgorithm(NaturalOrder())
GreedyColoringAlgorithm{:direct, NaturalOrder}(NaturalOrder(), false)

julia> ncolors(coloring(matrix, problem, algorithm))
487

julia> @benchmark coloring(matrix, problem, algorithm)
BenchmarkTools.Trial: 1 sample with 1 evaluation per sample.
 Single result which took 21.141 s (0.47% GC) to evaluate,
 with a memory estimate of 329.99 MiB, over 3668722 allocations.

chordal coloring

julia> algorithm = ChordalColoringAlgorithm()
ChordalColoringAlgorithm()

julia> ncolors(coloring(matrix, problem, algorithm))
201

julia> @benchmark coloring(matrix, problem, algorithm)
BenchmarkTools.Trial: 122 samples with 1 evaluation per sample.
 Range (min … max):  39.282 ms … 44.836 ms  ┊ GC (min … max): 0.00% … 1.65%
 Time  (median):     40.740 ms              ┊ GC (median):    1.13%
 Time  (mean ± σ):   40.990 ms ±  1.149 ms  ┊ GC (mean ± σ):  1.13% ± 0.83%

         ▂▂▂▆▂▂▂ ▂▆▂ █            ▂        ▄                   
  ██▆██▄▆█████████████▆▄█▆▆▄▆█▆▁▄▁█▆▄▄▆▁▆▆██▄▄▁▄▁▁▁▁▄▁▄▁▁▁▁▁▄ ▄
  39.3 ms         Histogram: frequency by time        44.2 ms <

 Memory estimate: 16.11 MiB, allocs estimate: 49.

The algorithm errors if the graph is not chordal.

julia> invp = invperm(perm);

julia> matrix[invp[1], invp[3]] = matrix[invp[3], invp[1]] = matrix[invp[4], invp[2]] = matrix[invp[2], invp[4]] = 0;

julia> dropzeros!(matrix);

julia> ncolors(coloring(matrix, problem, ChordalColoringAlgorithm()))
ERROR: Matrix does not have a chordal sparsity pattern.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] coloring(A::SparseMatrixCSC{Float64, Int64}, ::ColoringProblem{:symmetric, :column}, algo::ChordalColoringAlgorithm)
   @ SparseMatrixColorings ~/Source/git/SparseMatrixColorings.jl/src/chordal.jl:20
 [3] top-level scope
   @ REPL[18]:1

If you like the algorithm then I can add tests and so forth.

[samuelsonric]

Related?

[samuelsonric]

Oops! I'll rerun & update the benchmarks.

[gdalle]

Hi @samuelsonric, thank you for this suggestion!

Indeed, I had skipped the experimental section when reading Efficient Computation of Sparse Hessians Using Coloring and Automatic Differentiation (Gebremehdin et al., 2008), but such special cases deserve our attention. Banded matrices in particular are very relevant for SciML applications, and they generate are a special case of chordal graphs.

Here are a few remarks:

• Adding a new dependency to SparseMatrixColorings.jl is rather a big deal, since it affects many many users down the line. The only way we would be willing to consider it is in a package extension.
• New coloring algorithms which are not based on two-colored structures will not be compatible with our current decompression utilities. To test your code, we would need to solve Decompression from vector of colors? #67 first.
• From a user perspective, it might be rather confusing to see the number of coloring algorithms grow to include one for several families of graphs / matrices. What I had in mind to solve Optimized coloring for specific matrix types #65 was something like an OptimalColoringAlgorithm, which would work based on dispatch to associate matrices with known optimal colorings (like BandedMatrix from BandedMatrices.jl). The problem is that not every chordal graph can be detected just by dispatch on structured matrix families.
• Another alternative would be to encode the chordal aspect in the vertex ordering subroutine, instead of the global coloring algorithm. From what I understand, optimal coloring of a chordal graph can be achieved with a perfect elimination order. Is that what you're computing with CliqueTrees.permutation?

[samuelsonric]

Yes: CliqueTrees.permutation computes a peo!

[gdalle]

Yes: CliqueTrees.permutation computes a peo!

I had the vague notion that this was an NP-hard problem?

[samuelsonric]

Ok: I think a package extension makes sense. A chordal routine could easily be integrated into an OptimalColoringAlgorithm. Although chordality cannot be verified at compile time, it can be verified almost instantaneously at runtime using CliqueTrees.ischordal or CliqueTrees.isperfect.

[samuelsonric]

I had the vague notion that this was an NP-hard problem?

This is a bit subtle. CliqueTrees.permutation computes a vertex ordering in linear time. The ordering is a PEO if and only if the graph is chordal. Furthermore, whether the ordering is in fact perfect can also be determined in linear time.

This is the standard strategy for chordality testing. Computing a PEO for general graphs is NP hard!

[gdalle]

My suggestion would be to

• define and document a new struct PerfectEliminationOrder end in src/order.jl
• implement SparseMatrixColoring.vertices(g, ::PerfectEliminationOrder) in ext/SparseMatrixColoringsCliqueTrees.jl
• add tests in test/chordal.jl

This seems like a minimum-effort way to add this functionality without fully solving #67, because it would be inserted within the existing GreedyColoringAlgorithm. What do you think?

[samuelsonric]

Added a PerfectEliminationOrder algorithm. However, the number of colors generated by the greedy ordering algorithm seems to consistently be twice the chromatic number. Is this expected?

using BandedMatrices, BenchmarkTools, Random, SparseArrays, SparseMatrixColorings
n = 50 # vertices
m = 10 # bandwidth
matrix = sparse(brand(n, n, div(m, 2), 0))
matrix .+= matrix'
perm = randperm(n)
permute!(matrix, perm, perm)
problem = ColoringProblem(; structure=:symmetric, partition=:column)

LargestFirst

julia> algorithm = GreedyColoringAlgorithm(LargestFirst());

julia> ncolors(coloring(matrix, problem, algorithm))
13

PerfectEliminationOrder

julia> import CliqueTrees

julia> algorithm = GreedyColoringAlgorithm(PerfectEliminationOrder());

julia> ncolors(coloring(matrix, problem, algorithm))
11

julia> using SimpleGraphs, SimpleGraphAlgorithms

julia> chromatic_number(UndirectedGraph(matrix))
6

[samuelsonric]

Maybe I selected the wrong problem type?

decompression=:substitution got it

[samuelsonric]

Added tests!

[gdalle]

Thanks for this initiative! Here are a few remarks of mine ;)

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 100.00%. Comparing base (3cdcf96) to head (38ba304).
Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff            @@
##              main      #202   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files           13        14    +1     
  Lines         1628      1635    +7     
=========================================
+ Hits          1628      1635    +7     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[samuelsonric]

Should be everything.

[samuelsonric]

The docs should build successfully, now!

[gdalle]

Almost there!

[gdalle]

Wait, on second thought I'm not sure the order-based approach is right. It would get us an optimal distance-1 coloring if we applied a greedy distance-1 algorithm after, but we actually apply a greedy acyclic algorithm after. I struggle to see whether we have any optimality guarantees that way.

EDIT: actually I think it's fine, because Lemma 4.2 of the 2008 paper states

If a mapping is a distance-1 coloring for a chordal graph G, then is also an acyclic coloring for G.

[amontoison]

Thanks @samuelsonric for this PR!
I think we should add a small example in the documentation on how to use this new ordering with a banded matrix.