11.08.2017 16:03
FUNNYDMAN 0 2000 Источник


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).

Алгоритм 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. В примере наложено ограничение на количество итераций.

'''
    Pure Python/Numpy implementation of the Trust-Region Dogleg algorithm.
    Reference: https://optimization.mccormick.northwestern.edu/index.php/Trust-region_methods
'''

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

# Objective function    
def f(x):
    return (1-x[0])**2 + 100*(x[1]-x[0]**2)**2
    
# Gradient
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])
# Hessian
def hess(x):
    return np.array([[1200*x[0]**2 - 400*x[1]+2, -400*x[0]], [-400*x[0], 200]])

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

    # Compute the Newton point.
    # This is the optimum for the quadratic model function.
    # If it is inside the trust radius then return this point.
    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
    
    # Compute the Cauchy point.
    # This is the predicted optimum along the direction of steepest descent.
    pU = - (np.dot(gk, gk) / np.dot(gk, np.dot(Bk, gk))) * gk
    dot_pU = np.dot(pU, pU)
    norm_pU = sqrt(dot_pU)

    # If the Cauchy point is outside the trust region,
    # then return the point where the path intersects the boundary.
    if norm_pU >= trust_radius:
        return trust_radius * pU / norm_pU

    # Find the solution to the scalar quadratic equation.
    # Compute the intersection of the trust region boundary
    # and the line segment connecting the Cauchy and Newton points.
    # This requires solving a quadratic equation.
    # ||p_u + tau*(p_b - p_u)||**2 == trust_radius**2
    # Solve this for positive time t using the quadratic formula.
    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=100):
    xk = x0
    trust_radius = initial_trust_radius
    k = 0
    while True:
      
        gk = jac(xk)
        Bk = hess(xk)
        Hk = np.linalg.inv(Bk)
        
        pk = dogleg_method(Hk, gk, Bk, trust_radius)
       
        # Actual reduction.
        act_red = func(xk) - func(xk + pk)
        
        # Predicted reduction.
        pred_red = -(np.dot(gk, pk) + 0.5 * np.dot(pk, np.dot(Bk, pk)))
        
        # Rho.
        rhok = act_red / pred_red
        if pred_red == 0.0:
            rhok = 1e99
        else:
            rhok = act_red / pred_red
            
        # Calculate the Euclidean norm of pk.
        norm_pk = sqrt(np.dot(pk, pk))
        
        # Rho is close to zero or negative, therefore the trust region is shrunk.
        if rhok < 0.25:
            trust_radius = 0.25 * norm_pk
        else: 
        # Rho is close to one and pk has reached the boundary of the trust region, therefore the trust region is expanded.
            if rhok > 0.75 and norm_pk == trust_radius:
                trust_radius = min(2.0*trust_radius, max_trust_radius)
            else:
                trust_radius = trust_radius
        
        # Choose the position for the next iteration.
        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: {}".format(result))
print("Value of function at a point: {}".format(f(result)))


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

Комментарии (0)