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.

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
├ solver: none
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none

Variables

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

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

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:

Like YALMIP and CVX, JuMP is smart enough to not generate redundant constraints when declaring equality constraints between Symmetric or Hermitian matrices. In these cases @constraint(model, m == c) will not generate constraints for the lower diagonal and the imaginary part of the diagonal (in the complex case).

Experienced MATLAB users will probably be relieved to see that you must pass PSDCone() or HermitianPSDCone() to make a matrix positive semidefinite, because the >= ambiguity in YALMIP and CVX 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.9.0  -  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-04, max_scale = 1.0e+04
               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 =  657μs

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 Modeling with cones 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.43174978781546347, -0.43174978731622865)

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