python

Метод оптимизации Trust-Region DOGLEG. Пример реализации на Python

  • суббота, 12 августа 2017 г. в 03:12:14
https://habrahabr.ru/post/335224/
  • Программирование
  • Машинное обучение
  • Математика
  • Алгоритмы
  • Python




Trust-region метод (TRM) является одним из самых важных численных методов оптимизации в решении проблем нелинейного программирования (nonlinear programming problems). Метод базируется на определении региона вокруг лучшего решения, в котором квадратичная модель аппроксимирует целевую функцию.

Методы линейного поиска (line search) и методы trust-region генерируют шаги с помощью аппроксимации целевой функции квадратичной моделью, но использую они эту модель по-разному. Линейный поиск использует её для получения направления поиска и дальнейшего нахождения оптимального шага вдоль направления. Trust-region метод определяет область (регион) вокруг текущей итерации, в котором модель достаточно аппроксимирует целевую функцию. В целях повышения эффективности направление и длина шага выбираются одновременно.

Trust-region методы надежны и устойчивы, могут быть применены к плохо обусловленным задачам и имеют очень хорошие свойства сходимости. Хорошая сходимость обусловлена тем, что размер области TR (обычно определяется модулем радиус-вектора) на каждой итерации зависит от улучшений сделанных на предыдущих итерациях.

Если вычисления показывают достаточно хорошее приближение аппроксимирующей модели к целевой функции, тогда trust-region может быть увеличен. В противном случае, если аппроксимирующая модель работает недостаточно хорошо, trust-region должен быть сокращен.

В итоге получаем приблизительно такую картину:


Алгоритм


Шаг №1


Trust-region метод использует квадратичную модель. На каждой итерации шаг вычисляется путем решения следующей квадратичной задачи (subproblem):

$\min_{p\in R^n}\ m_k(p) = f_k + p^T g_k + \frac 1 2 p^T B_kp, \ \\ s.t. |p|\leq\Delta_k,$


где $f_k = f(x_k), g_k = \nabla f(x_k), B_k = \nabla^2f(x_k)$;
$\Delta_k>0$ — trust-region радиус;

Таким образом trust-region требует последовательного вычисления аппроксимаций квадратичной модели в которых целевая функция и условие (которое может быть записано $p^Tp\le\Delta_k^2$) тоже квадратично. Когда гессиан ($B_k$) положительно определен и $|B_k^{-1}\nabla f_k|\le \Delta_k$ решение легко определить — это безусловный минимум $p_k^B = -B_k^{-1}\nabla f_k$. В данном случае $p_k^B$ называют полным шагом. Решение не так очевидно в других случаях, однако поиск его не стоит слишком дорого. Нам необходимо найти лишь приблизительное решение, чтобы получить достаточную сходимость и хорошее поведение алгоритма на практике.

Существует несколько стратегий аппроксимации квадратичной модели, среди которых следующие:

Алгоритм Cauchy point

Концепция метода схожа с логикой работы алгоритма наискорейшего спуска (steepest descent). Точка Cauchy лежит на градиенте, который минимизирует квадратичную модель при условии наличия шага в trust-region. Посредством последовательного нахождения точек Cauchy, может быть найден локальный минимум. Метод имеет неэффективную сходимость подобно методу наискорейшего спуска.

Алгоритм Steihaug

Метод носит название его исследователя — Steihaug. Представляет из себя модифицированный метод сопряженных градиентов (conjugate gradient approach). Для каждого шага алгоритм определяет итеративную технику, которая практически соответствует стандартной процедуре метода сопряженных градиентов за исключением двух дополнительных условий окончания алгоритма: если следующий шаг находится вне trust-region, если встречается нулевое или отрицательное направление.

Алгоритм Dogleg

В статье будет подробно рассмотрен метод аппроксимации квадратичной модели dogleg, который является одним из самых распространенных и эффективных методов. Используется в случае, если матрица Гессе (или ее приближение) положительна определена.
Наиболее простая и интересная версия метода работает с многоугольником состоящим всего из двух отрезков, которые некоторым людям напоминают ногу собаки.

Шаг №2


Первая проблема, которая возникает при определении trust-region алгоритма — это выбор стратегии для поиска оптимального trust-region радиуса $\Delta_k$ на каждой итерации. Выбор основывается на сходстве функции $m_k$ и целевой функции $f$ на предыдущей итерации. После нахождения $p_k$ мы определяем следующее соотношение:

$\rho_k = \frac{actual\ reduction}{predicted \ reduction} = \frac {f(x_k) - f(x_k + p_k)}{m_k(0) - m_k(p_k)}$


Знаменатель всегда должен быть положительно определен, поскольку шаг $p_k$ получается путем минимизации квадратичной модели $m_k$ по региону, который включает в себя шаг $p=0$. Отношение используется для определения преемственности шага и последующего обновления trust-region радиуса.

Если $\rho_k < 0$ или $\rho_k \approx0$ мы уменьшаем размер trust-region области.

Если $\rho_k\approx1$, тогда модель хорошо соответствует целевой функции. В данном случае следует расширить trust-region на следующей итерации.

В других случаяx trust-region остается неизменным.

Шаг №3


Следующий алгоритм описывает процесс:

Определяем стартовую точку $x_0$, максимальный trust-region радиус $\stackrel{-}{\Delta}$, начальный trust-region радиус $\Delta_0 = \in(0, \stackrel{-}{\Delta})$ и константу $\eta\in [0, \frac{1}{4})$

for k = 0, 1, 2,… пока $x_k$ не оптимум.

Решаем:

$\min_{p\in R^n}\ m_k(p) = f_k + p^T g_k + \frac 1 2 p^T B_kp \ \\ s.t. |p|\leq\Delta_k$


где $p_k$-решение.
Вычисляем отношение:

$\rho_k = \frac {f(x_k) - f(x_k + p_k)}{m_k(0) - m_k(p_k)}$


Обновляем текущую точку:

$$display$$x_{k+1} = \begin{equation*} \begin{cases} x_k + p_k & if\ \rho_k>\eta,\\ x_k&\text{otherwise.} \end{cases} \end{equation*}$$display$$


Обновляем trust-region радиус:

$$display$$\Delta_{k+1}=\begin{equation*} \begin{cases} \frac{1}{4}\Delta_k & if\ \rho_k<\frac{1}{4},\\ min(2\Delta_k, \stackrel{-}{\Delta})&if\ \rho_k>\frac{3}{4}\ and\ |p_k|=\Delta_k,\\ \Delta_k&otherwise. \end{cases} \end{equation*}$$display$$


Алгоритм в развернутом виде



Заметьте, что радиус увеличивается только если $|p_k|$ достигает границы trust-region. Если шаг остается строго в регионе, тогда текущее значение не влияет на работу алгоритма и нет необходимости изменять значение радиуса на следующей итерации.

Алгоритм Dogleg


Метод начинается с проверки эффективности trust-region радиуса в решении $p^*(\Delta)$ квадратичной модели $m(p)$. Когда $B$ положительна определена, как уже было отмечено, оптимальным решением будет полный шаг $p^B = -B^{-1}g$. Когда эта точка может быть найдена, очевидно, она будет являться решением.

$p^*(\Delta) = p^B, \ \Delta \geq|p^B|.$



Когда $\Delta_k$ весьма мало, ограничение $|p|\leq\Delta_k$ гарантирует, что квадратный член в модели $m(p)$ оказывает небольшое влияние на решение. Реальное решение $p(\Delta)$ аппроксимируется также, как если бы мы оптимизировали линейную функцию $f + g^Tp$ при условии $|p|\leq\Delta$, тогда:

$p^*(\Delta) \approx -\Delta\frac{g}{|g|}$


Когда $\Delta$ достаточно мало.

Для средних значений $\Delta$ решение $p^*(\Delta)$, как правило, следует за криволинейной траекторией, как показано на картинке:



Dogleg метод аппроксимирует криволинейную траекторию $p^*(\Delta)$ линией состоящей из двух прямыx. Первый отрезок начинается с начала и простирается вдоль направления наискорейшего спуска (steepest descent direction) и определяется следующим образом:

$p^U = -\frac{g^Tg}{g^TBg}g$


Второй начинается с $p^U$ и продолжается до $p^B$.

Формально мы обозначим траекторию $\tilde p(\tau)$, где $\tau \in [0, 2]$;

$$display$$\tilde p(\tau) = \begin{equation*} \begin{cases} \tau p^U, &0\leq\tau\leq1,\\\ p^U + (\tau - 1)(p^B-p^U), &1\leq\tau\leq2. \end{cases} \end{equation*}$$display$$


Для поиска $\tau$ необходимо решить квадратное уравнение, следующего вида:

$|p^U + \tau*(p^B - p^U)|^2 = \Delta^2$


$(p^U)^2 + 2\tau(p^B-p^U)p^U + \tau^2(p^B-p^U)^2=\Delta^2 $


Находим дискриминант уравнения:

$D = 4(p^B-p^U)^2(p^U)^2 - 4(p^B-p^U)^2((p^U)^2 - \Delta^2)$


$\sqrt{D} = 2(p^B-p^U)\Delta$


Корень уравнения равен:

$\tau = \frac{-2(p^B-p^U)p^U + 2(p^B - p^U)\Delta}{2(p^B-p^U)^2} = \frac{\Delta-p^U}{p^B-p^U}$


Dogleg метод выбирает $p_k$, чтобы минимизировать модель $m$ вдоль этого пути. В действительности нет необходимости выполнять поиск, поскольку dogleg путь пересекает границу trust-region только один раз и точка пересечения может быть найдена аналитически.

Пример


С использованием алгоритма trust-region (dogleg) оптимизировать следующую функцию (функция Розенброка):

$f(x, y) = (1-x)^2 + 100(y-x^2)^2$


Находим градиент, якобиан и гессиан функции:

$\frac{\partial f}{\partial x} = -400(y-x^2)x - 2 + 2x$

$\frac{\partial f}{\partial y} = 200y - 200x^2$

$\frac{\partial^2 f}{\partial x \partial x} = 1200x^2 - 400y + 2$

$\frac{\partial^2 f}{\partial x \partial y} = \frac{\partial^2 f}{\partial y \partial x} = -400x$

$\frac{\partial^2 f}{\partial y \partial y} = 200$

Инициализируем необходимые переменные для работы алгоритма:
$\Delta_k$ = 1.0,
$\stackrel{-}{\Delta}$ = 100.0,
$x_k = x_0 = [5, 5]$,
$gtol = 0.15$ (необходимая точность),
$\eta =0.15$.

Итерация 1

Находим оптимальное решение квадратичной модели $p_k$:

Поскольку $p^U = [-1.4226, 0.1422]$ и $|p^U| > \Delta_k$

Следовательно:

$p_k = \frac{\Delta_kp^U}{|p^U|} = [-1.4226, 0.1422]$



Вычисляем $\rho_k$:

actual reduction = $f(x_k) - f(x_k + p_k) = 28038.11$

predicted reduction = $m_k(0) - m_k(p_k) = 26146.06$

$\rho_k = \frac{28038.11}{26146.06} = 1.0723$



$\Delta_k$ — остается неизменным.

Обновляем $x_k$:

$x_k = x_k + p_k = [4.004, 5.099]$



Итерация 2

Находим оптимальное решение квадратичной модели $p_k$:

Поскольку $p^U = [-1.0109 0.1261]$ и $|p^U| > \Delta_k$.

Следовательно:

$p_k = \frac{\Delta_kp^U}{|p^U|} = [-0.9923 0.1238]$



Вычисляем $\rho_k$:

actual reduction = $f(x_k) - f(x_k + p_k) = 10489.43$

predicted reduction = $m_k(0) - m_k(p_k) = 8996.73$

$\rho_k = \frac{10489.43}{8996.73} = 1.1659$



Т.к. $\rho_k > 0.75$ и $|p_k|=\Delta_k$:

$\Delta_k = min(2\Delta_k, \stackrel{-}{\Delta_k}) = 2.0$



Обновляем $x_k$:

$x_k = x_k + p_k = [ 3.01, 5.22]$



Итерация 3

Находим оптимальное решение квадратичной модели $p_k$:

$p_k = p^U + \tau * (p^B-p^U) = [-0.257, 1.983]$



где $\tau = 0.5058$.

Вычисляем $\rho_k$:

actual reduction = $f(x_k) - f(x_k + p_k) = 1470.62$

predicted reduction = $m_k(0) - m_k(p_k) = 1424.16$

$\rho_k = \frac{1470.62}{1424.16} = 1.0326$



$\Delta_k$ — остается прежним.

Обновляем $x_k$:

$x_k = x_k + p_k = [2.7551, 7.2066]$



Алгоритм продолжается до тех пока $|g_k|<$gtol или не произведено заданное количество итераций.

Таблица результатов работы алгоритма для функции Розенброка:


    
k $p_k$ $\rho_k$ $\Delta_k$ $x_k$
0 - - 1.0 [5, 5]
1 [-0.9950, 0.0994] 1.072 1.0 [4.0049, 5.0994]
2 [-0.9923, 0.1238] 1.1659 2.0 [3.0126, 5.2233]
3 [-0.2575, 1.9833] 1.0326 2.0 [2.7551, 7.2066]
4 [-0.0225, 0.2597] 1.0026 2.0 [ 2.7325, 7.4663]
5 [-0.3605, -1.9672] -0.4587 0.5 [2.7325, 7.4663]
6 [-0.0906, -0.4917] 0.9966 1.0 [ 2.6419, 6.9746]
7 [-0.1873, -0.9822] 0.8715 2.0 [ 2.4546, 5.9923]
8 [-0.1925, -0.9126] 1.2722 2.0 [ 2.2620, 5.0796]
9 [-0.1499, -0.6411] 1.3556 2.0 [ 2.1121, 4.4385]
10 [-0.2023, -0.8323] 1.0594 2.0 [ 1.9097, 3.6061]
11 [-0.0989, -0.3370] 1.2740 2.0 [ 1.8107, 3.2690]
12 [-0.2739, -0.9823] -0.7963 0.25495 [ 1.8107, 3.2690]
13 [-0.0707, -0.2449] 1.0811 0.5099 [ 1.7399, 3.0240]
14 [-0.1421, -0.4897] 0.8795 1.0198 [1.5978, 2.5343]
15 [-0.1254, -0.3821] 1.3122 1.0198 [ 1.4724, 2.1522]
16 [-0.1138, -0.3196] 1.3055 1.0198 [ 1.3585, 1.8326]
17 [-0.0997, -0.2580] 1.3025 1.0198 [ 1.2587, 1.5745]
18 [-0.0865, -0.2079] 1.2878 1.0198 [ 1.1722, 1.3666]
19 [-0.0689, -0.1541] 1.2780 1.0198 [ 1.1032, 1.2124]
20 [-0.0529, -0.1120] 1.2432 1.0198 [ 1.0503, 1.1004]
21 [-0.0322, -0.0649] 1.1971 1.0198 [1.0180, 1.0354]
22 [-0.0149, -0.0294] 1.1097 1.0198 [ 1.0031, 1.0060]
23 [-0.0001, -0.0002] 1.0012 1.0198 [ 1.00000024, 1.00000046]
24 [-2.37065e-07, -4.56344e-07] 1.0000 1.0198 [ 1.0, 1.0]


Аналитически находим минимум функции Розенброка, он достигается в точке $[1, 1]$. Таким образом можно убедиться в эффективности алгоритма.

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


Алгоритм реализован с использованием библиотеки numpy. В примере наложено ограничение на количество итераций.

#!/usr/bin/python
# -*- coding: utf-8 -*-
import numpy as np
import numpy.linalg as ln
from math import sqrt

def f(x):
    return (1-x[0])**2 + 100*(x[1]-x[0]**2)**2

def jac(x):
    return np.array([-400*(x[1] - x[0]**2)*x[0] - 2 + 2*x[0], 200*x[1] - 200*x[0]**2])

def hess(x):
    return np.array([[1200*x[0]**2 - 400*x[1]+2, -400*x[0]], [-400*x[0], 200]])


def dogleg_update(Hk, gk, Bk, trust_radius):

    # Calculate the full step and its norm.
    pB = -np.dot(Hk, gk)
    norm_pB = sqrt(np.dot(pB, pB))

    # Test if the full step is within the trust region.
    if norm_pB <= trust_radius:
        return pB

    # Calculate pU.
    pU = - (np.dot(gk, gk) / np.dot(gk, np.dot(Bk, gk))) * gk
    dot_pU = np.dot(pU, pU)
    norm_pU = sqrt(dot_pU)

    # Test if the step pU exits the trust region.
    if norm_pU >= trust_radius:
        return trust_radius * pU / norm_pU

    # Find the solution to the scalar quadratic equation.
    pB_pU = pB - pU
    dot_pB_pU = np.dot(pB_pU, pB_pU)
    dot_pU_pB_pU = np.dot(pU, pB_pU)
    fact = dot_pU_pB_pU**2 - dot_pB_pU * (dot_pU - trust_radius**2)
    tau = (-dot_pU_pB_pU + sqrt(fact)) / dot_pB_pU
    
    # Decide on which part of the trajectory to take.
    return pU + tau * pB_pU
    

def trust_region_dogleg(func, jac, hess, x0, initial_trust_radius=1.0,
                        max_trust_radius=100.0, eta=0.15, gtol=1e-4,
                        maxiter=400):
    xk = x0
    trust_radius = initial_trust_radius
    k = 0
    while True:
        
        # Initialize the search
        gk = jac(xk)
        Bk = hess(xk)
        Hk = np.linalg.inv(Bk)
        
        # Solve the sub-problem.
        pk = dogleg_update(Hk, gk, Bk, trust_radius)
       
        act_red = func(xk) - func(xk + pk)

        pred_red = -(np.dot(gk, pk) + 0.5 * np.dot(pk, np.dot(Bk, pk)))
        
        # Evaluate the ratio
        rhok = act_red / pred_red
        if pred_red == 0.0:
            rhok = 1e99
        else:
            rhok = act_red / pred_red
        norm_pk = sqrt(np.dot(pk, pk))
        
        # Update the trust radius according to the actual/predicted ratio
        if rhok < 0.25:
            trust_radius = 0.25 * norm_pk
        else: 
            if rhok > 0.75 and norm_pk == trust_radius:
                trust_radius = min(2.0*trust_radius, max_trust_radius)
            else:
                trust_radius = trust_radius
        
        if rhok > eta:
            xk = xk + pk
        else:
            xk = xk
            
        # Check if the gradient is small enough to stop
        if ln.norm(gk) < gtol:
            break
        
        # Check if we have looked at enough iterations
        if k >= maxiter:
            break
        k = k + 1
    return xk
    
result = trust_region_dogleg(f, jac, hess, [5, 5])
print("Result of trust region dogleg method:")
print(result)
print("Value of function at a point:")
print(f(result))


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

О всех неточностях сообщайте, пожалуйста, в личных сообщениях.