# SDP Sum of Squares optimization

I am trying to solve the following optimization problem:

``````using DynamicPolynomials
using SumOfSquares
using SDPAFamily

χ=1.0

@polyvar x1 x2 x3 x4

f = [2*x1*x4 - 2*x2*x3 - x1,
-2*x1*x3 - 2*x2*x4 - x2,
2*x1*x2,
x2*x2 - x1*x1]

model = SOSModel(() -> SDPAFamily.Optimizer{Float64}( # JuMP only supports Float64
variant = :sdpa_gmp, # use the arbitrary precision variant
params = SDPAFamily.STABLE_BUT_SLOW,
presolve=true,
#params = (  epsilonStar = 1e-10, # constraint tolerance
#            epsilonDash = 1e-10, # normalized duality gap tolerance
#            precision = 200) # arithmetric precision used in sdpa_gmp
))

X = monomials([x1, x2, x3, x4], 1:4);
@variable(model, V, Poly(X));

LHS = differentiate(V, x1)*f[1]+differentiate(V, x2)*f[2] + differentiate(V, x3)*f[3] +
+ differentiate(V, x4)*f[4] + x3^2 + x4^2 + differentiate(differentiate(V, x3),x3)/(4*χ) +
+ differentiate(differentiate(V, x4),x4)/(4*χ);

@variable(model, L);
@objective(model, Max, L);
@constraint(model, LHS >= L);
JuMP.optimize!(model)
``````

The problem is solved if `X` contains only low powers of monomials (from 1 to 5), but with an increase in this parameter, the solver ceases to produce even the solution that it found before. I tried a lot of different solvers, but without success.

Maybe someone can suggest what can be done to increase the number of monomials in `X`?

1 Like

what can be done to increase the number of monomials in `X` ?

What powers do you want to get up to?

I tried a lot of different solvers, but without success.

Which ones? Have you tried Mosek?

What powers do you want to get up to?

The bigger, the better. I would expect that as the power of `X` increases, the value of `L` will also increase. From an independent stochastic simulation, I know the “theoretical limit” of `L` and I want to understand how close I can get using Sum of Squares programming.

Which ones? Have you tried Mosek?

Yes, I tried Mosek, CDCS, Clarabel, COSMO, CSDP, Hypatia, ProxSDP, SCS, SeDuMi, SDPA-GMP. However, I always used the Float64 data type because JUMP does not support BigFloat. Can I solve this problem with julia but using the BigFloat type and how? Perhaps then I will be able to increase the degree of `X`.

@blegat might have some suggestions for numerical improvements. He is the SumOfSquares.jl expert.

What do you mean by “cease to produce any solution” ? Does it take too much time or does it finish but with an inconclusive status ?

It finishes with an inconclusive status. For example, the output for monomials of order 8 (Mosek solver):

``````Problem
Name                   :
Objective sense        : maximize
Type                   : CONIC (conic optimization problem)
Constraints            : 688
Affine conic cons.     : 0
Disjunctive cons.      : 0
Cones                  : 0
Scalar variables       : 495
Matrix variables       : 1
Integer variables      : 0

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 82
Eliminator terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 2                 time                   : 0.00
Lin. dep.  - tries                  : 1                 time                   : 0.00
Lin. dep.  - number                 : 17
Presolve terminated. Time: 0.00
GP based matrix reordering started.
GP based matrix reordering terminated.
Problem
Name                   :
Objective sense        : maximize
Type                   : CONIC (conic optimization problem)
Constraints            : 688
Affine conic cons.     : 0
Disjunctive cons.      : 0
Cones                  : 0
Scalar variables       : 495
Matrix variables       : 1
Integer variables      : 0

Optimizer  - solved problem         : the primal
Optimizer  - Constraints            : 462
Optimizer  - Cones                  : 1
Optimizer  - Scalar variables       : 319               conic                  : 319
Optimizer  - Semi-definite variables: 1                 scalarized             : 2145
Factor     - setup time             : 0.02              dense det. time        : 0.00
Factor     - ML order time          : 0.00              GP order time          : 0.00
Factor     - nonzeros before factor : 1.07e+05          after factor           : 1.07e+05
Factor     - dense dim.             : 0                 flops                  : 4.35e+07
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME
0   1.0e+00  1.0e+00  1.0e+00  0.00e+00   0.000000000e+00   0.000000000e+00   1.0e+00  0.02
1   2.8e-01  2.8e-01  1.3e-01  8.41e-01   5.658786273e-02   1.292251824e-01   2.8e-01  0.02
2   1.2e-01  1.2e-01  3.1e-02  8.56e-01   2.971659068e-01   3.385846790e-01   1.2e-01  0.03
3   2.2e-02  2.2e-02  1.8e-03  1.24e+00   4.004710586e-01   4.118891854e-01   2.2e-02  0.03
4   3.7e-03  3.7e-03  1.6e-04  1.13e+00   3.933816935e-01   3.947614139e-01   3.7e-03  0.03
5   1.0e-03  1.0e-03  3.4e-05  6.84e-01   3.845216114e-01   3.846743932e-01   1.0e-03  0.03
6   2.8e-04  2.8e-04  8.1e-06  2.92e-01   3.631527616e-01   3.628386926e-01   2.8e-04  0.05
7   9.4e-05  9.4e-05  2.5e-06  1.94e-01   3.486957620e-01   3.482766828e-01   9.4e-05  0.05
8   2.6e-05  2.6e-05  6.6e-07  7.92e-02   3.312750693e-01   3.307510053e-01   2.6e-05  0.05
9   7.5e-06  7.5e-06  1.7e-07  2.59e-01   3.203259811e-01   3.198739211e-01   7.5e-06  0.06
10  2.3e-06  2.3e-06  5.1e-08  3.58e-02   3.085929591e-01   3.081301685e-01   2.3e-06  0.06
11  7.2e-07  7.2e-07  1.4e-08  1.88e-01   3.004379714e-01   3.000531242e-01   7.2e-07  0.06
12  2.0e-07  2.0e-07  3.9e-09  7.59e-02   2.916395607e-01   2.912761518e-01   2.0e-07  0.08
13  5.8e-08  5.8e-08  1.0e-09  2.22e-01   2.855252086e-01   2.852327366e-01   5.8e-08  0.08
14  2.5e-08  1.6e-08  2.6e-10  8.99e-02   2.793028964e-01   2.790338667e-01   1.6e-08  0.08
15  8.6e-09  4.2e-09  5.8e-11  3.29e-01   2.749540167e-01   2.747579127e-01   4.2e-09  0.09
16  3.4e-09  1.6e-09  2.5e-11  -1.04e-02  2.703414885e-01   2.700904299e-01   1.6e-09  0.11
17  1.3e-09  5.8e-10  6.6e-12  5.23e-01   2.689030845e-01   2.687708780e-01   5.8e-10  0.11
18  1.2e-09  4.9e-10  5.5e-12  3.23e-01   2.684525884e-01   2.683226175e-01   4.8e-10  0.11
19  1.1e-09  3.3e-10  3.7e-12  2.58e-01   2.673859138e-01   2.672588002e-01   3.2e-10  0.12
20  1.1e-09  3.2e-10  3.6e-12  4.80e-01   2.673412672e-01   2.672160643e-01   3.1e-10  0.12
21  1.1e-09  3.2e-10  3.6e-12  4.80e-01   2.673412672e-01   2.672160643e-01   3.1e-10  0.14
22  1.1e-09  3.2e-10  3.6e-12  4.80e-01   2.673412672e-01   2.672160643e-01   3.1e-10  0.16
23  1.1e-09  3.2e-10  3.6e-12  1.00e+00   2.673411582e-01   2.672159702e-01   3.1e-10  0.17
24  1.1e-09  3.2e-10  3.6e-12  1.00e+00   2.673411037e-01   2.672159232e-01   3.1e-10  0.17
Optimizer terminated. Time: 0.20
``````

JuMP.primal_status(model)

``````UNKNOWN_RESULT_STATUS::ResultStatusCode = 8
``````

It looks quite close to feasible, maybe you can relax the feasibility tolerance : `MSK_DPAR_INTPNT_CO_TOL_PFEAS`
You could try to make the problem better conditioned by using other polynomial bases, see Constraints · SumOfSquares.

Thank you @blegat, I will try.

Hi, @blegat. Actually, your suggestions helped me a lot. I rescaled the initial system of equations and reduced the basis in the space of which the optimization takes place (based on some physical reasoning). Now the problem becomes:

``````using DynamicPolynomials
using SumOfSquares
using MosekTools

χ=1.0

@polyvar x1 x2 x3 x4

f = [2*x1*x4 - 2*x2*x3 - x1*χ^(1/3),
-2*x1*x3 - 2*x2*x4 - x2*χ^(1/3),
2*x1*x2,
x2*x2 - x1*x1]

model = SOSModel(Mosek.Optimizer)
set_optimizer_attribute(model, "MSK_DPAR_INTPNT_CO_TOL_PFEAS", 1e-3)

deg = 20

@polyvar n1 n2 J
arr = Polynomial{true, Float64}[]

for i in 0:ceil(deg/2)
for j in 0:ceil(deg/2)
for k in 0:ceil(deg/3)
if 1<=2*i+2*j+3*k<=deg
p = (n1^Int(i))*(n2^Int(j))*(J^Int(k))
push!(arr, p)
end
end
end
end

arr2 = subs(arr, n1=>x1*x1+x2*x2, n2=>x3*x3+x4*x4, J=>-4*x1*x2*x3+2*x1*x1*x4-2*x2*x2*x4)

basis = FixedPolynomialBasis(arr2)
@variable(model, V, Poly(basis))

LHS = differentiate(V, x1)*f[1]+differentiate(V, x2)*f[2] + differentiate(V, x3)*f[3] +
+ differentiate(V, x4)*f[4] + x3^2 + x4^2 + differentiate(differentiate(V, x3),x3)/(4) +
+ differentiate(differentiate(V, x4),x4)/(4);

@variable(model, L);
@objective(model, Max, L);
@constraint(model, LHS >= L);
JuMP.optimize!(model);

println(JuMP.primal_status(model))
println(objective_value(model)/χ^(2/3))
``````

The output for `deg=20` is following:

``````Problem
Name                   :
Objective sense        : maximize
Type                   : CONIC (conic optimization problem)
Constraints            : 11015
Affine conic cons.     : 0
Disjunctive cons.      : 0
Cones                  : 0
Scalar variables       : 192
Matrix variables       : 1
Integer variables      : 0

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 2
Eliminator terminated.
Eliminator - tries                  : 1                 time                   : 0.00
Lin. dep.  - tries                  : 1                 time                   : 0.00
Lin. dep.  - number                 : 152
Presolve terminated. Time: 0.00
GP based matrix reordering started.
GP based matrix reordering terminated.
Problem
Name                   :
Objective sense        : maximize
Type                   : CONIC (conic optimization problem)
Constraints            : 11015
Affine conic cons.     : 0
Disjunctive cons.      : 0
Cones                  : 0
Scalar variables       : 192
Matrix variables       : 1
Integer variables      : 0

Optimizer  - solved problem         : the primal
Optimizer  - Constraints            : 10565
Optimizer  - Cones                  : 1
Optimizer  - Scalar variables       : 191               conic                  : 191
Optimizer  - Semi-definite variables: 1                 scalarized             : 490545
Factor     - setup time             : 12.44             dense det. time        : 4.50
Factor     - ML order time          : 0.00              GP order time          : 0.00
Factor     - nonzeros before factor : 5.57e+07          after factor           : 5.57e+07
Factor     - dense dim.             : 0                 flops                  : 8.77e+11
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME
0   1.0e+00  1.0e+00  1.0e+00  0.00e+00   0.000000000e+00   0.000000000e+00   1.0e+00  13.11
1   4.8e-01  4.8e-01  3.3e-01  9.44e-01   4.550285069e-02   6.429942074e-02   4.8e-01  30.56
2   1.3e-01  1.3e-01  9.7e-02  5.14e-01   4.848911013e-01   1.575909731e-01   1.3e-01  49.72
3   4.0e-02  4.0e-02  5.1e-02  -6.23e-01  1.500290024e+00   1.720124574e-01   4.0e-02  70.61
4   1.4e-02  1.4e-02  2.6e-02  -6.48e-01  3.393740475e+00   2.310617171e-01   1.4e-02  90.11
5   4.7e-03  4.7e-03  1.2e-02  -6.54e-01  6.762067251e+00   2.970061982e-01   4.7e-03  108.58
6   1.4e-03  1.4e-03  5.2e-03  -6.23e-01  1.345494511e+01   3.677029505e-01   1.4e-03  128.56
7   5.3e-04  5.3e-04  2.4e-03  -4.72e-01  2.170525715e+01   4.316049177e-01   5.3e-04  147.66
8   3.2e-04  3.2e-04  1.3e-03  -1.14e-01  1.631440577e+01   4.815503283e-01   3.0e-04  166.56
9   5.2e-05  5.2e-05  2.4e-04  -1.70e-01  2.130817677e+01   5.856556728e-01   5.4e-05  187.61
10  2.1e-05  2.1e-05  8.5e-05  1.96e-01   1.664418491e+01   6.194935873e-01   2.1e-05  206.39
11  3.5e-06  3.5e-06  1.0e-05  3.19e-01   9.105502953e+00   6.693353513e-01   3.8e-06  261.75
12  7.9e-07  7.9e-07  1.3e-06  6.24e-01   3.315519459e+00   6.989735069e-01   8.5e-07  300.72
13  3.2e-07  1.4e-07  9.2e-08  8.62e-01   1.136843993e+00   7.058032411e-01   1.4e-07  322.11
14  7.1e-07  1.1e-07  6.1e-08  1.01e+00   1.036732868e+00   7.054957416e-01   1.1e-07  357.44
15  9.9e-07  9.6e-08  5.3e-08  1.00e+00   1.005906292e+00   7.055728351e-01   1.0e-07  377.91
16  9.8e-07  4.3e-08  1.6e-08  9.76e-01   8.441681925e-01   7.062274679e-01   4.5e-08  397.52
17  3.1e-07  1.2e-08  2.6e-09  9.83e-01   7.497089483e-01   7.073477438e-01   1.3e-08  417.08
18  5.7e-07  5.5e-09  7.7e-10  9.98e-01   7.267646158e-01   7.074757932e-01   6.0e-09  434.61
19  4.8e-07  1.4e-09  9.8e-11  9.97e-01   7.125862404e-01   7.076853996e-01   1.5e-09  452.86
20  1.1e-06  1.3e-09  8.4e-11  1.00e+00   7.121236355e-01   7.076887856e-01   1.4e-09  489.89
21  9.3e-07  6.5e-10  1.8e-11  1.00e+00   7.092900816e-01   7.076978265e-01   5.0e-10  509.56
22  4.5e-07  2.9e-10  2.6e-12  1.00e+00   7.081567134e-01   7.077133249e-01   1.4e-10  529.94
23  7.8e-07  1.9e-09  3.2e-13  1.00e+00   7.078313046e-01   7.077185736e-01   3.4e-11  548.22
24  7.8e-07  1.9e-09  3.2e-13  1.00e+00   7.078313046e-01   7.077185736e-01   3.4e-11  588.75
Optimizer terminated. Time: 612.20

UNKNOWN_RESULT_STATUS
0.7078313045648063
``````

and as I understand the result is again quite close to the feseability. Maybe you can tell me what other parameters of the solver should be paid attention to in order to try to improve its performance? I want to play around with the parameters a bit, but there are so many of them in the documentation that I can’t pick out the key ones.

You can see the parameters grouped by topic here:
https://docs.mosek.com/10.0/capi/param-groups.html#doc-param-groups

The likely ones relevant to you all start with `MSK_DPAR_INTPNT_CO_TOL_` under ‘Conic interior-point method’.

2 Likes

I just added an example on how to change the basis used for the Gram matrix in Mminimization of a univariate polynomial · SumOfSquares

1 Like