# Transitioning from MATLAB

*This tutorial was generated using Literate.jl.* *Download the source as a .jl file*.

YALMIP and CVX are two packages for mathematical optimization in MATLAB®. They are independently developed and are in no way affiliated with JuMP.

The purpose of this tutorial is to help new users to JuMP who have previously used YALMIP or CVX by comparing and contrasting their different features.

If you have not used Julia before, read the Getting started with Julia tutorial.

## Namespaces

Julia has namespaces, which MATLAB lacks. Therefore one needs to either use the command:

`using JuMP`

in order bring all names exported by JuMP into scope, or:

`import JuMP`

in order to merely make the JuMP package available. `import`

requires prefixing everything you use from JuMP with `JuMP.`

. In this tutorial we use the former.

## Models

YALMIP and CVX have a single, implicit optimization model that you build by defining variables and constraints.

In JuMP, we create an explicit model first, and then, when you declare variables, constraints, or the objective function, you specify to which model they are being added.

Create a new JuMP model with the command:

`model = Model()`

```
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.
```

## Variables

In most cases there is a direct translation between variable declarations. The following table shows some common examples:

JuMP | YALMIP | CVX |
---|---|---|

`@variable(model, x)` | `x = sdpvar` | `variable x` |

`@variable(model, x, Int)` | `x = intvar` | `variable x integer` |

`@variable(model, x, Bin)` | `x = binvar` | `variable x binary` |

`@variable(model, v[1:d])` | `v = sdpvar(d, 1)` | `variable v(d)` |

`@variable(model, m[1:d, 1:d])` | `m = sdpvar(d,d,'full')` | `variable m(d, d)` |

`@variable(model, m[1:d, 1:d] in ComplexPlane())` | `m = sdpvar(d,d,'full','complex')` | `variable m(d,d) complex` |

`@variable(model, m[1:d, 1:d], Symmetric)` | `m = sdpvar(d)` | `variable m(d,d) symmetric` |

`@variable(model, m[1:d, 1:d], Hermitian)` | `m = sdpvar(d,d,'hermitian','complex')` | `variable m(d,d) hermitian` |

Like CVX, but unlike YALMIP, JuMP can also constrain variables upon creation:

JuMP | CVX |
---|---|

`@variable(model, v[1:d] >= 0)` | `variable v(d) nonnegative` |

`@variable(model, m[1:d, 1:d], PSD)` | `variable m(d,d) semidefinite` |

`@variable(model, m[1:d, 1:d] in PSDCone())` | `variable m(d,d) semidefinite` |

`@variable(model, m[1:d, 1:d] in HermitianPSDCone())` | `variable m(d,d) complex semidefinite` |

JuMP can additionally set variable bounds, which may be handled more efficiently by a solver than an equivalent linear constraint. For example:

```
@variable(model, -1 <= x[i in 1:3] <= i)
upper_bound.(x)
```

```
3-element Vector{Float64}:
1.0
2.0
3.0
```

A more interesting case is when you want to declare, for example, `n`

real symmetric matrices. Both YALMIP and CVX allow you to put the matrices as the slices of a 3-dimensional array, via the commands `m = sdpvar(d, d, n)`

and `variable m(d, d, n) symmetric`

, respectively. With JuMP this is not possible. Instead, to achieve the same result one needs to declare a vector of `n`

matrices:

```
d, n = 3, 2
m = [@variable(model, [1:d, 1:d], Symmetric) for _ in 1:n]
```

```
2-element Vector{LinearAlgebra.Symmetric{VariableRef, Matrix{VariableRef}}}:
[_[4] _[5] _[7]; _[5] _[6] _[8]; _[7] _[8] _[9]]
[_[10] _[11] _[13]; _[11] _[12] _[14]; _[13] _[14] _[15]]
```

`m[1]`

```
3×3 LinearAlgebra.Symmetric{VariableRef, Matrix{VariableRef}}:
_[4] _[5] _[7]
_[5] _[6] _[8]
_[7] _[8] _[9]
```

`m[2]`

```
3×3 LinearAlgebra.Symmetric{VariableRef, Matrix{VariableRef}}:
_[10] _[11] _[13]
_[11] _[12] _[14]
_[13] _[14] _[15]
```

The analogous construct in MATLAB would be a cell array containing the optimization variables, which every discerning programmer avoids as cell arrays are rather slow. This is not a problem in Julia: a vector of matrices is almost as fast as a 3-dimensional array.

## Constraints

As in the case of variables, in most cases there is a direct translation between the packages:

JuMP | YALMIP | CVX |
---|---|---|

`@constraint(model, v == c)` | `v == c` | `v == c` |

`@constraint(model, v >= 0)` | `v >= 0` | `v >= 0` |

`@constraint(model, m >= 0, PSDCone())` | `m >= 0` | `m == semidefinite(length(m))` |

`@constraint(model, m >= 0, HermitianPSDCone())` | `m >= 0` | `m == hermitian_semidefinite(length(m))` |

`@constraint(model, [t; v] in SecondOrderCone())` | `cone(v, t)` | `{v, t} == lorentz(length(v))` |

`@constraint(model, [x, y, z] in MOI.ExponentialCone())` | `expcone([x, y, z])` | `{x, y, z} == exponential(1)` |

A subtlety appears when declaring equality constraints for matrices. In general, JuMP uses `@constraint(model, m .== c)`

, with the dot meaning broadcasting in Julia, except when `m`

is `Symmetric`

or `Hermitian`

: in this case `@constraint(model, m == c)`

is allowed, and is much better, as JuMP is smart enough to not generate redundant constraints for the lower diagonal and the imaginary part of the diagonal (in the complex case). Both YALMIP and CVX are also smart enough to do this and the syntax is always just `m == c`

.

Experienced YALMIP users will probably be relieved to see that you must pass `PSDCone()`

or `HermitianPSDCone()`

to make a matrix positive semidefinite, as the `>=`

ambiguity in YALMIP is common source of bugs.

## Setting the objective

Like CVX, but unlike YALMIP, JuMP has a specific command for setting an objective function:

`@objective(model, Min, sum(i * x[i] for i in 1:3))`

\[ x_{1} + 2 x_{2} + 3 x_{3} \]

Here the third argument is any expression you want to optimize, and `Min`

is an objective sense (the other possibility is `Max`

).

## Setting solver and options

In order to set an optimizer with JuMP, do:

```
import Clarabel
set_optimizer(model, Clarabel.Optimizer)
```

where "Clarabel" is an example solver. See the list of Supported solvers for other choices.

To configure the solver options you use the command:

`set_attribute(model, "verbose", true)`

where `verbose`

is an option specific to Clarabel.

A crucial difference is that with JuMP you must explicitly choose a solver before optimizing. Both YALMIP and CVX allow you to leave it empty and will try to guess an appropriate solver for the problem.

## Optimizing

Like YALMIP, but unlike CVX, with JuMP you need to explicitly start the optimization, with the command:

`optimize!(model)`

```
-------------------------------------------------------------
Clarabel.jl v0.7.1 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 15
constraints = 6
nnz(P) = 0
nnz(A) = 6
cones (total) = 1
: Nonnegative = 1, numel = 6
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-05, max_scale = 1.0e+05
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 1.0000e+01 -1.2500e+01 2.25e+00 0.00e+00 0.00e+00 1.00e+00 3.36e+00 ------
1 3.9744e+00 -5.5968e-01 4.53e+00 1.43e-16 1.27e-16 3.10e-01 6.92e-01 8.38e-01
2 1.1590e-01 -1.2437e-01 2.40e-01 4.88e-17 3.27e-17 2.81e-02 3.83e-02 9.73e-01
3 1.1746e-03 -1.2507e-03 2.43e-03 1.06e-16 7.36e-17 2.83e-04 3.87e-04 9.90e-01
4 1.1746e-05 -1.2507e-05 2.43e-05 1.44e-16 3.68e-17 2.83e-06 3.87e-06 9.90e-01
5 1.1746e-07 -1.2507e-07 2.43e-07 5.05e-15 4.78e-15 2.83e-08 3.87e-08 9.90e-01
6 1.1746e-09 -1.2507e-09 2.43e-09 1.59e-16 6.59e-17 2.83e-10 3.87e-10 9.90e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 1.25ms
```

The exclamation mark here is a Julia-ism that means the function is modifying its argument, `model`

.

## Querying solution status

After the optimization is done, you should check for the solution status to see what solution (if any) the solver found.

Like YALMIP and CVX, JuMP provides a solver-independent way to check it, via the command:

`is_solved_and_feasible(model)`

`true`

If the return value is `false`

, you should investigate with `termination_status`

, `primal_status`

, and `raw_status`

, See Solutions for more details on how to query and interpret solution statuses.

## Extracting variables

Like YALMIP, but unlike CVX, with JuMP you need to explicitly ask for the value of your variables after optimization is done, with the function call `value(x)`

to obtain the value of variable `x`

.

`value.(m[1][1, 1])`

`0.0`

A subtlety is that, unlike YALMIP, the function `value`

is only defined for scalars. For vectors and matrices you need to use Julia broadcasting: `value.(v)`

.

`value.(m[1])`

```
3×3 Matrix{Float64}:
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
```

There is also a specialized function for extracting the value of the objective, `objective_value(model)`

, which is useful if your objective doesn't have a convenient expression.

`objective_value(model)`

`-5.999999998825352`

## Dual variables

Like YALMIP and CVX, JuMP allows you to recover the dual variables. In order to do that, the simplest method is to name the constraint you're interested in, for example, `@constraint(model, bob, sum(v) == 1)`

and then, after the optimzation is done, call `dual(bob)`

. See Duality for more details.

## Reformulating problems

Perhaps the biggest difference between JuMP and YALMIP and CVX is how far the package is willing to go in reformulating the problems you give to it.

CVX is happy to reformulate anything it can, even using approximations if your solver cannot handle the problem.

YALMIP will only do exact reformulations, but is still fairly adventurous, for example, being willing to reformulate a nonlinear objective in terms of conic constraints.

JuMP does no such thing: it only reformulates objectives into objectives, and constraints into constraints, and is fairly conservative at that. As a result, you might need to do some reformulations manually, for which a good guide is the Tips and tricks tutorial.

## Vectorization

In MATLAB, it is absolutely essential to "vectorize" your code to obtain acceptable performance. This is because MATLAB is a slow interpreted language, which sends your commands to fast libraries. When you "vectorize" your code you are minimizing the MATLAB part of the work and sending it to the fast libraries instead.

There's no such duality with Julia.

Everything you write and most libraries you use will compile down to LLVM, so "vectorization" has no effect.

For example, if you are writing a linear program in MATLAB and instead of the usual `constraints = [v >= 0]`

you write:

```
for i = 1:n
constraints = [constraints, v(i) >= 0];
end
```

performance will be poor.

With Julia, on the other hand, there is hardly any difference between

`@constraint(model, v >= 0)`

and

```
for i in 1:n
@constraint(model, v[i] >= 0)
end
```

## Symmetric and Hermitian matrices

Julia has specialized support for symmetric and Hermitian matrices in the `LinearAlgebra`

package:

`import LinearAlgebra`

If you have a matrix that is numerically symmetric:

`x = [1 2; 2 3]`

```
2×2 Matrix{Int64}:
1 2
2 3
```

`LinearAlgebra.issymmetric(x)`

`true`

then you can wrap it in a `LinearAlgebra.Symmetric`

matrix to tell Julia's type system that the matrix is symmetric.

`LinearAlgebra.Symmetric(x)`

```
2×2 LinearAlgebra.Symmetric{Int64, Matrix{Int64}}:
1 2
2 3
```

Using a `Symmetric`

matrix lets Julia and JuMP use more efficient algorithms when they are working with symmetric matrices.

If you have a matrix that is nearly but not exactly symmetric:

```
x = [1.0 2.0; 2.001 3.0]
LinearAlgebra.issymmetric(x)
```

`false`

then you could, as you might do in MATLAB, make it numerically symmetric as follows:

`x_sym = 0.5 * (x + x')`

```
2×2 Matrix{Float64}:
1.0 2.0005
2.0005 3.0
```

In Julia, you can explicitly choose whether to use the lower or upper triangle of the matrix:

`x_sym = LinearAlgebra.Symmetric(x, :L)`

```
2×2 LinearAlgebra.Symmetric{Float64, Matrix{Float64}}:
1.0 2.001
2.001 3.0
```

`x_sym = LinearAlgebra.Symmetric(x, :U)`

```
2×2 LinearAlgebra.Symmetric{Float64, Matrix{Float64}}:
1.0 2.0
2.0 3.0
```

The same applies for Hermitian matrices, using `LinearAlgebra.Hermitian`

and `LinearAlgebra.ishermitian`

.

## Primal versus dual form

When you translate some optimization problems from YALMIP or CVX to JuMP, you might be surprised to see it get much faster or much slower, even if you're using exactly the same solver. The most likely reason is that YALMIP will always interpret the problem as the dual form, whereas CVX and JuMP will try to interpret the problem in the form most appropriate to the solver. If the problem is more naturally formulated in the primal form it is likely that YALMIP's performance will suffer, or if JuMP gets it wrong, its performance will suffer. It might be worth trying both primal and dual forms if you're having trouble, which can be done automatically with the package Dualization.jl.

For an in-depth explanation of this issue, see the Dualization tutorial.

## Rosetta stone

In this section, we show a complete example of the same optimization problem being solved with JuMP, YALMIP, and CVX. It is a semidefinite program that computes a lower bound on the random robustness of entanglement using the partial transposition criterion.

The code is complete, apart from the function that does partial transposition. With both YALMIP and CVX we use the function `PartialTranspose`

from QETLAB. With JuMP, we could use the function `Convex.partialtranspose`

from Convex.jl, but we reproduce it here for simplicity:

```
function partial_transpose(x::AbstractMatrix, sys::Int, dims::Vector)
@assert size(x, 1) == size(x, 2) == prod(dims)
@assert 1 <= sys <= length(dims)
n = length(dims)
s = n - sys + 1
p = collect(1:2n)
p[s], p[n+s] = n + s, s
r = reshape(x, (reverse(dims)..., reverse(dims)...))
return reshape(permutedims(r, p), size(x))
end
```

`partial_transpose (generic function with 1 method)`

### JuMP

The JuMP code to solve this problem is:

```
using JuMP
import Clarabel
import LinearAlgebra
function random_state_pure(d)
x = randn(Complex{Float64}, d)
y = x * x'
return LinearAlgebra.Hermitian(y / LinearAlgebra.tr(y))
end
function robustness_jump(d)
rho = random_state_pure(d^2)
id = LinearAlgebra.Hermitian(LinearAlgebra.I(d^2))
rhoT = LinearAlgebra.Hermitian(partial_transpose(rho, 1, [d, d]))
model = Model()
@variable(model, λ)
@constraint(model, PPT, rhoT + λ * id in HermitianPSDCone())
@objective(model, Min, λ)
set_optimizer(model, Clarabel.Optimizer)
set_attribute(model, "verbose", true)
optimize!(model)
if is_solved_and_feasible(model)
WT = dual(PPT)
return value(λ), real(LinearAlgebra.dot(WT, rhoT))
else
return "Something went wrong: $(raw_status(model))"
end
end
robustness_jump(3)
```

`(0.4317497878154882, -0.4317497873162501)`

### YALMIP

The corresponding YALMIP code is:

```
function robustness_yalmip(d)
rho = random_state_pure(d^2);
# PartialTranspose from https://github.com/nathanieljohnston/QETLAB
rhoT = PartialTranspose(rho, 1, [d d]);
lambda = sdpvar;
constraints = [(rhoT + lambda*eye(d^2) >= 0):'PPT'];
ops = sdpsettings(sdpsettings, 'verbose', 1, 'solver', 'sedumi');
sol = optimize(constraints, lambda, ops);
if sol.problem == 0
WT = dual(constraints('PPT'));
value(lambda)
real(WT(:).' * rhoT(:))
else
display(['Something went wrong: ', sol.info])
end
end
function rho = random_state_pure(d)
x = randn(d, 1) + 1i * randn(d, 1);
y = x * x';
rho = y / trace(y);
end
```

### CVX

The corresponding CVX code is:

```
function robustness_cvx(d)
rho = random_state_pure(d^2);
# PartialTranspose from https://github.com/nathanieljohnston/QETLAB
rhoT = PartialTranspose(rho, 1, [d d]);
cvx_begin
variable lambda
dual variable WT
WT : rhoT + lambda * eye(d^2) == hermitian_semidefinite(d^2)
minimise lambda
cvx_end
if strcmp(cvx_status, 'Solved')
lambda
real(WT(:)' * rhoT(:))
else
display('Something went wrong.')
end
end
function rho = random_state_pure(d)
x = randn(d, 1) + 1i * randn(d, 1);
y = x * x';
rho = y / trace(y);
end
```