Iterated greedy recoloring for row and column colorings

docs/src/dev.md

@@ -23,6 +23,7 @@ SparseMatrixColorings.bidirectional_pattern
 
 ```@docs
 SparseMatrixColorings.partial_distance2_coloring
+SparseMatrixColorings.iterated_greedy_recoloring!
 SparseMatrixColorings.star_coloring
 SparseMatrixColorings.acyclic_coloring
 SparseMatrixColorings.group_by_color

src/coloring.jl

@@ -78,6 +78,66 @@ function partial_distance2_coloring!(
     end
 end
 
+"""
+    iterated_greedy_recoloring!(
+        color::AbstractVector{<:Integer},
+        forbidden_colors::AbstractVector{<:Integer},
+        bg::BipartiteGraph,
+        ::Val{side},
+        iterations::Integer,
+    )
+
+Improve the distance-2 coloring `color` of the given `side` (`1` or `2`) in the bipartite graph `bg` in-place, with `iterations` passes of Culberson's iterated greedy algorithm.
+
+Each pass recolors the vertices greedily as in [`partial_distance2_coloring`](@ref), in an order where the vertices of each color class are consecutive.
+Since the vertices of a color class are independent, such a pass never increases the number of colors.
+The color classes are ordered by decreasing color in odd passes, and by decreasing size in even passes.
+
+# References
+
+> [_Iterated Greedy Graph Coloring_](https://doi.org/10.1090/dimacs/026/12), Culberson and Luo (1996)
+"""
+function iterated_greedy_recoloring!(
+    color::AbstractVector{<:Integer},
+    forbidden_colors::AbstractVector{<:Integer},
+    bg::BipartiteGraph{T},
+    ::Val{side},
+    iterations::Integer,
+) where {T,side}
+    n = length(color)
+    vertices_in_order = Vector{T}(undef, n)
+    class_size = Vector{T}(undef, n)
+    class_start = Vector{T}(undef, n)
+    for iteration in 1:iterations
+        k = maximum(color; init=zero(eltype(color)))
+        # order of the color classes
+        class_size[1:k] .= 0
+        for v in eachindex(color)
+            class_size[color[v]] += 1
+        end
+        classes = if isodd(iteration)
+            k:-1:1
+        else
+            sortperm(view(class_size, 1:k); rev=true)
+        end
+        # vertices sorted by class, in the order of the classes (counting sort)
+        position = 1
+        for c in classes
+            class_start[c] = position
+            position += class_size[c]
+        end
+        for v in eachindex(color)
+            c = color[v]
+            vertices_in_order[class_start[c]] = v
+            class_start[c] += 1
+        end
+        partial_distance2_coloring!(
+            color, forbidden_colors, bg, Val(side), vertices_in_order
+        )
+    end
+    return color
+end
+
 """
     star_coloring(
         g::AdjacencyGraph, vertices_in_order::AbstractVector, postprocessing::Bool;

src/interface.jl

@@ -69,12 +69,13 @@ It is passed as an argument to the main function [`coloring`](@ref).
 
 # Constructors
 
-    GreedyColoringAlgorithm{decompression}(order=NaturalOrder(); postprocessing=false, postprocessing_minimizes=:all_colors)
-    GreedyColoringAlgorithm(order=NaturalOrder(); postprocessing=false, postprocessing_minimizes=:all_colors, decompression=:direct)
+    GreedyColoringAlgorithm{decompression}(order=NaturalOrder(); postprocessing=false, postprocessing_minimizes=:all_colors, recoloring_iterations=0)
+    GreedyColoringAlgorithm(order=NaturalOrder(); postprocessing=false, postprocessing_minimizes=:all_colors, decompression=:direct, recoloring_iterations=0)
 
 - `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. This option does not affect row or column colorings.
 - `postprocessing_minimizes::Symbol`: which number of distinct colors is heuristically minimized by postprocessing, either `:all_colors`, `:row_colors` or `:column_colors`. This option only affects bidirectional colorings.
+- `recoloring_iterations::Integer`: number of passes of Culberson's iterated greedy recoloring applied to the coloring of each order (see `SparseMatrixColorings.iterated_greedy_recoloring!`). Each pass costs about as much as the initial greedy coloring and never increases the number of colors. This option only affects row and column colorings.
 - `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
@@ -100,20 +101,25 @@ struct GreedyColoringAlgorithm{decompression,N,O<:NTuple{N,AbstractOrder}} <:
     orders::O
     postprocessing::Bool
     postprocessing_minimizes::Symbol
+    recoloring_iterations::Int
 
     function GreedyColoringAlgorithm{decompression}(
         order_or_orders::Union{AbstractOrder,Tuple}=NaturalOrder();
         postprocessing::Bool=false,
         postprocessing_minimizes::Symbol=:all_colors,
+        recoloring_iterations::Integer=0,
     ) where {decompression}
         check_valid_algorithm(decompression)
+        if recoloring_iterations < 0
+            throw(ArgumentError("`recoloring_iterations` must be nonnegative"))
+        end
         if order_or_orders isa AbstractOrder
             orders = (order_or_orders,)
         else
             orders = order_or_orders
         end
         return new{decompression,length(orders),typeof(orders)}(
-            orders, postprocessing, postprocessing_minimizes
+            orders, postprocessing, postprocessing_minimizes, recoloring_iterations
         )
     end
 end
@@ -123,9 +129,10 @@ function GreedyColoringAlgorithm(
     postprocessing::Bool=false,
     decompression::Symbol=:direct,
     postprocessing_minimizes::Symbol=:all_colors,
+    recoloring_iterations::Integer=0,
 )
     return GreedyColoringAlgorithm{decompression}(
-        order_or_orders; postprocessing, postprocessing_minimizes
+        order_or_orders; postprocessing, postprocessing_minimizes, recoloring_iterations
     )
 end
 
@@ -247,7 +254,13 @@ function _coloring(
     bg = BipartiteGraph(A; symmetric_pattern)
     color_by_order = map(algo.orders) do order
         vertices_in_order = vertices(bg, Val(2), order)
-        return partial_distance2_coloring(bg, Val(2), vertices_in_order; forced_colors)
+        _color = partial_distance2_coloring(bg, Val(2), vertices_in_order; forced_colors)
+        if isnothing(forced_colors) && algo.recoloring_iterations > 0
+            iterated_greedy_recoloring!(
+                _color, similar(_color), bg, Val(2), algo.recoloring_iterations
+            )
+        end
+        return _color
     end
     color = argmin(maximum, color_by_order)
     if speed_setting isa WithResult
@@ -270,7 +283,13 @@ function _coloring(
     bg = BipartiteGraph(A; symmetric_pattern)
     color_by_order = map(algo.orders) do order
         vertices_in_order = vertices(bg, Val(1), order)
-        return partial_distance2_coloring(bg, Val(1), vertices_in_order; forced_colors)
+        _color = partial_distance2_coloring(bg, Val(1), vertices_in_order; forced_colors)
+        if isnothing(forced_colors) && algo.recoloring_iterations > 0
+            iterated_greedy_recoloring!(
+                _color, similar(_color), bg, Val(1), algo.recoloring_iterations
+            )
+        end
+        return _color
     end
     color = argmin(maximum, color_by_order)
     if speed_setting isa WithResult

test/order.jl

@@ -240,3 +240,29 @@ end
             ncolors(coloring(A, problem, algo))
     end
 end
+
+@testset "Iterated greedy recoloring" begin
+    rng = StableRNG(0)
+    @test_throws ArgumentError GreedyColoringAlgorithm(; recoloring_iterations=-1)
+    @testset "$partition" for partition in (:column, :row)
+        problem = ColoringProblem(; structure=:nonsymmetric, partition)
+        algo = GreedyColoringAlgorithm()
+        @testset "$((; m, n, p))" for (m, n, p) in ((100, 200, 0.02), (300, 300, 0.01))
+            A = sprand(rng, Bool, m, n, p)
+            k0 = ncolors(coloring(A, problem, algo))
+            k = map((1, 10, 50)) do recoloring_iterations
+                return ncolors(
+                    coloring(A, problem, GreedyColoringAlgorithm(; recoloring_iterations))
+                )
+            end
+            # a pass never increases the number of colors
+            @test k0 >= k[1] >= k[2] >= k[3]
+        end
+    end
+    # random pattern for which the recoloring helps
+    A = sprand(StableRNG(0), Bool, 2000, 2000, 0.003)
+    problem = ColoringProblem()
+    k0 = ncolors(coloring(A, problem, GreedyColoringAlgorithm()))
+    k = ncolors(coloring(A, problem, GreedyColoringAlgorithm(; recoloring_iterations=50)))
+    @test k < k0
+end

test/random.jl

@@ -49,6 +49,18 @@ end;
     end
 end;
 
+@testset "Column and row coloring with recoloring & decompression" begin
+    algo = GreedyColoringAlgorithm(; decompression=:direct, recoloring_iterations=10)
+    @testset "$partition" for partition in (:column, :row)
+        problem = ColoringProblem(; structure=:nonsymmetric, partition)
+        @testset "$((; m, n, p))" for (m, n, p) in asymmetric_params
+            # separate RNG to keep the instances of the other test sets
+            A0 = sprand(StableRNG(m + n), m, n, p)
+            test_coloring_decompression(A0, problem, algo)
+        end
+    end
+end;
+
 @testset "Symmetric coloring & direct decompression" begin
     problem = ColoringProblem(; structure=:symmetric, partition=:column)
     @testset for algo in (

test/type_stability.jl

@@ -26,6 +26,8 @@ end
     @testset "ADTypes" begin
         @test_opt column_coloring(A, GreedyColoringAlgorithm())
         @test_opt row_coloring(A, GreedyColoringAlgorithm())
+        @test_opt column_coloring(A, GreedyColoringAlgorithm(; recoloring_iterations=5))
+        @test_opt row_coloring(A, GreedyColoringAlgorithm(; recoloring_iterations=5))
         @test_opt symmetric_coloring(Symmetric(A), GreedyColoringAlgorithm())
 
         @inferred column_coloring(A, GreedyColoringAlgorithm())