Implement starting values for Constraint.QuadtoSOCBridge (#2115)

* Implement starting values for Constraint.QuadtoSOCBridge

* Fixes

* Fix format

* Remove show

* Add removed line back ??

* Update src/attributes.jl

Co-authored-by: Mathieu Besançon <[email protected]>

* Update src/attributes.jl

Co-authored-by: Mathieu Besançon <[email protected]>

* Update src/attributes.jl

Co-authored-by: Mathieu Besançon <[email protected]>

* Remove docstring notes

* Add test with copy_to

---------

Co-authored-by: Mathieu Besançon <[email protected]>

Author: blegat

src/Bridges/Constraint/bridges/quad_to_soc.jl

@@ -9,21 +9,25 @@
 
 `QuadtoSOCBridge` converts quadratic inequalities
 ```math
-\\frac{1}{2}x^T Q x + a^T x + b \\le 0
+\\frac{1}{2}x^T Q x + a^T x \\le ub
 ```
 into [`MOI.RotatedSecondOrderCone`](@ref) constraints, but it only applies when
 ``Q`` is positive definite.
 
 This is because, if `Q` is positive definite, there exists `U` such that
 ``Q = U^T U``, and so the inequality can then be rewritten as;
 ```math
-\\|U x\\|_2^2 \\le 2 (-a^T x - b)
+\\|U x\\|_2^2 \\le 2 (-a^T x + ub)
 ```
 
-Therefore, `QuadtoSOCBridge` implements the following reformulation:
+Therefore, `QuadtoSOCBridge` implements the following reformulations:
 
-  * ``\\frac{1}{2}x^T Q x + a^T x + b \\le 0`` into
-    ``(1, -a^T x - b, Ux) \\in RotatedSecondOrderCone``
+  * ``\\frac{1}{2}x^T Q x + a^T x \\le ub`` into
+    ``(1, -a^T x + ub, Ux) \\in RotatedSecondOrderCone``
+    where ``Q = U^T U``
+  * ``\\frac{1}{2}x^T Q x + a^T x \\ge lb`` into
+    ``(1, a^T x - lb, Ux) \\in RotatedSecondOrderCone``
+    where ``-Q = U^T U``
 
 ## Source node
 
@@ -50,6 +54,7 @@ struct QuadtoSOCBridge{T} <: AbstractBridge
     dimension::Int  # dimension of the SOC constraint
     less_than::Bool # whether the constraint was ≤ or ≥
     set_constant::T # the constant that was on the set
+    index_to_variable_map::Vector{MOI.VariableIndex}
 end
 
 const QuadtoSOC{T,OT<:MOI.ModelLike} =
@@ -84,17 +89,17 @@ function bridge_constraint(
         )
     end
     # Construct the VectorAffineFunction. We're aiming for:
-    #  |          1 |
-    #  | -a^T x - b | ∈ RotatedSecondOrderCone()
-    #  |     Lx + 0 |
+    #  |          1  |
+    #  | -a^T x + ub | ∈ RotatedSecondOrderCone()
+    #  |     Ux + 0  |
     # Start with the -a^T x terms...
     vector_terms = MOI.VectorAffineTerm{T}[
         MOI.VectorAffineTerm(
             2,
             MOI.ScalarAffineTerm(scale * term.coefficient, term.variable),
         ) for term in func.affine_terms
     ]
-    # Add the Lx terms...
+    # Note that `L = U'` so we add the terms of L'x
     L, p = SparseArrays.sparse(F.L), F.p
     I, J, V = SparseArrays.findnz(L)
     for i in 1:length(V)
@@ -108,31 +113,32 @@ function bridge_constraint(
             MOI.VectorAffineTerm(J[i] + 2, MOI.ScalarAffineTerm(V[i], xi)),
         )
     end
-    # This is the [1, b, 0] vector...
+    # This is the [1, ub, 0] vector...
     set_constant = MOI.constant(set)
-    vector_constant = vcat(
-        one(T),
-        scale * (func.constant - set_constant),
-        zeros(T, size(L, 1)),
-    )
+    MOI.throw_if_scalar_and_constant_not_zero(func, typeof(set))
+    vector_constant = vcat(one(T), -scale * set_constant, zeros(T, size(L, 1)))
     f = MOI.VectorAffineFunction(vector_terms, vector_constant)
     dimension = MOI.output_dimension(f)
     soc = MOI.add_constraint(model, f, MOI.RotatedSecondOrderCone(dimension))
-    return QuadtoSOCBridge(soc, dimension, less_than, set_constant)
+    return QuadtoSOCBridge(
+        soc,
+        dimension,
+        less_than,
+        set_constant,
+        index_to_variable_map,
+    )
 end
 
 function _matrix_from_quadratic_terms(
     terms::Vector{MOI.ScalarQuadraticTerm{T}},
 ) where {T}
     variable_to_index_map = Dict{MOI.VariableIndex,Int}()
-    index_to_variable_map = Dict{Int,MOI.VariableIndex}()
-    n = 0
+    index_to_variable_map = MOI.VariableIndex[]
     for term in terms
         for variable in (term.variable_1, term.variable_2)
             if !(variable in keys(variable_to_index_map))
-                n += 1
-                variable_to_index_map[variable] = n
-                index_to_variable_map[n] = variable
+                push!(index_to_variable_map, variable)
+                variable_to_index_map[variable] = length(index_to_variable_map)
             end
         end
     end
@@ -150,6 +156,7 @@ function _matrix_from_quadratic_terms(
         end
     end
     # Duplicate terms are summed together in `sparse`
+    n = length(index_to_variable_map)
     return SparseArrays.sparse(I, J, V, n, n), index_to_variable_map
 end
 
@@ -209,13 +216,32 @@ function MOI.delete(model::MOI.ModelLike, bridge::QuadtoSOCBridge)
     return
 end
 
-# Attributes, Bridge acting as a constraint
+function MOI.supports(
+    model::MOI.ModelLike,
+    attr::Union{MOI.ConstraintPrimalStart,MOI.ConstraintDualStart},
+    ::Type{QuadtoSOCBridge{T}},
+) where {T}
+    F, S = MOI.VectorAffineFunction{T}, MOI.RotatedSecondOrderCone
+    return MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex) &&
+           MOI.supports(model, attr, MOI.ConstraintIndex{F,S})
+end
+
 function MOI.get(
     model::MOI.ModelLike,
-    attr::MOI.ConstraintPrimal,
+    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintPrimalStart},
     bridge::QuadtoSOCBridge,
 )
+    # The constraint primal is x'Qx/2 + a'x
+    # If `less_than` then `Q = U'U` and we have the value of
+    # `Ux` and `-a'x + ub`, so we get it with
+    # `(Ux)'Ux / 2 - (-a'x + ub) + ub`
+    # Otherwise, `Q = -U'U` and we have the value of
+    # `Ux` and `a'x - ub`, so we get it with
+    # `-(Ux)'Ux / 2 + (a'x - lb) + ub`
     soc = MOI.get(model, attr, bridge.soc)
+    if soc === nothing
+        return nothing
+    end
     output = sum(soc[i]^2 for i in 3:bridge.dimension)
     output /= 2
     output -= soc[1] * soc[2]
@@ -226,6 +252,68 @@ function MOI.get(
     return output
 end
 
+function _primal_start_or_error(model, attr, v)
+    var_attr = MOI.VariablePrimalStart()
+    value = MOI.get(model, MOI.VariablePrimalStart(), v)
+    if isnothing(value)
+        error(
+            "In order to set the `$attr`, the",
+            "`MOI.Bridges.Constraint.QuadtoSOCBridge` needs to get the ",
+            "`$var_attr` but it is not set. Set the `$var_attr` first before ",
+            "setting the `$attr` in order to fix this.",
+        )
+    end
+    return value
+end
+
+function MOI.set(
+    model::MOI.ModelLike,
+    attr::MOI.ConstraintPrimalStart,
+    bridge::QuadtoSOCBridge{T},
+    value,
+) where {T}
+    # `value` represent `x'Qx/2 + a'x + ε` where `ε` is
+    # the difference between the value of the slack variable and the value of the function.
+    # That is, if `less_than`, we set
+    # |       1       |
+    # | -a'x - ε + ub |
+    # |     U * x     |
+    # which is obtained as
+    # |          1          |
+    # | x'Qx/2 - value + ub |
+    # |        U * x        |
+    # Otherwise, we set
+    # |       1       |
+    # |  a'x + ε - lb |
+    # |     U * x     |
+    # which is obtained as
+    # |          1          |
+    # | value - x'Qx/2 - lb |
+    # |        U * x        |
+    # where we compute `x'Qx/2` and `U * x` using the starting values of the variable.
+    soc = MOI.get(model, MOI.ConstraintFunction(), bridge.soc)
+    Ux = MOI.Utilities.eval_variables(MOI.Utilities.eachscalar(soc)[3:end]) do v
+        return _primal_start_or_error(model, attr, v)
+    end
+    if bridge.less_than
+        s2 = Ux'Ux / 2 - value + bridge.set_constant
+    else
+        s2 = Ux'Ux / 2 + value - bridge.set_constant
+    end
+    MOI.set(model, attr, bridge.soc, [1; s2; Ux])
+    return
+end
+
+function MOI.set(
+    model::MOI.ModelLike,
+    attr::MOI.ConstraintPrimalStart,
+    bridge::QuadtoSOCBridge,
+    ::Nothing,
+)
+    MOI.set(model, attr, bridge.soc, nothing)
+    return
+end
+
 # Lemma: If (1, s, x), (v, u, y) in RotatedSecondOrderCone and
 #        (1, s, x) ⋅ (v, u, y) = 0, then we have
 #        y = -u*x, v = u*||x||_2^2/2.
@@ -261,13 +349,48 @@ end
 # is exactly the same. Q.E.D.
 function MOI.get(
     model::MOI.ModelLike,
-    attr::MOI.ConstraintDual,
+    attr::Union{MOI.ConstraintDual,MOI.ConstraintDualStart},
     bridge::QuadtoSOCBridge,
 )
-    λ = MOI.get(model, attr, bridge.soc)[2]
+    dual = MOI.get(model, attr, bridge.soc)
+    if dual === nothing
+        return nothing
+    end
+    λ = dual[2]
     return bridge.less_than ? -λ : λ
 end
 
+# Let `(v, u, y)` be the dual of the RSOC and `λ` be the dual of the quadratic.
+# From same reasoning as above, we know that `u` is `-λ`.
+# From the Lemma above, we have ``
+# `y = -u * x`, `v = u*||x||_2^2/2`
+function MOI.set(
+    model::MOI.ModelLike,
+    attr::MOI.ConstraintDualStart,
+    bridge::QuadtoSOCBridge{T},
+    λ,
+) where {T}
+    u = bridge.less_than ? -λ : λ
+    x = T[
+        _primal_start_or_error(model, attr, xi) for
+        xi in bridge.index_to_variable_map
+    ]
+    v = u * sum(x .^ 2) / 2
+    y = -u * x
+    MOI.set(model, attr, bridge.soc, [v; u; y])
+    return
+end
+
+function MOI.set(
+    model::MOI.ModelLike,
+    attr::MOI.ConstraintDualStart,
+    bridge::QuadtoSOCBridge,
+    ::Nothing,
+)
+    MOI.set(model, attr, bridge.soc, nothing)
+    return
+end
+
 function MOI.get(
     ::MOI.ModelLike,
     ::MOI.ConstraintSet,

test/Bridges/Constraint/quad_to_soc.jl

@@ -106,7 +106,8 @@ function test_quadratic_constraints_with_2_variables()
 end
 
 function test_qcp()
-    mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
+    model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
+    mock = MOI.Utilities.MockOptimizer(model)
     config = MOI.Test.Config()
     bridged_mock = MOI.Bridges.Constraint.QuadtoSOC{Float64}(mock)
     MOI.Utilities.set_mock_optimize!(
@@ -159,6 +160,15 @@ function test_qcp()
     MOI.Test.test_quadratic_constraint_LessThan(bridged_mock, config)
     MOI.empty!(bridged_mock)
     MOI.Test.test_quadratic_constraint_GreaterThan(bridged_mock, config)
+    F = MOI.ScalarQuadraticFunction{Float64}
+    S = MOI.GreaterThan{Float64}
+    ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{F,S}()))
+    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
+        err = ErrorException(
+            "In order to set the `$attr`, the`MOI.Bridges.Constraint.QuadtoSOCBridge` needs to get the `MathOptInterface.VariablePrimalStart()` but it is not set. Set the `MathOptInterface.VariablePrimalStart()` first before setting the `$attr` in order to fix this.",
+        )
+        @test_throws err MOI.set(bridged_mock, attr, ci, 0.0)
+    end
     return
 end
 
@@ -207,6 +217,7 @@ function test_runtests()
         variables: x, y
         [1.0, 0.0, 2.0 * x, 2.0 * y] in RotatedSecondOrderCone(4)
         """,
+        variable_start = 0.0,
     )
     MOI.Bridges.runtests(
         MOI.Bridges.Constraint.QuadtoSOCBridge,
@@ -219,9 +230,54 @@ function test_runtests()
         [1.0, 0.0, 2.0 * x, 2.0 * y] in RotatedSecondOrderCone(4)
         """,
     )
+    MOI.Bridges.runtests(
+        MOI.Bridges.Constraint.QuadtoSOCBridge,
+        """
+        variables: x, y
+        2.0 * x * x + 0.5 * y * y + 2.0 * x <= 3.0
+        """,
+        """
+        variables: x, y
+        [1.0, -2.0 * x + 3.0, 2.0 * x, 1.0 * y] in RotatedSecondOrderCone(4)
+        """,
+    )
+    MOI.Bridges.runtests(
+        MOI.Bridges.Constraint.QuadtoSOCBridge,
+        """
+        variables: x, y, z
+        -2.0 * x * x + -1.0 * y * y + -2.0 * x * y + 3.0 * z >= 4.0
+        """,
+        """
+        variables: x, y, z
+        [1.0, 3.0 * z + -4.0, 2.0 * x + 1.0 * y, 1.0 * y] in RotatedSecondOrderCone(4)
+        """,
+    )
     return
 end
 
+"""
+    test_copy_to_start()
+
+Tests that `copy_to` copies the variable starting values first and the
+constraint primal and dual starting values second as is needed by the
+[`MOI.Bridges.Constraint.QuadtoSOCBridge`](@ref).
+"""
+function test_copy_to_start()
+    src = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
+    x = MOI.add_variable(src)
+    c = MOI.add_constraint(src, 1.0x * x + 2.0x, MOI.LessThan(-1.0))
+    MOI.set(src, MOI.ConstraintPrimalStart(), c, 1.0)
+    MOI.set(src, MOI.ConstraintDualStart(), c, -1.0)
+    MOI.set(src, MOI.VariablePrimalStart(), x, -1.0)
+
+    model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
+    dest = MOI.Bridges.Constraint.QuadtoSOC{Float64}(model)
+    index_map = MOI.copy_to(dest, src)
+    @test MOI.get(dest, MOI.ConstraintPrimalStart(), index_map[c]) == 1.0
+    @test MOI.get(dest, MOI.ConstraintDualStart(), index_map[c]) == -1.0
+    @test MOI.get(dest, MOI.VariablePrimalStart(), index_map[x]) == -1.0
+end
+
 end  # module
 
 TestConstraintQuadToSOC.runtests()