Пользовательские гессианы

В этом руководстве объясняется, как написать пользовательский оператор (см. раздел Пользовательские операторы) с гессиановой матрицей, явным образом предоставленной пользователем.

Более сложный пример см. в разделе Nested optimization problems.

В этом руководстве используются следующие пакеты:

using JuMP
import Ipopt

Пример функции Розенброка

В качестве простого примера рассмотрим функцию Розенброка:

rosenbrock(x...) = (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2

rosenbrock (generic function with 1 method)

с градиентным вектором:

function ∇rosenbrock(g::AbstractVector, x...)
    g[1] = 400 * x[1]^3 - 400 * x[1] * x[2] + 2 * x[1] - 2
    g[2] = 200 * (x[2] - x[1]^2)
    return
end

∇rosenbrock (generic function with 1 method)

и гессиановой матрицей:

function ∇²rosenbrock(H::AbstractMatrix, x...)
    H[1, 1] = 1200 * x[1]^2 - 400 * x[2] + 2
    # H[1, 2] = -400 * x[1] <-- не требуется, так как гессиан симметричен
    H[2, 1] = -400 * x[1]
    H[2, 2] = 200.0
    return
end

∇²rosenbrock (generic function with 1 method)

Предполагается, что гессианова матрица H инициализирована нулями, и поскольку она симметрична, достаточно заполнить ненулевыми значениями элементы нижнего треугольника.

Тип передаваемой матрицы H зависит от системы автоматического дифференцирования, поэтому первый аргумент функции гессиана должен поддерживать тип AbstractMatrix (он может быть отличен от Matrix{Float64}). Однако изначально предполагается только то, что H поддерживает size(H) и setindex!.

Теперь, когда у нас есть функция, ее градиент и гессиан, мы можем построить модель JuMP, добавить оператор и использовать его в макросе:

model = Model(Ipopt.Optimizer)
@variable(model, x[1:2])
@operator(model, op_rosenbrock, 2, rosenbrock, ∇rosenbrock, ∇²rosenbrock)
@objective(model, Min, op_rosenbrock(x[1], x[2]))
optimize!(model)
@assert is_solved_and_feasible(model)
solution_summary(model; verbose = true)

* Solver : Ipopt

* Status
  Result count       : 1
  Termination status : LOCALLY_SOLVED
  Message from the solver:
  "Solve_Succeeded"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 8.14943e-28
  Dual objective value : 0.00000e+00
  Primal solution :
    x[1] : 1.00000e+00
    x[2] : 1.00000e+00
  Dual solution :

* Work counters
  Solve time (sec)   : 4.06561e-02
  Barrier iterations : 14