How to become an author
Search
Write a publication
Pull to refresh

Метод BFGS или один из самых эффективных методов оптимизации. Пример реализации на Python

Reading time4 min
Views51K
Python*Algorithms*Mathematics*Machine learning*Programming*


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

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

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


Алгоритм



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

Шаг №1
Инициализируем начальную точку $x_0$;
Задаем точность поиска $ε$ > 0;
Определяем начальное приближение $H_0 = B_0^{-1}$, где $B_0^{-1}$ — обратный гессиан функции;

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

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

$p_k = -H_k* \nabla f_k$



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

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


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

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


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


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

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

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

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


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


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

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

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


где $ρ_k$

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


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

Замечание


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

$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}$



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

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





Пример


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

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

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


Итерация №0:

$x_0 = (1, 1)$


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

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


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

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


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


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

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

$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}$


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

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

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

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


Тогда

$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$

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

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


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

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


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


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

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


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



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

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


Итерация 1:

$x_1 = (-2.571, 2.786)$


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


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

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



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

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


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


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

$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$

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


Тогда:

$x_2 = (-3.571, 0.786)$


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


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

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


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


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


Итерация 2:

$x_2 = (-4, 1)$


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

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


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

Еще о методе


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

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

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

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


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

'''
    Pure Python/Numpy implementation of the BFGS algorithm.
    Reference: https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm
'''

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


# Objective function
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):
    """
    Minimize a function func using the BFGS algorithm.
    
    Parameters
    ----------
    func : f(x)
        Function to minimise.
    x0 : ndarray
        Initial guess.
    fprime : fprime(x)
        The gradient of `func`.
    """
    
    if maxiter is None:
        maxiter = len(x0) * 200

    # initial values
    k = 0
    gfk = fprime(x0)
    N = len(x0)
    # Set the Identity matrix I.
    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)
        
        # Line search constants for the Wolfe conditions.
        # Repeating the line search
        
        # 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)


result, k = bfgs_method(f, f1, np.array([1, 1]))

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


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

С вами был FUNNYDMAN. Всем пока!
Tags:
Hubs:
Total votes 27: ↑24 and ↓3+21
Comments4
22
Karma
0
Rating
@FUNNYDMAN

software engineer

Send message

Articles

Your account

Sections

Information

Services

Support
© 2006–2024, Habr