Coloring in JuMP with start values

Specific Domains Optimization (Mathematical)
, ,
1

Hi there,
I am tinkering with an optimal coloring algorithm for sparse autodiff. Here is the most basic constraint programming formulation I could find for partial distance-2 column coloring.

Problem statement: color all columns of a matrix such that if two columns have a non-zero coefficient in the same row, they have different colors.

using JuMP
using HiGHS
using SparseArrays
using Random

function partial_distance2_column_coloring(A)
    m, n = size(A)
    model = Model(optimizer_with_attributes(HiGHS.Optimizer, "time_limit" => 1.0))
    # one color per column
    @variable(model, 1 <= color[j=1:n] <= n, Int)
    @variable(model, ncolors, Int)
    # count and minimize number of distinct colors
    @constraint(model, [ncolors; color] in MOI.CountDistinct(n + 1))
    @objective(model, Min, ncolors)
    # start values don't work??
    for j in 1:n
        set_start_value(color[j], j)
    end
    set_start_value(ncolors, n)
    # all neighbors of a row must have distinct colors
    for i in 1:m
        neighbors = [j for j in 1:n if A[i, j]]
        @constraint(model, color[neighbors] in MOI.AllDifferent(length(neighbors)))
    end
    optimize!(model)
    @assert primal_status(model) == MOI.FEASIBLE_POINT
    return round.(Int, value.(color))
end

A = sprand(MersenneTwister(0), Bool, 100, 100, 0.05)  # too large to solve in 1s
partial_distance2_column_coloring(A)

Ideally, I’d like the solver to return a feasible solution no matter when it is cut off.
My problem is that the trivial start values I provide don’t seem to work: HiGHS tell me they are infeasible, even though they assign a different color to each column. As a result, it doesn’t find a feasible solution:

[more log]
Attempting to find feasible solution by solving LP for user-supplied values of discrete variables
Coefficient ranges:
  Matrix [1e+00, 1e+02]
  Cost   [1e+00, 1e+00]
  Bound  [1e+00, 1e+02]
  RHS    [1e+00, 1e+02]
Presolving model
Problem status detected on presolve: Infeasible
[more log]
ERROR: AssertionError: primal_status(model) == MOI.FEASIBLE_POINT

Did I miss something obvious in the constraints? Or do I misunderstand the point of start_value?

2 Likes
2

If I run your problem with smaller matrices, like 40x40, I do get a feasible solution but with a warning that the initial point is infeasible. I also see that for larger matrices, the number of integer variables and non-zero values in the HiGHS formulation quickly explodes, which is probably due to new variables and constraints introduced by MOI.AllDifferent and MOI.CountDistinct. I am not sure how JuMP maps the provided initial points to the new variables created by those bridges.

Perhaps one option would be to write out those constraints by hand and then set the start values yourself so you know they are correct?

1 Like
3

That’s what I tried to express: I do know that these starting values are feasible. When you have a coloring problem and use one color per vertex, you can’t really go wrong. In fact, when I run fix(color[i], i), it works fine.

4

The logging from HiGHS is a bit misleading.

  • You have set the primal value for the variables you have defined.
  • But JuMP reformulates CountDistinct and AllDifferent into a larger MILP reformulation
  • The reformulation introduces additional binary variables and constraints
  • It seems that the new variables have a default start of 0.0
  • Rather than looking for a solution using only the start values you have set, it looks for a solution using all of the discrete start values
  • The reformulated fixed problem is infeasible, even though your full problem with the start you provided is feasible.

As an action item, I’ll see if we can improve how we pass starting solutions to HiGHS.

For your problem, you can just manually do the MILP reformulation that JuMP is doing under the hood. Something like:

function partial_distance2_column_x_coloring_2(A)
    m, n = size(A)
    model = Model(HiGHS.Optimizer)
    @variables(model, begin
        x_color[j in 1:n, c in 1:n], Bin, (start = (j == c))
        x_is_color_chosen[1:n], Bin, (start = 1)
    end)
    @constraints(model, begin
        # Each node had one exactly one x_color
        [j in 1:n], sum(x_color[j, :]) == 1
        # x_is_color_chosen = 1 if any nodes are x_colored
        [c in 1:n], sum(x_color[:, c]) <= n * x_is_color_chosen[c]
        # Only one neighbor can have the same x_color
        [i in 1:m, c in 1:n], sum(x_color[j, c] for j in 1:n if A[i, j]) <= 1
        # Symmetry reduction
        [c in 2:n], x_is_color_chosen[c-1] >= x_is_color_chosen[c]
    end)
    fix(x_color[1, 1], 1) # More symmetry
    @objective(model, Min, sum(x_is_color_chosen))
    optimize!(model)
    @assert primal_status(model) == MOI.FEASIBLE_POINT
    return round.(Int, value.(x_color))
end

partial_distance2_column_x_coloring_2(A)
2 Likes
5

Okay, we can fix this in HiGHS.jl: Use setSparseSolution · Issue #299 · jump-dev/HiGHS.jl · GitHub.

(At some point a new function was added to the C API that didn’t exist when we originally wrote the wrapper.)

2 Likes
6

Thanks for the investigation! I’m not crazy, yay!!!

7

And yes it would be very nice to fix this because ultimately I do want to use a constraint programming solver, not a MILP solver (was just using HiGHS for an MWE). From what I understand, CP is much better for coloring tasks?

8

Fixed by Use Highs_setSparseSolution for primal starts by odow · Pull Request #300 · jump-dev/HiGHS.jl · GitHub just in the process of tagging a new release.

The answer for what works best almost certainly depends on the input data set.

1 Like