[TambyVanderpooten] round solutions to Int during the algorithm

src/algorithms/TambyVanderpooten.jl

@@ -33,32 +33,42 @@ This algorithm is restricted to problems with:
 struct TambyVanderpooten <: AbstractAlgorithm end
 
 """
-    const _TAMBY_VANDERPOOTEN_ATOL = 1e-3
+    _round_objective(y::Float64, constant::Float64)::Int
 
-This constant is important.
+Return `y - constant`, rounded to the nearest integer.
 
-When we have a solution, we need to decide if the `Y` vector is at the edge of
-one of the bounds. The paper assumes exact arithmetic. Solvers have tolerances.
+The paper assumes exact arithmetic, and the algorithm relies on exact
+comparisons, like whether `y[k] == u[k]` or `y[k] < u[k]`, to decide if a point
+is on the edge of a search zone or strictly inside it. Solvers have tolerances,
+so the same nondominated point can be returned with different noise, for
+example, as `1234.00005` and `1233.99998`. If we compare these values directly,
+the algorithm can revisit the same point and cycle.
 
-In theory the algorithm works only for problems with integer-valued objectives.
-So we could make this as big as `0.5` or something. We previously used `1e-6`,
-but this caused problems when the solver's feasibility tolerance meant that the
-objective value was something like `1234.00005`.
+The algorithm works only for problems with integer-valued objectives (we use an
+ε of `1.0` for the strict inequalities), so we round the non-constant part of
+every objective value to the nearest integer and store it as an `Int`. This
+ensures that all comparisons are exact.
 
-In future, we could consider how to rewrite the algorithm to make it more robust
-to numerical issues like this.
+All objective values in the algorithm (`yI`, `yN`, `U_N`, `V`, and the keys of
+`solutions`) are stored without the objective constant. We still return the
+original floating point objective values to the user. This does not change
+the algorithm because the comparisons are elementwise and `_h` uses only
+differences between points. It also means that the right-hand side of the
+ε-constraints, which have their constant moved to the set, is the stored value.
 """
-const _TAMBY_VANDERPOOTEN_ATOL = 1e-3
+function _round_objective(y::Float64, constant::Float64)
+    return round(Int, y - constant)
+end
 
 # This is Algorithm 1 from the paper.
 function _update_search_region(
-    U_N::Dict{Vector{Float64},Vector{Vector{Vector{Float64}}}},
-    y::Vector{Float64},
-    yN::Vector{Float64},
+    U_N::Dict{Vector{Int},Vector{Vector{Vector{Int}}}},
+    y::Vector{Int},
+    yN::Vector{Int},
 )
     p = length(y)
-    bounds_to_remove = Vector{Float64}[]
-    bounds_to_add = Dict{Vector{Float64},Vector{Vector{Vector{Float64}}}}()
+    bounds_to_remove = Vector{Int}[]
+    bounds_to_add = Dict{Vector{Int},Vector{Vector{Vector{Int}}}}()
     for (u, u_n) in U_N
         if _is_less(y, u, 0) # k=0 here because we want to check every element
             push!(bounds_to_remove, u)
@@ -102,33 +112,34 @@ function _is_less(y, u, k)
     return true
 end
 
-function _get_child(u::Vector{Float64}, y::Vector{Float64}, k::Int)
+function _get_child(u::Vector{Int}, y::Vector{Int}, k::Int)
     @assert length(u) == length(y)
     child = copy(u)
     child[k] = y[k]
     return child
 end
 
-# This is function h(u, k) in the paper.
-function _h(k::Int, u::Vector{Float64}, yI::Vector{Float64})
+# This is function h(u, k) in the paper. We use Float64 because the product can
+# overflow an Int.
+function _h(k::Int, u::Vector{Int}, yI::Vector{Int})
     h = 1.0
     for (i, (u_i, yI_i)) in enumerate(zip(u, yI))
         if i != k
-            h *= u_i - yI_i
+            h *= Float64(u_i - yI_i)
         end
     end
     return h
 end
 
 # This is Problem (6) in the paper. The upper bound M is our nadir point yN.
 function _select_search_zone(
-    U_N::Dict{Vector{Float64},Vector{Vector{Vector{Float64}}}},
-    yI::Vector{Float64},
-    yN::Vector{Float64},
+    U_N::Dict{Vector{Int},Vector{Vector{Vector{Int}}}},
+    yI::Vector{Int},
+    yN::Vector{Int},
 )
     k_star, u_star, v_star = 1, yN, -Inf
     for k in 1:length(yI), u in keys(U_N)
-        if !isapprox(u[k], yN[k]; atol = _TAMBY_VANDERPOOTEN_ATOL)
+        if u[k] != yN[k]
             v = _h(k, u, yI)
             if v > v_star
                 k_star, u_star, v_star = k, u, v
@@ -144,37 +155,39 @@ function minimize_multiobjective!(
     inner::MOI.ModelLike,
     f::MOI.AbstractVectorFunction,
 )
-    solutions = Dict{Vector{Float64},Dict{MOI.VariableIndex,Float64}}()
+    solutions = Dict{Vector{Int},SolutionPoint}()
     status = _minimize_multiobjective!(
         algorithm,
         model,
         inner,
         f,
         solutions,
     )::MOI.TerminationStatusCode
-    return status, SolutionPoint[SolutionPoint(X, Y) for (Y, X) in solutions]
+    return status, collect(values(solutions))
 end
 
 # To ensure the type stability of `.inner` and `.f`, we use a function barrier.
 #
 # We also pass in `solutions` to simplify the function and to ensure that we can
 # return at any point from the this inner function with a partial list of primal
-# solutions.
+# solutions. The keys of `solutions` are the rounded objective values that we
+# use in the algorithm, and the values are the original solutions.
 function _minimize_multiobjective!(
     algorithm::TambyVanderpooten,
     model::Optimizer,
     inner::MOI.ModelLike,
     f::MOI.AbstractVectorFunction,
-    solutions::Dict{Vector{Float64},Dict{MOI.VariableIndex,Float64}},
+    solutions::Dict{Vector{Int},SolutionPoint},
 )
     @assert MOI.get(inner, MOI.ObjectiveSense()) == MOI.MIN_SENSE
     warm_start_supported =
         MOI.supports(inner, MOI.VariablePrimalStart(), MOI.VariableIndex)
     variables = MOI.get(inner, MOI.ListOfVariableIndices())
     n = MOI.output_dimension(f)
     # Compute the ideal (yI) and nadir (yN) points
-    yI, yN = zeros(n), zeros(n)
+    yI, yN = zeros(Int, n), zeros(Int, n)
     scalars = MOI.Utilities.scalarize(f)
+    obj_constants = MOI.constant.(scalars)
     for (i, f_i) in enumerate(scalars)
         MOI.set(inner, MOI.ObjectiveFunction{typeof(f_i)}(), f_i)
         MOI.set(inner, MOI.ObjectiveSense(), MOI.MIN_SENSE)
@@ -185,7 +198,8 @@ function _minimize_multiobjective!(
         end
         _, y_i = _compute_point(model, variables, f_i)
         _log_subproblem_solve(model, variables)
-        model.ideal_point[i] = yI[i] = y_i
+        model.ideal_point[i] = y_i
+        yI[i] = _round_objective(y_i, obj_constants[i])
         MOI.set(inner, MOI.ObjectiveSense(), MOI.MAX_SENSE)
         optimize_inner!(model)
         status = MOI.get(inner, MOI.TerminationStatus())
@@ -197,7 +211,7 @@ function _minimize_multiobjective!(
         _log_subproblem_solve(model, variables)
         # This is not really the nadir point, it's an upper bound on the
         # nadir point because we later add epsilon constraints with a -1.
-        yN[i] = y_i + 1
+        yN[i] = _round_objective(y_i, obj_constants[i]) + 1
     end
     MOI.set(inner, MOI.ObjectiveSense(), MOI.MIN_SENSE)
     # Instead of adding and deleting constraints during the algorithm, we add
@@ -206,15 +220,14 @@ function _minimize_multiobjective!(
         MOI.Utilities.normalize_and_add_constraint(
             inner,
             f_i,
-            MOI.LessThan(yN[i]),
+            MOI.LessThan(yN[i] + obj_constants[i]),
         ) for (i, f_i) in enumerate(scalars)
     ]
-    obj_constants = MOI.constant.(scalars)
-    U_N = Dict{Vector{Float64},Vector{Vector{Vector{Float64}}}}(
+    U_N = Dict{Vector{Int},Vector{Vector{Vector{Int}}}}(
         # The nadir point, except for the ideal point in position k
         yN => [[_get_child(yN, yI, k)] for k in 1:n],
     )
-    V = [Tuple{Vector{Float64},Vector{Float64}}[] for k in 1:n]
+    V = [Tuple{Vector{Int},Vector{Int}}[] for k in 1:n]
     status = MOI.OPTIMAL
     sum_f = MOI.Utilities.operate(+, Float64, scalars...)
     while !isempty(U_N)
@@ -230,11 +243,11 @@ function _minimize_multiobjective!(
         # computing the nadir point.
         for i in 1:n
             u_i = ifelse(i == k, yN[i], u[i])
-            set = MOI.LessThan(u_i - 1.0 - obj_constants[i])
+            set = MOI.LessThan(Float64(u_i - 1))
             MOI.set(inner, MOI.ConstraintSet(), ε_constraints[i], set)
         end
         if warm_start_supported && !isempty(solutions)
-            for (x_i, start_i) in solutions[last(U_N[u][k])]
+            for (x_i, start_i) in solutions[last(U_N[u][k])].x
                 MOI.set(inner, MOI.VariablePrimalStart(), x_i, start_i)
             end
         end
@@ -244,23 +257,23 @@ function _minimize_multiobjective!(
             break
         end
         # Now solve problem Π²(k, u)
-        y_k = MOI.get(inner, MOI.ObjectiveValue())::Float64
-        set = MOI.LessThan(y_k - obj_constants[k])
+        obj = MOI.get(inner, MOI.ObjectiveValue())::Float64
+        y_k = _round_objective(obj, obj_constants[k])
+        set = MOI.LessThan(Float64(y_k))
         MOI.set(inner, MOI.ConstraintSet(), ε_constraints[k], set)
         MOI.set(inner, MOI.ObjectiveFunction{typeof(sum_f)}(), sum_f)
         optimize_inner!(model)
         status = MOI.get(inner, MOI.TerminationStatus())
         if !_is_scalar_status_optimal(status)
             break
         end
-        X, Y = _compute_point(model, variables, f)
-        _log_subproblem_solve(model, Y)
+        X, Y_float = _compute_point(model, variables, f)
+        _log_subproblem_solve(model, Y_float)
+        Y = _round_objective.(Y_float, obj_constants)
         push!(V[k], (u, Y))
-        # We want `if !(Y in U_N[u][k])` but this tests exact equality. We
-        # want an approximate comparison.
-        if all(!isapprox(Y; atol = _TAMBY_VANDERPOOTEN_ATOL), U_N[u][k])
+        if !(Y in U_N[u][k])
             _update_search_region(U_N, Y, yN)
-            solutions[Y] = X
+            solutions[Y] = SolutionPoint(X, Y_float)
         end
         _clean_search_region(U_N, yI, V, k)
     end
@@ -288,7 +301,7 @@ end
 
 function _clean_search_region_inner(u′, U_N, yI, V, k)
     for (u′_k, yI_k) in zip(u′, yI)
-        if isapprox(u′_k, yI_k; atol = _TAMBY_VANDERPOOTEN_ATOL)
+        if u′_k == yI_k
             return true
         end
     end
@@ -303,7 +316,7 @@ function _clean_search_region_inner(u′, U_N, yI, V, k)
 end
 
 function _comparison_line_16(u′, u, y_k, k)
-    if !≈(y_k[k], u′[k]; atol = _TAMBY_VANDERPOOTEN_ATOL)
+    if y_k[k] != u′[k]
         return false
     end
     for (i, (u′_i, u_i)) in enumerate(zip(u′, u))

test/algorithms/test_TambyVanderpooten.jl

@@ -173,6 +173,68 @@ function test_solve_failures()
     return
 end
 
+function _solve_noisy_mock(mock, counter::Ref{Int})
+    highs = HiGHS.Optimizer()
+    MOI.set(highs, MOI.Silent(), true)
+    index_map = MOI.copy_to(highs, mock)
+    MOI.optimize!(highs)
+    x = [index_map[xi] for xi in MOI.get(mock, MOI.ListOfVariableIndices())]
+    # Perturb the solution by a deterministic amount that is within the
+    # tolerances that a typical MIP solver might return.
+    noise() = 1e-6 * sin(counter[] += 1)
+    x_primal = [xi + noise() for xi in MOI.get(highs, MOI.VariablePrimal(), x)]
+    status = MOI.get(highs, MOI.TerminationStatus())
+    MOI.Utilities.mock_optimize!(mock, status, x_primal)
+    obj = MOI.get(highs, MOI.ObjectiveValue())
+    MOI.set(mock, MOI.ObjectiveValue(), obj + noise())
+    return
+end
+
+function test_noisy_solutions()
+    C = [
+        [3 13 3 20 9; 4 13 3 4 9; 18 9 3 9 17; 1 16 16 11 3; 18 4 7 13 11],
+        [1 12 4 6 20; 18 5 13 20 7; 9 3 3 2 15; 11 18 1 7 13; 6 7 3 15 18],
+        [13 13 13 14 15; 16 3 20 7 13; 12 17 9 8 3; 6 10 4 18 5; 4 9 7 14 11],
+    ]
+    n = size(C[1], 1)
+    solutions = Dict{Bool,Set{Vector{Int}}}()
+    for noisy in (false, true)
+        counter = Ref(0)
+        model = MOA.Optimizer() do
+            mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
+            MOI.Utilities.set_mock_optimize!(
+                mock,
+                m -> noisy ? _solve_noisy_mock(m, counter) : _solve_mock(m),
+            )
+            return mock
+        end
+        MOI.set(model, MOA.Algorithm(), MOA.TambyVanderpooten())
+        MOI.set(model, MOI.Silent(), true)
+        MOI.set(model, MOI.TimeLimitSec(), 60.0)
+        x = reshape(MOI.add_variables(model, n * n), n, n)
+        MOI.add_constraint.(model, x, MOI.ZeroOne())
+        for i in 1:n
+            MOI.add_constraint(model, sum(1.0 .* x[i, :]), MOI.EqualTo(1.0))
+            MOI.add_constraint(model, sum(1.0 .* x[:, i]), MOI.EqualTo(1.0))
+        end
+        f_k = [sum(1.0 * C[k] .* x) for k in 1:3]
+        f = MOI.Utilities.operate(vcat, Float64, f_k...)
+        MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+        MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
+        MOI.optimize!(model)
+        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
+        N = MOI.get(model, MOI.ResultCount())
+        Y_float = [MOI.get(model, MOI.ObjectiveValue(i)) for i in 1:N]
+        # The algorithm rounds internally, but we return the solver's values.
+        @test any(y -> any(!isinteger, y), Y_float) == noisy
+        Y = [round.(Int, y) for y in Y_float]
+        @test allunique(Y)
+        solutions[noisy] = Set(Y)
+    end
+    @test solutions[false] == solutions[true]
+    return
+end
+
 end  # module TestTambyVanderpooten
 
 TestTambyVanderpooten.run_tests()