Contributed by: Benoît Legat
Consider the polynomial optimization problem (a variation from Lasserre2009; Example 2.2) of minimizing the polynomial x3−x2+2xy−y2+y3 over the polyhedron defined by the inequalities x≥0,y≥0 and x+y≥1.
using DynamicPolynomials
@polyvar x 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
Basic semialgebraic Set defined by no equality 3 inequalities x ≥ 0 y ≥ 0 -2 + y^2 + x^2 ≥ 0
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 for a list of the available choices.
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)
#12 (generic function with 1 method)
A Sum-of-Squares certificate that p≥α over the domain S, ensures that α is a lower bound to the polynomial optimization problem.
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)
Iter: 18 Ap: 9.59e-01 Pobj: -3.5512000e-03 Ad: 9.59e-01 Dobj: -3.5512523e-03
Success: SDP solved
Primal objective value: -3.5512000e-03
Dual objective value: -3.5512523e-03
Relative primal infeasibility: 1.38e-13
Relative dual infeasibility: 5.36e-10
Real Relative Gap: -5.20e-08
XZ Relative Gap: 1.22e-09
DIMACS error measures: 1.49e-13 0.00e+00 1.44e-09 0.00e+00 -5.20e-08 1.22e-09
------------------------------------------------------------------
SCS v3.2.7 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 11, constraints m: 20
cones: z: primal zero / dual free vars: 1
l: linear vars: 1
s: psd vars: 18, ssize: 3
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
max_iters: 100000, normalize: 1, rho_x: 1.00e-06
acceleration_lookback: 10, acceleration_interval: 10
compiled with openmp parallelization enabled
lin-sys: sparse-direct-amd-qdldl
nnz(A): 22, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 2.57e+01 1.01e+00 2.00e+06 -1.00e+06 1.00e-01 1.44e-04
25| 1.30e+11 4.45e+10 8.00e+18 -4.00e+18 1.00e-01 3.16e-04
------------------------------------------------------------------
status: unbounded
timings: total: 3.17e-04s = setup: 5.67e-05s + solve: 2.60e-04s
lin-sys: 1.44e-05s, cones: 1.63e-04s, accel: 1.46e-06s
------------------------------------------------------------------
objective = -inf
------------------------------------------------------------------
We can see that the problem is infeasible, meaning that no lower bound was found.
solution_summary(model)
* Solver : Dual model with SCS attached * Status Result count : 1 Termination status : INFEASIBLE Message from the solver: "unbounded" * Candidate solution (result #1) Primal status : INFEASIBLE_POINT Dual status : INFEASIBILITY_CERTIFICATE Objective value : NaN Dual objective value : -1.00000e+00 * Work counters Solve time (sec) : 3.16804e-04
We now define the Schmüdgen's certificate:
import StarAlgebras as SA
import MultivariateBases as MB
const SOS = SumOfSquares
const SOSC = SOS.Certificate
struct Schmüdgen{IC<:SOSC.AbstractIdealCertificate,CT<:SOS.SOSLikeCone,B<:SA.AbstractBasis} <: SOSC.AbstractPreorderCertificate
ideal_certificate::IC
cone::CT
basis::B
maxdegree::Int
end
SOSC.cone(certificate::Schmüdgen) = certificate.cone
function SOSC.preprocessed_domain(::Schmüdgen, domain::BasicSemialgebraicSet, p)
return SOSC.with_variables(domain, p)
end
function SOSC.preorder_indices(::Schmüdgen, domain::SOSC.WithVariables)
n = length(domain.inner.p)
if n >= Sys.WORD_SIZE
error("There are $(2^n - 1) products in Schmüdgen's certificate, they cannot even be indexed with `$Int`.")
end
return map(SOSC.PreorderIndex, 1:(2^n-1))
end
function SOSC.multiplier_basis(certificate::Schmüdgen, index::SOSC.PreorderIndex, domain::SOSC.WithVariables)
q = SOSC.generator(certificate, index, domain)
return SOSC.maxdegree_gram_basis(certificate.basis, variables(domain), SOSC.multiplier_maxdegree(certificate.maxdegree, q))
end
function SOSC.multiplier_basis_type(::Type{Schmüdgen{IC, CT, BT}}, ::Type{M}) where {IC,CT,BT,M}
return MB.explicit_basis_type(BT)
end
function SOSC.generator(::Schmüdgen, index::SOSC.PreorderIndex, domain::SOSC.WithVariables)
I = [i for i in eachindex(domain.inner.p) if !iszero(index.value & (1 << (i - 1)))]
return prod([domain.inner.p[i] for i in eachindex(domain.inner.p) if !iszero(index.value & (1 << (i - 1)))])
end
SOSC.ideal_certificate(certificate::Schmüdgen) = certificate.ideal_certificate
SOSC.ideal_certificate(::Type{<:Schmüdgen{IC}}) where {IC} = IC
SOS.matrix_cone_type(::Type{<:Schmüdgen{IC, CT}}) where {IC, CT} = SOS.matrix_cone_type(CT)
Let's try again with our the Schmüdgen certificate:
model = SOSModel(solver)
@variable(model, α)
@objective(model, Max, α)
basis = MB.FullBasis{MB.Monomial,typeof(x * y)}()
ideal_certificate = SOSC.Newton(SOSCone(), basis, basis, tuple())
certificate = Schmüdgen(ideal_certificate, SOSCone(), basis, maxdegree(p))
@constraint(model, c, p >= α, domain = S, certificate = certificate)
optimize!(model)
solution_summary(model)
------------------------------------------------------------------
SCS v3.2.7 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 15, constraints m: 49
cones: z: primal zero / dual free vars: 1
l: linear vars: 3
s: psd vars: 45, ssize: 5
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
max_iters: 100000, normalize: 1, rho_x: 1.00e-06
acceleration_lookback: 10, acceleration_interval: 10
compiled with openmp parallelization enabled
lin-sys: sparse-direct-amd-qdldl
nnz(A): 67, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 8.53e+00 5.66e-01 3.99e+01 -2.00e+01 1.00e-01 1.69e-04
250| 1.40e-03 5.14e-05 7.85e-06 8.28e-01 1.00e-01 3.95e-03
------------------------------------------------------------------
status: solved
timings: total: 3.96e-03s = setup: 7.40e-05s + solve: 3.88e-03s
lin-sys: 2.46e-04s, cones: 3.40e-03s, accel: 3.92e-05s
------------------------------------------------------------------
objective = 0.828440
------------------------------------------------------------------
* Solver : Dual model with SCS attached * Status Result count : 1 Termination status : OPTIMAL Message from the solver: "solved" * Candidate solution (result #1) Primal status : FEASIBLE_POINT Dual status : FEASIBLE_POINT Objective value : 8.28436e-01 Dual objective value : 8.28444e-01 * Work counters Solve time (sec) : 3.95599e-03
This notebook was generated using Literate.jl.