# SDP Sum of Squares optimization

**URL:** <https://discourse.julialang.org/t/sdp-sum-of-squares-optimization/88487>\
**Category:** Optimization (Mathematical)\
**Tags:** question, jump\
**Created:** [October 9, 2022, 4:14pm UTC](https://discourse.julialang.org/t/sdp-sum-of-squares-optimization/88487 "2022-10-09T16:14:12Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![parf](https://avatars.discourse-cdn.com/v4/letter/p/779978/32.png) [@parf](https://discourse.julialang.org/u/parf)\
**Post date:** [October 9, 2022, 4:14pm UTC](https://discourse.julialang.org/t/sdp-sum-of-squares-optimization/88487/1 "2022-10-09T16:14:12Z")

</div>

I am trying to solve the following optimization problem:

```julia
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`?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [October 9, 2022, 8:35pm UTC](https://discourse.julialang.org/t/sdp-sum-of-squares-optimization/88487/2 "2022-10-09T20:35:19Z")

</div>

> 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?

---

<div class="post-metadata">

**Author:** ![parf](https://avatars.discourse-cdn.com/v4/letter/p/779978/32.png) [@parf](https://discourse.julialang.org/u/parf)\
**Post date:** [October 9, 2022, 9:19pm UTC](https://discourse.julialang.org/t/sdp-sum-of-squares-optimization/88487/3 "2022-10-09T21:19:02Z")

</div>

> 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`.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [October 10, 2022, 11:08pm UTC](https://discourse.julialang.org/t/sdp-sum-of-squares-optimization/88487/4 "2022-10-10T23:08:47Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![blegat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/blegat/32/217090_2.png) [@blegat](https://discourse.julialang.org/u/blegat)\
**Post date:** [October 11, 2022, 1:38pm UTC](https://discourse.julialang.org/t/sdp-sum-of-squares-optimization/88487/5 "2022-10-11T13:38:29Z")

</div>

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 ?

---

<div class="post-metadata">

**Author:** ![parf](https://avatars.discourse-cdn.com/v4/letter/p/779978/32.png) [@parf](https://discourse.julialang.org/u/parf)\
**Post date:** [October 11, 2022, 4:43pm UTC](https://discourse.julialang.org/t/sdp-sum-of-squares-optimization/88487/6 "2022-10-11T16:43:40Z")

</div>

> [@blegat](#):
>
> 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):

```julia
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 - threads : 4               
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)

```julia
UNKNOWN_RESULT_STATUS::ResultStatusCode = 8

```

---

<div class="post-metadata">

**Author:** ![blegat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/blegat/32/217090_2.png) [@blegat](https://discourse.julialang.org/u/blegat)\
**Post date:** [October 12, 2022, 11:28am UTC](https://discourse.julialang.org/t/sdp-sum-of-squares-optimization/88487/7 "2022-10-12T11:28:22Z")

</div>

It looks quite close to feasible, maybe you can relax the feasibility tolerance : [`MSK_DPAR_INTPNT_CO_TOL_PFEAS`](https://docs.mosek.com/10.0/capi/parameters.html#mosek.dparam.intpnt_co_tol_pfeas)  
You could try to make the problem better conditioned by using other polynomial bases, see [Constraints · SumOfSquares](https://jump.dev/SumOfSquares.jl/stable/constraints/#Changing-the-polynomial-basis).

---

<div class="post-metadata">

**Author:** ![parf](https://avatars.discourse-cdn.com/v4/letter/p/779978/32.png) [@parf](https://discourse.julialang.org/u/parf)\
**Post date:** [October 13, 2022, 7:31am UTC](https://discourse.julialang.org/t/sdp-sum-of-squares-optimization/88487/8 "2022-10-13T07:31:22Z")

</div>

Thank you @blegat, I will try.

---

<div class="post-metadata">

**Author:** ![parf](https://avatars.discourse-cdn.com/v4/letter/p/779978/32.png) [@parf](https://discourse.julialang.org/u/parf)\
**Post date:** [October 14, 2022, 11:28am UTC](https://discourse.julialang.org/t/sdp-sum-of-squares-optimization/88487/9 "2022-10-14T11:28:47Z")

</div>

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:

```julia
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:

```julia
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 - threads : 4               
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.

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [October 16, 2022, 11:45pm UTC](https://discourse.julialang.org/t/sdp-sum-of-squares-optimization/88487/10 "2022-10-16T23:45:56Z")

</div>

> [@parf](#):
>
> 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](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’.

---

<div class="post-metadata">

**Author:** ![blegat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/blegat/32/217090_2.png) [@blegat](https://discourse.julialang.org/u/blegat)\
**Post date:** [October 27, 2022, 1:45pm UTC](https://discourse.julialang.org/t/sdp-sum-of-squares-optimization/88487/11 "2022-10-27T13:45:08Z")

</div>

I just added an example on how to change the basis used for the Gram matrix in [Mminimization of a univariate polynomial · SumOfSquares](https://jump.dev/SumOfSquares.jl/latest/generated/Getting%20started/univariate/)