Polynomial Optimization
Contributed by: Benoît Legat
Introduction
Consider the polynomial optimization problem (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
using SumOfSquares
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)Local search
A local solver only uses the local information given by the the value, gradient and Hessian of the objective function and constraints at a given solution. When it converges, it therefore only guarantees that the found solution is a local minimum. In this example, the optimal solutions are $(x, y) = (1, 0)$ and $(x, y) = (0, 1)$ with objective value $0$ but Ipopt only finds the local minimum $(1/2, 1/2)$ with objective value $1/4$.
import Ipopt
model = Model(Ipopt.Optimizer)
@variable(model, a >= 0)
@variable(model, b >= 0)
@constraint(model, a + b >= 1)
@objective(model, Min, a^3 - a^2 + 2a*b - b^2 + b^3)
optimize!(model)Success: SDP solved
Primal objective value: 0.0000000e+00
Dual objective value: 0.0000000e+00
Relative primal infeasibility: 5.23e-17
Relative dual infeasibility: 5.00e-11
Real Relative Gap: 0.00e+00
XZ Relative Gap: 1.00e-10
DIMACS error measures: 7.85e-17 0.00e+00 5.00e-11 0.00e+00 0.00e+00 1.00e-10
This is Ipopt version 3.14.19, running with linear solver MUMPS 5.9.0.
Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 2
Number of nonzeros in Lagrangian Hessian.............: 3
Total number of variables............................: 2
variables with only lower bounds: 2
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 1
inequality constraints with only lower bounds: 1
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.9999940e-06 9.80e-01 1.33e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 5.3716721e-05 9.40e-01 1.18e+00 -1.0 3.97e-01 - 3.23e-02 5.03e-02h 1
2 2.8887051e-01 0.00e+00 4.87e+01 -1.0 4.95e-01 2.0 1.00e+00 1.00e+00f 1
3 2.8997185e-01 0.00e+00 1.00e-06 -1.0 1.33e-03 - 1.00e+00 1.00e+00f 1
4 2.8022658e-01 0.00e+00 3.97e-01 -2.5 1.19e-02 1.5 1.00e+00 1.00e+00f 1
5 2.5585775e-01 0.00e+00 3.45e-01 -2.5 3.10e-02 1.0 1.00e+00 1.00e+00f 1
6 2.5292585e-01 0.00e+00 1.16e-01 -2.5 1.03e-02 0.6 1.00e+00 3.73e-01f 2
7 2.5283150e-01 0.00e+00 1.54e-04 -2.5 1.25e-04 0.1 1.00e+00 1.00e+00f 1
8 2.5021392e-01 0.00e+00 1.14e-02 -3.8 3.48e-03 0.5 1.00e+00 1.00e+00f 1
9 2.5014901e-01 0.00e+00 9.49e-05 -3.8 8.65e-05 0.0 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 2.5000198e-01 0.00e+00 5.74e-04 -5.7 1.96e-04 0.5 1.00e+00 1.00e+00f 1
11 2.5000184e-01 0.00e+00 1.48e-06 -5.7 1.90e-07 0.9 1.00e+00 1.00e+00f 1
12 2.5000000e-01 0.00e+00 6.39e-06 -8.6 2.46e-06 0.4 1.00e+00 1.00e+00f 1
13 2.5000000e-01 0.00e+00 1.67e-10 -8.6 2.41e-11 0.8 1.00e+00 1.00e+00f 1
Number of Iterations....: 13
(scaled) (unscaled)
Objective...............: 2.4999999500590336e-01 2.4999999500590336e-01
Dual infeasibility......: 1.6744738928764491e-10 1.6744738928764491e-10
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 2.5059035596801710e-09 2.5059035596801710e-09
Overall NLP error.......: 2.5059035596801710e-09 2.5059035596801710e-09
Number of objective function evaluations = 15
Number of objective gradient evaluations = 14
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 15
Number of equality constraint Jacobian evaluations = 0
Number of inequality constraint Jacobian evaluations = 1
Number of Lagrangian Hessian evaluations = 13
Total seconds in IPOPT = 0.004
EXIT: Optimal Solution Found.As we can see below, the termination status is LOCALLY_SOLVED and not of OPTIMAL because Ipopt only guarantees local optimality.
solution_summary(model)solution_summary(; result = 1, verbose = false)
├ solver_name : Ipopt
├ Termination
│ ├ termination_status : LOCALLY_SOLVED
│ ├ result_count : 1
│ └ raw_status : Solve_Succeeded
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ └ objective_value : 2.50000e-01
└ Work counters
├ solve_time (sec) : 4.51207e-03
└ barrier_iterations : 13Indeed, the solution found is not globally optimal:
value(a), value(b)(0.49999999667061656, 0.4999999966705879)Note that the problem can be written equivalently as follows using registered functions. The difference is that the gradient and Hessian will be computed via the Symbolic Differentiation provided by MultivariatePolynomials instead of JuMP's Automatic Differentiation:
f(a, b) = p(x => a, y => b)
∇p = differentiate(p, [x, y])
function ∇f(g, a, b)
for i in eachindex(g)
g[i] = ∇p[i](x => a, y => b)
end
end
∇²p = differentiate(∇p, [x, y])
function ∇²f(H, a, b)
for j in axes(∇²p, 2)
for i in j:size(∇²p, 1)
H[i, j] = ∇²p[i, j](x => a, y => b)
end
end
end
using Ipopt
gmodel = Model(Ipopt.Optimizer)
@variable(gmodel, α >= 0)
@variable(gmodel, β >= 0)
@constraint(gmodel, α + β >= 1)
register(gmodel, :f, 2, f, ∇f, ∇²f)
@NLobjective(gmodel, Min, f(α, β))
optimize!(gmodel)This is Ipopt version 3.14.19, running with linear solver MUMPS 5.9.0.
Number of nonzeros in equality constraint Jacobian...: 0
Number of nonzeros in inequality constraint Jacobian.: 2
Number of nonzeros in Lagrangian Hessian.............: 3
Total number of variables............................: 2
variables with only lower bounds: 2
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 1
inequality constraints with only lower bounds: 1
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.9999940e-06 9.80e-01 1.33e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 5.3716721e-05 9.40e-01 1.18e+00 -1.0 3.97e-01 - 3.23e-02 5.03e-02h 1
2 2.8887051e-01 0.00e+00 4.87e+01 -1.0 4.95e-01 2.0 1.00e+00 1.00e+00f 1
3 2.8997185e-01 0.00e+00 1.00e-06 -1.0 1.33e-03 - 1.00e+00 1.00e+00f 1
4 2.8022658e-01 0.00e+00 3.97e-01 -2.5 1.19e-02 1.5 1.00e+00 1.00e+00f 1
5 2.5585775e-01 0.00e+00 3.45e-01 -2.5 3.10e-02 1.0 1.00e+00 1.00e+00f 1
6 2.5292585e-01 0.00e+00 1.16e-01 -2.5 1.03e-02 0.6 1.00e+00 3.73e-01f 2
7 2.5283150e-01 0.00e+00 1.54e-04 -2.5 1.25e-04 0.1 1.00e+00 1.00e+00f 1
8 2.5021392e-01 0.00e+00 1.14e-02 -3.8 3.48e-03 0.5 1.00e+00 1.00e+00f 1
9 2.5014901e-01 0.00e+00 9.49e-05 -3.8 8.65e-05 0.0 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 2.5000198e-01 0.00e+00 5.74e-04 -5.7 1.96e-04 0.5 1.00e+00 1.00e+00f 1
11 2.5000184e-01 0.00e+00 1.48e-06 -5.7 1.90e-07 0.9 1.00e+00 1.00e+00f 1
12 2.5000000e-01 0.00e+00 6.39e-06 -8.6 2.46e-06 0.4 1.00e+00 1.00e+00f 1
13 2.5000000e-01 0.00e+00 1.67e-10 -8.6 2.41e-11 0.8 1.00e+00 1.00e+00f 1
Number of Iterations....: 13
(scaled) (unscaled)
Objective...............: 2.4999999500590339e-01 2.4999999500590339e-01
Dual infeasibility......: 1.6744738928764491e-10 1.6744738928764491e-10
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 2.5059035596801710e-09 2.5059035596801710e-09
Overall NLP error.......: 2.5059035596801710e-09 2.5059035596801710e-09
Number of objective function evaluations = 15
Number of objective gradient evaluations = 14
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 15
Number of equality constraint Jacobian evaluations = 0
Number of inequality constraint Jacobian evaluations = 1
Number of Lagrangian Hessian evaluations = 13
Total seconds in IPOPT = 0.052
EXIT: Optimal Solution Found.Even if we have the algebraic expressions of gradient and Hessian, Ipopt is not using these symbolic expressions but only local information hence it can still only provide local guarantees:
solution_summary(gmodel)solution_summary(; result = 1, verbose = false)
├ solver_name : Ipopt
├ Termination
│ ├ termination_status : LOCALLY_SOLVED
│ ├ result_count : 1
│ └ raw_status : Solve_Succeeded
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ ├ objective_value : 2.50000e-01
│ └ dual_objective_value : 7.50000e-01
└ Work counters
├ solve_time (sec) : 5.23570e-02
└ barrier_iterations : 13and the same solution is found:
value(α), value(β)(0.49999999667061656, 0.4999999966705879)Sum-of-Squares approach
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 SCS
scs = SCS.Optimizer
import Dualization
dual_scs = Dualization.dual_optimizer(scs)#7 (generic function with 1 method)The Sum-of-Squares approach can be applied to the same model by simply changing its optimizer to SumOfSquares.Optimizer. This optimizer computes a lower bound using the Sum-of-Squares relaxation detailed in the section "How it works" below. The SDP solver solving this relaxation is given as argument to its constructor:
set_optimizer(model, () -> SumOfSquares.Optimizer(dual_scs))
optimize!(model)------------------------------------------------------------------
SCS v3.3.1 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 45, constraints m: 70
cones: z: primal zero / dual free vars: 31
s: psd vars: 39, ssize: 4
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: 5
compiled with openmp parallelization enabled
lin-sys: sparse-direct-amd-qdldl
nnz(A): 120, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 6.61e+00 1.01e+00 2.04e+01 -9.77e+00 1.00e-01 2.42e-04
75| 6.15e-05 2.12e-06 3.95e-06 -9.08e-06 1.00e-01 2.13e-03
------------------------------------------------------------------
status: solved
timings: total: 2.14e-03s = setup: 1.24e-04s + solve: 2.01e-03s
lin-sys: 9.77e-05s, cones: 1.73e-03s, accel: 6.53e-05s
------------------------------------------------------------------
objective = -0.000009
------------------------------------------------------------------The termination status is now OPTIMAL: the relaxation was solved to optimality so its objective value, queried with objective_bound, is a global lower bound of 0 for the polynomial problem. The result_count is however zero: no candidate solution could be recovered from the relaxation (we detail why at the end of the section "How it works"). Combining this lower bound with the solution found by Ipopt, we know at this point that the optimal value is in the interval $[0, 1/4]$.
solution_summary(model)solution_summary(; result = 1, verbose = false)
├ solver_name : SumOfSquares
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count : 0
│ ├ raw_status : solved
│ └ objective_bound : -7.10150e-06
├ Solution (result = 1)
│ ├ primal_status : NO_SOLUTION
│ └ dual_status : NO_SOLUTION
└ Work counters
└ solve_time (sec) : 1.67704e+00SAGE approach
The Sum-of-Squares certificate is not the only nonnegativity certificate that can be used in such relaxation. The SAGE certificate leads to a relative entropy program, solved with PolyJuMP.SAGE.Optimizer as follows. The candidate solution is recovered from the dual of the relaxation, which is not unique, so it depends on the solver used. We use the interior-point solver Clarabel that converges to the center of the optimal face which gives a reproducible candidate.
import Clarabel
set_optimizer(model, () -> PolyJuMP.SAGE.Optimizer(Clarabel.Optimizer))
optimize!(model)-------------------------------------------------------------
Clarabel.jl v0.11.1 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 1546
constraints = 1680
nnz(P) = 0
nnz(A) = 3957
cones (total) = 482
: Zero = 1, numel = 135
: Nonnegative = 1, numel = 105
: Exponential = 480, numel = (3,3,3,3,...,3)
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 0.0000e+00 -0.0000e+00 0.00e+00 9.53e-01 1.10e+00 1.00e+00 1.00e+00 ------
1 5.2933e-01 4.6570e-01 6.36e-02 4.72e-01 3.74e-01 3.09e-01 4.32e-01 6.43e-01
2 -5.7978e-01 -5.9803e-01 1.82e-02 1.55e-01 8.29e-02 4.72e-02 1.06e-01 9.08e-01
3 -4.8531e-01 -4.9374e-01 8.43e-03 6.07e-02 2.98e-02 9.54e-03 4.04e-02 8.29e-01
4 -1.8918e-01 -1.9163e-01 2.45e-03 1.86e-02 9.16e-03 2.57e-03 1.35e-02 7.02e-01
5 -1.3440e-01 -1.3619e-01 1.78e-03 1.33e-02 6.33e-03 1.75e-03 9.73e-03 4.06e-01
6 -3.6672e-02 -3.7176e-02 5.03e-04 3.88e-03 1.65e-03 5.13e-04 2.77e-03 7.92e-01
7 -3.6632e-02 -3.7116e-02 4.84e-04 3.65e-03 1.55e-03 4.73e-04 2.62e-03 1.33e-01
8 -1.4602e-02 -1.4802e-02 1.99e-04 1.41e-03 5.37e-04 1.67e-04 9.88e-04 7.92e-01
9 -5.8164e-03 -5.8958e-03 7.95e-05 5.53e-04 1.86e-04 6.43e-05 3.71e-04 7.92e-01
10 -1.9749e-03 -2.0000e-03 2.51e-05 1.72e-04 5.30e-05 1.97e-05 1.13e-04 7.92e-01
11 -1.2236e-03 -1.2381e-03 1.45e-05 9.81e-05 2.84e-05 1.09e-05 6.35e-05 6.34e-01
12 -3.0527e-04 -3.0886e-04 3.59e-06 2.52e-05 6.59e-06 2.93e-06 1.57e-05 7.92e-01
13 -1.3031e-04 -1.3168e-04 1.37e-06 9.53e-06 2.33e-06 1.10e-06 5.87e-06 7.92e-01
14 -3.6616e-05 -3.7028e-05 4.12e-07 3.02e-06 6.67e-07 3.71e-07 1.77e-06 7.92e-01
15 -3.0885e-05 -3.1202e-05 3.17e-07 2.26e-06 4.97e-07 2.71e-07 1.34e-06 4.06e-01
16 -9.6001e-06 -9.6958e-06 9.57e-08 7.04e-07 1.44e-07 8.71e-08 4.04e-07 7.92e-01
17 -3.2585e-06 -3.2874e-06 2.89e-08 2.16e-07 4.16e-08 2.72e-08 1.22e-07 7.92e-01
18 -1.2189e-06 -1.2296e-06 1.07e-08 8.40e-08 1.51e-08 1.10e-08 4.51e-08 7.92e-01
19 -3.9392e-07 -3.9716e-07 3.24e-09 2.63e-08 4.47e-09 3.57e-09 1.35e-08 7.92e-01
20 -3.2395e-07 -3.2645e-07 2.50e-09 2.02e-08 3.38e-09 2.74e-09 1.02e-08 4.06e-01
21 -1.1097e-07 -1.1170e-07 7.25e-10 6.46e-09 1.02e-09 9.50e-10 3.08e-09 7.92e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 23.8msThe SAGE relaxation also certifies (up to the tolerance of the solver) the lower bound 0. This time, a candidate solution is recovered from the dual of the relaxation. It is feasible but its objective value $27/32 \approx 0.84$ does not close the gap with the lower bound; the interval is still $[0, 1/4]$.
solution_summary(model)solution_summary(; result = 1, verbose = false)
├ solver_name : PolyJuMP.SAGE
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count : 1
│ ├ raw_status : SOLVED
│ └ objective_bound : 1.10970e-07
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : NO_SOLUTION
│ └ objective_value : 8.40551e-01
└ Work counters
└ solve_time (sec) : 1.57799e+00The recovered candidate is the following:
value(a), value(b)(0.7490523459480604, 0.7490496247425655)How it works
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 and finds zero.
model = SOSModel(dual_scs)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c3, p >= α, domain = S)
optimize!(model)------------------------------------------------------------------
SCS v3.3.1 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 10, constraints m: 25
cones: z: primal zero / dual free vars: 1
s: psd vars: 24, ssize: 4
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: 5
compiled with openmp parallelization enabled
lin-sys: sparse-direct-amd-qdldl
nnz(A): 37, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.12e+01 1.01e+00 3.16e+01 -1.54e+01 1.00e-01 1.96e-04
75| 1.26e-06 3.20e-07 1.87e-07 2.28e-07 1.00e-01 1.58e-03
------------------------------------------------------------------
status: solved
timings: total: 1.58e-03s = setup: 8.76e-05s + solve: 1.50e-03s
lin-sys: 3.85e-05s, cones: 1.26e-03s, accel: 1.01e-04s
------------------------------------------------------------------
objective = 0.000000
------------------------------------------------------------------This time, the termination status is OPTIMAL but this does not necessarily mean that we found the optimal solution to the polynomial optimization problem. This only means that SCS found an optimal solution to the Sum-of-Squares relaxation.
solution_summary(model)solution_summary(; result = 1, verbose = false)
├ solver_name : Dual model with SCS attached
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count : 1
│ └ raw_status : solved
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ ├ objective_value : 3.21788e-07
│ └ dual_objective_value : 1.34456e-07
└ Work counters
└ solve_time (sec) : 1.58375e-03The feasibility of the primal solution guarantees that the objective value 0 is a lower bound to the polynomial optimization problem. The optimality means that it's the best lower bound we can get (at this degree of the hierarcy). 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.
ν3 = moment_matrix(c3)
atomic_measure(ν3, 1e-3) # Returns nothing as the dual is not atomicFortunately, 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 a different effect in the choice of the maximum total degree of some $s_i$.
model = SOSModel(dual_scs)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c4, p >= α, domain = S, maxdegree = 4)
optimize!(model)------------------------------------------------------------------
SCS v3.3.1 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 10, constraints m: 25
cones: z: primal zero / dual free vars: 1
s: psd vars: 24, ssize: 4
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: 5
compiled with openmp parallelization enabled
lin-sys: sparse-direct-amd-qdldl
nnz(A): 37, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.12e+01 1.01e+00 3.16e+01 -1.54e+01 1.00e-01 3.23e-04
75| 1.26e-06 3.20e-07 1.87e-07 2.28e-07 1.00e-01 1.72e-03
------------------------------------------------------------------
status: solved
timings: total: 1.72e-03s = setup: 8.51e-05s + solve: 1.64e-03s
lin-sys: 3.91e-05s, cones: 1.26e-03s, accel: 1.04e-04s
------------------------------------------------------------------
objective = 0.000000
------------------------------------------------------------------We can see that the basis of the moment matrix didn't increase:
moment_matrix(c4)MomentMatrix with row/column basis:
SubBasis{Monomial}([1, y, x])
And entries in a 3×3 SymMatrix{Float64}:
1.0000002372798955 0.4999998349645504 0.49999983496455164
0.4999998349645504 0.4999997865612889 -2.3527197411372402e-7
0.49999983496455164 -2.3527197411372402e-7 0.4999997865612893This is because of the Newton polytope reduction that determined that gram matrix will be zero for these degrees so it reduced the problem back to the equivalent of maxdegree 3 Let's turn this off with newton_polytope = nothing
function sos(solver, deg)
model = SOSModel(solver)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c, p >= α, domain = S, maxdegree = deg, newton_polytope = nothing)
optimize!(model)
return model
end
dual_model4 = sos(dual_scs, 4)------------------------------------------------------------------
SCS v3.3.1 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 15, constraints m: 40
cones: z: primal zero / dual free vars: 1
s: psd vars: 39, ssize: 4
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: 5
compiled with openmp parallelization enabled
lin-sys: sparse-direct-amd-qdldl
nnz(A): 52, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 6.69e+00 4.94e-01 2.05e+01 -9.83e+00 1.00e-01 2.40e-04
75| 3.62e-05 1.05e-06 3.95e-06 3.77e-06 1.00e-01 2.05e-03
------------------------------------------------------------------
status: solved
timings: total: 2.05e-03s = setup: 1.15e-04s + solve: 1.93e-03s
lin-sys: 5.27e-05s, cones: 1.71e-03s, accel: 5.56e-05s
------------------------------------------------------------------
objective = 0.000004
------------------------------------------------------------------We see that the lower bound is still 0:
solution_summary(dual_model4)solution_summary(; result = 1, verbose = false)
├ solver_name : Dual model with SCS attached
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count : 1
│ └ raw_status : solved
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ ├ objective_value : 1.79206e-06
│ └ dual_objective_value : 5.73971e-06
└ Work counters
└ solve_time (sec) : 2.04907e-03Let's now look at which solution we can extract from the moment matrix:
dual_ν4 = moment_matrix(dual_model4[:c])MomentMatrix with row/column basis:
SubBasis{Monomial}([1, y, x, y^2, x*y, x^2])
And entries in a 6×6 SymMatrix{Float64}:
1.0000029341830468 0.49999783146691473 … 0.4999343565938622
0.49999783146691473 0.499935123900179 3.6886780619883884e-5
0.49999783146692595 5.907240172334005e-5 0.4998796305415934
0.4999343565938504 0.4998796305415815 0.8694642575517366
5.526500460034848e-5 3.9926933686792856e-5 -1.07208915800904
0.4999343565938622 3.6886780619883884e-5 … 5.313963225889048Looking at the singular values, 4 seems to be a reasonable rank:
using LinearAlgebra
svdvals(Matrix(dual_ν4.Q))6-element Vector{Float64}:
6.719973081607136
4.506850936016416
1.416446754861985
0.43752408381934577
0.41646903469355584
2.6543439393011307e-18The atoms we extract lie on the line x + y = 1 but they are not the global minimizers. Which atoms we get depends on the solver: with SCS v2.6, it is (0.5, 0.5), the solution found by Ipopt.
atomic_measure(dual_ν4, FixedRank(4))Atomic measure on the variables x, y with 2 atoms:
at [0.8306194493257181, 0.16938064701165073] with weight 1.0513391752139143
at [0.1693806470117425, 0.8306194493256541] with weight 1.0513391752137435This process is quite sensitive numerically so let's try to solve it without dualization as well:
model4 = sos(scs, 4)------------------------------------------------------------------
SCS v3.3.1 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem: variables n: 40, constraints m: 54
cones: z: primal zero / dual free vars: 15
s: psd vars: 39, ssize: 4
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: 5
compiled with openmp parallelization enabled
lin-sys: sparse-direct-amd-qdldl
nnz(A): 91, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res | gap | obj | scale | time (s)
------------------------------------------------------------------
0| 1.56e+01 1.35e+00 2.66e+01 -1.31e+01 1.00e-01 2.41e-04
150| 3.21e-05 2.00e-05 2.21e-05 2.04e-05 1.00e-01 4.15e-03
------------------------------------------------------------------
status: solved
timings: total: 4.16e-03s = setup: 1.26e-04s + solve: 4.03e-03s
lin-sys: 1.57e-04s, cones: 3.67e-03s, accel: 5.33e-05s
------------------------------------------------------------------
objective = 0.000020
------------------------------------------------------------------We see that the lower bound is again 0:
solution_summary(model4)solution_summary(; result = 1, verbose = false)
├ solver_name : SCS
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count : 1
│ └ raw_status : solved
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ ├ objective_value : -3.14947e-05
│ └ dual_objective_value : -9.35254e-06
└ Work counters
└ solve_time (sec) : 4.15621e-03The moment matrix is the following
ν4 = moment_matrix(model4[:c])MomentMatrix with row/column basis:
SubBasis{Monomial}([1, y, x, y^2, x*y, x^2])
And entries in a 6×6 SymMatrix{Float64}:
0.99999537176974 0.5000009724207274 … 0.5000139771043142
0.5000009724207274 0.5000114500907771 -8.924967759520866e-6
0.5000009724197981 -7.172829545643228e-6 0.5000096640434137
0.5000139771041751 0.5000096640448924 0.2494845238631467
-1.9315969249553834e-5 4.556688537533686e-6 0.23193932626130517
0.5000139771043142 -8.924967759520866e-6 … 0.7925650285131572Looking at the singular values, 3 seems to be a reasonable rank:
svdvals(Matrix(ν4.Q))6-element Vector{Float64}:
2.187886714968942
1.0220315022890947
0.5704369436083492
0.03321004928362795
0.02106762527528648
1.3284570130763252e-16This 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)$.
atomic_measure(ν4, FixedRank(3))Atomic measure on the variables x, y with 2 atoms:
at [0.9949406027939827, 0.005057726814789662] with weight 0.5513776437204988
at [-0.01616766382274894, 1.0161659934321372] with weight 0.5235644210787144A deeper look into atom extraction
The moment matrix is transformed 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 the lowest singular values. When this system of equation has an infinite number of solutions, atomic_measure 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 atomic_measure returned nothing.
ν3 = moment_matrix(c3)
SumOfSquares.MultivariateMoments.compute_support!(ν3, LeadingRelativeRankTol(1e-3))With maxdegree = 4, we obtain the system
\[\begin{aligned} x + y & = 1\\ y^2 & = y\\ xy & = 0\\ -y + y^2 - x*y & = 0 y^2 - 2y + 1 & = x^2 \end{aligned}\]
ν4 = moment_matrix(model4[:c])
SumOfSquares.MultivariateMoments.compute_support!(ν4, FixedRank(3))This system can be reduced to the equivalent system
\[\begin{aligned} x + y & = 1\\ y^2 & = y \end{aligned}\]
which has the solutions $(0, 1)$ and $(1, 0)$.
SemialgebraicSets.compute_gröbner_basis!(ideal(ν4.support))
ν4.supportAlgebraic Set defined by 2 equalities
-0.999998329608769 + 0.9999999999993909*y + x = 0
0.005139489993259282 - 1.0212237202469265*y + y^2 = 0
The solutions of this system then give the minimizers
collect(ν4.support)2-element Vector{Vector{Float64}}:
[0.9949406027939821, 0.005057726814790055]
[-0.016167663822748528, 1.0161659934321365]The function atomic_measure then reuses 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 how the function obtained the result $0.5 \delta(x-1, y) + 0.5 \delta(x, y-1)$ given in the previous section.
HomotopyContinuation
As discussed in the previous section, the atom extraction relies on the solution of a system of algebraic equations. The atomic_measure function takes an optional algebraic_solver argument that is used to solve this system of equation. If no solver is provided, the default solver of SemialgebraicSets.jl is used which currently computes the Gröbner basis, then the multiplication matrices and then the Schur decomposition of a random combination of these matrices. As the system of equations is obtained from a numerical solution and is represented using floating point coefficients, homotopy continuation is recommended as it is more numerically robust than Gröbner basis computation. The following uses homotopy continuation to solve the system of equations. Note that we fix the seed to be reproducible, it's not needed if reproducibility is not a concern.
using HomotopyContinuation
algebraic_solver = SemialgebraicSetsHCSolver(; excess_residual_tol = 1e-1, real_tol = 1e-1, compile = false, seed = UInt32(1))
atomic_measure(ν4, FixedRank(3), Echelon(), algebraic_solver)Atomic measure on the variables x, y with 2 atoms:
at [0.9895441237196457, 0.01826716876552563] with weight 0.5537973229430497
at [-0.000865044495158338, 1.0039503412191175] with weight 0.5342399591578915As the system has 3 equations for 2 variables and the coefficients of the equations are to be treated with tolerance since they originate from the solution of an SDP, we need to set excess_residual_tol and real_tol to a high tolerance otherwise, HomotopyContinuation would consider that there is no solution. Indeed, as the system is overdetermined (it has more equations than variables) HomotopyContinuation expects to have excess solution hence it filters out excess solution among the solution found. It determines which solution are in excess by comparing the infinity norm of the residuals of the equations at the solution with excess_residual_tol. It also filters out solution for which the absolute value of the imaginary part of one of the entry is larger than real_tol and strips out the imaginary part. We also fix the seed of the random start system used by HomotopyContinuation so that the result is reproducible. The raw solutions obtained by HomotopyContinuation can be obtained as follows:
F = HomotopyContinuation.System(ν4.support)
res = HomotopyContinuation.solve(F; algebraic_solver.options...)
path_results(res)4-element Vector{HomotopyContinuation.PathResult}:
PathResult:
• return_code → :excess_solution
• solution → ComplexF64[0.9895441237196457 - 0.02002483585466767im, 0.01826716876552563 + 0.031162379307783587im]
• accuracy → 0.065622
• residual → 0.03434
• condition_jacobian → 21.064
• steps → 31 / 0
• extended_precision → false
• path_number → 1
PathResult:
• return_code → :excess_solution
• solution → ComplexF64[0.5834989994939692 + 1.5957268401913185im, 1.246548093297424 - 1.146125200534452im]
• accuracy → 4.5533
• residual → 0.82456
• condition_jacobian → 2.8019
• steps → 37 / 0
• extended_precision → false
• path_number → 2
PathResult:
• return_code → :excess_solution
• solution → ComplexF64[-0.000865044495158338 - 0.04618109702403455im, 1.0039503412191175 + 0.030109508203704394im]
• accuracy → 0.078938
• residual → 0.047018
• condition_jacobian → 13.579
• steps → 22 / 0
• extended_precision → false
• path_number → 3
PathResult:
• return_code → :excess_solution
• solution → ComplexF64[134.33166499917485 - 8.700542331968991im, 135.0283275064215 - 11.315735713724424im]
• accuracy → 1298.0
• residual → 68.526
• condition_jacobian → 574.21
• steps → 41 / 0
• extended_precision → false
• path_number → 4
The printed residual above shows why 1e-1 allows to filter how the 2 actual solutions from the 2 excess solutions.
This page was generated using Literate.jl.