# Polynomial Optimization

## Introduction

Consider the polynomial optimization problem 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$.


```julia
using DynamicPolynomials
@polyvar x y
p = x^3 - x^2 + 2x*y -y^2 + y^3
using SemialgebraicSets
S = @set x >= 0 && y >= 0 && x + y >= 1
p(x=>1, y=>0), p(x=>1//2, y=>1//2), p(x=>0, y=>1)
```




    (0, 1//4, 0)



The optimal solutions are $(x, y) = (1, 0)$ and $(x, y) = (0, 1)$ with objective value $0$ but [Ipopt](https://github.com/JuliaOpt/Ipopt.jl/) only finds the local minimum $(1/2, 1/2)$ with objective value $1/4$.


```julia
using JuMP
using Ipopt
model = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
@variable(model, a >= 0)
@variable(model, b >= 0)
@constraint(model, a + b >= 1)
@NLobjective(model, Min, a^3 - a^2 + 2a*b - b^2 + b^3)
optimize!(model)

@show termination_status(model)
@show value(a)
@show value(b)
@show objective_value(model)
```

    
    ******************************************************************************
    This program contains Ipopt, a library for large-scale nonlinear optimization.
     Ipopt is released as open source code under the Eclipse Public License (EPL).
             For more information visit http://projects.coin-or.org/Ipopt
    ******************************************************************************
    
    termination_status(model) = LOCALLY_SOLVED::TerminationStatusCode = 4
    value(a) = 0.4999999966705758
    value(b) = 0.49999999667062867
    objective_value(model) = 0.24999999500590342





    0.24999999500590342



Note that the problem can be written equivalently as follows using [registered functions](http://www.juliaopt.org/JuMP.jl/v0.19.0/nlp/#User-defined-Functions-1).


```julia
using JuMP
using Ipopt
model = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
@variable(model, a >= 0)
@variable(model, b >= 0)
@constraint(model, a + b >= 1)
peval(a, b) = p(x=>a, y=>b)
register(model, :peval, 2, peval, autodiff=true)
@NLobjective(model, Min, peval(a, b))
optimize!(model)

@show termination_status(model)
@show value(a)
@show value(b)
@show objective_value(model)
```

    termination_status(model) = LOCALLY_SOLVED::TerminationStatusCode = 4
    value(a) = 0.499999995002878
    value(b) = 0.4999999950028773
    objective_value(model) = 0.24999999250431656





    0.24999999250431656



## Sum-of-Squares approach

We will now see how to find the optimal solution using Sum of Squares Programming.


```julia
using SumOfSquares
```

We first need to pick an SDP solver, see [here](http://www.juliaopt.org/JuMP.jl/dev/installation/#Getting-Solvers-1) for a list of the available choices.


```julia
using CSDP
optimizer = optimizer_with_attributes(CSDP.Optimizer, "printlevel" => 0);
```


```julia
using MosekTools
optimizer = optimizer_with_attributes(Mosek.Optimizer, "QUIET" => 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 program searches for the largest upper bound and finds zero.


```julia
model = SOSModel(optimizer)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c3, p >= α, domain = S)

optimize!(model)
@show termination_status(model)
@show objective_value(model)
```

    termination_status(model) = OPTIMAL::TerminationStatusCode = 1
    objective_value(model) = -2.0092710828478744e-10





    -2.0092710828478744e-10



Using the solution $(1/2, 1/2)$ found by Ipopt of objective value $1/4$
and this certificate of lower bound $0$ we know that the optimal objective value is in the interval $[0, 1/4]$
but we still do not know what it is (if we consider that we did not try the solutions $(1, 0)$ and $(0, 1)$ as done in the introduction).
If the dual of the constraint `c3` was atomic, its atoms would have given optimal solutions of objective value $0$ but that is not the case.


```julia
using MultivariateMoments
ν3 = moment_matrix(c3)
extractatoms(ν3, 1e-3) # Returns nothing as the dual is not atomic
```

Fortunately, there is a hierarchy of programs with increasingly better bounds that can be solved until we get one with atom dual variables.
This comes from the way the Sum-of-Squares constraint with domain `S` is formulated.
The polynomial $p - \alpha$ is guaranteed to be nonnegative over the domain `S` if there exists Sum-of-Squares polynomials $s_0$, $s_1$, $s_2$, $s_3$ such that
$$ p - \alpha = s_0 + s_1 x + s_2 y + s_3 (x + y - 1). $$
Indeed, in the domain `S`, $x$, $y$ and $x + y - 1$ are nonnegative so the right-hand side is a sum of squares hence is nonnegative.
Once the degrees of $s_1$, $s_2$ and $s_3$ have been decided, the degree needed for $s_0$ will be determined but we have a freedom in choosing the degrees of $s_1$, $s_2$ and $s_3$.
By default, they are chosen so that the degrees of $s_1 x$, $s_2 y$ and $s_3 (x + y - 1)$ match those of $p - \alpha$ but this can be overwritten using the $maxdegree$ keyword argument.

### The maxdegree keyword argument

The maximum total degree (i.e. maximum sum of the exponents of $x$ and $y$) of the monomials of $p$ is 3 so the constraint in the program above is equivalent to `@constraint(model, p >= α, domain = S, maxdegree = 3)`.
That is, since $x$, $y$ and $x + y - 1$ have total degree 1, the sum of squares polynomials $s_1$, $s_2$ and $s_3$ have been chosen with maximum total degree $2$.
Since these polynomials are sums of squares, their degree must be even so the next maximum total degree to try is 4.
For this reason, the keywords `maxdegree = 4` and `maxdegree = 5` have the same effect in this example.
In general, if the polynomials in the domain are not all odd or all even, each value of `maxdegree` has different effect in the choice of the maximum total degree of $s_i$.


```julia
model = SOSModel(optimizer)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c5, p >= α, domain = S, maxdegree = 5)

optimize!(model)
@show termination_status(model)
@show objective_value(model)
```

    termination_status(model) = OPTIMAL::TerminationStatusCode = 1
    objective_value(model) = 8.700006062680539e-8





    8.700006062680539e-8



This time, the dual variable is atomic as it is the moments of the measure
$$0.5 \delta(x-1, y) + 0.5 \delta(x, y-1)$$
where $\delta(x, y)$ is the dirac measure centered at $(0, 0)$.
Therefore the program provides both a certificate that $0$ is a lower bound and a certificate that it is also an upper bound since it is attained at the global minimizers $(1, 0)$ and $(0, 1)$.


```julia
using MultivariateMoments
ν5 = moment_matrix(c5)
extractatoms(ν5, 1e-3)
```




    Atomic measure on the variables x, y with 2 atoms:
     at [-0.000242397, 1.00024] with weight 0.4936307922624579
     at [0.999972, 2.76679e-5] with weight 0.5063248377957446



## A deeper look into atom extraction

The `extractatoms` function requires a `ranktol` argument that we have set to `1e-3` in the preceding section.
This argument is used to transform the dual variable into a system of polynomials equations whose solutions give the atoms.
This transformation uses the SVD decomposition of the moment matrix and discards the equations corresponding to a singular value lower than `ranktol`.
When this system of equation has an infinite number of solutions, `extractatoms` concludes that the measure is not atomic.
For instance, with `maxdegree = 3`, we obtain the system
$$x + y = 1$$
which contains a whole line of solution.
This explains `extractatoms` returned `nothing`.


```julia
ν3 = moment_matrix(c3)
MultivariateMoments.computesupport!(ν3, 1e-3)
```




    Algebraic Set defined by 1 equality
     -x - 0.9999999999999998*y + 1.0000000000569824 == 0




With `maxdegree = 5`, we obtain the system
\begin{align}
  x + y & = 1\\
  y^2 & = y\\
  xy & = 0\\
  x^2 + y & = 1
\end{align}


```julia
ν5 = moment_matrix(c5)
MultivariateMoments.computesupport!(ν5, 1e-3)
```




    Algebraic Set defined by 4 equalities
     -x - 0.9999999962812296*y + 1.0000000580930948 == 0
     -y^2 + 1.0002701267987304*y - 2.7674560303336815e-5 == 0
     -x*y - 1.293349230893601e-6*y + 3.3938404378219955e-5 == 0
     -x^2 - 1.0002601608230963*y + 1.0002330609868753 == 0




This system can be reduced to the equivalent system
\begin{align}
  x + y & = 1\\
  y^2 & = y
\end{align}
which has the solutions $(0, 1)$ and $(1, 0)$.


```julia
SemialgebraicSets.computegröbnerbasis!(ideal(ν5.support))
ν5.support
```




    Algebraic Set defined by 2 equalities
     x + 0.9999999962812296*y - 1.0000000580930948 == 0
     y^2 - 1.0002701267987304*y + 2.7674560303336815e-5 == 0




The function `extractatoms` then reuse the matrix of moments to find the weights $1/2$, $1/2$ corresponding to the diracs centered respectively at $(0, 1)$ and $(1, 0)$.
This details the how the function obtained the result
$$0.5 \delta(x-1, y) + 0.5 \delta(x, y-1)$$
given in the previous section.

b937aa0413a5dc3bd1b897d04643af46690e63ee
ipynb
Jupyter Notebook
examples/Polynomial_Optimization.ipynb
stephentu/SumOfSquares.jl
e781643d15fd5cf49b5b4e2cd45b484e0cc45b77
MIT
examples/Polynomial_Optimization.ipynb
stephentu/SumOfSquares.jl
e781643d15fd5cf49b5b4e2cd45b484e0cc45b77
MIT
examples/Polynomial_Optimization.ipynb
stephentu/SumOfSquares.jl
e781643d15fd5cf49b5b4e2cd45b484e0cc45b77
MIT
Qwen/Qwen-72B
1. YES
2. YES