Add support for setting MOI.VariablePrimalStart (#588)

Author: odow

docs/src/manual/advanced.md

@@ -34,6 +34,31 @@ p.constraints[1].dual
 constraint.dual
 ```
 
+## Warmstarting
+
+If you're solving the same problem many times with different values of a
+parameter, Convex.jl can initialize many solvers with the solution to the
+previous problem, which sometimes speeds up the solution time. This is called a
+**warm start**.
+
+To use this feature, pass the optional argument `warmstart=true` to the `solve!`
+method.
+
+```julia
+using Convex, SCS
+n = 1_000
+y = rand(n)
+x = Variable(n)
+lambda = Variable(Positive())
+fix!(lambda, 100)
+problem = minimize(sumsquares(y - x) + lambda * sumsquares(x - 10))
+@time solve!(problem, SCS.Optimizer)
+# Now warmstart. If the solver takes advantage of warmstarts, this run will be
+# faster
+fix!(lambda, 105)
+@time solve!(problem, SCS.Optimizer; warmstart = true)
+```
+
 ## Fixing and freeing variables
 
 Convex.jl allows you to fix a variable `x` to a value by calling the `fix!`

src/solution.jl

@@ -7,11 +7,56 @@ function add_variables!(model, var::AbstractVariable)
     return MOI.add_variables(model, length(var))
 end
 
+function _primal_start(
+    ::Type{T},
+    x::Vector{MOI.VariableIndex},
+    start::AbstractArray,
+) where {T}
+    return convert(Vector{T}, real.(reshape(start, length(x))))
+end
+
+function _primal_start(
+    ::Type{T},
+    x::Tuple{Vector{MOI.VariableIndex},Vector{MOI.VariableIndex}},
+    start::AbstractArray,
+) where {T}
+    y = reshape(start, length(x[1]))
+    return convert(Vector{T}, real.(y)), convert(Vector{T}, imag.(y))
+end
+
+function _primal_start(::Type{T}, x::Vector, start::Number) where {T}
+    return _primal_start(T, x, fill(start, length(x)))
+end
+
+function _primal_start(::Type{T}, x::Tuple, start::Number) where {T}
+    return _primal_start(T, x, fill(start, length(x[1])))
+end
+
+function _add_variable_primal_start(context::Convex.Context{T}) where {T}
+    attr = MOI.VariablePrimalStart()
+    for (id, moi_indices) in context.var_id_to_moi_indices
+        x = context.id_to_variables[id]
+        if x.value === nothing
+            continue
+        elseif moi_indices isa Tuple  # Variable is complex
+            start_re, start_im = _primal_start(T, moi_indices, x.value)
+            MOI.set(context.model, attr, moi_indices[1], start_re)
+            MOI.set(context.model, attr, moi_indices[2], start_im)
+        else
+            @assert moi_indices isa Vector{MOI.VariableIndex}
+            start = _primal_start(T, moi_indices, x.value)
+            MOI.set(context.model, attr, moi_indices, start)
+        end
+    end
+    return
+end
+
 """
     solve!(
         problem::Problem,
         optimizer_factory;
         silent_solver = false,
+        warmstart::Bool = true,
     )
 
 Solves the problem, populating `problem.optval` with the optimal value, as well
@@ -22,15 +67,27 @@ Optional keyword arguments:
 
  * `silent_solver`: whether the solver should be silent (and not emit output or
    logs) during the solution process.
+ * `warmstart` (default: `false`): whether the solver should start the
+   optimization from a previous optimal value (according to the current primal
+   value of the variables in the problem, which can be set by [`set_value!`](@ref).
 """
-function solve!(p::Problem, optimizer_factory; silent_solver = false)
+function solve!(
+    p::Problem,
+    optimizer_factory;
+    silent_solver = false,
+    warmstart::Bool = false,
+)
     if problem_vexity(p) in (ConcaveVexity(), NotDcp())
         throw(DCPViolationError())
     end
     context = Context(p, optimizer_factory)
     if silent_solver
         MOI.set(context.model, MOI.Silent(), true)
     end
+    attr = MOI.VariablePrimalStart()
+    if warmstart && MOI.supports(context.model, attr, MOI.VariableIndex)
+        _add_variable_primal_start(context)
+    end
     if context.detected_infeasible_during_formulation[]
         p.status = MOI.INFEASIBLE
     else

test/test_problem_depot.jl

@@ -7,7 +7,6 @@ import Clarabel
 # import ECOS
 import GLPK
 import MathOptInterface as MOI
-# import SCS
 
 function runtests()
     for name in names(@__MODULE__; all = true)
@@ -40,29 +39,22 @@ end
 test_problem_depot_load_Float64() = _test_problem_depot_load(Float64)
 test_problem_depot_load_BigFloat() = _test_problem_depot_load(BigFloat)
 
-# function test_SCS_with_warmstarts()
-#     Convex.ProblemDepot.run_tests(
-#         exclude = [
-#             r"mip",
-#             # TODO(odow): investigate
-#             r"sdp_lieb_ando",
-#             # Tolerance issue with SCS 3.0
-#             r"sdp_Real_Variables_with_complex_equality_constraints",
-#         ],
-#     ) do p
-#         return solve!(
-#             p,
-#             MOI.OptimizerWithAttributes(
-#                 SCS.Optimizer,
-#                 "verbose" => 0,
-#                 "eps_rel" => 1e-6,
-#                 "eps_abs" => 1e-6,
-#             );
-#             warmstart = true,
-#         )
-#     end
-#     return
-# end
+function test_Clarabel_warmstarts()
+    # `sdp_quantum_relative_entropy3_lowrank` failed on CI for Ubuntu with
+    #   Expression: ≈(p.optval, evaluate(quantum_relative_entropy(B, A)), atol = atol, rtol = rtol)
+    #    Evaluated: -4.887297347885561e-6 ≈ Inf (atol=0.001, rtol=0.0)
+    Convex.ProblemDepot.run_tests(;
+        exclude = [r"mip", r"sdp_quantum_relative_entropy3_lowrank"],
+    ) do p
+        return solve!(
+            p,
+            Clarabel.Optimizer;
+            silent_solver = true,
+            warmstart = true,
+        )
+    end
+    return
+end
 
 function test_Clarabel()
     # `sdp_quantum_relative_entropy3_lowrank` failed on CI for Ubuntu with

test/test_utilities.jl

@@ -1246,6 +1246,44 @@ function test_set_value_complex()
     return
 end
 
+function test_variable_primal_start()
+    x1 = Variable()
+    x2 = Variable(2)
+    x2.value = [1.0, 2.0]
+    x3 = Variable(3, 3)
+    x3.value = 3.0
+    x4 = ComplexVariable()
+    x4.value = 1 + 2im
+    x5 = ComplexVariable(2, 2)
+    x5.value = [1 2; 3im 4]
+    problem = satisfy([x1 >= 0, x2 <= 0, x3 <= 1, x4 == 5, x5 == 0])
+    context = Convex.Context(
+        problem,
+        () ->
+            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
+    )
+    Convex._add_variable_primal_start(context)
+    for (key, query) in (
+        x1 => [nothing],
+        x2 => [1.0, 2.0],
+        x3 => fill(3.0, 9),
+        x4 => ([1.0], [2.0]),
+        x5 => ([1.0, 0.0, 2.0, 4.0], [0.0, 3.0, 0.0, 0.0]),
+    )
+        variable = context.var_id_to_moi_indices[key.id_hash]
+        start = if variable isa Tuple
+            (
+                MOI.get(context.model, MOI.VariablePrimalStart(), variable[1]),
+                MOI.get(context.model, MOI.VariablePrimalStart(), variable[2]),
+            )
+        else
+            MOI.get(context.model, MOI.VariablePrimalStart(), variable)
+        end
+        @test start == query
+    end
+    return
+end
+
 end  # TestUtilities
 
 TestUtilities.runtests()