Certificate

Contributed by: Benoît Legat

Introduction

Consider the polynomial optimization problem (a variation from (Lasserre, 2009; 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$.

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 Clarabel 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 Clarabel solves the dual problem.

import Clarabel
using Dualization
solver = dual_optimizer(Clarabel.Optimizer)

#7 (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)

-------------------------------------------------------------
           Clarabel.jl v0.11.1  -  Clever Acronym

                   (c) Paul Goulart
                University of Oxford, 2022
-------------------------------------------------------------

problem:
  variables     = 11
  constraints   = 20
  nnz(P)        = 0
  nnz(A)        = 22
  cones (total) = 5
    : Zero        = 1,  numel = 1
    : Nonnegative = 1,  numel = 1
    : PSDTriangle = 3,  numel = (6,6,6)

settings:
  linear algebra: direct / qdldl, precision: 64 bit (1 thread)
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
  static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step
---------------------------------------------------------------------------------------------
  0  -4.9489e-01  -4.9489e-01  9.61e-09  7.64e-01  6.51e-01  1.00e+00  3.40e+00   ------
  1  -2.2460e+01  -1.5003e+01  4.97e-01  6.98e-02  6.09e-02  8.46e+00  4.09e-01  8.87e-01
  2  -2.8137e+02  -1.9564e+02  4.38e-01  5.09e-03  4.78e-03  8.67e+01  4.59e-02  9.53e-01
  3  -3.5303e+03  -2.5774e+03  3.70e-01  4.40e-04  3.97e-04  9.54e+02  4.22e-03  9.13e-01
  4  -3.3724e+04  -2.6710e+04  2.63e-01  8.10e-05  5.59e-05  7.02e+03  7.25e-04  9.69e-01
  5  -8.0048e+04  -5.2489e+04  5.25e-01  2.33e-05  2.02e-05  2.76e+04  2.41e-04  9.09e-01
  6  -8.0840e+05  -5.1919e+05  5.57e-01  2.27e-06  2.00e-06  2.89e+05  2.37e-05  9.43e-01
  7  -4.5856e+06  -3.1944e+06  4.35e-01  5.75e-07  4.21e-07  1.39e+06  5.68e-06  9.90e-01
  8  -2.4991e+07  -1.7798e+07  4.04e-01  1.17e-07  8.16e-08  7.19e+06  1.11e-06  8.19e-01
  9  -2.6734e+07  -7.8819e+06  2.39e+00  3.63e-08  4.30e-08  1.89e+07  4.92e-07  9.27e-01
 10  -4.1111e+07  -1.4996e+07  1.74e+00  2.67e-08  3.02e-08  2.61e+07  3.57e-07  3.24e-01
 11  -4.5308e+07  -1.1821e+07  2.83e+00  2.20e-08  2.64e-08  3.35e+07  2.98e-07  3.55e-01
 12  -1.2910e+08  -4.2257e+07  2.06e+00  8.44e-09  9.66e-09  8.68e+07  1.13e-07  6.58e-01
 13  -7.7076e+08  -5.5086e+08  3.99e-01  5.61e-09  3.08e-09  2.20e+08  7.48e-08  9.28e-01
 14  -4.1029e+09  -2.8929e+09  4.18e-01  1.25e-09  6.46e-10  1.21e+09  1.61e-08  7.95e-01
---------------------------------------------------------------------------------------------
Terminated with status = dual infeasible
solve time = 2.50ms

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

solution_summary(model)

solution_summary(; result = 1, verbose = false)
├ solver_name          : Dual model with Clarabel attached
├ Termination
│ ├ termination_status : INFEASIBLE
│ ├ result_count       : 1
│ └ raw_status         : DUAL_INFEASIBLE
└ Solution (result = 1)
  ├ primal_status        : INFEASIBLE_POINT
  ├ dual_status          : INFEASIBILITY_CERTIFICATE
  ├ objective_value      : -2.89286e+09
  └ dual_objective_value : -4.10292e+09

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}(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)

solution_summary(; result = 1, verbose = false)
├ solver_name          : Dual model with Clarabel attached
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 1
│ └ raw_status         : SOLVED
└ Solution (result = 1)
  ├ primal_status        : FEASIBLE_POINT
  ├ dual_status          : FEASIBLE_POINT
  ├ objective_value      : 8.28427e-01
  └ dual_objective_value : 8.28427e-01

This page was generated using Literate.jl.