# Extracting minimizers

[L09] Laurent, Monique. Sums of squares, moment matrices and optimization over polynomials. Emerging applications of algebraic geometry (2009): 157-270.

## Introduction

Consider the polynomial optimization problem [L09, Example 6.23] of minimizing the linear function $-x_1 - x_2$ over the basic semialgebraic set defined by the inequalities $x_2 \le 2x_1^4 - 8x_1^3 + 8x_1^2 + 2$, $x_2 \le 4x_1^4 - 32x_1^3 + 88x_1^2 - 96x_1 + 36$ and the box constraints $0 \le x_1 \le 3$ and $0 \le x_2 \le 4$, World Scientific, 2009.

using DynamicPolynomials
@polyvar x[1:2]
p = -sum(x)
using SumOfSquares
K = @set x >= 0 && x <= 3 && x >= 0 && x <= 4 && x <= 2x^4 - 8x^3 + 8x^2 + 2 && x <= 4x^4 - 32x^3 + 88x^2 - 96x + 36
Basic semialgebraic Set defined by no equality
6 inequalities
x ≥ 0
3 - x ≥ 0
x ≥ 0
4 - x ≥ 0
2 - x + 8*x^2 - 8*x^3 + 2*x^4 ≥ 0
36 - x - 96*x + 88*x^2 - 32*x^3 + 4*x^4 ≥ 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.

import CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
MathOptInterface.OptimizerWithAttributes(CSDP.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.Silent() => true])

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 function searches for the largest lower bound and finds zero using the dth level of the hierarchy.

function solve(d)
model = SOSModel(solver)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c, p >= α, domain = K, maxdegree = d)
optimize!(model)
println(solution_summary(model))
return model
end
solve (generic function with 1 method)

The first level of the hierarchy gives a lower bound of -7

model4 = solve(4)
A JuMP Model
Maximization problem with:
Variable: 1
Objective function type: VariableRef
Vector{AffExpr}-in-SumOfSquares.SOSPolynomialSet{BasicSemialgebraicSet{Int64, DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Int64}, FullSpace}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, SumOfSquares.Certificate.Putinar{SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}}}: 1 constraint
Model mode: AUTOMATIC
CachingOptimizer state: ATTACHED_OPTIMIZER
Solver name: CSDP
Names registered in the model: c, α

The second level improves the lower bound

model5 = solve(5)
A JuMP Model
Maximization problem with:
Variable: 1
Objective function type: VariableRef
Vector{AffExpr}-in-SumOfSquares.SOSPolynomialSet{BasicSemialgebraicSet{Int64, DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Int64}, FullSpace}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, SumOfSquares.Certificate.Putinar{SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}}}: 1 constraint
Model mode: AUTOMATIC
CachingOptimizer state: ATTACHED_OPTIMIZER
Solver name: CSDP
Names registered in the model: c, α

The third level finds the optimal objective value as lower bound...

model7 = solve(7)
A JuMP Model
Maximization problem with:
Variable: 1
Objective function type: VariableRef
Vector{AffExpr}-in-SumOfSquares.SOSPolynomialSet{BasicSemialgebraicSet{Int64, DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Int64}, FullSpace}, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, SumOfSquares.Certificate.Putinar{SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}}}: 1 constraint
Model mode: AUTOMATIC
CachingOptimizer state: ATTACHED_OPTIMIZER
Solver name: CSDP
Names registered in the model: c, α

...and proves it by exhibiting the minimizer.

ν7 = moment_matrix(model7[:c])
η = atomic_measure(ν7, 1e-3) # Returns nothing as the dual is not atomic
Atomic measure on the variables x, x with 1 atoms:
at [2.32952019617943, 3.1784930775648848] with weight 0.9999999990880203

We can indeed verify that the objective value at x_opt is equal to the lower bound.

x_opt = η.atoms.center
p(x_opt)
-5.508013273744314