10.3. Метод Ньютона

Приблизим локально значение функции \(f\) в точке \(\mathbf{x}_k + \mathbf{d}\) квадратичной функцией \(m_k(\mathbf{d})\), используя разложение Тейлора

\[f(\mathbf{x}_k + \mathbf{d}) \approx f(\mathbf{x}_k) + \nabla f(\mathbf{x}_k)^\top \mathbf{d} + \frac{1}{2} \mathbf{d}^\top \nabla^2 f(\mathbf{x}_k) \mathbf{d} = m_k(\mathbf{d}).\]

Предположим, что гессиан \(\nabla^2 f(\mathbf{x}_k)\) положительно определён, тогда в качестве направления убывания возьмём вектор \(\mathbf{p}\), который минимизирует модельную функцию \(m_k(\mathbf{d})\). Для этого найдём ноль градиента функции \(m_k\)

\[\mathbf{0} = \nabla m_k(\mathbf{p}) = \nabla f(\mathbf{x}_k) + \nabla^2 f(\mathbf{x}_k) \mathbf{d},\]

откуда

(10.7)\[\mathbf{d} = - (\nabla^2 f(\mathbf{x}_k))^{-1} \nabla f(\mathbf{x}_k).\]

Данный вектор будем использовать в качестве направления убывания на \(k\)-ом шаге оптимизации. Направление вида (10.7) называют Ньютоновским. В сочетании с линейным поиском получаем метод Ньютона для задачи оптимизации. Однако, прежде сделаем модификацию.

10.3.1. Модификация метода

Отметим предположения, при которых (10.7) справедливо. Во-первых, размер направления \(\|\mathbf{d}\|\) не должен быть слишком большим, во-вторых, мы предполагали положительную определённость гессиана.

На практике размер шага \(\|\mathbf{d}\|\) подбирается линейным поиском \(\mathbf{x}_k + \alpha \mathbf{d}\), начиная с \(\alpha = 1\). Вдали от минимума подходящие шаги могу оказаться \(\alpha < 1\), однако, вблизи минимума полный шаг чаще подходит.

Направление \(\mathbf{d}\) (10.7) может быть не убывающим вдали от минимума, поскольку гессиан может быть не положительно определён. Чтобы гарантировать положительную определённость, используют модифицированное разложение Холецкого. Мы уже упоминали о Разложении Холецкого. Его же модифицированная версия принимает на вход произвольную матрицу \(\mathbf{A}\) и, пытаясь выполнить стандартное разложение, корректирует элементы \(\mathbf{A}\) так, чтобы итоговая матрица \(\mathbf{A}^{(m)}\) была положительно определена и не слишком сильно отличалась от исходной:

\[\mathbf{A} \approx \mathbf{A}^{(m)} = \mathbf{L}\mathbf{L}^\top.\]

Наша модификация алгоритма Ньютона заключается в использовании на \(k\)-ом шаге оптимизации не исходного гессиана \(\mathbf{H}\), а его приближения \(\mathbf{H}^{(m)}\), полученного из модифицированного разложения Холецкого. Вдали от минимума гессиан и приближение могут отличаться, однако вблизи минимума разложение практически не внесёт поправок.

В принципе, мы могли бы проверять на итерации, является ли гессиан положительно определённым и, если нет, то использовать \(\mathbf{H}^{(m)}\). Однако, процедура проверки строится на попытке выполнить стандартное разложение Холецкого, поэтому большой выгоды от этого нет.

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

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

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

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

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

Функция 10.10 (newton)

Метод Ньютона (шаблон)

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

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

- `f::Function`: по вектору x возвращает значение функции, `f(x)`: `Vector` → `Real`;
- `∇f::Function`: по вектору x возвращает градиент `f`, `∇f(x)`: `Vector` → `Vector`;
- `hess::Function`: по вектору x возвращает гессиан `f`, `hess(x)`: `Vector` → `Matrix`.

Возвращает `NewtonResult` - результат оптимизации.
"""
function newtontemplate(
    f::Function,
    ∇f::Function,
    hess::Function,
    x0::AbstractVector;
    gtol::Real=1e-5,
    maxiter::Integer=200,
)
    x = float.(x0)  # копирование и приведение к флоат-числам

    # проверяем сходимость по норме-2 градиента
    g = ∇f(x)
    # ... && return NewtonResult(true, x, 0)

    for i in 1:maxiter
        # выбор направления d из модифицированного Гессиана Hᵐ
        # см. mcholesky!, LinearAlgebra.ldiv!
        # ...
        # d = ...

        # выбор шага вдоль d
        # α = ...
        # ... && return NewtonResult(false, x, i)  # α не найдено
        
        # совершаем шаг
        # x = ...

        # проверяем сходимость
        g = ∇f(x)
        # ... && return NewtonResult(true, x, i)
    end
    return NewtonResult(false, x, maxiter)
end