Метод BFGS, итерационный метод численной оптимизации, назван в честь его исследователей: Broyden, Fletcher, Goldfarb, Shanno. Относится к классу так называемых квазиньютоновских методов. В отличие от ньютоновских методов в квазиньютоновских не вычисляется напрямую гессиан функции, т.е. нет необходимости находить частные производные второго порядка. Вместо этого гессиан вычисляется приближенно, исходя из сделанных до этого шагов.

Существует несколько модификаций метода:
L-BFGS (ограниченное использование памяти) — используется в случае большого количества неизвестных.
L-BFGS-B — модификация с ограниченным использованием памяти в многомерном кубе.

Метод эффективен и устойчив, поэтому зачастую применяется в функциях оптимизации. Например в SciPy, популярной библиотеки для языка python, в функции optimize по умолчанию применяется BFGS, L-BFGS-B.

Алгоритм

Пусть задана некоторая функция $inline$f(x, y)$inline$ и мы решаем задачу оптимизации: $inline$min f(x, y)$inline$.
Где в общем случае $inline$f(x, y)$inline$ является не выпуклой функцией, которая имеет непрерывные вторые производные.

Шаг №1
Инициализируем начальную точку $inline$x_0$inline$;
Задаем точность поиска $inline$ε$inline$ > 0;
Определяем начальное приближение гессиана функции $inline$H_0$inline$;

Каким нужно выбрать начальное приближение $inline$H_0$inline$?
К сожалению не существует общей формулы, которая хорошо бы работала во всех случаях. В качестве начального приближения можно взять гессиан функции, вычисленный в начальной точке $inline$x_0$inline$. Иначе можно использовать хорошо обусловленную, невырожденную матрицу, на практике часто берут единичную матрицу.

Шаг №2
Находим точку, в направлении которой будем производить поиск, она определяется следующим образом:

$$display$$p_k = -H_k* nabla f_k$$display$$

Шаг №3
Вычисляем $inline$x_{k+1}$inline$ через рекуррентное соотношение:

$$display$$x_{k+1} = x_k + α_k * p_k$$display$$

Коэффициент $inline$α_k$inline$ находим используя линейный поиск (linear search), где $inline$α_k$inline$ удовлетворяет условиям Вольфе (Wolfe conditions ^[1]):

$$display$$f(x_k + α_k * p_k) leq f(x_k) + c_1 * α_k * nabla f_k^T *p_k$$display$$

$$display$$nabla f(x_k + α_k * p_k)^T * p_k geq c_2 * nabla f_k^T * p_k$$display$$

Константы $inline$с_1$inline$ и $inline$с_2$inline$ выбирают следующим образом: $inline$0 leq c_1 leq c_2 leq 1$inline$. В большинстве реализаций: $inline$c_1 = 0.0001$inline$ и $inline$с_2 = 0.9$inline$.

Фактически мы находим такое $inline$α_k$inline$ при котором значение функции $inline$f(x_k + α_k * p_k)$inline$ минимально.

Шаг №4
Определяем вектора:

$$display$$s_k = x_{k+1} - x_k$$display$$

$$display$$y_k = nabla f_{k+1} - nabla f_k$$display$$

$inline$s_k$inline$ — шаг алгоритма на итерации, $inline$y_k$inline$ — изменение градиента на итерации.

Шаг №5
Обновляем гессиан функции, согласно следующей формуле:

$$display$$H_{k+1} = (I - ρ_k * s_k * y_k^T)H_k(I - ρ_k * y_k * s_k^T) + ρ * s_k * s_k^T$$display$$

где $inline$ρ_k$inline$

$$display$$ρ_k = frac {1}{y_k^T s_k}$$display$$

$inline$I$inline$ — единичная матрица.

Замечание

Выражение вида $inline$y_k * s_k^T$inline$ является внешним произведением (outer product ^[2]) двух векторов.
Пусть определены два вектора $inline$U$inline$ и $inline$V$inline$, тогда их внешнее произведение эквивалентно матричному произведению $inline$UV^T$inline$. Например, для векторов на плоскости:

$$display$$UV^T = begin{pmatrix} U_1 \ U_2 end{pmatrix}begin{pmatrix} V_1 & V_2 end{pmatrix} = begin{pmatrix} U_1V_1 & U_1V_2 \ U_2V_1 & U_2V_2 end{pmatrix}$$display$$

Шаг №6
Алгоритм продолжает выполнятся до тех пор пока истинно неравенство: $inline$|nabla f_k| > ε$inline$.

Алгоритм схематически

Пример

Найти экстремум следующей функции: $inline$f(x, y) = x^2 - xy + y^2 + 9x - 6y + 20$inline$
В качестве начальной точки возьмем $inline$x_0 = (1, 1)$inline$;
Необходимая точность $inline$ε = 0.001$inline$;

Находим градиент функции $inline$f(x, y)$inline$:

$$display$$nabla f = begin{pmatrix} 2x - y + 9 \ -x + 2y - 6 end{pmatrix} $$display$$

Итерация №0:

$$display$$x_0 = (1, 1)$$display$$

Находим градиент в точке $inline$x_0$inline$:

$$display$$nabla f(x_0) = begin{pmatrix} 10 \ -5 end{pmatrix}$$display$$

Проверка на окончание поиска:

$$display$$|nabla f(x_0)| = sqrt {10^2 + (-5)^2} = 11.18$$display$$

$$display$$|nabla f(x_0)| = 11.18 > 0.001$$display$$

Неравенство не выполняется, поэтому продолжаем процесс итераций.

Находим точку в направлении которой будем производить поиск:

$$display$$p_0 = -H_0 * nabla f(x_0) = -begin{pmatrix} 1 & 0 \ 0 & 1 end{pmatrix} begin{pmatrix} 10 \ -5 end{pmatrix} = begin{pmatrix} -10 \ 5 end{pmatrix}$$display$$

где $inline$H_0 = I $inline$ — единичная матрица.

Ищем такое $inline$α_0$inline$$inline$geq$inline$0, чтобы $inline$f(x_0 + α_0*p_0) = min(f(x_0 + α_0*p_0))$inline$

Аналитически задачу можно свести к нахождению экстремума функции от одной переменной:

$$display$$x_0 + α_0 * p_0 = (1, 1) + α_0(-10, 5) = (1 - 10α_0, 1 + 5α_0)$$display$$

Тогда

$inline$f(α_0) = (1 - 10α_0)^2 - (1 - 10α_0)(1 + 5α_0) + (1 + 5α_0)^2 + 9(1 - 10α_0) - 6(1 + 5α_0) + 20$inline$

Упростив выражение, находим производную:

$$display$$frac{partial f}{partial x} = 350α_0 - 125 = 0 Rightarrow α_0 = 0.357$$display$$

Находим следующую точку $inline$x_1$inline$:

$$display$$x_1 = x_0 + α_0*p_0 = (-2.571, 2.786)$$display$$

$$display$$s_0 = x_1 - x_0 = (-2.571, 2.786) - (1, 1) = (-3.571, 1.786)$$display$$

Вычисляем значение градиента в т. $inline$x_1$inline$:

$$display$$nabla f(x_1) = begin{pmatrix} 1.071 \ 2.143 end{pmatrix}$$display$$

$$display$$y_0 = nabla f(x_1) - nabla f(x_0) = (1.071, 2.143) - (10, -5) = (-8.929, 7.143)$$display$$

Находим приближение гессиана:

$$display$$H_1 = begin{pmatrix} 0.694 & 0.367 \ 0.367 & 0.709 end{pmatrix}$$display$$

Итерация 1:

$$display$$x_1 = (-2.571, 2.786)$$display$$

$$display$$nabla f(x_1) = begin{pmatrix} 1.071 \ 2.143 end{pmatrix}$$display$$

Проверка на окончание поиска:

$$display$$|nabla f(x_1)| = 2.396 > 0.001$$display$$

Неравенство не выполняется, поэтому продолжаем процесс итераций:

$$display$$H_1 = begin{pmatrix} 0.694 & 0.367 \ 0.367 & 0.709 end{pmatrix}$$display$$

$$display$$p_1 = -H_1nabla f(x_1) = -begin{pmatrix} 0.694 & 0.367 \ 0.367 & 0.709end{pmatrix}begin{pmatrix} 1.071 \ 2.143 end{pmatrix} = -begin{pmatrix} 1.531 \ 1.913end{pmatrix}$$display$$

Находим $inline$α_1$inline$ > 0:

$inline$f(α_1) = -2.296α_1 + (-1.913α_1 + 2.786)^2 - (-1.913α_1 + 2.786)(-1.531α_1 - 2.571) + (-1.531α_1 - 2.571)^2 - 19.86$inline$

$$display$$frac{partial f}{partial α_1} = 6.15α_1 - 5.74 = 0 Rightarrow α_1 = 0.933 $$display$$

Тогда:

$$display$$x_2 = (-3.571, 0.786)$$display$$

$$display$$s_1 = x_2 - x_1 = (-4, 1) - (-2.571, 2.786) = (-1.429, -1.786)$$display$$

Вычисляем значение градиента в точке $inline$x_2$inline$:

$$display$$nabla f(x_2) = begin{pmatrix} 0 \ 0 end{pmatrix}$$display$$

$$display$$g_1 = nabla f(x_2) - nabla f (x_1) = (0, 0) - (1.071, 2.143) = (-1.071, -2.143)$$display$$

$$display$$H_2 = begin{pmatrix} 0.667 & 0.333 \ 0.333 & 0.667 end{pmatrix}$$display$$

Итерация 2:

$$display$$x_2 = (-4, 1)$$display$$

Проверка на окончание поиска:

$$display$$|nabla f(x_2)| = 0 < 0.001$$display$$

Неравенство выполняется следовательно поиска заканчивается. Аналитически находим экстремум функции он достигается в точке $inline$x_2$inline$.

Еще о методе

Каждая итерация может быть совершена со стоимостью $inline$O(n^2)$inline$ ( плюс стоимость вычисления функции и оценки градиента). Здесь нет $inline$O(n^3)$inline$ операций таких, как решение линейных систем или сложных математических операций. Алгоритм устойчив и имеет сверхлинейную сходимость, чего достаточно для большинства практических задач. Даже если методы Ньютона сходятся гораздо быстрее (квадратично), стоимость каждой итерации выше, поскольку необходимо решать линейные системы. Неоспоримое преимущество алгоритма, конечно, состоит в том, что нет необходимости вычислять вторые производные.

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

Формула BFGS имеет самокорректирующиеся свойства. Если матрица $inline$H$inline$ не верно оценивает кривизну функции и если эта плохая оценка замедляет алгоритм, тогда апроксимация гессиана стремится исправить ситуацию за несколько шагов. Самокорректирующие свойства алгоритма работают только в том случае, если реализован соответствующий линейный поиск (соблюдены условия Вольфе).

Пример реализации на Python

Алгоритмы реализован с использованием библиотек numpy и scipy. Линейный поиск был реализован посредством использования функция line_search() из модуля scipy.optimize. В примере наложено ограничение на количество итераций.

#!/usr/bin/python
# -*- coding: utf-8 -*-
import numpy as np
import numpy.linalg as ln
import scipy as sp
import scipy.optimize


# Objective function
# f(x, y) = x^2 - xy + y^2 + 9x - 6y + 20
def f(x):
    return x[0]**2 - x[0]*x[1] + x[1]**2 + 9*x[0] - 6*x[1] + 20


# Derivative

def f1(x):
    return np.array([2 * x[0] - x[1] + 9, -x[0] + 2*x[1] - 6])


def bfgs_method(f, fprime, x0, maxiter=None, epsi=10e-3):
    
    if maxiter is None:
        maxiter = len(x0) * 50

    # initial values
    k = 0
    gfk = fprime(x0)
    N = len(x0)
    I = np.eye(N, dtype=int)
    Hk = I
    xk = x0
   
    while ln.norm(gfk) > epsi and k < maxiter:
        
        # pk - direction of search
        
        pk = -np.dot(Hk, gfk)

        # Repeating the linesearch
        # line_search returns not only alpha
        # but only this value is interesting for us

        line_search = sp.optimize.line_search(f, f1, xk, pk)
        alpha_k = line_search[0]
        
        xkp1 = xk + alpha_k * pk
        sk = xkp1 - xk
        xk = xkp1
        
        gfkp1 = fprime(xkp1)
        yk = gfkp1 - gfk
        gfk = gfkp1
        
        k += 1
        
        ro = 1.0 / (np.dot(yk, sk))
        A1 = I - ro * sk[:, np.newaxis] * yk[np.newaxis, :]
        A2 = I - ro * yk[:, np.newaxis] * sk[np.newaxis, :]
        Hk = np.dot(A1, np.dot(Hk, A2)) + (ro * sk[:, np.newaxis] *
                                                 sk[np.newaxis, :])
                                                 
    return (xk, k)


x0 = np.array([1, 1])

result, k = bfgs_method(f, f1, x0)

print('Result of BFGS method:')
print('Final Result (best point): %s' % (result))
print('Iteration Count: %s' % (k))

Спасибо за интерес проявленный к моей статье. Надеюсь она была Вам полезна и Вы узнали много нового.

Лучший способ отблагодарить автора — это оценить статью. Таким образом Вы даете понять, что мной не зря проделана работа. Если заметили неточности, пожалуйста, сообщите.

С вами был FYNNYDMAN. Всем пока!

Автор: FUNNYDMAN

Источник ^[3]