Метод оптимизации Trust-Region DOGLEG. Пример реализации на 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).

    Алгоритм 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)))
    


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

    Подробнее
    Реклама
    Комментарии 0

    Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.