Pre-processing time and Mosek Solve inordinately long

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.

I’ll leave @blegat to answer this.

Do you have a reproducible example?

1 Like

The code finished running overnight. I am not sure if the Mosek sparsity error is being thrown by Mosek or by JuMP now, as the problem sent from Julia to Mosek completes in ~200 seconds, and the problem is much smaller than the one YALMIP sends to Moesk:

MOSEK warning 705 (MSK_RES_WRN_ZEROS_IN_SPARSE_ROW): 1 (nearly) zero elements are specified in sparse row ''(3234) of matrix 'A'.
MOSEK warning 705 (MSK_RES_WRN_ZEROS_IN_SPARSE_ROW): 3 (nearly) zero elements are specified in sparse row ''(3417) of matrix 'A'.
MOSEK warning 705 (MSK_RES_WRN_ZEROS_IN_SPARSE_ROW): 3 (nearly) zero elements are specified in sparse row ''(3483) of matrix 'A'.
MOSEK warning 705 (MSK_RES_WRN_ZEROS_IN_SPARSE_ROW): 1 (nearly) zero elements are specified in sparse row ''(3506) of matrix 'A'.
MOSEK warning 705 (MSK_RES_WRN_ZEROS_IN_SPARSE_ROW): 1 (nearly) zero elements are specified in sparse row ''(11043) of matrix 'A'.
MOSEK warning 705 (MSK_RES_WRN_ZEROS_IN_SPARSE_ROW): 3 (nearly) zero elements are specified in sparse row ''(11226) of matrix 'A'.
MOSEK warning 705 (MSK_RES_WRN_ZEROS_IN_SPARSE_ROW): 3 (nearly) zero elements are specified in sparse row ''(11292) of matrix 'A'.
MOSEK warning 705 (MSK_RES_WRN_ZEROS_IN_SPARSE_ROW): 1 (nearly) zero elements are specified in sparse row ''(11315) of matrix 'A'.
Problem
  Name                   :
  Objective sense        : minimize
  Type                   : CONIC (conic optimization problem)
  Constraints            : 62270
  Affine conic cons.     : 0
  Disjunctive cons.      : 0
  Cones                  : 0
  Scalar variables       : 1462
  Matrix variables       : 176 (scalarized: 213296)
  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.01
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.02
GP based matrix reordering started.
GP based matrix reordering terminated.
Optimizer  - threads                : 4
Optimizer  - solved problem         : the primal
Optimizer  - Constraints            : 62270
Optimizer  - Cones                  : 1
Optimizer  - Scalar variables       : 1463              conic                  : 1337
Optimizer  - Semi-definite variables: 176               scalarized             : 213296
Factor     - setup time             : 59.33
Factor     - dense det. time        : 23.47             GP order time          : 0.09
Factor     - nonzeros before factor : 1.97e+08          after factor           : 2.05e+08
Factor     - dense dim.             : 2                 flops                  : 9.97e+11
Factor     - GP saved nzs           : 5.71e+07          GP saved flops         : 5.17e+11
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME
0   4.3e+00  1.0e+00  1.0e+00  0.00e+00   0.000000000e+00   0.000000000e+00   1.0e+00  59.54
1   3.9e+00  9.0e-01  6.3e-01  1.12e+00   0.000000000e+00   -2.141610855e-01  9.0e-01  68.35
2   3.5e+00  8.1e-01  2.4e-01  4.67e+00   0.000000000e+00   -2.010974046e-01  8.1e-01  76.20
3   1.8e+00  4.2e-01  4.6e-02  2.58e+00   0.000000000e+00   -4.223159703e-02  4.2e-01  84.57
4   1.1e+00  2.5e-01  2.0e-02  2.30e+00   0.000000000e+00   -1.713102369e-02  2.5e-01  92.68
5   8.2e-01  1.9e-01  1.4e-02  1.17e+00   0.000000000e+00   -1.351085923e-02  1.9e-01  100.81
6   5.2e-01  1.2e-01  5.6e-03  1.73e+00   0.000000000e+00   -7.143249986e-03  1.2e-01  108.99
7   9.7e-02  2.3e-02  1.9e-04  1.63e+00   0.000000000e+00   -1.138514772e-03  2.3e-02  118.59
8   5.9e-02  1.4e-02  8.7e-05  1.48e+00   0.000000000e+00   -6.181437437e-04  1.4e-02  126.64
9   1.8e-02  4.2e-03  1.3e-05  1.30e+00   0.000000000e+00   -1.670972967e-04  4.2e-03  135.40
10  3.1e-03  7.1e-04  6.7e-07  1.19e+00   0.000000000e+00   -2.575307935e-05  7.1e-04  145.17
11  1.3e-03  3.0e-04  1.7e-07  1.26e+00   0.000000000e+00   -9.680049110e-06  3.0e-04  153.14
12  1.9e-04  4.5e-05  9.4e-09  1.13e+00   0.000000000e+00   -1.362385208e-06  4.5e-05  163.19
13  2.6e-05  6.1e-06  4.7e-10  1.06e+00   0.000000000e+00   -1.814125895e-07  6.1e-06  172.67
14  2.3e-05  5.3e-06  3.7e-10  1.01e+00   0.000000000e+00   -1.550888941e-07  5.3e-06  180.94
15  9.9e-06  2.3e-06  1.1e-10  1.01e+00   0.000000000e+00   -6.764456253e-08  2.3e-06  189.29
16  1.4e-06  3.2e-07  5.6e-12  1.00e+00   0.000000000e+00   -9.607204827e-09  3.2e-07  199.25
Optimizer terminated. Time: 199.44

47416.878854 seconds (266.72 G allocations: 26.844 TiB, 42.96% gc time, 0.01% compilation time: 100% of which was recompilation)

@odow The Julia code for a representative example of the non-expanded version is on a Github shared with @blegat. I can share it privately, and I can tweak it to give an example for the expanded version. It looks like the compute time is focused on how JuMP is sending the problem to Mosek, and that appears to take the hours of preparation time. So basically I am wondering if there is a way to reduce that preparation time by “helping” the code out a priori.

(As an aside, I would venture a guess that this does not have to do with symmetries, as the only symmetries I am imposing are sign symmetries, and the same method of sign symmetry identification is implemented in YALMIP.)

Can you push a small reproducible script to the private repository so that I can take a look ?

1 Like

I pushed a minimal example, let me know if it runs on your machine or if you run into any errors.

I ran your example, doing:

@time MOI.Utilities.attach_optimizer(model.moi_backend)

in order to isolate the time spent sending the problem to MosekTools. It took 13s, I guess you decreased the size of the problem when building the MWE. Among these 13 s, almost all the time is spent SumOfSquares.jl/src/Certificate/newton_polytope.jl at fd7336e413e58a7a5f6546c13ae281cf2b22bb05 · jump-dev/SumOfSquares.jl · GitHub
That would actually have been a good motivation for my ISMP talk in July ^^
Can you confirm that this remains the bottleneck for larger instances as well ?
You can see this for instance with @profview MOI.Utilities.attach_optimizer(model.moi_backend) in VS code

Apologies on the late reply, I was at a conference and did not check the email this account is attached to. I was a bit surprised when you said the example I gave you ran so quickly, so I checked the tail split subfolder and found some bugs which I have now fixed. My guess is you ran the command on the example in the main folder? All is now contained in the tail_split subfolder and should run cleanly. Regardless, now we can compare the two examples.

I am running your suggested @time function after setting up the problem now on what I have edited in the Github repo under the tail split subfolder. Thanks for making me aware of this particular timer, I did not know about it. It is now running for a ridiculous amount of time again, as it should, and I will edit this post with the results tomorrow morning, when it will finish running.
The point of me replying now is to let you know the Github has been corrected, and you can run the command on the subfolder for when you may want to poke around the compiled model.

Edit: Running @time on the attach_optimizer command only produced the following output

51562.822447 seconds (266.70 G allocations: 26.843 TiB, 41.04% gc time, 0.00% compilation time)

I downloaded VS code to see if @profview would run, but that is not a command recognized by Julia with the package Profile. Not sure if I should be doing something more specific here.

@profile is running for now in the VS code, will see what output that gives in another 13 hours

Profiling such a long example isn’t always helpful. Can’t you make a smaller problem that demonstrates the same problem?

Itd be most helpful if you could post a public reproducible example.