Merge pull request #109 from JuliaOpt/bl/psd2x2_variable_bridge

Add PSD2x2 variable bridge

Author: blegat

src/SumOfSquares.jl

@@ -35,6 +35,7 @@ const MOIB = MOI.Bridges
 
 # Variable Bridges
 include("variable_bridge.jl")
+include("psd2x2_variable_bridge.jl")
 include("scaled_diagonally_dominant_variable_bridge.jl")
 include("generic_variable_bridge.jl")
 include("copositive_inner_variable_bridge.jl")

src/attributes.jl

@@ -7,33 +7,46 @@ A constraint attribute for the monomials indexing the
 struct CertificateMonomials <: MOI.AbstractConstraintAttribute end
 
 """
-    struct GramMatrixAttribute <: MOI.AbstractConstraintAttribute end
+    GramMatrixAttribute(N)
+    GramMatrixAttribute()
 
 A constraint attribute for the [`GramMatrix`](@ref) of a constraint, that is,
 the positive semidefinite matrix `Q` indexed by the monomials in the vector `X`
 such that ``X^\\top Q X`` is the sum-of-squares certificate of the constraint.
-The
 """
-struct GramMatrixAttribute <: MOI.AbstractConstraintAttribute end
+struct GramMatrixAttribute <: MOI.AbstractConstraintAttribute
+    N::Int
+end
+GramMatrixAttribute() = GramMatrixAttribute(1)
 
 """
-    struct MomentMatrixAttribute <: MOI.AbstractConstraintAttribute end
+    MomentMatrixAttribute(N)
+    MomentMatrixAttribute()
 
 A constraint attribute fot the `MomentMatrix` of a constraint.
 """
-struct MomentMatrixAttribute <: MOI.AbstractConstraintAttribute end
+struct MomentMatrixAttribute <: MOI.AbstractConstraintAttribute
+    N::Int
+end
+MomentMatrixAttribute() = MomentMatrixAttribute(1)
 
 """
-    struct LagrangianMultipliers <: MOI.AbstractConstraintAttribute end
+    LagrangianMultipliers(N)
+    LagrangianMultipliers()
 
-A constraint attribute fot the `LagrangianMultipliers` assiciated to the
+A constraint attribute fot the `LagrangianMultipliers` associated to the
 inequalities of the domain of a constraint. There is one multiplier per
 inequality and no multiplier for equalities as the equalities are handled by
 reducing the polynomials over the ideal they generate instead of explicitely
 creating multipliers.
 """
-struct LagrangianMultipliers <: MOI.AbstractConstraintAttribute end
+struct LagrangianMultipliers <: MOI.AbstractConstraintAttribute
+    N::Int
+end
+LagrangianMultipliers() = LagrangianMultipliers(1)
 
+# Needs to declare it set by optimize that it is not queried in the Caching
+# optimize, even of `CertificateMonomials` which is set befor optimize.
 function MOI.is_set_by_optimize(::Union{CertificateMonomials,
                                         GramMatrixAttribute,
                                         MomentMatrixAttribute,

src/copositive_inner.jl

@@ -44,9 +44,11 @@ which is the same as, using the `domain` keyword,
 @constraint(model, x^2 - 2x*y + y^2 in SOSCone(), domain = @set x*y ≥ 0)
 ```
 
-For consistency with its equivalent forms, the [`GramMatrixAttribute`](@ref) for
-this constraint is given by the gram matrix in the `psd_inner` cone, i.e. which
-should be equal to `Q - Λ`.
+As an important difference with its equivalent forms, the
+[`GramMatrixAttribute`](@ref) for the copositive constraint is given by matrix
+`Q` while for the equivalent form using the  domain` keyword, the value
+of the attribute would correspond to the the gram matrix in the `psd_inner`
+cone, i.e. which should be equal to `Q - Λ`.
 
 [BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R.
 *Semidefinite Optimization and Convex Algebraic Geometry*.

src/copositive_inner_variable_bridge.jl

@@ -40,7 +40,7 @@ function MOI.delete(model::MOI.ModelLike, bridge::CopositiveInnerVariableBridge)
 end
 
 function MOI.get(model::MOI.ModelLike,
-                 attr::Union{MomentMatrixAttribute, GramMatrixAttribute},
+                 attr::Union{MOI.ConstraintDual, MOI.ConstraintPrimal},
                  bridge::CopositiveInnerVariableBridge)
     return MOI.get(model, attr, bridge.variable_bridge)
 end

src/generic_variable_bridge.jl

@@ -42,12 +42,13 @@ function MOI.delete(model::MOI.ModelLike, bridge::GenericVariableBridge)
     end
 end
 
-function MOI.get(model::MOI.ModelLike, ::MomentMatrixAttribute,
+function MOI.get(model::MOI.ModelLike,
+                 attr::MOI.ConstraintPrimal,
                  bridge::GenericVariableBridge)
-    return MOI.get(model, MOI.ConstraintDual(), bridge.gram_constraint)
+    return MOI.get(model, MOI.VariablePrimal(attr.N), bridge.gram_matrix)
 end
-
-function MOI.get(model::MOI.ModelLike, ::GramMatrixAttribute,
+function MOI.get(model::MOI.ModelLike,
+                 attr::MOI.ConstraintDual,
                  bridge::GenericVariableBridge)
-    return MOI.get(model, MOI.VariablePrimal(), bridge.gram_matrix)
+    return MOI.get(model, attr, bridge.gram_constraint)
 end

src/psd2x2_variable_bridge.jl

@@ -0,0 +1,67 @@
+struct PositiveSemidefinite2x2VariableBridge{T} <: AbstractVariableBridge
+    variables::Vector{MOI.VariableIndex}
+    rsoc::MOI.ConstraintIndex{MOI.VectorOfVariables,
+                              MOI.RotatedSecondOrderCone}
+end
+
+function add_variable_bridge(
+    ::Type{PositiveSemidefinite2x2VariableBridge{T}}, model::MOI.ModelLike,
+    s::PositiveSemidefinite2x2ConeTriangle) where {T}
+    x = MOI.add_variables(model, 3)
+    rsoc = MOI.add_constraint(model, MOI.VectorOfVariables(x), MOI.RotatedSecondOrderCone(3))
+    Q12 = MOIU.operate(/, T, MOI.SingleVariable(x[3]), convert(T, √2))
+    g = typeof(Q12)[MOI.SingleVariable(x[1]), Q12, MOI.SingleVariable(x[2])]
+    return g, PositiveSemidefinite2x2VariableBridge{T}(x, rsoc)
+end
+
+function MOIB.added_constraint_types(
+    ::Type{PositiveSemidefinite2x2VariableBridge{T}}) where {T}
+    return [(MOI.VectorOfVariables, MOI.RotatedSecondOrderCone)]
+end
+
+function variable_bridge_type(::Type{PositiveSemidefinite2x2ConeTriangle},
+                              T::Type)
+    return PositiveSemidefinite2x2VariableBridge{T}
+end
+
+
+# Attributes, VariableBridge acting as an model
+function MOI.get(::PositiveSemidefinite2x2VariableBridge,
+                 ::MOI.NumberOfVariables)
+    return 3
+end
+function MOI.get(
+    ::PositiveSemidefinite2x2VariableBridge,
+    ::MOI.NumberOfConstraints{MOI.VectorOfVariables,
+                              MOI.RotatedSecondOrderCone})
+    return 1
+end
+function MOI.get(
+    bridge::PositiveSemidefinite2x2VariableBridge,
+    ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables,
+                                  MOI.RotatedSecondOrderCone})
+    return [bridge.rsoc]
+end
+
+# Indices
+function MOI.delete(model::MOI.ModelLike,
+                    bridge::PositiveSemidefinite2x2VariableBridge)
+    MOI.delete(model, bridge.rsoc)
+    for vi in bridge.variables
+        MOI.delete(model, vi)
+    end
+end
+
+function MOI.get(model::MOI.ModelLike,
+                 attr::MOI.ConstraintPrimal,
+                 bridge::PositiveSemidefinite2x2VariableBridge)
+    value = MOI.get(model, attr, bridge.rsoc)
+    return [value[1], value[3] / √2, value[2]]
+end
+function MOI.get(model::MOI.ModelLike,
+                 attr::MOI.ConstraintDual,
+                 bridge::PositiveSemidefinite2x2VariableBridge)
+    dual = MOI.get(model, attr, bridge.rsoc)
+    # / 2 (because of different scalar product) * √2 (because of A^{-*} = 1 / √2
+    return [dual[1], dual[3] / √2, dual[2]]
+end

src/scaled_diagonally_dominant_variable_bridge.jl

@@ -1,13 +1,10 @@
-struct ScaledDiagonallyDominantVariableBridge{T} <: AbstractVariableBridge
-    side_dimension::Int
-    variables::Vector{NTuple{3, MOI.VariableIndex}}
-    psd2x2::Vector{MOI.ConstraintIndex{MOI.VectorOfVariables,
-                                       PositiveSemidefinite2x2ConeTriangle}}
+struct ScaledDiagonallyDominantVariableBridge{T, VBT} <: AbstractVariableBridge
+    psd2x2::Vector{VBT}
 end
 
 function add_variable_bridge(
-    ::Type{ScaledDiagonallyDominantVariableBridge{T}}, model::MOI.ModelLike,
-    s::ScaledDiagonallyDominantConeTriangle) where {T}
+    ::Type{ScaledDiagonallyDominantVariableBridge{T, VBT}}, model::MOI.ModelLike,
+    s::ScaledDiagonallyDominantConeTriangle) where {T, VBT}
     # `p.Q` is SDD iff it is the sum of psd matrices Mij that are zero except
     # for entries ii, ij and jj [Lemma 9, AM17].
     n = s.side_dimension
@@ -17,62 +14,53 @@ function add_variable_bridge(
     F = MOI.ScalarAffineFunction{T}
     g = F[zero(F) for i in 1:MOI.dimension(s)]
     num_off_diag = MOI.dimension(s) - n
-    variables = Vector{NTuple{3, MOI.VariableIndex}}(undef, num_off_diag)
-    psd2x2 = Vector{MOI.ConstraintIndex{
-        MOI.VectorOfVariables, PositiveSemidefinite2x2ConeTriangle}}(
-            undef, num_off_diag)
+    psd2x2 = Vector{VBT}(undef, num_off_diag)
     diag_idx(i) = div(i * (i + 1), 2)
     k = 0
     k2x2 = 0
     for j in 1:n
         for i in 1:(j-1)
             k += 1
             k2x2 += 1
-            vii, vij, vjj = MOI.add_variables(model, 3)
-            variables[k2x2] = (vii, vij, vjj)
-            Mii = MOI.SingleVariable(vii)
-            MOIU.operate!(+, T, g[diag_idx(i)], Mii)
-            Mij = MOI.SingleVariable(vij)
-            MOIU.operate!(+, T, g[k], Mij)
-            Mjj = MOI.SingleVariable(vjj)
-            MOIU.operate!(+, T, g[diag_idx(j)], Mjj)
             # PSD constraints on 2x2 matrices are SOC representable
-            psd2x2[k2x2] = MOI.add_constraint(
-                model, MOI.VectorOfVariables([vii, vij, vjj]),
-                PositiveSemidefinite2x2ConeTriangle())
+            x, psd2x2[k2x2] = add_variable_bridge(VBT, model, PositiveSemidefinite2x2ConeTriangle())
+            MOIU.operate!(+, T, g[diag_idx(i)], x[1])
+            MOIU.operate!(+, T, g[k], x[2])
+            MOIU.operate!(+, T, g[diag_idx(j)], x[3])
         end
         k += 1 # diagonal entry `(i, i)`
     end
-    return g, ScaledDiagonallyDominantVariableBridge{T}(n, variables, psd2x2)
+    return g, ScaledDiagonallyDominantVariableBridge{T, VBT}(psd2x2)
 end
 
 function MOIB.added_constraint_types(
-    ::Type{ScaledDiagonallyDominantVariableBridge{T}}) where {T}
-    return [(MOI.VectorOfVariables, PositiveSemidefinite2x2ConeTriangle)]
+    ::Type{ScaledDiagonallyDominantVariableBridge{T, VBT}}) where {T, VBT}
+    added = Tuple{DataType, DataType}[]
+    return append!(added, MOIB.added_constraint_types(VBT))
 end
 
 function variable_bridge_type(::Type{ScaledDiagonallyDominantConeTriangle},
                               T::Type)
-    return ScaledDiagonallyDominantVariableBridge{T}
+    VBT = variable_bridge_type(PositiveSemidefinite2x2ConeTriangle, T)
+    return ScaledDiagonallyDominantVariableBridge{T, VBT}
 end
 
 
 # Attributes, VariableBridge acting as an model
-function MOI.get(bridge::ScaledDiagonallyDominantVariableBridge,
-                 ::MOI.NumberOfVariables)
-    return 3 * length(bridge.variables)
+function MOI.get(bridge::ScaledDiagonallyDominantVariableBridge, attr::MOI.NumberOfVariables)
+    return reduce(+, MOI.get.(bridge.psd2x2, attr), init=0)
 end
-function MOI.get(
-    bridge::ScaledDiagonallyDominantVariableBridge,
-    ::MOI.NumberOfConstraints{MOI.VectorOfVariables,
-                              ScaledDiagonallyDominantConeTriangle})
-    return length(bridge.psd2x2)
+function MOI.get(bridge::ScaledDiagonallyDominantVariableBridge,
+                 attr::MOI.NumberOfConstraints)
+    return reduce(+, MOI.get.(bridge.psd2x2, attr), init=0)
 end
-function MOI.get(
-    bridge::ScaledDiagonallyDominantVariableBridge,
-    ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables,
-                                  ScaledDiagonallyDominantConeTriangle})
-    return bridge.psd2x2
+function MOI.get(bridge::ScaledDiagonallyDominantVariableBridge,
+                 attr::MOI.ListOfConstraintIndices{F, S}) where {F, S}
+    list = MOI.ConstraintIndex{F, S}[]
+    for variable_bridge in bridge.psd2x2
+        append!(list, MOI.get(variable_bridge, attr))
+    end
+    return list
 end
 
 # Indices
@@ -81,9 +69,4 @@ function MOI.delete(model::MOI.ModelLike,
     for ci in bridge.psd2x2
         MOI.delete(model, ci)
     end
-    for (vii, vij, vjj) in bridge.variables
-        MOI.delete(model, vii)
-        MOI.delete(model, vij)
-        MOI.delete(model, vjj)
-    end
 end

src/sos_polynomial_bridge.jl

@@ -90,7 +90,6 @@ function MOI.get(bridge::SOSPolynomialBridge,
     return MOI.get(bridge.variable_bridge, attr)
 end
 
-
 # Indices
 function MOI.delete(model::MOI.ModelLike, bridge::SOSPolynomialBridge)
     # First delete the constraints in which the Gram matrix appears
@@ -144,12 +143,13 @@ end
 function MOI.get(model::MOI.ModelLike,
                  attr::GramMatrixAttribute,
                  bridge::SOSPolynomialBridge)
-    return build_gram_matrix(MOI.get(model, attr, bridge.variable_bridge),
+    return build_gram_matrix(MOI.get(model, MOI.ConstraintPrimal(attr.N), bridge.variable_bridge),
                              bridge.certificate_monomials)
 end
 function MOI.get(model::MOI.ModelLike,
                  attr::MomentMatrixAttribute,
                  bridge::SOSPolynomialBridge)
-    return build_moment_matrix(MOI.get(model, attr, bridge.variable_bridge),
+    return build_moment_matrix(MOI.get(model, MOI.ConstraintDual(attr.N),
+                                       bridge.variable_bridge),
                                bridge.certificate_monomials)
 end

src/sos_polynomial_in_semialgebraic_set_bridge.jl

@@ -130,7 +130,7 @@ end
 function MOI.get(model::MOI.ModelLike, ::LagrangianMultipliers,
                  bridge::SOSPolynomialInSemialgebraicSetBridge)
     @assert eachindex(bridge.lagrangian_bridges) == eachindex(bridge.lagrangian_monomials)
-    map(i -> build_gram_matrix(MOI.get(model, GramMatrixAttribute(),
+    map(i -> build_gram_matrix(MOI.get(model, MOI.ConstraintPrimal(),
                                        bridge.lagrangian_bridges[i]),
                                bridge.lagrangian_monomials[i]),
         eachindex(bridge.lagrangian_bridges))

test/Mock/BPT12e399.jl

@@ -1,11 +1,11 @@
 config = MOI.Test.TestConfig()
-optimize!_max(mock) = MOIU.mock_optimize!(mock, [ 6.0, 1.0, -3.0, 9.0],
+optimize!_max(mock) = MOIU.mock_optimize!(mock, [ 6.0, 1.0, 9.0, -3.0*√2],
     (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[3, 1, 1/3]],
-    (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => [[3, 1/3, √2]]
+    (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone) => [[3, 1/3, √2]]
 )
-optimize!_min(mock) = MOIU.mock_optimize!(mock, [-6.0, 1.0,  3.0, 9.0],
+optimize!_min(mock) = MOIU.mock_optimize!(mock, [-6.0, 1.0, 9.0, 3.0*√2],
     (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[3, -1, 1/3]],
-    (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => [[3, 1/3, -√2]]
+    (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone) => [[3, 1/3, -√2]]
 )
 mock = bridged_mock(optimize!_max, optimize!_min)
 Tests.BPT12e399_test(mock, config)

test/Mock/bivariate_quadratic.jl

@@ -1,6 +1,6 @@
 config = MOI.Test.TestConfig()
-optimize!(mock) = MOIU.mock_optimize!(mock, [2.0, 1.0, 1.0, 1.0],
-    (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => [[1.0, 1.0, -√2]],
+optimize!(mock) = MOIU.mock_optimize!(mock, [2.0, 1.0, 1.0, √2],
+    (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone) => [[1.0, 1.0, -√2]],
     (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[1.0, -1.0, 1.0]])
 mock = bridged_mock(optimize!)
 Tests.sos_bivariate_quadratic_test(mock, config)

test/Mock/concave_then_convex_cubic.jl

@@ -9,15 +9,15 @@ Tests.dsos_concave_then_convex_cubic_test(mock, config)
 optimize!(mock) = MOIU.mock_optimize!(mock, [1.0; zeros(5); 3.0; zeros(3); 3.0; zeros(7); 3.0; zeros(3); 3.0; zeros(5)],
     (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) => zeros(4),
     (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [zeros(7), zeros(7)],
-    (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => [zeros(3) for i in 1:8],
+    (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone) => [zeros(3) for i in 1:8],
 )
 mock = bridged_mock(optimize!)
 Tests.sdsos_concave_then_convex_cubic_test(mock, config)
 optimize!(mock) = MOIU.mock_optimize!(mock, [1.0; zeros(5); 6.0; zeros(8); 6.0; zeros(6)],
     (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) => zeros(4),
     (MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle) => [zeros(6), zeros(6)],
     (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [zeros(7), zeros(7)],
-    (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => [zeros(3), zeros(3)],
+    (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone) => [zeros(3), zeros(3)],
 )
 mock = bridged_mock(optimize!)
 Tests.sos_concave_then_convex_cubic_test(mock, config)

test/Mock/lyapunov_switched_system.jl

@@ -12,7 +12,7 @@ config = MOI.Test.TestConfig(atol=1e-5, rtol=1e-5)
 #                 0   1+α]
 # P - A2'P*A*2 = [1+α 0
 #                 0   β]
-optimize!(mock) = MOIU.mock_optimize!(mock, [α, 0.0, α, 1 + α, 0.0, β, 1 + α, β, 0.0])
+optimize!(mock) = MOIU.mock_optimize!(mock, [α, α, 0.0, β, 1 + α, 0.0, 1 + α, β, 0.0])
 mock = bridged_mock(optimize!)
 Tests.quadratic_feasible_lyapunov_switched_system_test(mock, config)
 # TODO quadratic infeasible and quartic

test/Tests/bivariate_quadratic.jl

@@ -19,6 +19,8 @@ function bivariate_quadratic_test(optimizer,
     @objective(model, Max, α + 1)
     optimize!(model)
 
+    @test certificate_monomials(cref) == [x, 1]
+
     @test termination_status(model) == MOI.OPTIMAL
     @test objective_value(model) ≈ 3.0 atol=atol rtol=rtol
 

test/Tests/utilities.jl

@@ -136,6 +136,7 @@ end
 # Constraint dual values for inner bridged model
 function inner_inspect(model, atol=1e-4)
     inner = _inner(backend(model))
+    @show MOI.get(inner, MOI.NumberOfVariables())
     for (F, S) in MOI.get(inner, MOI.ListOfConstraints())
         @show F
         @show S