Merge pull request #108 from JuliaOpt/bl/moments

Support MomentsAttribute

Author: blegat

Project.toml

@@ -21,7 +21,7 @@ JuMP = "0.19"
 MathOptInterface = "~0.8.2"
 MultivariateMoments = "≥ 0.2.1"
 MultivariatePolynomials = "~0.2.7"
-PolyJuMP = "0.3"
+PolyJuMP = "0.3.1"
 SemialgebraicSets = "0.1"
 julia = "1"
 

src/attributes.jl

@@ -10,7 +10,7 @@ struct CertificateMonomials <: MOI.AbstractConstraintAttribute end
     struct GramMatrixAttribute <: MOI.AbstractConstraintAttribute end
 
 A constraint attribute for the [`GramMatrix`](@ref) of a constraint, that is,
-the positive semidefinte matrix `Q` indexed by the monomials in the vector `X`
+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
 """

src/sos_polynomial_bridge.jl

@@ -7,6 +7,8 @@ struct SOSPolynomialBridge{T, F <: MOI.AbstractVectorFunction,
     variable_bridge::VBS
     certificate_monomials::MVT
     zero_constraint::MOI.ConstraintIndex{F, PolyJuMP.ZeroPolynomialSet{DT, BT, MT, MVT}}
+    domain::DT
+    monomials::MVT
 end
 
 function SOSPolynomialBridge{T, F, DT, VBS, MCT, BT, MT, MVT}(
@@ -30,7 +32,7 @@ function SOSPolynomialBridge{T, F, DT, VBS, MCT, BT, MT, MVT}(
     coefs = MOIU.vectorize(MP.coefficients(q))
     zero_constraint = MOI.add_constraint(model, coefs, set)
     return SOSPolynomialBridge{T, F, DT, VBS, MCT, BT, MT, MVT}(
-        variable_bridge, X, zero_constraint)
+        variable_bridge, X, zero_constraint, s.domain, s.monomials)
 end
 
 function MOI.supports_constraint(::Type{SOSPolynomialBridge{T}},
@@ -116,16 +118,25 @@ function MOI.get(::MOI.ModelLike,
                  ::SOSPolynomialBridge)
     throw(ValueNotSupported())
 end
-function MOI.get(::MOI.ModelLike,
-                 ::MOI.ConstraintDual,
-                 ::SOSPolynomialBridge)
-    throw(DualNotSupported())
+
+function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual,
+                 bridge::SOSPolynomialBridge)
+    dual = MOI.get(model, attr, bridge.zero_constraint)
+    set = MOI.get(model, MOI.ConstraintSet(), bridge.zero_constraint)
+    μ = measure(dual, set.monomials)
+    I = ideal(bridge.domain)
+    return [dot(rem(mono, I), μ) for mono in bridge.monomials]
 end
 function MOI.get(model::MOI.ModelLike,
                  attr::MOI.ConstraintDual,
                  bridge::SOSPolynomialBridge{T, <:MOI.AbstractVectorFunction, FullSpace}) where {T}
     return MOI.get(model, attr, bridge.zero_constraint)
 end
+function MOI.get(model::MOI.ModelLike, attr::PolyJuMP.MomentsAttribute,
+                 bridge::SOSPolynomialBridge)
+    return MOI.get(model, attr, bridge.zero_constraint)
+end
+
 function MOI.get(::MOI.ModelLike, ::CertificateMonomials,
                  bridge::SOSPolynomialBridge)
     return bridge.certificate_monomials

src/sos_polynomial_in_semialgebraic_set_bridge.jl

@@ -26,6 +26,7 @@ struct SOSPolynomialInSemialgebraicSetBridge{
     lagrangian_monomials::Vector{MVT}
     lagrangian_bridges::Vector{VBS}
     constraint::MOI.ConstraintIndex{F, SOSPolynomialSet{DT, CT, BT, MT, MVT, NPT}}
+    monomials::MVT
 end
 
 function SOSPolynomialInSemialgebraicSetBridge{T, F, DT, CT, VBS, BT, MT, MVT, NPT}(
@@ -55,7 +56,7 @@ function SOSPolynomialInSemialgebraicSetBridge{T, F, DT, CT, VBS, BT, MT, MVT, N
                                     new_set)
 
     return SOSPolynomialInSemialgebraicSetBridge{
-        T, F, DT, CT, VBS, BT, MT, MVT, NPT}(λ_monos, λ_bridges, constraint)
+        T, F, DT, CT, VBS, BT, MT, MVT, NPT}(λ_monos, λ_bridges, constraint, set.monomials)
 end
 
 function MOI.supports_constraint(::Type{SOSPolynomialInSemialgebraicSetBridge{T}},
@@ -107,10 +108,17 @@ function MOI.get(::MOI.ModelLike,
                  ::SOSPolynomialInSemialgebraicSetBridge)
     throw(ValueNotSupported())
 end
-function MOI.get(::MOI.ModelLike,
-                 ::MOI.ConstraintDual,
-                 ::SOSPolynomialInSemialgebraicSetBridge)
-    throw(DualNotSupported())
+
+function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual,
+                 bridge::SOSPolynomialInSemialgebraicSetBridge)
+    dual = MOI.get(model, attr, bridge.constraint)
+    set = MOI.get(model, MOI.ConstraintSet(), bridge.constraint)
+    μ = measure(dual, set.monomials)
+    return [dot(mono, μ) for mono in bridge.monomials]
+end
+function MOI.get(model::MOI.ModelLike, attr::PolyJuMP.MomentsAttribute,
+                 bridge::SOSPolynomialInSemialgebraicSetBridge)
+    return MOI.get(model, attr, bridge.constraint)
 end
 
 function MOI.get(model::MOI.ModelLike,

test/BPT12e399.jl

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

test/Mock/BPT12e399.jl

@@ -1,10 +1,23 @@
 config = MOI.Test.TestConfig()
-optimize!(mock) = MOIU.mock_optimize!(mock, [1.0; zeros(5); 6.0; zeros(12); 6.0; zeros(10)])
+optimize!(mock) = MOIU.mock_optimize!(mock, [1.0; zeros(5); 6.0; zeros(12); 6.0; zeros(10)],
+    (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => zeros(26),
+    (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) => zeros(4),
+    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [zeros(7), zeros(7)],
+)
 mock = bridged_mock(optimize!)
 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)])
+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],
+)
 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)])
+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)],
+)
 mock = bridged_mock(optimize!)
 Tests.sos_concave_then_convex_cubic_test(mock, config)

test/Mock/concave_then_convex_cubic.jl

@@ -18,3 +18,6 @@ end
 @testset "Lyapunov Switched System" begin
     include("lyapunov_switched_system.jl")
 end
+@testset "BPT12e399" begin
+    include("BPT12e399.jl")
+end

test/Mock/mock_tests.jl

@@ -0,0 +1,94 @@
+# Adapted from Example 3.99 of [BPT12].
+#
+# [BPT12] Blekherman, G.; Parrilo, P. & Thomas, R.
+# Semidefinite Optimization and Convex Algebraic Geometry
+# Society for Industrial and Applied Mathematics, 2012
+
+using Test
+using SumOfSquares
+using DynamicPolynomials
+
+function BPT12e399_test(optimizer, config::MOIT.TestConfig)
+    atol = config.atol
+    rtol = config.rtol
+
+    model = _model(optimizer)
+
+    @variable(model, α)
+
+    @polyvar x y
+    cref = @constraint(model, 10 - (x^2 + α*y) in SOSCone(),
+                       domain = @set x^2 + y^2 == 1)
+
+    @objective(model, Max, α)
+    JuMP.optimize!(model)
+
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ 6.0 atol=atol rtol=rtol
+
+    @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
+    @test JuMP.value(α) ≈ 6 atol=1e-6
+
+    @test_throws SumOfSquares.ValueNotSupported value(cref)
+    p = gram_matrix(cref)
+    @test getmat(p) ≈ [1 -3; -3 9] atol=atol rtol=rtol
+    @test p.x == [y, 1]
+
+    @test dual_status(model) == MOI.FEASIBLE_POINT
+    μ = dual(cref)
+    @test μ isa AbstractMeasure{Float64}
+    @test length(moments(μ)) == 3
+    @test moment_value(moments(μ)[1]) ≈ -8/3 atol=atol rtol=rtol
+    @test monomial(moments(μ)[1]) == x^2
+    @test moment_value(moments(μ)[2]) ≈ 1 atol=atol rtol=rtol
+    @test monomial(moments(μ)[2]) == y
+    @test moment_value(moments(μ)[3]) ≈ 1/3 atol=atol rtol=rtol
+    @test monomial(moments(μ)[3]) == 1
+
+    μ = moments(cref)
+    @test μ isa AbstractMeasure{Float64}
+    @test length(moments(μ)) == 3
+    @test moment_value(moments(μ)[1]) ≈ 3 atol=atol rtol=rtol
+    @test monomial(moments(μ)[1]) == y^2
+    @test moment_value(moments(μ)[2]) ≈ 1 atol=atol rtol=rtol
+    @test monomial(moments(μ)[2]) == y
+    @test moment_value(moments(μ)[3]) ≈ 1/3 atol=atol rtol=rtol
+    @test monomial(moments(μ)[3]) == 1
+
+    @objective(model, Min, α)
+    JuMP.optimize!(model)
+
+    @test termination_status(model) == MOI.OPTIMAL
+    @test objective_value(model) ≈ -6.0 atol=atol rtol=rtol
+
+    @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
+    @test JuMP.value(α) ≈ -6 atol=1e-6
+
+    @test_throws SumOfSquares.ValueNotSupported value(cref)
+    p = gram_matrix(cref)
+    @test getmat(p) ≈ [1 3; 3 9] atol=atol rtol=rtol
+    @test p.x == [y, 1]
+
+    @test dual_status(model) == MOI.FEASIBLE_POINT
+    μ = dual(cref)
+    @test μ isa AbstractMeasure{Float64}
+    @test length(moments(μ)) == 3
+    @test moment_value(moments(μ)[1]) ≈ -8/3 atol=atol rtol=rtol
+    @test monomial(moments(μ)[1]) == x^2
+    @test moment_value(moments(μ)[2]) ≈ -1 atol=atol rtol=rtol
+    @test monomial(moments(μ)[2]) == y
+    @test moment_value(moments(μ)[3]) ≈ 1/3 atol=atol rtol=rtol
+    @test monomial(moments(μ)[3]) == 1
+
+    μ = moments(cref)
+    @test μ isa AbstractMeasure{Float64}
+    @test length(moments(μ)) == 3
+    @test moment_value(moments(μ)[1]) ≈ 3 atol=atol rtol=rtol
+    @test monomial(moments(μ)[1]) == y^2
+    @test moment_value(moments(μ)[2]) ≈ -1 atol=atol rtol=rtol
+    @test monomial(moments(μ)[2]) == y
+    @test moment_value(moments(μ)[3]) ≈ 1/3 atol=atol rtol=rtol
+    @test monomial(moments(μ)[3]) == 1
+end
+
+sd_tests["BPT12e399"] = BPT12e399_test

test/Tests/BPT12e399.jl

@@ -2,33 +2,16 @@ module Tests
 
 include("utilities.jl")
 
-include("term_test.jl")
-include("bivariate_quadratic_test.jl")
-include("concave_then_convex_cubic_test.jl")
-include("horn_test.jl")
-
-const linear_tests = Dict(
-    "dsos_term" => dsos_term_test,
-    "dsos_bivariate_quadratic" => dsos_bivariate_quadratic_test,
-    "dsos_concave_then_convex_cubic" => dsos_concave_then_convex_cubic_test,
-    "dsos_horn" => dsos_horn_test
-)
-
-const soc_tests = Dict(
-    "sdsos_term" => sdsos_term_test,
-    "sdsos_bivariate_quadratic" => sdsos_bivariate_quadratic_test,
-    "sdsos_concave_then_convex_cubic" => sdsos_concave_then_convex_cubic_test,
-    "sdsos_horn" => sdsos_horn_test
-)
-
-const sd_tests = Dict(
-    "sos_term" => sos_term_test,
-    "sos_bivariate_quadratic" => sos_bivariate_quadratic_test,
-    "sos_concave_then_convex_cubic" => sos_concave_then_convex_cubic_test,
-    "sos_horn" => sos_horn_test,
-)
-
+const linear_tests = Dict{String, Function}()
+const soc_tests = Dict{String, Function}()
+const sd_tests = Dict{String, Function}()
+
+include("term.jl")
+include("bivariate_quadratic.jl")
+include("horn.jl")
+include("concave_then_convex_cubic.jl")
 include("lyapunov_switched_system.jl")
+include("BPT12e399.jl")
 
 @test_suite linear
 @test_suite soc

test/Tests/Tests.jl

@@ -30,14 +30,17 @@ function bivariate_quadratic_test(optimizer,
     @test getmat(p) ≈ ones(2, 2) atol=atol rtol=rtol
     @test p.x == [x, 1]
 
+    a = moment_value.(moments(dual(cref)))
+    @test a[2] ≈ -1.0 atol=atol rtol=rtol
+    @test a[1] + a[3] ≈ 2.0 atol=atol rtol=rtol
+
     @test dual_status(model) == MOI.FEASIBLE_POINT
-    μ = dual(cref)
-    @test μ isa AbstractMeasure{Float64}
-    @test length(moments(μ)) == 3
-    a = moment_value.(moments(μ))
-    @test a[2] ≈ -1.0  atol=atol rtol=rtol
-    @test a[1] + a[3] ≈ 2.0  atol=atol rtol=rtol
-    @test monomial.(moments(μ)) == [x^2, x, 1]
+    for μ in [dual(cref), moments(cref)]
+        @test μ isa AbstractMeasure{Float64}
+        @test length(moments(μ)) == 3
+        @test a ≈ moment_value.(moments(μ)) atol=atol rtol=rtol
+        @test monomial.(moments(μ)) == [x^2, x, 1]
+    end
 
     ν = moment_matrix(cref)
     @test getmat(ν) ≈ [a[1] a[2]
@@ -60,5 +63,8 @@ function bivariate_quadratic_test(optimizer,
           0)))
 end
 sos_bivariate_quadratic_test(optimizer, config)   = bivariate_quadratic_test(optimizer, config, SOSCone())
+sd_tests["sos_bivariate_quadratic"] = sos_bivariate_quadratic_test
 sdsos_bivariate_quadratic_test(optimizer, config) = bivariate_quadratic_test(optimizer, config, SDSOSCone())
+soc_tests["sdsos_bivariate_quadratic"] = sdsos_bivariate_quadratic_test
 dsos_bivariate_quadratic_test(optimizer, config)  = bivariate_quadratic_test(optimizer, config, DSOSCone())
+linear_tests["dsos_bivariate_quadratic"] = dsos_bivariate_quadratic_test

test/Tests/bivariate_quadratic_test.jl renamed to test/Tests/bivariate_quadratic.jl

@@ -29,19 +29,30 @@ function concave_then_convex_cubic_test(optimizer, config::MOIT.TestConfig,
     @test_throws SumOfSquares.ValueNotSupported value(cref_convex)
     @test_throws SumOfSquares.ValueNotSupported value(cref_concave)
 
-    @test_throws SumOfSquares.DualNotSupported dual(cref_convex)
-    @test_throws SumOfSquares.DualNotSupported dual(cref_concave)
+    # The monomials contain the variables created for the Hessian so we cannot
+    # check them easily
+    for μ in [dual(cref_convex), dual(cref_concave)]
+        @test μ isa AbstractMeasure{Float64}
+        @test length(moments(μ)) == 2
+    end
+    for μ in [moments(cref_convex), moments(cref_concave)]
+        @test μ isa AbstractMeasure{Float64}
+        @test length(moments(μ)) == 7
+    end
 end
 
 function sos_concave_then_convex_cubic_test(optimizer, config)
     concave_then_convex_cubic_test(optimizer, config,
                                    MOI.PositiveSemidefiniteConeTriangle)
 end
+sd_tests["sos_concave_then_convex_cubic"] = sos_concave_then_convex_cubic_test
 function sdsos_concave_then_convex_cubic_test(optimizer, config)
     concave_then_convex_cubic_test(optimizer, config,
                                    SumOfSquares.ScaledDiagonallyDominantConeTriangle)
 end
+soc_tests["sdsos_concave_then_convex_cubic"] = sdsos_concave_then_convex_cubic_test
 function dsos_concave_then_convex_cubic_test(optimizer, config)
     concave_then_convex_cubic_test(optimizer, config,
                                    SumOfSquares.DiagonallyDominantConeTriangle)
 end
+linear_tests["dsos_concave_then_convex_cubic"] = dsos_concave_then_convex_cubic_test

test/Tests/concave_then_convex_cubic_test.jl renamed to test/Tests/concave_then_convex_cubic.jl

@@ -59,5 +59,8 @@ function horn_test(optimizer,
     end
 end
 sos_horn_test(optimizer, config)   = horn_test(optimizer, config, SOSCone())
+sd_tests["sos_horn"] = sos_horn_test
 sdsos_horn_test(optimizer, config) = horn_test(optimizer, config, SDSOSCone())
+soc_tests["sdsos_horn"] = sdsos_horn_test
 dsos_horn_test(optimizer, config)  = horn_test(optimizer, config, DSOSCone())
+linear_tests["dsos_horn"] = dsos_horn_test