I have converted a YALMIP sum-of-squares optimization code to a Julia sum-of-squares optimization code. The base level of the code, from which I want to expand, seems to produce similar-sized problems in Julia and YALMIP. I print the Julia output here. (YALMIP is giving me some version issues, so I will have to try to add that later when I can restart Matlab.)
Julia
MOSEK warning 705 (MSK_RES_WRN_ZEROS_IN_SPARSE_ROW): 1 (nearly) zero elements are specified in sparse row ''(568) of matrix 'A'.
MOSEK warning 705 (MSK_RES_WRN_ZEROS_IN_SPARSE_ROW): 3 (nearly) zero elements are specified in sparse row ''(625) of matrix 'A'.
MOSEK warning 705 (MSK_RES_WRN_ZEROS_IN_SPARSE_ROW): 3 (nearly) zero elements are specified in sparse row ''(658) of matrix 'A'.
MOSEK warning 705 (MSK_RES_WRN_ZEROS_IN_SPARSE_ROW): 1 (nearly) zero elements are specified in sparse row ''(675) of matrix 'A'.
MOSEK warning 705 (MSK_RES_WRN_ZEROS_IN_SPARSE_ROW): 1 (nearly) zero elements are specified in sparse row ''(1316) of matrix 'A'.
MOSEK warning 705 (MSK_RES_WRN_ZEROS_IN_SPARSE_ROW): 3 (nearly) zero elements are specified in sparse row ''(1373) of matrix 'A'.
MOSEK warning 705 (MSK_RES_WRN_ZEROS_IN_SPARSE_ROW): 3 (nearly) zero elements are specified in sparse row ''(1406) of matrix 'A'.
MOSEK warning 705 (MSK_RES_WRN_ZEROS_IN_SPARSE_ROW): 1 (nearly) zero elements are specified in sparse row ''(1423) of matrix 'A'.
Problem
Name :
Objective sense : minimize
Type : CONIC (conic optimization problem)
Constraints : 3601
Affine conic cons. : 0
Disjunctive cons. : 0
Cones : 0
Scalar variables : 471
Matrix variables : 37 (scalarized: 11856)
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. - primal attempts : 1 successes : 1
Lin. dep. - dual attempts : 0 successes : 0
Lin. dep. - primal deps. : 0 dual deps. : 0
Presolve terminated. Time: 0.01
GP based matrix reordering started.
GP based matrix reordering terminated.
Optimizer - threads : 4
Optimizer - solved problem : the primal
Optimizer - Constraints : 3590
Optimizer - Cones : 1
Optimizer - Scalar variables : 460 conic : 455
Optimizer - Semi-definite variables: 37 scalarized : 11856
Factor - setup time : 0.31
Factor - dense det. time : 0.05 GP order time : 0.14
Factor - nonzeros before factor : 1.03e+06 after factor : 1.46e+06
Factor - dense dim. : 2 flops : 7.57e+08
Factor - GP saved nzs : 2.28e+05 GP saved flops : 2.26e+08
The size of the Julia model before calling optimize!(model) is given here
A JuMP Model
Feasibility problem with:
Variables: 458
`Vector{AffExpr}`-in-`SumOfSquares.SOSPolynomialSet{FullSpace, Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}, MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}, SumOfSquares.Certificate.Sparsity.Ideal{SignSymmetry, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}}}`: 15 constraints
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Mosek
The main difference is not the size of the problem being passed, but rather these strange sparse matrix warnings Mosek is giving us. Note that all of the constraints I impose are sum-of-squares constraints, there are no other types of constraints being considered. It did not produce different results compared to YALMIP for a large sweep parameter search, so I did not worry about it.
Until, I made a small expansion to the code.
The way I build off of the code is I (approximately) triple 4 of the constraints by splitting one of the simple variables into three separate positive variables. The new constraints are created using a minor modification of the original constraints.
If I run the expanded code with YALMIP, it pre-processes quite quickly (a few minutes), and Mosek takes a while to run but it is still solvable in about 35-40 minutes.
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 154407
Cones : 0
Scalar variables : 218
Matrix variables : 281
Integer variables : 0
Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries : 1 time : 0.00
Lin. dep. - tries : 1 time : 0.02
Lin. dep. - number : 0
Presolve terminated. Time: 0.04
GP based matrix reordering started.
GP based matrix reordering terminated.
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 154407
Cones : 0
Scalar variables : 218
Matrix variables : 281
Integer variables : 0
Optimizer - threads : 16
Optimizer - solved problem : the primal
Optimizer - Constraints : 154407
Optimizer - Cones : 0
Optimizer - Scalar variables : 218 conic : 0
Optimizer - Semi-definite variables: 281 scalarized : 198203
Factor - setup time : 370.44 dense det. time : 0.00
Factor - ML order time : 28.86 GP order time : 281.91
Factor - nonzeros before factor : 7.40e+08 after factor : 7.54e+08
Factor - dense dim. : 0 flops : 5.04e+12
Factor - GP saved nzs : 3.96e+09 GP saved flops : 2.64e+14
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 5.3e+00 1.0e+00 2.7e+01 0.00e+00 2.614865243e+01 0.000000000e+00 1.0e+00 370.85
1 4.8e+00 9.4e-01 1.7e+01 2.14e+00 1.771499467e+01 0.000000000e+00 9.2e-01 438.12
2 3.4e+00 6.6e-01 2.8e+00 3.24e+00 4.116435798e+00 0.000000000e+00 6.4e-01 508.09
3 2.7e+00 5.3e-01 2.0e+00 8.15e+00 1.016905384e+00 0.000000000e+00 5.1e-01 582.14
4 1.8e+00 3.6e-01 8.8e-01 2.92e+00 4.125881262e-01 0.000000000e+00 3.5e-01 657.40
5 8.3e-01 1.6e-01 2.1e-01 2.28e+00 1.209489166e-01 0.000000000e+00 1.6e-01 737.18
6 7.4e-01 1.4e-01 1.8e-01 1.11e+00 1.105135213e-01 0.000000000e+00 1.4e-01 811.60
7 5.7e-01 1.1e-01 1.0e-01 1.53e+00 7.389093983e-02 0.000000000e+00 1.1e-01 886.61
8 3.1e-01 6.1e-02 3.5e-02 1.59e+00 3.337636852e-02 0.000000000e+00 5.9e-02 970.10
9 4.7e-02 9.1e-03 7.0e-04 1.54e+00 3.772519840e-03 0.000000000e+00 8.9e-03 1061.62
10 2.8e-02 5.4e-03 2.9e-04 1.37e+00 2.053881229e-03 0.000000000e+00 5.3e-03 1140.67
11 2.0e-02 3.9e-03 1.7e-04 1.24e+00 1.417854136e-03 0.000000000e+00 3.8e-03 1215.24
12 9.7e-03 1.9e-03 5.0e-05 1.19e+00 6.417097031e-04 0.000000000e+00 1.8e-03 1292.67
13 3.9e-03 7.5e-04 1.1e-05 1.13e+00 2.463192535e-04 0.000000000e+00 7.4e-04 1369.12
14 9.7e-04 1.9e-04 1.1e-06 1.07e+00 6.051571352e-05 0.000000000e+00 1.9e-04 1448.05
15 3.1e-04 6.1e-05 1.8e-07 1.04e+00 1.886065009e-05 0.000000000e+00 5.9e-05 1525.38
16 1.3e-04 2.5e-05 4.1e-08 1.13e+00 7.240006424e-06 0.000000000e+00 2.4e-05 1602.11
17 1.7e-05 3.3e-06 1.3e-09 1.17e+00 8.431872660e-07 0.000000000e+00 3.2e-06 1685.79
18 1.1e-05 2.2e-06 7.2e-10 1.08e+00 5.471010230e-07 0.000000000e+00 2.1e-06 1775.25
19 1.9e-06 3.7e-07 5.5e-11 1.05e+00 9.248212079e-08 0.000000000e+00 3.6e-07 1872.91
20 1.4e-07 2.7e-08 1.1e-12 1.02e+00 6.554579871e-09 0.000000000e+00 2.6e-08 1951.36
21 1.2e-12 1.0e-12 2.2e-20 1.00e+00 5.811184064e-14 0.000000000e+00 2.3e-13 2038.72
Optimizer terminated. Time: 2039.75
Interior-point solution summary
Problem status : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal. obj: 5.8111840641e-14 nrm: 5e-08 Viol. con: 4e-15 var: 0e+00 barvar: 0e+00
Dual. obj: 0.0000000000e+00 nrm: 2e+00 Viol. con: 0e+00 var: 7e-18 barvar: 9e-15
Optimizer summary
Optimizer - time: 2039.75
Interior-point - iterations : 21 time: 2039.71
Basis identification - time: 0.00
Primal - iterations : 0 time: 0.00
Dual - iterations : 0 time: 0.00
Clean primal - iterations : 0 time: 0.00
Clean dual - iterations : 0 time: 0.00
Simplex - time: 0.00
Primal simplex - iterations : 0 time: 0.00
Dual simplex - iterations : 0 time: 0.00
Mixed integer - relaxations: 0 time: 0.00
Elapsed time is 2043.013840 seconds.
>>
The Julia code takes an entire day to run (when I am not getting weird licensing issues). It takes a little over an hour to set up the problem for Mosek to solve. Here is the setup for this Julia model.
A JuMP Model
Feasibility problem with:
Variables: 1336
`Vector{AffExpr}`-in-`SumOfSquares.SOSPolynomialSet{FullSpace, Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}, MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}, SumOfSquares.Certificate.Sparsity.Ideal{SignSymmetry, SumOfSquares.Certificate.Newton{SOSCone, MonomialBasis, SumOfSquares.Certificate.NewtonFilter{SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}}}}}}`: 30 constraints
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Mosek
I want to post this before it finishes running (it has been running for hours now), but I will update this with Julia’s solved output once it is done. Mosek also gives sparse matrix errors here.
I verified (multiple times) that the constraints I have in both systems are the same. All constraints are of the form
model = SOSModel(Mosek.Optimizer)
@polyvar(a[1:m]
V_coeffs = @variable(model, 1:length([1:length(basis2_90V) + length(basis3_90V)])
poly = sum(V_coeffs.*vcat(basis2_90V, basis3_90V));
@constraint(model, poly >= 0, sparsity=SignSymmetry())
Is there a better way to formulate these constraints to send the problem to Mosek in better shape?
@blegat This is a variation of the same code you see in the Github I sent you a while back. We are expanding the q variable in that code into three variables, specifically, q = q_1 + q_2 + q_3, and we have to enforce that each of these are positive so we define them to be q_i = (q_{sqrt,i})^2 (we do NOT enforce the positivity of q_i as a constraint). This changes the maximum degree of the polynomials we construct, but in a sense only superficially; it should not change the way the Gram matrices are formed under the hood.