Вычисление гессианов

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

Это руководство для опытных пользователей, предполагающее взаимодействие с низкоуровневым нелинейным интерфейсом MathOptInterface.

JuMP по умолчанию экспортирует символ `MOI` как псевдоним для пакета

MathOptInterface.jl. Мы рекомендуем сделать это более явным в вашем коде, добавив следующие строки: julia import MathOptInterface as MOI

Для нелинейной программы:

гессиан лагранжиана вычисляется следующим образом:

где — прямая точка, — скаляр (обычно ), а — вектор весов, соответствующий лагранжевой двойственной задаче для ограничений.

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

using JuMP
import Ipopt
import LinearAlgebra
import Random
import SparseArrays

Базовая модель

Для демонстрации взаимодействия с низкоуровневым нелинейным интерфейсом понадобится пример модели. Точная форма модели не важна; мы используем модель из руководства The Rosenbrock function с некоторыми дополнительными ограничениями для демонстрации различных функций низкоуровневого интерфейса.

model = Model(Ipopt.Optimizer)
set_silent(model)
@variable(model, x[i = 1:2], start = -i)
@constraint(model, g_1, x[1]^2 <= 1)
@constraint(model, g_2, (x[1] + x[2])^2 <= 2)
@objective(model, Min, (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2)
optimize!(model)
@assert is_solved_and_feasible(model)

Аналитическое решение

Приложив немного усилий, можно аналитически вывести правильный гессиан:

function analytic_hessian(x, σ, μ)
    g_1_H = [2.0 0.0; 0.0 0.0]
    g_2_H = [2.0 2.0; 2.0 2.0]
    f_H = zeros(2, 2)
    f_H[1, 1] = 2.0 + 1200.0 * x[1]^2 - 400.0 * x[2]
    f_H[1, 2] = f_H[2, 1] = -400.0 * x[1]
    f_H[2, 2] = 200.0
    return σ * f_H + μ' * [g_1_H, g_2_H]
end

analytic_hessian (generic function with 1 method)

Вот различные точки:

analytic_hessian([1, 1], 0, [0, 0])

2×2 Matrix{Float64}:
 0.0  0.0
 0.0  0.0

analytic_hessian([1, 1], 0, [1, 0])

2×2 Matrix{Float64}:
 2.0  0.0
 0.0  0.0

analytic_hessian([1, 1], 0, [0, 1])

2×2 Matrix{Float64}:
 2.0  2.0
 2.0  2.0

analytic_hessian([1, 1], 1, [0, 0])

2×2 Matrix{Float64}:
  802.0  -400.0
 -400.0   200.0

Создание нелинейной модели

JuMP возлагает задачу автоматического дифференцирования на подмодуль MOI.Nonlinear. Следовательно, чтобы вычислить гессиан лагранжиана, необходимо создать объект MOI.Nonlinear.Model:

rows = Any[]
nlp = MOI.Nonlinear.Model()
for (F, S) in list_of_constraint_types(model)
    if F <: VariableRef
        continue  # Пропускаем границы переменных
    end
    for ci in all_constraints(model, F, S)
        push!(rows, ci)
        object = constraint_object(ci)
        MOI.Nonlinear.add_constraint(nlp, object.func, object.set)
    end
end
MOI.Nonlinear.set_objective(nlp, objective_function(model))
nlp

A Nonlinear.Model with:
 1 objective
 0 parameters
 0 expressions
 2 constraints

Важно сохранить индексы ограничений в векторе rows, чтобы знать порядок ограничений в нелинейной модели.

Далее необходимо преобразовать модель в MOI.Nonlinear.Evaluator, указав бэкенд автоматического дифференцирования. В данном случае мы используем MOI.Nonlinear.SparseReverseMode:

evaluator = MOI.Nonlinear.Evaluator(
    nlp,
    MOI.Nonlinear.SparseReverseMode(),
    index.(all_variables(model)),
)

Nonlinear.Evaluator with available features:

  * :Grad

  * :Jac

  * :JacVec

  * :Hess

  * :HessVec

  * :ExprGraph

Прежде чем использовать средство оценки для вычислений, его необходимо инициализировать. Чтобы узнать, что можно запросить, используйте MOI.features_available:

MOI.features_available(evaluator)

6-element Vector{Symbol}:
 :Grad
 :Jac
 :JacVec
 :Hess
 :HessVec
 :ExprGraph

Подробности см. в документации по MOI, но для получения гессиановой матрицы необходимо инициализировать :Hess:

MOI.initialize(evaluator, [:Hess])

MOI представляет гессиан как разреженную матрицу. Получаем схему разреженности следующим образом.

hessian_sparsity = MOI.hessian_lagrangian_structure(evaluator)

7-element Vector{Tuple{Int64, Int64}}:
 (1, 1)
 (2, 2)
 (2, 1)
 (1, 1)
 (1, 1)
 (2, 2)
 (2, 1)

У схемы разреженности есть несколько важных свойств:

Каждый элемент (i, j) означает структурный ненулевой элемент в строке i и столбце j.
Список может содержать повторяющиеся элементы; в этом случае следует сложить значения.
Список не обязательно должен быть отсортирован.
Список может содержать любое сочетание индексов нижнего и верхнего треугольников. Этот формат соответствует разреженной триплетной форме SparseArray в Julia, поэтому разреженное представление гессиана можно преобразовать в массив Julia типа SparseArray следующим образом.

I = [i for (i, _) in hessian_sparsity]
J = [j for (_, j) in hessian_sparsity]
V = zeros(length(hessian_sparsity))
n = num_variables(model)
H = SparseArrays.sparse(I, J, V, n, n)

2×2 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
 0.0   ⋅
 0.0  0.0

Местоположение нулей, конечно, не представляет особого интереса. На самом деле нам нужно вычислить значение гессиановой матрицы в некоторой точке.

MOI.eval_hessian_lagrangian(evaluator, V, ones(n), 1.0, ones(length(rows)))
H = SparseArrays.sparse(I, J, V, n, n)

2×2 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
  806.0     ⋅
 -398.0  202.0

На практике зачастую требуется вычислить значение гессиана в точке оптимального решения.

Сначала вычислим прямое решение. Для этого понадобится вектор переменных в том порядке, в котором они были переданы в решатель:

x = all_variables(model)

2-element Vector{VariableRef}:
 x[1]
 x[2]

Здесь x[1] — это переменная, соответствующая столбцу 1, и т. д. Вот оптимальное прямое решение:

x_optimal = value.(x)

2-element Vector{Float64}:
 0.7903587565231842
 0.6238546272155127

Далее необходимо получить оптимальное двойственное решение, связанное с нелинейными ограничениями (важно использовать порядок ограничений, в котором они были добавлены в nlp):

y_optimal = dual.(rows)

2-element Vector{Float64}:
 -8.038451738599348e-8
 -0.05744089305771262

Теперь можно вычислить гессиан в оптимальной прямо-двойственной точке:

MOI.eval_hessian_lagrangian(evaluator, V, x_optimal, 1.0, y_optimal)
H = SparseArrays.sparse(I, J, V, n, n)

2×2 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
  501.944     ⋅
 -316.258  199.885

Однако этот гессиан не совсем верен, так как он несимметричен. Это можно исправить, заполнив соответствующие недиагональные члены:

function fill_off_diagonal(H)
    ret = H + H'
    row_vals = SparseArrays.rowvals(ret)
    non_zeros = SparseArrays.nonzeros(ret)
    for col in 1:size(ret, 2)
        for i in SparseArrays.nzrange(ret, col)
            if col == row_vals[i]
                non_zeros[i] /= 2
            end
        end
    end
    return ret
end

fill_off_diagonal(H)

2×2 SparseMatrixCSC{Float64, Int64} with 4 stored entries:
  501.944  -316.258
 -316.258   199.885

Соберем все вместе:

function compute_optimal_hessian(model::Model)
    rows = Any[]
    nlp = MOI.Nonlinear.Model()
    for (F, S) in list_of_constraint_types(model)
        for ci in all_constraints(model, F, S)
            push!(rows, ci)
            object = constraint_object(ci)
            MOI.Nonlinear.add_constraint(nlp, object.func, object.set)
        end
    end
    MOI.Nonlinear.set_objective(nlp, objective_function(model))
    x = all_variables(model)
    backend = MOI.Nonlinear.SparseReverseMode()
    evaluator = MOI.Nonlinear.Evaluator(nlp, backend, index.(x))
    MOI.initialize(evaluator, [:Hess])
    hessian_sparsity = MOI.hessian_lagrangian_structure(evaluator)
    I = [i for (i, _) in hessian_sparsity]
    J = [j for (_, j) in hessian_sparsity]
    V = zeros(length(hessian_sparsity))
    MOI.eval_hessian_lagrangian(evaluator, V, value.(x), 1.0, dual.(rows))
    H = SparseArrays.sparse(I, J, V, length(x), length(x))
    return Matrix(fill_off_diagonal(H))
end

H_star = compute_optimal_hessian(model)

2×2 Matrix{Float64}:
  501.944  -316.258
 -316.258   199.885

Сравним полученное решение с аналитическим:

analytic_hessian(value.(x), 1.0, dual.([g_1, g_2]))

2×2 Matrix{Float64}:
  501.944  -316.258
 -316.258   199.885

Если посмотреть на собственные значения гессиана:

LinearAlgebra.eigvals(H_star)

2-element Vector{Float64}:
   0.4443995924983142
 701.3843426037456

то мы увидим, что все они положительные. Следовательно, гессиан является положительно определенным, и решение, найденное решателем Ipopt, является локально минимальным.

Якобианы

Помимо гессиана, можно также запрашивать другие компоненты нелинейной модели. Например, якобиан ограничений можно запросить с помощью MOI.jacobian_structure и MOI.eval_constraint_jacobian.

function compute_optimal_jacobian(model::Model)
    rows = Any[]
    nlp = MOI.Nonlinear.Model()
    for (F, S) in list_of_constraint_types(model)
        for ci in all_constraints(model, F, S)
            if !(F <: VariableRef)
                push!(rows, ci)
                object = constraint_object(ci)
                MOI.Nonlinear.add_constraint(nlp, object.func, object.set)
            end
        end
    end
    MOI.Nonlinear.set_objective(nlp, objective_function(model))
    x = all_variables(model)
    backend = MOI.Nonlinear.SparseReverseMode()
    evaluator = MOI.Nonlinear.Evaluator(nlp, backend, index.(x))
    # Инициализация якобиана
    MOI.initialize(evaluator, [:Jac])
    # Запрос структуры якобиана
    sparsity = MOI.jacobian_structure(evaluator)
    I, J, V = first.(sparsity), last.(sparsity), zeros(length(sparsity))
    # Запрос значений якобиана
    MOI.eval_constraint_jacobian(evaluator, V, value.(x))
    return SparseArrays.sparse(I, J, V, length(rows), length(x))
end

compute_optimal_jacobian(model)

2×2 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
 1.58072   ⋅
 2.82843  2.82843

Сравним с аналитическим решением:

y = value.(x)
[2y[1] 0; 2y[1]+2y[2] 2y[1]+2y[2]]

2×2 Matrix{Float64}:
 1.58072  0.0
 2.82843  2.82843