10.4. Метод BFGS

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

Квазиньютоновские методы оптимизации также используют ньютоновское направление убывания. Однако, в отличие от метода Ньютона, эти методы не вычисляют гессиан на итерации, а обновляют его по некоторому правилу (например, используя матрицу Гессе с предыдущей итерации). Методы отличаются между собой формулой обновления, ниже мы рассмотрим широко-используемый метод Бройдена – Флетчера – Гольдфарба – Шанно (Broyden – Fletcher – Goldfarb – Shanno) или, более кратко, метод BFGS.

10.4.1. Метод

Обозначим приближение гессиана на \(k\)-ой итерации через \(\mathbf{B}_k\). После того, как направление вычислено как направление Ньютона

(10.8)\[\mathbf{d}_k = -\mathbf{B}_k^{-1} \nabla f(\mathbf{x}_k),\]

находится новое приближение решения

\[\mathbf{x}_{k+1} = \mathbf{x}_k + \alpha \mathbf{d}_k,\]

где \(\alpha\) подобрано линейным поиском.

Поскольку \(\mathbf{x}_{k+1}\) известно, можно учесть кривизну функции по векторам

\[\mathbf{s}_k = \mathbf{x}_{k+1} - \mathbf{x}_k ,\quad \mathbf{y}_k = \nabla f(\mathbf{x}_{k+1}) - \nabla f(\mathbf{x}_k),\]

используя их для нового приближения гессиана \(\mathbf{B}_{k+1}\).

В методе BFGS используется следующее правило обновления

Определение 10.11 (Правило BFGS)

(10.9)\[\mathbf{B}_{k+1} = \mathbf{B}_k + \frac{\mathbf{y}_k\mathbf{y}_k^\top}{\mathbf{y}_k^\top\mathbf{s}_k} - \frac{ \mathbf{B}_k \mathbf{s}_k \mathbf{s}_k^\top \mathbf{B}_k^\top }{ \mathbf{s}_k^\top\mathbf{B}_k\mathbf{s}_k } .\]

Однако, \(\mathbf{B}_k\) требует решение линейной системы (10.8). Чтобы этого не делать, можно аналитически обратить (10.9) и работать с обратным приближением гессиана \(\mathbf{B}^{-1}_k\), тогда (10.8) уже будет представлен не систему уравнений, а матрично-векторное умножение.

Определение 10.12 (Правило BFGS для обратного гессиана)

\[\mathbf{B}^{-1}_{k+1} = \mathbf{B}^{-1}_k + \frac{ (\mathbf{s}_k^\top \mathbf{y}_k + \mathbf{y}_k^\top \mathbf{B}^{-1}_k \mathbf{y}_k) }{ (\mathbf{s}_k^\top \mathbf{y}_k)^2 } (\mathbf{s}_k \mathbf{s}_k^\top) - \frac{ \mathbf{B}^{-1}_k \mathbf{y}_k \mathbf{s}_k^\top + \mathbf{s}_k \mathbf{y}_k^\top \mathbf{B}^{-1}_k }{ \mathbf{s}^\top_k \mathbf{y}_k } .\]

В качестве начального приближения \(\mathbf{B}_0\) используют обычно истинный гессиан или его модифицированное разложение Холецкого. Также используют диагональные матрицы, например, единичную или построенную по разностной схеме из градиента в начальной точке.

Для начального приближения обратного гессиана можно также использовать диагональные матрицы. Если требуется истинный обратный гессиан, можно перед запуском алгоритма высчитать его один раз. Кроме того, для устойчивости алгоритма, можно сделать и модифицированное разложение Холецкого для \(\mathbf{B}^{-1}_0\).

10.4.2. Реализация

Модифицированное разложение Холецкого приведено в Приложении.

Ниже дан шаблон метода BFGS с обратным гессианом, который нужно будет завершить в домашнем задании.

Структура данных для результата алгоритма

struct BFGSResult{T<:Real}
    converged::Bool      # метод сошёлся
    argument::Vector{T}  # найденный минимум
    iters::Int           # число прошедших итераций
end

Функция 10.13 (bfgs)

Метод BFGS (шаблон)

"""
    bfgs(f, ∇f, x0[, invH0; maxiter=200, gtol=1e-5])

Поиск минимума методом BFGS функции `f` с градиентом `∇f`, начиная с `x0`.
Оптимизация длится не более `maxiter` итераций, при этом, если норма градиента не превышает `gtol`, завершается досрочно.

- `f::Function`: по вектору x возвращает значение функции, `f(x)`: `Vector` → `Real`;
- `∇f::Function`: по вектору x возвращает градиент `f`, `∇f(x)`: `Vector` → `Vector`;
- `invH0::Matrix`: приближение обратного гессиана в `x0`, по умолчанию -- единичная матрица.
"""
function bfgs(f, ∇f, x0, invH0=diagm(0=>ones(length(x)));
    maxiter=200,
    gtol=1e-5,
)
    x = float.(x0)
    
    # проверяем сходимость по норме-2 градиента
    g = ∇f(x)
    # ... && BFGSResult(true, x, 0)

    # обратный квази-гессиан
    # invB = ...

    for i in 1:maxiter
        # выбор направления
        # d = ...
        
        # подбор шага вдоль d
        # α = ...
        # ... && return BFGSResult(false, x, i)  # α не найдено
        
        # совершение шага
        # x = ...
        
        # проверка сходимости
        g = ∇f(x)
        # ... && return BFGSResult(true, x, i)

        # обновление обратного квази-гессиана для следующей итерации
        # invB = ...
    end
    return BFGSResult(false, x, maxiter)
end