Документация Engee

Переход с MATLAB

YALMIP и CVX — это два пакета для математической оптимизации в MATLAB®. Они разрабатываются независимо и никак не связаны с JuMP.

Цель данного руководства — помочь новым пользователям JuMP, которые ранее использовали YALMIP или CVX, путем сравнения и сопоставления возможностей.

Если у вас нет опыта работы с Julia, прочитайте руководство Getting started with Julia.

Пространства имен

В Julia есть пространства имен, которых нет в MATLAB. Поэтому нужно использовать либо команду:

using JuMP

чтобы ввести в область все имена, экспортируемые JuMP, либо команду:

import JuMP

чтобы просто сделать доступным пакет JuMP. import требует добавления префикса JuMP. ко всем используемым компонентам из JuMP. В данном руководстве применяется первый вариант.

Модели

В YALMIP и CVX имеется единая неявная модель оптимизации, которая создается путем определения переменных и ограничений.

В JuMP сначала создается явная модель, а затем при объявлении переменных, ограничений или целевой функции указывается, в какую модель они добавляются.

Модель JuMP создается с помощью следующей команды:

model = Model()
A JuMP Model
├ solver: none
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none

Переменные

В большинстве случаев объявления переменных преобразуются напрямую. В следующей таблице представлены некоторые стандартные примеры:

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

Как и в CVX, но в отличие от YALMIP, в JuMP также можно ограничивать переменные при создании:

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 можно дополнительно устанавливать границы переменных, которые могут обрабатываться решателем более эффективно, чем эквивалентное линейное ограничение. Например:

@variable(model, -1 <= x[i in 1:3] <= i)
upper_bound.(x)
3-element Vector{Float64}:
 1.0
 2.0
 3.0

Более интересный случай возникает, когда нужно объявить, например, n вещественных симметричных матриц. И YALMIP, и CVX позволяют добавлять матрицы в виде срезов трехмерного массива с помощью команд m = sdpvar(d, d, n) и variable m(d, d, n) symmetric соответственно. В JuMP это невозможно. Вместо этого для достижения того же результата необходимо объявить вектор из n матриц:

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]

Аналогичной конструкцией в MATLAB был бы массив ячеек, содержащий переменные оптимизации, чего любой грамотный программист избегает, так как массивы ячеек обрабатываются довольно медленно. В Julia это не проблема: вектор матриц обрабатывается почти так же быстро, как трехмерный массив.

Ограничения

Как и в случае с переменными, в большинстве случаев существует прямое преобразование между пакетами:

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)

Как и в YALMIP и CVX, в JuMP есть механизм, позволяющий не создавать избыточные ограничения при объявлении ограничений равенства между матрицами типа Symmetric или Hermitian. В этих случаях @constraint(model, m == c) не будет генерировать ограничения для нижней диагонали и мнимой части диагонали (в комплексном случае).

Опытные пользователи MATLAB, вероятно, будут рады узнать, что для того, чтобы сделать матрицу положительно полуопределенной, необходимо передать PSDCone() или HermitianPSDCone(), так как неоднозначность >= в YALMIP и CVX часто вызывает ошибки.

Задание целевой функции

Как и в CVX, но в отличие от YALMIP, в JuMP есть специальная команда для задания целевой функции:

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

Здесь третий аргумент — это любое выражение, которое нужно оптимизировать, а Min — это целевое назначение (другой вариант — Max).

Задание решателя и параметров

Задать оптимизатор в JuMP можно следующим образом:

import Clarabel
set_optimizer(model, Clarabel.Optimizer)

где Clarabel — это пример решателя. Другие варианты см. в списке в разделе Supported solvers.

Параметры решателя настраиваются с помощью следующей команды:

set_attribute(model, "verbose", true)

где verbose — это параметр, относящийся к Clarabel.

Принципиальное отличие заключается в том, что в JuMP необходимо явно выбрать решатель перед оптимизацией. И YALMIP, и CVX позволяют оставить этот параметр пустым, чтобы они сами могли попытаться подобрать подходящий решатель для задачи.

Оптимизация

Как и в YALMIP, но в отличие от CVX, в JuMP необходимо явным образом запустить оптимизацию с помощью следующей команды:

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  9.46e-17  9.32e-17  2.83e-04  3.87e-04  9.90e-01
  4   1.1746e-05  -1.2507e-05  2.43e-05  1.72e-16  7.55e-17  2.83e-06  3.87e-06  9.90e-01
  5   1.1746e-07  -1.2507e-07  2.43e-07  5.02e-15  4.74e-15  2.83e-08  3.87e-08  9.90e-01
  6   1.1746e-09  -1.2507e-09  2.43e-09  1.18e-16  9.32e-17  2.83e-10  3.87e-10  9.90e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 2.28ms

Восклицательный знак здесь — это характерный для Julia символ, означающий, что функция изменяет свой аргумент, в данном случае model.

Запрос состояния решения

По завершении оптимизации следует проверить состояние решения, чтобы узнать, какое решение (если таковое имеется) нашел решатель.

Как и в YALMIP и CVX, в JuMP доступен независимый от решателя способ такой проверки с помощью следующей команды:

is_solved_and_feasible(model)
true

Если возвращается значение false, следует изучить состояние с помощью termination_status, primal_status и raw_status. Дополнительные сведения о запросе и интерпретации состояний решений см. в разделе Решения.

Извлечение переменных

Как и в YALMIP, но в отличие от CVX, в JuMP необходимо явным образом запрашивать значения переменных после завершения оптимизации с помощью вызова функции value(x), возвращающей значение переменной x.

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

Тонкость в том, что, в отличие от YALMIP, функция value определена только для скалярных значений. Для векторов и матриц необходимо использовать трансляцию Julia: 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

Существует также специализированная функция для извлечения значения целевой функции, objective_value(model), которая полезна, если у целевой функции нет удобного выражения.

objective_value(model)
-5.999999998825352

Двойственные переменные

Так же как YALMIP и CVX, JuMP позволяет восстанавливать двойные переменные. Для этого проще всего присвоить имя интересующему вас ограничению, например @constraint(model, bob, sum(v) == 1), а затем после завершения оптимизации вызвать dual(bob). Дополнительные сведения см. в разделе Duality.

Переформулировка задач

Возможно, самое большое различие между JuMP с одной стороны и YALMIP и CVX с другой заключается в том, насколько сильно пакет может переформулировать предоставленные задачи.

CVX легко переформулирует все, что можно, даже используя аппроксимации, если решатель не справляется с задачей.

YALMIP выполняет только точные переформулировки, но все же делает это достаточно свободно, например может переформулировать нелинейную задачу в терминах конических ограничений.

В JuMP такое невозможно: целевые функции переформулируются только в целевые функции, а ограничения — в ограничения, причем довольно осторожно. В результате некоторые переформулировки может потребоваться произвести вручную. Для этого можно прибегнуть к руководству Советы.

Векторизация

В MATLAB код совершенно необходимо «векторизовать», чтобы получить приемлемую производительность. Причина в том, что MATLAB — медленный интерпретируемый язык, который отправляет ваши команды быстрым библиотекам. Когда вы «векторизуете» код, то сводите к минимуму часть работы, приходящуюся на MATLAB, и вместо этого перекладываете задачи на быстрые библиотеки.

В Julia такой двойственности нет.

Весь написанный вами код и большинство используемых библиотек будут скомпилированы в LLVM, поэтому «векторизация» не имеет никакого эффекта.

Например, если вы создаете линейную программу в MATLAB и вместо обычного выражения constraints = [v >= 0] пишете следующее:

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

производительность будет низкой.

В Julia же практически нет разницы между

@constraint(model, v >= 0)

и

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

Симметричные и эрмитовы матрицы

В Julia поддержка симметричных и эрмитовых матриц реализована особо в пакете LinearAlgebra:

import LinearAlgebra

Если имеется численно симметричная матрица:

x = [1 2; 2 3]
2×2 Matrix{Int64}:
 1  2
 2  3
LinearAlgebra.issymmetric(x)
true

то ее можно заключить в матрицу LinearAlgebra.Symmetric, чтобы сообщить системе типов Julia, что матрица симметричная.

LinearAlgebra.Symmetric(x)
2×2 LinearAlgebra.Symmetric{Int64, Matrix{Int64}}:
 1  2
 2  3

Использование матрицы типа Symmetric позволяет Julia и JuMP применять более эффективные алгоритмы при работе с симметричными матрицами.

Если есть матрица, которая почти, но не совсем симметрична:

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

то, как и в MATLAB, ее можно сделать численно симметричной следующим образом:

x_sym = 0.5 * (x + x')
2×2 Matrix{Float64}:
 1.0     2.0005
 2.0005  3.0

В Julia можно явным образом выбрать, следует ли использовать нижний или верхний треугольник матрицы:

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

То же самое касается эрмитовых матриц, для которых используются LinearAlgebra.Hermitian и LinearAlgebra.ishermitian.

Прямая и двойственная формы

При переносе некоторых задач оптимизации из YALMIP или CVX в JuMP к вашему удивлению может оказаться, что они стали работать намного быстрее или намного медленнее, даже если применяется тот же самый решатель. Скорее всего, причина заключается в том, что YALMIP всегда интерпретирует задачу как имеющую двойственную форму, тогда как CVX и JuMP пытаются интерпретировать ее в наиболее подходящей для решателя форме. Если задача более естественно формулируется в прямой форме, то, скорее всего, производительность в YALMIP будет ниже, а если в JuMP задача будет распознана неправильно, то и это скажется на производительности. Если возникают проблемы, возможно, стоит попробовать как прямую, так и двойственную формы, что можно сделать автоматически с помощью пакета Dualization.jl.

Подробное объяснение этой проблемы см. в руководстве Dualization.

Розеттский камень

В этом разделе показан пример полного решения одной и той же задачи оптимизации с помощью JuMP, YALMIP и CVX. Это полуопределенная программа, которая вычисляет нижнюю границу случайной устойчивости запутанности, используя критерий частной транспозиции.

Код является полным, за исключением функции, которая выполняет частную транспозицию. И в YALMIP, и в CVX используется функция PartialTranspose из QETLAB. В JuMP можно было бы использовать функцию Convex.partialtranspose из Convex.jl, но для простоты мы воспроизведем ее здесь:

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

Код JuMP для решения этой задачи:

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.4317497878154874, -0.4317497873162498)

YALMIP

Соответствующий код YALMIP:

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

Соответствующий код CVX:

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