YALMIP vs JuMP

raktim · November 6, 2019, 4:05am

I compared JuMP and YALMIP’s performance in parsing convex problem. I noticed a difference in the problem size of the final optimization problem that is passed to Mosek. The problem from JuMP was quite large, than YALMIP.

YALMIP generated problem has 5 constraints, 0 scalar variables, 2 matrix variables
JuMP generated problem has 18 constraints, 5 scalar variables, 2 matrix variables

Both give the exact same answer. Any comments on why this is the case? Am I formulating the problem incorrectly in JuMP? Any help in this regard will be greatly appreciated. Details of the test problem are below.

Thanks.

Problem Statement

Design of LQR for a system \dot{x} = Ax + Bu, with u=Kx. This is solved using a SDP problem. I am not getting into the details of LQR theory here. The codes for YALMIP and JuMP are as following.

YALMIP (MATLAB)

clc; clear;

A = [0 1; -1 -1];
B = [0;1];
Q = [1 0; 0 0.1];
R = 1;

Y = sdpvar(2,2);
W = sdpvar(1,2)
Z = zeros(2,1);

M11 = (AY + BW) + (AY + BW)‘;
M = [M11 Y W’;
Y -inv(Q) Z;
W Z’ -inv( R );]
constr = [M <= 0, Y>=0];
optimize(constr,-trace(Y));
disp(trace(value(Y)));

JuMP Code

using JuMP, MosekTools, LinearAlgebra
A = [0. 1; -1 -1];
B = [0.;1];

model = Model(with_optimizer(Mosek.Optimizer));
@variable(model, Y[1:2, 1:2], Symmetric);
@variable(model, W[1:2]);

Q = [1 0; 0 0.1];
R = 1;
Z = zeros(2,1);

M11 = (AY + BW’) + (AY + BW’)‘;
M = [M11 Y W;
Y -inv(Q) Z;
W’ Z’ -inv(R);]

@SDconstraint(model,Symmetric(M) <= zeros(size(M)));
@SDconstraint(model,Symmetric(Y) >= zeros(size(Y)));
@objective(model,Max,tr(Y));
println(“Calling optimize!(…)”);
JuMP.optimize!(model);
println(tr(value.(Y)));

MOSEK Output from YALMIP

MOSEK Version 9.0.89 (Build date: 2019-5-24 09:53:18)
Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com
Platform: MACOSX/64-X86

Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 5
Cones : 0
Scalar variables : 0
Matrix variables : 2
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.00
Lin. dep. - number : 0
Presolve terminated. Time: 0.01
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 5
Cones : 0
Scalar variables : 0
Matrix variables : 2
Integer variables : 0

Optimizer - threads : 8
Optimizer - solved problem : the primal
Optimizer - Constraints : 5
Optimizer - Cones : 1
Optimizer - Scalar variables : 3 conic : 3
Optimizer - Semi-definite variables: 1 scalarized : 15
Factor - setup time : 0.00 dense det. time : 0.00
Factor - ML order time : 0.00 GP order time : 0.00
Factor - nonzeros before factor : 15 after factor : 15
Factor - dense dim. : 0 flops : 6.14e+02
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 1.7e+00 9.0e+00 1.3e+01 0.00e+00 1.200000000e+01 0.000000000e+00 1.0e+00 0.02
1 3.6e-01 1.9e+00 4.2e+00 -7.61e-01 7.493854220e+00 1.637944828e+00 2.1e-01 0.06
2 9.2e-02 4.9e-01 6.4e-01 4.11e-01 5.431599917e+00 3.638549110e+00 5.4e-02 0.06
3 2.0e-02 1.1e-01 9.4e-02 2.68e-01 6.802314950e+00 6.369265942e+00 1.2e-02 0.06
4 2.4e-04 1.3e-03 1.2e-04 9.08e-01 6.662079753e+00 6.655739510e+00 1.4e-04 0.06
5 8.9e-06 4.7e-05 8.0e-07 1.00e+00 6.660644505e+00 6.660411659e+00 5.2e-06 0.07
6 8.9e-07 4.7e-06 2.6e-08 1.00e+00 6.660618400e+00 6.660594963e+00 5.2e-07 0.07
7 8.7e-08 4.6e-07 7.8e-10 1.00e+00 6.660615958e+00 6.660613675e+00 5.1e-08 0.07
8 8.7e-09 4.6e-08 2.5e-11 1.00e+00 6.660615311e+00 6.660615082e+00 5.1e-09 0.07
9 3.6e-10 4.4e-09 2.1e-13 1.00e+00 6.660615168e+00 6.660615158e+00 2.1e-10 0.07
Optimizer terminated. Time: 0.10

Interior-point solution summary
Problem status : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal. obj: 6.6606151676e+00 nrm: 4e+00 Viol. con: 4e-09 barvar: 0e+00
Dual. obj: 6.6606151582e+00 nrm: 1e+01 Viol. con: 0e+00 barvar: 2e-08
Optimizer summary
Optimizer - time: 0.10
Interior-point - iterations : 9 time: 0.07
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

6.660615158242211

MOSEK Output from JuMP

Problem
Name :
Objective sense : max
Type : CONIC (conic optimization problem)
Constraints : 18
Cones : 0
Scalar variables : 5
Matrix variables : 2
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 started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries : 2 time : 0.00
Lin. dep. - tries : 1 time : 0.00
Lin. dep. - number : 0
Presolve terminated. Time: 0.00
Problem
Name :
Objective sense : max
Type : CONIC (conic optimization problem)
Constraints : 18
Cones : 0
Scalar variables : 5
Matrix variables : 2
Integer variables : 0

Optimizer - threads : 8
Optimizer - solved problem : the primal
Optimizer - Constraints : 18
Optimizer - Cones : 2
Optimizer - Scalar variables : 9 conic : 9
Optimizer - Semi-definite variables: 1 scalarized : 15
Factor - setup time : 0.00 dense det. time : 0.00
Factor - ML order time : 0.00 GP order time : 0.00
Factor - nonzeros before factor : 153 after factor : 153
Factor - dense dim. : 0 flops : 2.76e+03
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 9.0e+00 1.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.00
1 1.9e+00 2.1e-01 3.3e-01 -7.58e-01 1.494060954e+00 -2.347484915e-01 2.1e-01 0.00
2 4.8e-01 5.4e-02 4.8e-02 3.97e-01 3.216540796e+00 2.621094335e+00 5.4e-02 0.00
3 1.3e-01 1.5e-02 1.1e-02 1.56e-01 6.236366207e+00 5.794301541e+00 1.5e-02 0.00
4 2.8e-03 3.1e-04 3.0e-05 8.65e-01 6.647109017e+00 6.639514147e+00 3.1e-04 0.00
5 6.2e-05 6.9e-06 9.9e-08 1.00e+00 6.660351189e+00 6.660190874e+00 6.9e-06 0.00
6 5.8e-06 6.5e-07 2.8e-09 1.00e+00 6.660595719e+00 6.660580869e+00 6.5e-07 0.00
7 4.9e-07 5.4e-08 6.8e-11 1.00e+00 6.660614264e+00 6.660613014e+00 5.4e-08 0.00
8 4.8e-08 5.3e-09 2.1e-12 1.00e+00 6.660615113e+00 6.660614990e+00 5.3e-09 0.00
9 6.2e-10 6.1e-10 3.1e-15 1.00e+00 6.660615160e+00 6.660615159e+00 6.9e-11 0.00
Optimizer terminated. Time: 0.01

6.660615160346445

Per · November 6, 2019, 7:00am

Looking at the timing for YALMIP:

Optimizer terminated. Time: 0.10

comparted to the timing for JuMP:

Optimizer terminated. Time: 0.01

I would say JuMP did better, even if the number of variables and constraints was higher.

Amin_Yahyaabadi · November 6, 2019, 7:05am

@raktim Nice investigation!
By the way, If you put your codes inside ``` it is easier to read.

blegat · November 6, 2019, 8:29am

You need to replace

@variable(model, Y[1:2, 1:2], Symmetric)
@SDconstraint(model,Symmetric(Y) >= zeros(size(Y)))

By

@variable(model, Y[1:2, 1:2], PSD)

raktim · November 6, 2019, 2:53pm

Improved things a bit. Now we have 15 constraints, 2 scalar variables, and 2 matrix variables. I am more worried about the increased number of constraints.

For the problem I am trying to solve, which is a bit more complex:

JuMP results in:
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 40604
Cones : 0
Scalar variables : 105
Matrix variables : 4
Integer variables : 0

Optimal cost: 1.358
Time: 1087.15 s

YALMIP results in
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 66
Cones : 0
Scalar variables : 1
Matrix variables : 3
Integer variables : 0

Optimal cost: 1.358
Time: 0.55 s

Both return the same answer. But why does JuMP impose so many more constraints. This is affecting solution time significantly. Perhaps YALMIP is more efficient in implementing SD constraints, or I am missing something here.

odow · November 6, 2019, 3:38pm

There are a few reasons why this could be happening.

Please read Please read: make it easier to help you and provide a working example of the code you are trying to improve. It is impossible to offer advice without knowing what you’ve tried!

Make sure to include what versions of Julia/JuMP you are using.

raktim · November 7, 2019, 12:41am

Hi odow,
I posted the code (both YALMIP and JuMP) in the initial posting.

Versions:
JuMP v0.20.1
MathOptInterface v0.9.7
MosekTools v0.9.1
Julia 1.2

raktim · November 7, 2019, 12:45am

Julia code link:

MATLAB YALMIP code link:

The issue is JuMP seems to be imposing more constraints than YALMIP.

odow · November 7, 2019, 8:39pm

You should take a read of the relevant documentation:
https://www.juliaopt.org/JuMP.jl/v0.20.0/constraints/#Semidefinite-constraints-1
https://www.juliaopt.org/JuMP.jl/v0.20.0/constraints/#JuMP.PSDCone

You probably want

a = @constraint(model, Symmetric(-M) in PSDCone())

instead of

b = @SDconstraint(model,Symmetric(M) <= zeros(size(M)))

If you follow the example in the PSDCone documentation, you will see that the @SDconstraint approach is adding many more constraints (as you guess).

julia> jump_function(constraint_object(a))
15-element Array{GenericAffExpr{Float64,VariableRef},1}:
 -2 Y[1,2]                       
 -Y[2,2] + Y[1,1] + Y[1,2] - W[1]
 2 Y[1,2] + 2 Y[2,2] - 2 W[2]    
 -Y[1,1]                         
 -Y[1,2]                         
 1                               
 -Y[1,2]                         
 -Y[2,2]                         
 0                               
 10                              
 -W[1]                           
 -W[2]                           
 0                               
 0                               
 1                               

julia> jump_function(constraint_object(b))
25-element Array{GenericAffExpr{Float64,VariableRef},1}:
 -2 Y[1,2]                       
 -Y[2,2] + Y[1,1] + Y[1,2] - W[1]
 -Y[1,1]                         
 -Y[1,2]                         
 -W[1]                           
 -Y[2,2] + Y[1,1] + Y[1,2] - W[1]
 2 Y[1,2] + 2 Y[2,2] - 2 W[2]    
 -Y[1,2]                         
 -Y[2,2]                         
 -W[2]                           
 -Y[1,1]                         
 -Y[1,2]                         
 1                               
 0                               
 0                               
 -Y[1,2]                         
 -Y[2,2]                         
 0                               
 10                              
 0                               
 -W[1]                           
 -W[2]                           
 0                               
 0                               
 1                               
``

blegat · November 7, 2019, 9:51pm

Conic solvers usually don’t support both conic variables and conic constraints.
They either have the primal form as interface (variables in cones and equality constraints):

min c' x
A x = b
x in K

or the dual form as interface (free variables and and affine constraints in cones):

min c' x
A x + b in K
x free

Mosek for instance uses the primal form.

When you model, your model is sometimes closer to the primal or sometimes to the dual form.
If your model is closer to the primal form, you should give it as is to Mosek.
Otherwise, you should give the dual form to Mosek.
As you have seen, modeling the primal form is error prone as writing @constraint(model, X >= zeros(X)) creates an affine constraints hence need to create slack variables in primal form!
So the user should be more disciplined and create make sure the constraint on the variables are not interpreted as affine constraints.
The modeling language should also allow creating constraints on variables (which is not the case in YALMIP and that can prevent you to create the best model sometimes, see @mforets presentation at JuMP-dev)
For these reasons, YALMIP always use the dual form for modeling and dualize the problem when the solver needs the primal form.

In JuMP, we just give the problem to the solver as it is modeled by the user and if you want to dual to be given instead you have to say it explicitly using Dualization.jl.
In your case, you can do

using JuMP, MosekTools, LinearAlgebra
A = [0. 1; -1 -1];
B = [0.;1];

model = Model(with_optimizer(() -> Dualization.DualOptimizer(Mosek.Optimizer())))
@variable(model, Y[1:2, 1:2], PSD)
@variable(model, W[1:2])

Q = [1 0; 0 0.1]
R = 1
Z = zeros(2,1)

M11 = (A*Y + B*W') + (A*Y + B*W')'
M = [M11 Y W;
     Y -inv(Q) Z;
     W' Z' -inv(R);]

@constraint(model, Symmetric(-M) in PSDCone())
@objective(model, Max, tr(Y))
JuMP.optimize!(model)
println(tr(value.(Y)))

Note that here, doing @constraint(model, Y >= zeros(size(Y))) and @constraint(model, M <= 0) instead would lead the same result as you are modeling in dual form since Mosek is in primal form and you use Dualization.jl.

This gives

Problem
  Name                   :                 
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 5               
  Cones                  : 0               
  Scalar variables       : 0               
  Matrix variables       : 2               
  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.00            
Lin. dep.  - number                 : 0               
Presolve terminated. Time: 0.00    
Problem
  Name                   :                 
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 5               
  Cones                  : 0               
  Scalar variables       : 0               
  Matrix variables       : 2               
  Integer variables      : 0               

Optimizer  - threads                : 4               
Optimizer  - solved problem         : the primal      
Optimizer  - Constraints            : 5
Optimizer  - Cones                  : 1
Optimizer  - Scalar variables       : 3                 conic                  : 3               
Optimizer  - Semi-definite variables: 1                 scalarized             : 15              
Factor     - setup time             : 0.00              dense det. time        : 0.00            
Factor     - ML order time          : 0.00              GP order time          : 0.00            
Factor     - nonzeros before factor : 15                after factor           : 15              
Factor     - dense dim.             : 0                 flops                  : 6.14e+02        
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  
0   1.7e+00  9.0e+00  1.3e+01  0.00e+00   1.200000000e+01   0.000000000e+00   1.0e+00  0.00  
1   3.6e-01  1.9e+00  4.2e+00  -7.61e-01  7.493854220e+00   1.637944828e+00   2.1e-01  0.00  
2   9.2e-02  4.9e-01  6.4e-01  4.11e-01   5.431599917e+00   3.638549110e+00   5.4e-02  0.00  
3   2.0e-02  1.1e-01  9.4e-02  2.68e-01   6.802314950e+00   6.369265942e+00   1.2e-02  0.00  
4   2.4e-04  1.3e-03  1.2e-04  9.08e-01   6.662079753e+00   6.655739510e+00   1.4e-04  0.00  
5   8.9e-06  4.7e-05  8.0e-07  1.00e+00   6.660644505e+00   6.660411659e+00   5.2e-06  0.00  
6   8.9e-07  4.7e-06  2.6e-08  1.00e+00   6.660618400e+00   6.660594963e+00   5.2e-07  0.00  
7   8.7e-08  4.6e-07  7.8e-10  1.00e+00   6.660615958e+00   6.660613675e+00   5.1e-08  0.00  
8   8.7e-09  4.6e-08  2.5e-11  1.00e+00   6.660615311e+00   6.660615082e+00   5.1e-09  0.00  
9   3.6e-10  2.5e-09  2.1e-13  1.00e+00   6.660615168e+00   6.660615158e+00   2.1e-10  0.00  
Optimizer terminated. Time: 0.00    

6.660615158242325

as with YALMIP

raktim · November 7, 2019, 10:14pm

This makes things clear! Thank you for taking time to address this.

odow · November 7, 2019, 11:53pm

I feel like this is great content for the JuMP documentation.

raktim · November 8, 2019, 12:34am

I think so too! Coming from cvx/YALMIP, I find it hard to transfer to JuMP. I am documenting these, especially for folks coming from control/estimation and are used to cvx/YALMIP. Perhaps it can be included with JuMP examples in the future.

With Dualization, things are better. However I am still running into difficulties in implementing constraints such as s\ge0 where s\in\mathcal{R}^n.

The documentation for Dualization is probably aimed at developers at this time.
I have to spend some time on this to figure these basic things.

It may be better to dualize the problem just before calling \texttt{optimize!(...)}. But then I don’t know how to recover the primal variables from the dual solution.
Where can I get information about this?

I am going to mark this topic as resolved.

blegat · November 8, 2019, 7:19am

recover the primal variables from the dual solution.

Dualization dualize both the model and the results so you can just get the primal solution and Dualization take care of asking the dual of the corresponding constraint instead.