Contributed by: Benoît Legat


Consider the polynomial optimization problem (a variation from [L09, Example 2.2]) of minimizing the polynomial $x^3 - x^2 + 2xy - y^2 + y^3$ over the polyhedron defined by the inequalities $x \ge 0, y \ge 0$ and $x + y \geq 1$.

[L09] Lasserre, J. B. Moments, positive polynomials and their applications. World Scientific, 2009.

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)
#11 (generic function with 1 method)

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 lower bound.

model = SOSModel(solver)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c, p >= α, domain = S)
Success: SDP solved
Primal objective value: 0.0000000e+00
Dual objective value: 0.0000000e+00
Relative primal infeasibility: 4.15e-17
Relative dual infeasibility: 5.00e-11
Real Relative Gap: 0.00e+00
XZ Relative Gap: 3.63e-10
DIMACS error measures: 6.34e-17 0.00e+00 5.00e-11 0.00e+00 0.00e+00 3.63e-10
	       SCS v3.2.3 - 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
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.38e-04
    25| 1.30e+11  4.45e+10  8.00e+18 -4.00e+18  1.00e-01  3.02e-04
status:  unbounded
timings: total: 3.03e-04s = setup: 7.95e-05s + solve: 2.24e-04s
	 lin-sys: 1.52e-05s, cones: 1.72e-04s, accel: 1.70e-06s
objective = -inf

We can see that the problem is infeasible, meaning that no lower bound was found.

* Solver : Dual model with SCS attached

* Status
  Result count       : 1
  Termination status : INFEASIBLE
  Message from the solver:

* Candidate solution (result #1)
  Primal status      : INFEASIBLE_POINT
  Objective value    : NaN
  Dual objective value : -1.00000e+00

* Work counters
  Solve time (sec)   : 3.03212e-04

We now define the Schmüdgen's certificate:

import MultivariateBases as MB
const SOS = SumOfSquares
const SOSC = SOS.Certificate
struct Schmüdgen{IC <: SOSC.AbstractIdealCertificate, CT <: SOS.SOSLikeCone, BT <: SOS.AbstractPolynomialBasis} <: SOSC.AbstractPreorderCertificate

SOSC.cone(certificate::Schmüdgen) = certificate.cone

function SOSC.preprocessed_domain(::Schmüdgen, domain::BasicSemialgebraicSet, p)
    return SOSC.with_variables(domain, p)

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`.")
    return map(SOSC.PreorderIndex, 1:(2^n-1))

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))
function SOSC.multiplier_basis_type(::Type{Schmüdgen{IC, CT, BT}}) where {IC, CT, BT}
    return BT

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)))])

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, α)
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)
* Solver : Dual model with SCS attached

* Status
  Result count       : 1
  Termination status : OPTIMAL
  Message from the solver:

* 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)   : 4.71127e-03

This page was generated using Literate.jl.