Certificate

Contributed by: Benoît Legat

Introduction

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)
optimize!(model)

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.

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.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
    ideal_certificate::IC
    cone::CT
    basis::Type{BT}
    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}}) where {IC, CT, BT}
    return 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, α)
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)
solution_summary(model)

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

This page was generated using Literate.jl.