Метод оптимизации Нелдера — Мида. Пример реализации на Python



    Метод Нелдера — Мида — метод оптимизации (поиска минимума) функции от нескольких переменных. Простой и в тоже время эффективный метод, позволяющий оптимизировать функции без использования градиентов. Метод надежен и, как правило, показывает хорошие результаты, хотя и отсутствует теория сходимости. Может использоваться в функции optimize из модуля scipy.optimize популярной библиотеки для языка python, которая используется для математических расчетов.

    Алгоритм заключается в формировании симплекса (simplex) и последующего его деформирования в направлении минимума, посредством трех операций:

    1) Отражение (reflection);
    2) Растяжение (expansion);
    3) Сжатие (contract);

    Симплекс представляет из себя геометрическую фигуру, являющуюся n — мерным обобщением треугольника. Для одномерного пространства — это отрезок, для двумерного — треугольник. Таким образом n — мерный симплекс имеет n + 1 вершину.

    Алгоритм


    1) Пусть $f(x, y)$ функция, которую необходимо оптимизировать. На первом шаге выбираем три случайные точки (об этом чуть позже) и формируем симплекс (треугольник). Вычисляем значение функции в каждой точке: $f(V_1)$, $f(V_2)$, $f(V_3)$.

    Сортируем точки по значениям функции $f(x, y)$ в этих точках, таким образом получаем двойное неравенство: $f(V_2) \leq f(V_1) \leq f(V_3).$

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

    b = $V_2$, g = $V_1$, w = $V_3$, где best, good, worst — соответственно.



    2) На следующем шаге находим середину отрезка, точками которого являются g и b. Т.к. координаты середины отрезка равны полусумме координат его концов, получаем:

    $mid = \left ( \frac {x_1 + x_2} 2; \frac {y_1 + y_2} 2 \right)$


    В более общем виде можно записать так:

    $mid = \frac 1 n\sum_{i=1}^n x_i$



    3) Применяем операцию отражения:
    Находим точку $x_r$, следующим образом:

    $x_r = mid + α(mid - w)$


    Т.е. фактически отражаем точку w относительно mid. В качестве коэффициента берут как правило 1. Проверяем нашу точку: если $f(x_r) < f(g)$, то это хорошая точка. А теперь попробуем расстояние увеличить в 2 раза, вдруг нам повезет и мы найдем точку еще лучше.



    4) Применяем операцию растяжения:
    Находим точку $x_e$ следующим образом:

    $x_e = mid + γ(x_r - mid)$


    В качестве γ принимаем γ = 2, т.е. расстояние увеличиваем в 2 раза.

    Проверяем точку $x_e$:

    Если $f(x_e) < f(b)$, то нам повезло и мы нашли точку лучше, чем есть на данный момент, если бы этого не произошло, мы бы остановились на точке $x_r$.

    Далее заменяем точку w на $x_e$, в итоге получаем:



    5) Если же нам совсем не повезло и мы не нашли хороших точек, пробуем операцию сжатия.
    Как следует из названия операции мы будем уменьшать наш отрезок и искать хорошие точки внутри треугольника.

    Пробуем найти хорошую точку $x_c$:

    $x_c = mid + β(w - mid)$


    Коэффициент β принимаем равным 0.5, т.е. точка $x_c$ на середине отрезка wmid.



    Существует еще одна операция — shrink (сокращение). В данном случае, мы переопределяем весь симплекс. Оставляем только «лучшую» точку, остальные определяем следующим образом:

    $x_j = b + δ(x_j - b)$



    Коэффициент δ берут равным 0.5.

    По существу передвигаем точки по направлению к текущей «лучшей» точке. Преобразование выглядит следующим образом:



    Необходимо отметить, что данная операция дорого обходится, поскольку необходимо заменять точки в симплексе. К счастью было установлено, при проведении большого количества экспериментов, что shrink — трансформация редко случается на практике.

    Алгоритм заканчивается, когда:

    1) Было выполнено необходимое количество итераций.
    2) Площадь симплекса достигла определенной величины.
    3) Текущее лучшее решение достигло необходимой точности.

    Как и в большинстве эвристических методов, не существует идеального способа выбора инициализирующих точек. Как уже было сказано, можно брать случайные точки, находящиеся недалеко друг от друга для формирования симплекса; но есть решение и получше, которое используется в реализации алгоритма в MATHLAB:

    Выбор первой точки $V_1$ поручаем пользователю, если он имеет некоторое представление о возможном хорошем решении, в противном случае выбирается случайным образом. Остальные точки выбираются исходя из $V_1$, на небольшом расстоянии вдоль направления каждого измерения:

    $V_{i + 1} = V_i + h(V_1, i)*U_i$


    где $U_i$ — единичный вектор.
    $h(V_1, i)$ определяется таким образом:
    $h(V_1, i)$ = 0.05, если коэффициент при $U_i$ в определении $V_1$ не нулевой.
    $h(V_1, i)$ = 0.00025, если коэффициент при $U_i$ в определении нулевой.

    Пример:


    Найти экстремум следующей функции: $f(x, y) = x^2 + xy + y^2 - 6x - 9y$

    В качестве начальных возьмем точки:

    $V_1(0, 0), V_2(1, 0), V_3(0, 1)$


    Вычислим значение функции в каждой точке:

    $f(V_1) = f(0, 0) = 0$


    $f(V_2) = f(1, 0) = -5$


    $f(V_3) = f(0, 1) = -8$



    Переобозначим точки следующим образом:

    $b = V_3(0, 1), g = V_2(1, 0), w = V_1(0, 0)$





    Находим середину отрезка bg:

    $mid = \frac{b+g}2 = \left(\frac 1 2; \frac 1 2 \right)$


    Находим точку $x_r$ (операция отражения):

    $x_r = mid + α(mid - w),$


    если α=1, тогда:

    $x_r = 2*mid - w = 2 \left(\frac 1 2; \frac 1 2 \right) - \left(0, 0 \right) = (1, 1)$





    Проверяем точку $x_r$:

    $f(x_r) = -12$, т.к. $f(x_r) < f(b)$ пробуем увеличить отрезок (операция растяжения).

    $x_e = mid + γ(x_r - mid), $


    если γ = 2, тогда:

    $x_e = 2x_r - mid$


    $x_e = 2(1, 1) - \left(\frac 1 2, \frac 1 2 \right) = (1.5, 1.5)$





    Проверяем значение функции в точке $x_e$:

    $f(x_e) = f(1.5, 1.5) = -15.75$



    Оказалось, что точка $x_e $ «лучше» точки b. Следовательно мы получаем новые вершины:

    $V_1(1.5, 1.5), V_2(1, 0), V_3(0, 1)$


    И алгоритм начинается сначала.

    Таблица значений для 10 итераций:
    Best Good Worst
    $f(0, 1) = -8$ $f(1.0, 0) = -5$ $f(0, 0) = 0$
    $f(1.5, 1.5) = -15.75$ $f(0, 1) = -8$ $f(1.0, 0) = -5$
    $f(0.25, 3.75) = -20.187$ $f(1.5, 1.5) = -15.75$ $f(0, 1) = -8$
    $f(0.25, 3.75) = -20.187$ $f(1.75, 4.25) = -20.1875$ $f(1.5, 1.5) = -15.75$
    $f(1.125, 3.375) = -20.671$ $f(1.75, 4.25) = -20.1875$ $f(0.25, 3.75) = -20.1875$
    $f(1.140, 3.796) = -20.9638$ $f(1.125, 3.375) = -20.6718$ $f(1.75, 4.25) = -20.1875$
    $f(1.140, 3.796) = -20.9638$ $f(1.287, 3.751) = -20.8668$ $f(1.125, 3.375) = -20.6718$
    $f(1.140, 3.796) = -20.9638$ $f(1.236, 3.874) = -20.9521$ $f(1.287, 3.751) = -20.8668$
    $f(0.990, 4.002) = -20.9951$ $f(1.140, 3.796) = -20.9638$ $f(1.2365, 3.874) = -20.9520$
    $f(0.990, 4.002) = -20.9951$ $f(0.895, 3.925) = -20.9855$ $f(1.140, 3.796) = -20.9638$




    Аналитически находим экстремум функции, он достигается в точке $f(1, 4) = -21$.
    После 10 итераций мы получаем достаточно точное приближение: $f(0.990, 4.002) = -20.999916$

    Еще о методе:


    Алгоритм Нелдера — Мида в основном используется для выбора параметра в машинном обучении. В сущности, симплекс-метод используется для оптимизации параметров модели. Это связано с тем, что данный метод оптимизирует целевую функцию довольно быстро и эффективно (особенно там, где не используется shrink — модификация).

    С другой стороны, в силу отсутствия теории сходимости, на практике метод может приводить к неверному ответу даже на гладких (непрерывно дифференцируемых) функциях. Также возможна ситуация, когда рабочий симплекс находится далеко от оптимальной точки, а алгоритм производит большое число итераций, при этом мало изменяя значения функции. Эвристический метод решения этой проблемы заключается в запуске алгоритма несколько раз и ограничении числа итераций.

    Реализация на языке программирования python:


    Создаем вспомогательный класс Vector и перегружаем операторы для возможности производить с векторами базовые операции. Я намерено не использовал вспомогательные библиотеки для реализации алгоритма, т.к. в таком случае зачастую снижается восприятие.

    #!/usr/bin/python
    # -*- coding: utf-8 -*-
    class Vector(object):
        def __init__(self, x, y):
            """ Create a vector, example: v = Vector(1,2) """
            self.x = x
            self.y = y
    
        def __repr__(self):
            return "({0}, {1})".format(self.x, self.y)
    
        def __add__(self, other):
            x = self.x + other.x
            y = self.y + other.y
            return Vector(x, y)
    
        def __sub__(self, other):
            x = self.x - other.x
            y = self.y - other.y
            return Vector(x, y)
    
        def __rmul__(self, other):
            x = self.x * other
            y = self.y * other
            return Vector(x, y)
    
        def __truediv__(self, other):
            x = self.x / other
            y = self.y / other
            return Vector(x, y)
    
        def c(self):
            return (self.x, self.y)
            
    # objective function
    def f(point):
        x, y = point
        return x**2 + x*y + y**2 - 6*x - 9*y
    
    def nelder_mead(alpha=1, beta=0.5, gamma=2, maxiter=10):
        
        # initialization
        v1 = Vector(0, 0)
        v2 = Vector(1.0, 0)
        v3 = Vector(0, 1)
    
        for i in range(maxiter):
            adict = {v1:f(v1.c()), v2:f(v2.c()), v3:f(v3.c())}
            points = sorted(adict.items(), key=lambda x: x[1])
            
            b = points[0][0]
            g = points[1][0]
            w = points[2][0]
            
            
            mid = (g + b)/2
    
            # reflection
            xr = mid + alpha * (mid - w)
            if f(xr.c()) < f(g.c()):
                w = xr
            else:
                if f(xr.c()) < f(w.c()):
                    w = xr
                c = (w + mid)/2
                if f(c.c()) < f(w.c()):
                    w = c
            if f(xr.c()) < f(b.c()):
    
                # expansion
                xe = mid + gamma * (xr - mid)
                if f(xe.c()) < f(xr.c()):
                    w = xe
                else:
                    w = xr
            if f(xr.c()) > f(g.c()):
                
                # contraction
                xc = mid + beta * (w - mid)
                if f(xc.c()) < f(w.c()):
                    w = xc
    
            # update points
            v1 = w
            v2 = g
            v3 = b
        return b
    
    print("Result of Nelder-Mead algorithm: ")
    xk = nelder_mead()
    print("Best poits is: %s"%(xk))
    


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

    Подробнее
    Реклама
    Комментарии 17
    • –4
      Хмм, а чем Ваша реализация лучше, чем вот эта: https://docs.scipy.org/doc/scipy-0.19.0/reference/optimize.minimize-neldermead.html#optimize-minimize-neldermead?
      • +5
        Спасибо за Ваш комментарий. Цель моей статьи не показать алгоритм лучше, а рассказать о методе и представить возможную реализацию.
        • 0
          Как забавно получилось, вы молча добавили первый абзац к своей публикации после моего второго комментария, а я в результате получил кучу минусов.

          И, да, вы ошиблись, утверждая, что обсуждаемый метод:

          Используется по умолчанию в функции optimize из модуля scipy.optimize

          Нас самом деле по умолчанию используются: one of BFGS, L-BFGS-B, SLSQP, depending if the problem has constraints or bounds.

          • 0
            Спасибо за ваш ответ, Andy_U. Первый абзац был изначально, и я его никак не видоизменял. Да, вы правы. По умолчанию используется BFGS, L-BFGS-B, SLSQP. Обязательно внесу поправки в статью. Спасибо за уточнение. Всего доброго!
            • 0
              И вам спасибо за совет не читать статьи по диагонали.
      • –3
        Может быть, тогда стоило начать статью с какой-нибудь подобной фразы?
        • +1
          Может у кого-нибудь есть реализация алгоритма на других языках программирования? С радостью расширю статью, естественно указав автора реализации

          А каков смысл?
          Деформируемый многогранник Нелдера-Мида — метод старый, реализаций немеряно. У меня со студенчества (это середина 80-х) остались реализации на Алголе для БЭСМ-6 и Фортране для ЕС и СМ-4. Не думаю, что расширение ими статьи будет кому-то полезно:)
          • +1

            Писал этой весной на Scilab, затем переносил в ANSYS APDL — сначала не увидел встроенных процедур, а когда их абсолютно случайно увидел их в СТАРОМ 11-м Хэлпе и прогнал примеры — не понравилась точность встроенных в ANSYS оптимизаторов. Ничётак получилось, несколько процентов форы отжал у ANSYS =)
            И автор прав, требуется гонять алгоритм несколько раз до победного и правильно выбирать длину вектора старта и вектора корректировки на очередном витке цикла, иначе велик шанс не дойти до точки минимума. Хотя сам алгоритм работает шустро.

          • 0
            Интересно, спасибо. Странно, правда, что в этом методе не используются значения функции в точках, — все равно же их считать приходится. Ну то есть те же производные по направлениям, соединяющим точки. Ясно же, что если в одном направлении функция уменьшается на 8, а в другом только на 5, то логичнее сделать шаг в первом направлении в 8/5 раза больше, чем во втором. Быстрее сходиться должно. Думаю, такой алгоритм тоже известен.
            • +1
              Здравствуйте, dmagin! Спасибо за Ваш комментарий. Под ваше описание подходит метод покоординатного спуска (градиентные методы).
              • 0

                Тут смысл именно в том, что это derivative free метод поиска (т. н. direct search). Функции могут быть шумные и негладкие, да вообще любые. Кстати, более продвинутые методы direct search также давно существуют, например, Generalized Pattern Search, Mesh Adaptive Search, но они и более медленные.

                • 0
                  Да, спасибо, я уже заглянул в вики, чтобы понять, в чем фишка. В статье об этом не упомянуто, а пример приведен на гладкой функции, что и вызвало вопрос.
              • 0
                Я когда-то на Delphi писал для диплома )
                Еще на Фортране у меня был, откуда-то качал из математических библиотек.
                Нелдер-Мид хорош, когда функция недифференцируема.
                • +1
                  Спасибо за статью, все очень доходчиво и понятно.
                  • 0

                    Офтопик, но в каком пакете Вы делали геометрические картинки?

                    • 0
                      Здравствуйте, Arastas. Я использовал сайт: http://yotx.ru/
                    • +1
                      Писала этот алгоритм (он еще в некоторых источниках называется «метод деформируемого многогранника») на c++ для оптимизации функции с неограниченным количеством параметров.

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