Variable-wise bounds for putinar (#282)

* Variable-wise bounds for putinar

* Fixes

* Fix

* Fixes

* Exclude sparsity tests

* Fixes

* Fix format

* Fix

* Fixes

* Fix

* Update solver failing tests

* Fix format

* Reenable CSDP tests

They are fixed by jump-dev/CSDP.jl#89

* Fixes

* Fix certificate example

CSDP fails because an empty constraint is added

* Remove polynomial optimization tutorial

* Fix doc

Author: blegat

docs/Project.toml

@@ -7,6 +7,7 @@ Cyclotomics = "da8f5974-afbb-4dc8-91d8-516d5257c83b"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Dualization = "191a621a-6537-11e9-281d-650236a99e60"
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 GroupsCore = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
@@ -21,6 +22,7 @@ MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 PermutationGroups = "8bc5a954-2dfc-11e9-10e6-cd969bffa420"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1"
 SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2"

docs/src/tutorials/Extension/certificate.jl

@@ -17,26 +17,40 @@
 using Test #src
 using DynamicPolynomials
 @polyvar x y
-p = x^3 - x^2 + 2x*y -y^2 + y^3 + x^3 * y
+p = x^3 - x^2 + 2x*y - y^2 + y^3 + x^3 * y
 using SumOfSquares
 S = @set x >= 0 && y >= 0 && x^2 + y^2 >= 2
 
 # We will now see how to find the optimal solution using Sum of Squares Programming.
 # We first need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v1.8/installation/#Supported-solvers) for a list of the available choices.
-
-import CSDP
-solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
+# Note that SumOfSquares generates a *standard form* SDP (i.e., SDP variables
+# and equality constraints) while SCS expects a *geometric form* SDP (i.e.,
+# free variables and symmetric matrices depending affinely on these variables
+# constrained to belong to the PSD cone).
+# JuMP will transform the standard from to the geometric form will create the PSD
+# variables as free variables and then constrain then to be PSD.
+# While this will work, since the dual of a standard from is in in geometric form,
+# dualizing the problem will generate a smaller SDP.
+# We use therefore `Dualization.dual_optimizer` so that SCS solves the dual problem.
+
+import SCS
+using Dualization
+solver = dual_optimizer(SCS.Optimizer)
 
 # A Sum-of-Squares certificate that $p \ge \alpha$ over the domain `S`, ensures that $\alpha$ is a lower bound to the polynomial optimization problem.
-# The following program searches for the largest upper bound and finds zero.
+# The following program searches for the largest lower bound.
 
 model = SOSModel(solver)
 @variable(model, α)
 @objective(model, Max, α)
 @constraint(model, c, p >= α, domain = S)
 optimize!(model)
-@show termination_status(model)
-@show objective_value(model)
+
+# We can see that the problem is infeasible, meaning that no lower bound was found.
+
+@test termination_status(model) == MOI.INFEASIBLE #src
+@test dual_status(model) == MOI.INFEASIBILITY_CERTIFICATE #src
+solution_summary(model)
 
 # We now define the Schmüdgen's certificate:
 
@@ -91,6 +105,8 @@ ideal_certificate = SOSC.Newton(SOSCone(), MB.MonomialBasis, tuple())
 certificate = Schmüdgen(ideal_certificate, SOSCone(), MB.MonomialBasis, maxdegree(p))
 @constraint(model, c, p >= α, domain = S, certificate = certificate)
 optimize!(model)
-@show termination_status(model)
-@show objective_value(model)
+@test termination_status(model) == MOI.OPTIMAL #src
+@test primal_status(model) == MOI.FEASIBLE_POINT #src
+@test objective_value(model) ≈ 0.8284 rtol=1e-3 #src
 @test length(lagrangian_multipliers(c)) == 7 #src
+solution_summary(model)

docs/src/tutorials/Polynomial Optimization/polynomial_optimization.jl

@@ -1,13 +1,14 @@
-struct EmptyBridge{T} <: MOI.Bridges.Constraint.AbstractBridge end
+struct EmptyBridge{T,F<:MOI.AbstractVectorFunction} <:
+       MOI.Bridges.Constraint.AbstractBridge end
 
 function MOI.Bridges.Constraint.bridge_constraint(
-    ::Type{EmptyBridge{T}},
+    ::Type{EmptyBridge{T,F}},
     model::MOI.ModelLike,
-    f::MOI.AbstractVectorFunction,
+    f::F,
     s::SOS.EmptyCone,
-) where {T}
+) where {T,F}
     @assert MOI.output_dimension(f) == MOI.dimension(s)
-    return EmptyBridge{T}()
+    return EmptyBridge{T,F}()
 end
 
 function MOI.supports_constraint(
@@ -25,14 +26,28 @@ function MOI.Bridges.added_constraint_types(::Type{<:EmptyBridge})
 end
 function MOI.Bridges.Constraint.concrete_bridge_type(
     ::Type{<:EmptyBridge{T}},
-    ::Type{<:MOI.AbstractVectorFunction},
+    ::Type{F},
     ::Type{SOS.EmptyCone},
-) where {T}
-    return EmptyBridge{T}
+) where {T,F<:MOI.AbstractVectorFunction}
+    return EmptyBridge{T,F}
 end
 
 # Indices
-function MOI.delete(model::MOI.ModelLike, bridge::EmptyBridge) end
+function MOI.delete(::MOI.ModelLike, ::EmptyBridge) end
+
+function _empty_function(::Type{MOI.VectorOfVariables})
+    return MOI.VectorOfVariables(MOI.VariableIndex[])
+end
+function _empty_function(::Type{F}) where {F}
+    return MOI.Utilities.zero_with_output_dimension(F, 0)
+end
+function MOI.get(
+    ::MOI.ModelLike,
+    ::MOI.ConstraintFunction,
+    ::EmptyBridge{T,F},
+) where {T,F}
+    return _empty_function(F)
+end
 
 function MOI.get(
     ::MOI.ModelLike,

src/Bridges/Constraint/empty.jl

@@ -161,10 +161,16 @@ function MOI.get(
 end
 
 # Indices
-_delete_variables(model, Q::Vector{MOI.VariableIndex}) = MOI.delete(model, Q)
+function _delete_variables(model, Q::Vector{MOI.VariableIndex})
+    if !isempty(Q)
+        # FIXME Since there is not variables in the list, we cannot
+        # identify the `EmptyBridge` to delete
+        MOI.delete(model, Q)
+    end
+end
 function _delete_variables(model, Qs::Vector{Vector{MOI.VariableIndex}})
     for Q in Qs
-        MOI.delete(model, Q)
+        _delete_variables(model, Q)
     end
 end
 function MOI.delete(model::MOI.ModelLike, bridge::SOSPolynomialBridge)

src/Bridges/Constraint/sos_polynomial.jl

@@ -8,8 +8,8 @@
 
 A matrix is SDD iff it is the sum of psd matrices Mij that are zero except
 for entries ii, ij and jj [Lemma 9, AM17]. This bridge substitute the
-constrained variables in [`SumOfSquares.ScaledDiagonallyDominantConeTriangle`](@ref)
-into a sum of constrained variables in [`SumOfSquares.PositiveSemidefinite2x2ConeTriangle`](@ref).
+constrained variables in [`SOS.ScaledDiagonallyDominantConeTriangle`](@ref)
+into a sum of constrained variables in [`SOS.PositiveSemidefinite2x2ConeTriangle`](@ref).
 
 [AM17] Ahmadi, A. A. & Majumdar, A.
 *DSOS and SDSOS Optimization: More Tractable Alternatives to Sum of Squares and Semidefinite Optimization*