Revisit the post-processing

[amontoison]

close #264

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 100.00%. Comparing base (a314d12) to head (262a09a).

Additional details and impacted files

@@            Coverage Diff             @@
##              main      #267    +/-   ##
==========================================
  Coverage   100.00%   100.00%            
==========================================
  Files           20        20            
  Lines         2027      2128   +101     
==========================================
+ Hits          2027      2128   +101     

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

[gdalle]

Also, how does this all generalize to the acyclic bicoloring?

[amontoison]

Also, how does this all generalize to the acyclic bicoloring?

It is exactly the same in terms of handling trivial edges.
The only difference is that we need to check the non-leaf nodes in non-trivial trees and mark them as needed for decompression, but this process is fully deterministic, just like checking hubs in non-trivial stars.

It was just late yesterday, and I didn’t find the energy to finalize the PR.
I will continue it tonight.

[amontoison]

@gdalle I have an issue in my new function postprocess_star_bicoloring, I need to check if the hub h is in the column partition (1 <= h <= n) or the row partition (n+1 <= h <= n+m) but I didn't find a way to obtain n or m from the arguments.
Because we pass the AdjacencyGraph of the augmented matrix, we can only recover n+m.

[gdalle]

We could store this information in the adjacency graph, to know whether it comes from a symmetric matrix or not, and in the latter case have access to the row/column split

[gdalle]

Or we pass an additional argument to the coloring

[amontoison]

We could store this information in the adjacency graph, to know whether it comes from a symmetric matrix or not, and in the latter case have access to the row/column split

I prefer this approach because we can know that we have "augmented adjacency matrix" this way, and thus know that we are in bicoloring context for the function postprocess!.

I will not need to pass the boolean parameter bicoloring everywhere.

[gdalle]

I suggest replacing AdjacencyGraph with:

struct AdjacencyGraph{T<:Integer,bicoloring}
    S::SparsityPatternCSC{T}
    edge_to_index::Vector{T}
    original_size::Tuple{Int,Int}
end

We no longer need has_diagonal because we have bicoloring => !has_diagonal

[amontoison]

I suggest replacing AdjacencyGraph with:

struct AdjacencyGraph{T<:Integer,bicoloring}
    S::SparsityPatternCSC{T}
    edge_to_index::Vector{T}
    original_size::Tuple{Int,Int}
end

Ok for me 👍

We no longer need has_diagonal because we have bicoloring => !has_diagonal

Yes but we don't have !has_diagonal => bicoloring. Is it fine for you that the user can't specify anymore has_diagonal = false in a non-bicoloring case?

I think adding a field bicoloring::Bool in AdjacencyGraph can handle more corner cases.

[gdalle]

Yes but we don't have !has_diagonal => bicoloring. Is it fine for you that the user can't specify anymore has_diagonal = false in a non-bicoloring case?

has_diagonal was never part of our public API to begin with, so the user never gets to specify it. I think the cases where someone wants to color a Hessian whose entire diagonal is zero are sufficiently rare for us to ignore them.
In addition, they would lead to type instabilities because they require a check on the values of the matrix before creating the AdjacencyGraph object.

[amontoison]

has_diagonal was never part of our public API to begin with, so the user never gets to specify it. I think the cases where someone wants to color a Hessian whose entire diagonal is zero are sufficiently rare for us to ignore them.

Fair point 👍
I updated the PR.

[gdalle]

Aside from my remarks above, do you have any ideas for test cases where our heuristics actually produce the result we want?

[amontoison]

Aside from my remarks above, do you have any ideas for test cases where our heuristics actually produce the result we want?

@gdalle I added the heuristic that I quickly described on Messenger in the following commit: 5fa3f7a
I only need to fuse the post-processing functions together and it should be good.

[amontoison]

@gdalle It is finally ready for a review!!

[amontoison]

@gdalle I don't know what to do for the code coverage of the last two missing lines.
One solution is to run many tests with the acyclic coloring + random order and find one case for the arrowhead matrix that works.

[gdalle]

Isn't it a bit worrying that across all of our randomized tests we don't hit this case even once? I think that should tell us to pause and think about why

[amontoison]

I did some tests and it is because he always got a par for the number columns / rows to neutralize per color.
Using more rectangular random problems could help.

It is harder to get trivial edges in acyclic coloring than star coloring.

[amontoison]

@gdalle I was able to run the code in the two lines not covered by the unit tests with a random order and the following matrices:

6×12 SparseMatrixCSC{Float64, Int64} with 11 stored entries:
  ⋅   1.0   ⋅    ⋅    ⋅   1.0   ⋅    ⋅   1.0   ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅ 
  ⋅   1.0  1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅   1.0
  ⋅    ⋅    ⋅   1.0   ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅ 

6×12 SparseMatrixCSC{Float64, Int64} with 11 stored entries:
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅   1.0  1.0   ⋅ 
  ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅ 
  ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅   1.0   ⋅   1.0   ⋅   1.0   ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅

3×6 SparseMatrixCSC{Float64, Int64} with 6 stored entries:
 1.0   ⋅   1.0  1.0   ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅   1.0   ⋅ 
  ⋅    ⋅    ⋅   1.0   ⋅   1.0

3×6 SparseMatrixCSC{Float64, Int64} with 5 stored entries:
 1.0   ⋅    ⋅    ⋅    ⋅    ⋅ 
  ⋅   1.0  1.0   ⋅    ⋅   1.0
  ⋅   1.0   ⋅    ⋅    ⋅    ⋅ 

3×6 SparseMatrixCSC{Float64, Int64} with 7 stored entries:
  ⋅    ⋅    ⋅    ⋅   1.0   ⋅ 
  ⋅    ⋅   1.0   ⋅    ⋅   1.0
 1.0  1.0  1.0   ⋅   1.0   ⋅ 

3×5 SparseMatrixCSC{Float64, Int64} with 5 stored entries:
  ⋅   1.0  1.0  1.0   ⋅ 
 1.0   ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅   1.0   ⋅    ⋅ 

3×4 SparseMatrixCSC{Float64, Int64} with 5 stored entries:
  ⋅   1.0  1.0   ⋅ 
  ⋅    ⋅    ⋅   1.0
 1.0   ⋅   1.0   ⋅

3×4 SparseMatrixCSC{Float64, Int64} with 6 stored entries:
  ⋅    ⋅    ⋅   1.0
 1.0  1.0   ⋅    ⋅ 
 1.0   ⋅   1.0  1.0

3×4 SparseMatrixCSC{Float64, Int64} with 4 stored entries:
  ⋅    ⋅   1.0  1.0
  ⋅    ⋅   1.0   ⋅ 
 1.0   ⋅    ⋅    ⋅ 

3×4 SparseMatrixCSC{Float64, Int64} with 6 stored entries:
 1.0  1.0   ⋅   1.0
 1.0   ⋅    ⋅    ⋅ 
  ⋅   1.0  1.0   ⋅ 

3×3 SparseMatrixCSC{Float64, Int64} with 4 stored entries:
  ⋅   1.0   ⋅ 
 1.0   ⋅   1.0
  ⋅    ⋅   1.0

3×3 SparseMatrixCSC{Float64, Int64} with 4 stored entries:
 1.0   ⋅   1.0
 1.0   ⋅    ⋅ 
  ⋅   1.0   ⋅ 

3×3 SparseMatrixCSC{Float64, Int64} with 4 stored entries:
 1.0  1.0   ⋅ 
  ⋅    ⋅   1.0
  ⋅   1.0   ⋅

I will add the 3x3 matrix in the tests and we should be good.

[gdalle]

Thanks!
I'm still on vacation but I'll review when I get back to work

[amontoison]

For after the holidays 🏖️ 🎅

I followed your suggestion Guillaume, and tried to check the impact of postprocessing on a collection of “real” problems.
I downloaded all 846 matrices from the SuiteSparse Matrix Collection where the number of rows and columns is between 1000 and 10000, and I checked the number of problems on which postprocessing results in a gain in total colors.

It also helped me detect that I had forgotten to split unused row and column colors at the end of postprocessing (for the bicoloring), which drastically changed the results. See this commit: e999c65.
Initially, I got worse results than the dummy heuristic we currently have on main!

Before the PR (on the branch main), the results are:

star_bicoloring | random | 18 / 846
star_bicoloring | natural | 794 / 846
star_bicoloring | largest_first | 532 / 846
star_bicoloring | smallest_last | 508 / 846
star_bicoloring | indidence_degree | 214 / 846
star_bicoloring | dynnamic_largest_first | 305 / 846

acyclic_bicoloring | random | 14 / 846
acyclic_bicoloring | natural | 76 / 846
acyclic_bicoloring | largest_first | 52 / 846
acyclic_bicoloring | smallest_last | 41 / 846
acyclic_bicoloring | indidence_degree | 9 / 846
acyclic_bicoloring | dynnamic_largest_first | 65 / 846

and with this PR:

star_bicoloring | random | :all_colors | 29 / 846
star_bicoloring | random | :row_colors | 29 / 846
star_bicoloring | random | :column_colors | 29 / 846

star_bicoloring | natural | :all_colors | 837 / 846
star_bicoloring | natural | :row_colors | 837 / 846
star_bicoloring | natural | :column_colors | 837 / 846

star_bicoloring | largest_first | :all_colors | 605 / 846
star_bicoloring | largest_first | :row_colors | 605 / 846
star_bicoloring | largest_first | :column_colors | 605 / 846

star_bicoloring | smallest_last | :all_colors | 594 / 846
star_bicoloring | smallest_last | :row_colors | 594 / 846
star_bicoloring | smallest_last | :column_colors | 594 / 846

star_bicoloring | indidence_degree | :all_colors | 523 / 846
star_bicoloring | indidence_degree | :row_colors | 523 / 846
star_bicoloring | indidence_degree | :column_colors | 523 / 846

star_bicoloring | dynamic_largest_first | :all_colors | 415 / 846
star_bicoloring | dynamic_largest_first | :row_colors | 415 / 846
star_bicoloring | dynamic_largest_first | :column_colors | 415 / 846

acyclic_bicoloring | random | :all_colors | 28 / 846
acyclic_bicoloring | random | :row_colors | 28 / 846
acyclic_bicoloring | random | :column_colors | 28 / 846

acyclic_bicoloring | natural | :all_colors | 141 / 846
acyclic_bicoloring | natural | :row_colors | 141 / 846
acyclic_bicoloring | natural | :column_colors | 141 / 846

acyclic_bicoloring | largest_first | :all_colors | 99 / 846
acyclic_bicoloring | largest_first | :row_colors | 99 / 846
acyclic_bicoloring | largest_first | :column_colors | 99 / 846

acyclic_bicoloring | smallest_last | :all_colors | 84 / 846
acyclic_bicoloring | smallest_last | :row_colors | 84 / 846
acyclic_bicoloring | smallest_last | :column_colors | 84 / 846

acyclic_bicoloring | indidence_degree | :all_colors | 68 / 846
acyclic_bicoloring | indidence_degree | :row_colors | 68 / 846
acyclic_bicoloring | indidence_degree | :column_colors | 68 / 846

acyclic_bicoloring | dynamic_largest_first | :all_colors | 197 / 846
acyclic_bicoloring | dynamic_largest_first | :row_colors | 197 / 846
acyclic_bicoloring | dynamic_largest_first | :column_colors | 197 / 846

The number of problems showing an improvement is independent of postprocessing_minimizes (on the subset of 846 matrices).
We can observe that, depending on the ordering, the number of problems showing an improvement after postprocessing varies greatly.
Postprocessing has more impact on star bicoloring than on acyclic bicoloring.
I explain this by the fact that acyclic bicoloring allows more complex graph structures (trees) than star bicoloring (stars), which probably limits the number of trivial edges.
I will perform a similar run for symmetric matrices to compare star coloring and acyclic coloring.

[amontoison]

The code that I used is:

using SparseMatrixColorings
using SuiteSparseMatrixCollection
using MatrixMarket
using StableRNGs

ssmc = ssmc_db()
matrices = ssmc[(1000 .≤ ssmc.nrows) .& (ssmc.nrows .≤ 10000) .& (1000 .≤ ssmc.ncols) .& (ssmc.ncols .≤ 10000), :]
paths = fetch_ssmc(matrices, format="MM")
nmatrices = length(paths)

stats = Dict{String, Dict{Symbol, Dict{String, Dict{Symbol, Tuple{Int,Int,Int}}}}}()
verbose = false
legacy = false

decompressions = [("star_bicoloring", :direct),
                  ("acyclic_bicoloring", :substitution)]

orders = [("random", RandomOrder(StableRNG(0), 0)),
          ("natural", NaturalOrder()),
          ("largest_first", LargestFirst()),
          ("smallest_last", SmallestLast()),
          ("indidence_degree", IncidenceDegree()),
          ("dynamic_largest_first" , DynamicLargestFirst())]

variants_pp = legacy ? (:none, :default) : (:none, :all_colors, :row_colors, :column_colors)

for i = 1:nmatrices
  name = matrices.name[i]
  stats[name] = Dict{Symbol, Dict{String, Dict{Symbol, Tuple{Int,Int,Int}}}}()
  path = joinpath(paths[i], "$(name).mtx")
  A = MatrixMarket.mmread(path) 
  for decompression in (:direct, :substitution)
    stats[name][decompression] = Dict{String, Dict{Symbol, Tuple{Int,Int,Int}}}()
    for (ordering, order) in orders
      stats[name][decompression][ordering] = Dict{Symbol, Tuple{Int,Int,Int}}()
      for postprocessing_minimizes in variants_pp
        postprocessing = (postprocessing_minimizes == :none) ? false : true
        verbose && println("matrix = $name | decompression = $decompression | order = $ordering | postprocessing = $postprocessing")
        problem = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional)
        if legacy
          algo = GreedyColoringAlgorithm(order; decompression, postprocessing)
        else
          algo = GreedyColoringAlgorithm(order; decompression, postprocessing, postprocessing_minimizes)
        end
        result = coloring(A, problem, algo)
        ncolors_row = row_groups(result) |> length
        ncolors_col = column_groups(result) |> length
        verbose && println("rcolors: $ncolors_row | ccolors: $ncolors_col")
        stats[name][decompression][ordering][postprocessing_minimizes] = (ncolors_row, ncolors_col, ncolors_row + ncolors_col)
      end
      verbose && println()
    end
    verbose && println("-----------------------------------------------")
  end
end

if legacy
  for (bicoloring, decompression) in decompressions
    for (ordering, order) in orders
      gain = 0
      for i = 1:nmatrices
        name = matrices.name[i]
        (r0,c0,s0) = stats[name][decompression][ordering][:none]
        (r1,c1,s1) = stats[name][decompression][ordering][:default]
        (s1 < s0) && (gain += 1)
      end
      println("$bicoloring | $ordering | $gain / $nmatrices")
    end
  end
else
  for (bicoloring, decompression) in decompressions
    for (ordering, order) in orders
      gain1 = 0  # all_colors
      gain2 = 0  # row_colors
      gain3 = 0  # column_colors
      for i = 1:nmatrices
        name = matrices.name[i]
        (r0,c0,s0) = stats[name][decompression][ordering][:none]
        (r1,c1,s1) = stats[name][decompression][ordering][:all_colors]
        (r2,c2,s2) = stats[name][decompression][ordering][:row_colors]
        (r3,c3,s3) = stats[name][decompression][ordering][:column_colors]
        (s1 < s0) && (gain1 += 1)
        (s2 < s0) && (gain2 += 1)
        (s3 < s0) && (gain3 += 1) 
      end
      println("$bicoloring | $ordering | :all_colors | $gain1 / $nmatrices")
      println("$bicoloring | $ordering | :row_colors | $gain2 / $nmatrices")
      println("$bicoloring | $ordering | :column_colors | $gain3 / $nmatrices")
    end
  end
end

The variable legacy is used to specify if we run the code on the branch main or not.

[gdalle]

Thanks for running the tests! If I understand correctly, the takeaways are:

1. The new postprocessing reduces the number of colors for a larger set of instances
2. It is unclear whether the postprocessing_minimizes argument actually changes anything
3. The order has a big impact, which to me suggests that we should try using it not just when we initially assign colors, but also when we go back through the list of vertices to mark them as essential or not (e.g. here)

[amontoison]

1. The new postprocessing reduces the number of colors for a larger set of instances

Yes, but only for bicoloring, I didn't see any difference with the new post-processing for star and acyclic colorings.

2. It is unclear whether the postprocessing_minimizes argument actually changes anything

The total number of colors doesn't change with postprocessing_minimizes but we remove more row or column colors depending on it.

3. The order has a big impact, which to me suggests that we should try using it not just when we initially assign colors, but also when we go back through the list of vertices to mark them as essential or not (e.g. here)

It doesn't change anything at the line that you specified.
Whatever is the order of the vertices, we just mark necessary or not (independent or the ordering).

What is important is the order of treatment of the trivial trees or stars (last step of post-processing).

[gdalle]

What is important is the order of treatment of the trivial trees or stars (last step of post-processing).

Sorry, yes, replace the line I pointed to with e.g. these ones. Perhaps there we could permute stars or trees according to the same order used for coloring, or its reverse order?

[amontoison]

What is important is the order of treatment of the trivial trees or stars (last step of post-processing).

Sorry, yes, replace the line I pointed to with e.g. these ones. Perhaps there we could permute stars or trees according to the same order used for coloring, or its reverse order?

For the acyclic coloring, we store ane iterate over the list of trees, so I don't have any idea about how we could use the ordering of the vertices for that.

For the star coloring, we iterate over the stars multiple times. We iterate over the star of each vertex and not the stars directly (we can't because we only encode the hub and not the spokes in our implementation of the stars).
But what we want in pratice is to iterate over the stars and I don't know how we can permute them based on the ordering of the vertices.

It is an opened question and the treatment of the trivial stars / trees is a standalone new graph problem from what I understand. We have something to dig here if we want.

We should not forget to explain it in our paper when we will get the review.

[amontoison]

@gdalle Can we merge this PR?

[gdalle]

I'll take a cursory look today or tomorrow

[gdalle]

(Check your DMs)