Pull to refresh

Алгоритм Левенберга — Марквардта для нелинейного метода наименьших квадратов и его реализация на Python

Reading time 9 min
Views 64K



Нахождение экстремума(минимума или максимума) целевой функции является важной задачей в математике и её приложениях(в частности, в машинном обучении есть задача curve-fitting). Наверняка каждый слышал о методе наискорейшего спуска (МНС) и методе Ньютона (МН). К сожалению, эти методы имеют ряд существенных недостатков, в частности — метод наискорейшего спуска может очень долго сходиться в конце оптимизации, а метод Ньютона требует вычисления вторых производных, для чего требуется очень много вычислений.



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



\frac{1}{2}\sum \limits_{i=1}^{N}(y_i'-y_i)^2 = \frac{1}{2}\sum \limits_{i=1}^{N}r_i^2 \tag{1}


Алгоритм Левенберга — Марквардта является нелинейным методом наименьших квадратов. Статья содержит:


  • объяснение алгоритма
  • объяснение методов: наискорейшего спуска, Ньтона, Гаусса-Ньютона
  • приведена реализация на Python с исходниками на github
  • сравнение методов


В коде использованы дополнительные библиотеки, такие как numpy, matplotlib. Если у вас их нет — очень рекомендую установить их из пакета Anaconda for Python



Зависимости
Алгоритм Левенберга — Марквардта опирается на методы, приведённые в блок-схеме




Поэтому, сначала, необходимо изучить их. Этим и займёмся

Определения


  • F=F(x) — наша целевая функция. Мы будем минимизировать F. В этом случае, F является функцией потерь
  • \nabla f(x) — градиент функции f в точке x
  • x^*x, при котором F(x^*) является локальным минимумом, т.е. если существует проколотая окрестность \dot{U}(x), такая что \forall x \in \dot{U}(x)\ \ \ F(x^*)\le F(x)
  • x^+ — глобальный минимум, если \forall x \in R^n\ \ \ F(x^+)\le F(x), т.е. F не имеет значений меньших, чем F(x^+)
  • J_f=J_f(x) = Jматрица Якоби для функции f: R^n \to R^m в точке x. Т.е. это таблица всех частных производных первого порядка. По сути, это аналог градиента для f:R^n \to R, так как в этом случае мы имеем дело с отображением из n-мерного вектора в m-мерный, поэтому нельзя просто посчитать первые производные по одному измерению, как это происходит в градиенте. Подробнее
  • H_f = H_f(x) = Hматрица Гессе (матрица вторых производных). Необходима для квадратичной аппроксимации

Выбор функции



В математической оптимизации есть функции, на которых тестируются новые методы.
Одна из таких функция — Функция Розенброка. В случае функции двух переменных она определяется как



f(x,y) = (a-x)^2 + b(y-x^2)^2


Я принял a=0.5,\ b=0.5. Т.е. функция имеет вид:



f(x,y) = \frac{1}{2}(1-x)^2 + \frac{1}{2}(y-x^2)^2


Будем рассматривать поведение функции на интервале -2 \le x,y \le 2




Эта функция определена неотрицательно, имеет минимум z = 0$ в точке $(x=1,y=1)
В коде проще инкапсулировать все данные о функции в один класс и брать класс той функции, которая потребуется. Результат зависит от начальной точки оптимизации. Выберем её как (x=-2,y=-2). Как видно из графика, в этой точке функция принимает наибольшее значение на интервале.



functions.py


class Rosenbrock:

    initialPoint = (-2, -2)
    camera = (41, 75)
    interval = [(-2, 2), (-2, 2)]

   """
   Целевая функция
   """

   @staticmethod
   def function(x):

       return 0.5*(1-x[0])**2 + 0.5*(x[1]-x[0]**2)**2

   """
   Для нелинейного МНК - возвращает вектор-функцию r
   """
   @staticmethod
   def function_array(x):
       return np.array([1 - x[0] , x[1] - x[0] ** 2]).reshape((2,1))

   @staticmethod
   def gradient(x):
       return np.array([-(1-x[0]) - (x[1]-x[0]**2)*2*x[0], (x[1] - x[0]**2)])

   @staticmethod
   def hesse(x):
       return np.array(((1 -2*x[1] + 6*x[0]**2, -2*x[0]), (-2 * x[0], 1)))

   @staticmethod
   def jacobi(x):
       return np.array([ [-1, 0], [-2*x[0], 1]])

   """
   Векторизация для отрисовки поверхности
   Детали: http://www.mathworks.com/help/matlab/matlab_prog/vectorization.html
   """
   @staticmethod
   def getZMeshGrid(X, Y):
       return 0.5*(1-X)**2 + 0.5*(Y - X**2)**2

Метод наискорейшего спуска(Steepest Descent)



Сам метод крайне прост. Принимаем F(x)=f(x), т.е. целевая функция совпадает с заданной.
Нужно найти d_{нс}(x) — направление наискорейшего спуска функции f в точке x.
f(x) может быть линейно аппроксимирована в точке x:



f(x+d) \approx f(x) + \nabla f(x)^Td, \ d \in R^n , ||d|| \to 0 \tag{2}


\lim_{d \to 0}f(x)-f(x+d) = - \nabla f(x)^Td  \stackrel{(3.a)} = - || \nabla f(x)^T|| \ ||d|| cos \theta ,\tag{3}


где \theta — угол между вектором d \ и \nabla f(x)^T. (3.a) следует из скалярного произведения



Так как мы минимизируем f(x), то чем больше разница в (3), тем лучше. При \theta = \pi выражение будет максимально( cos \theta = -1, норма вектора всегда неотрицательна), а \theta = \pi будет только если вектора d \ и \  \nabla f(x)^T будут противоположны, поэтому



d_{нс} =  -\nabla f(x)^T


Направление у нас верное, но делая шаг длиной ||d_{нс}|| можно уйти не туда. Делаем шаг меньше:



d_{нс} =  -\alpha \nabla f(x)^T,  0 < \alpha < 1


Теоретически, чем меньше шаг, тем лучше. Но тогда пострадает скорость сходимости. Рекомендуемое значение \alpha = 0.05



В коде это выглядит так: сначала базовый класс-оптимизатор. Передаём всё, что понадобится в дальнейшем( матрицы Гессе, Якоби, сейчас не нужны, но понадобятся для других методов)


class Optimizer:
    def __init__(self, function, initialPoint, gradient=None, jacobi=None, hesse=None,
                 interval=None, epsilon=1e-7, function_array=None, metaclass=ABCMeta):
        self.function_array = function_array
        self.epsilon = epsilon
        self.interval = interval
        self.function = function
        self.gradient = gradient
        self.hesse = hesse
        self.jacobi = jacobi
        self.name = self.__class__.__name__.replace('Optimizer', '')
        self.x = initialPoint
        self.y = self.function(initialPoint)

   "Возвращает следующую координату по ходу оптимизационного процесса"
   @abstractmethod
   def next_point(self):
       pass

   """
   Движемся к следующей точке
   """

   def move_next(self, nextX):
       nextY = self.function(nextX)
       self.y = nextY
       self.x = nextX
       return self.x, self.y
 


Код самого оптимизатора:
class SteepestDescentOptimizer(Optimizer):
    ...
    def next_point(self):
        nextX = self.x - self.learningRate * self.gradient(self.x)
        return self.move_next(nextX)
 


Результат оптимизации:




Итерация X Y Z
25 0.383 -0.409 0.334
75 0.693 0.32 0.058
532 0.996 0.990 10^{-6}

Бросается в глаза: как быстро шла оптимизация в 0-25 итерациях, в 25-75 уже медленне, а в конце потребовалось 457 итераций, чтобы приблизиться к нулю вплотную. Такое поведение очень свойственно для МНС: очень хорошая скорость сходимости вначале, плохая в конце.



Метод Ньютона



Сам Метод Ньютона ищет корень уравнения, т.е. такой x, что f(x)=0. Это не совсем то, что нам нужно, т.к. функция может иметь экстремум не обязательно в нуле.



А есть ещё Метод Ньютона для оптимизации. Когда говорят о МН в контексте оптимизации — имеют в виду его. Я сам, учась в институте, спутал по глупости эти методы и не мог понять фразу «Метод Ньютона имеет недостаток — необходимость считать вторые производные».



Рассмотрим для f(x): R \to R
Принимаем F(x)=f(x), т.е. целевая функция совпадает с заданной.



Разлагаем f(x) в ряд Тейлора, только в отличии от МНС нам нужно квадратичное приближение:



f(x+d) \approx f(x) + f^{'}(x)d + \frac{1}{2} f^{''}(x)d^2, \ d \in R^n , ||d|| \to 0 \tag{4}


Несложно показать, что если f{'}(x) \ne 0, то функция не может иметь экстремум в x. Точка x^*$ в которой $f{'}(x) = 0 называется стационарной.



Продифференцируем обе части по d. Наша цель, чтобы f(x+d)^{'}=0, поэтому решаем уравнение:



0 = f(x+d)^{'} = f^{'}(x) + f^{''}(x)d  \\
d_{н}=- \frac{f^{'}(x)}{f^{''}(x)}


d_н — это направление экстремума, но оно может быть как максимумом, так и минимумом. Чтобы узнать — является ли точка x+d_н минимумом — нужно проанализировать вторую производную. Если f^{''}(x)>0$, то $f(x+d_н) является локальным минимумом, если f^{''}(x)<0 — максимумом.



В многомерном случае первая производная заменяется на градиент, вторая — на матрицу Гессе. Делить матрицы нельзя, вместо этого умножают на обратную(соблюдая сторону, т.к. коммутативность отсутствует):



f(x): R^n \to R \\
H(x)d_{н}=- \nabla f(x)\\
d_{н}=- H^{-1}(x)\nabla f(x)


Аналогично одномерному случаю — нужно проверить, правильно ли мы идём? Если матрица Гессе положительно определена, значит направление верное, иначе используем МНС.



В коде:


def is_pos_def(x):
    return np.all(np.linalg.eigvals(x) > 0)

class NewtonOptimizer(Optimizer):
    def next_point(self):
        hesse = self.hesse(self.x)
        # if Hessian matrix if positive - Ok, otherwise we are going in wrong direction, changing to gradient descent
        if is_pos_def(hesse):
            hesseInverse = np.linalg.inv(hesse)
            nextX = self.x - self.learningRate * np.dot(hesseInverse, self.gradient(self.x))
        else:
            nextX = self.x - self.learningRate * self.gradient(self.x)

   return self.move_next(nextX)
 


Результат:




Итерация X Y Z
25 -1.49 0.63 4.36
75 0.31 -0.04 0.244
179 0.995 -0.991 10^{-6}

Сравните с МНС. Там был очень сильный спуск до 25 итерации( практически упали с горы), но потом сходимость сильно замедлилась. В МН, напротив, мы сначала медленно спускаемся с горы, но затем движемся быстрее. У МНС ушло с 25 по 532 итерацию, чтобы дойти до нуля с z=0.334. МН же оптимизировал 4.36 за 154 последних итераций.



Это частое поведение: МН обладает квадратичной скоростью сходимости, если начинать с точки, близкой к локальному экстремуму. МНС же хорошо работает далеко от экстремума.



МН использует информацию кривизны, что было видно на рисунке выше(плавный спуск с горки).
Ещё пример, демонстрирующий эту идею: на рисунке ниже красный вектор — это направление МНС, а зелёный — МН




[Нелинейный vs линейный] метод наименьших квадратов



В МНК у нас есть модель y=f(\beta_1,..\beta_n;x), имеющая n параметров, которые настраиваются так, чтобы минимизировать



\frac{1}{2}\sum \limits_{i=1}^{N}(y_i'-y_i)^2 = \frac{1}{2}\sum \limits_{i=1}^{N}r_i^2


, где y_i'i-е наблюдение.



В линейном МНК у нас есть $m$ уравнений, каждое из которых мы можем представить как линейное уравнение



x_i\beta_1+x_i\beta_2+..x_i\beta_n = y_i


Для линейного МНК решение единственно. Существуют мощные методы, такие как QR декомпозиция, SVD декомпозиция, способные найти решение для линейного МНК за 1 приближённое решение матричного уравнения Ax=b.



В нелинейном МНК параметр \beta_i может сам быть представлен функцией, например \beta_i^2. Так же, может быть произведение параметров, например



\beta_1 \beta_2


и т.д.
Здесь же приходится находить решение итеративно, причём решение зависит от выбора начальной точки.



Методы ниже имеют дело как раз с нелинейным случаем. Но, сперва, рассмотрим нелиненый МНК в контексте нашей задачи — минимизации функции



f:R^2 \to R \\
F(x_1,x_2) = f(x_1,x_2) = \frac{1}{2}(1-x_1)^2 + \frac{1}{2}(x_2-x_1^2)^2 =  \frac{1}{2}r_1^2(x_1,x_2) + \frac{1}{2}r_2^2(x_1,x_2)


Ничего не напоминает? Это как раз форма МНК! Введём вектор-функцию r



r: R^2 \to R^2 \\
r= \left [ \begin{matrix} 1-x_1 \\ x_2-x_1^2 \end{matrix} \right ]


и будем подбирать x_1,x_2 так, чтобы решить систему уравнений(хотя бы приближённо):



\begin{cases} 
r_1 = 1-x_1 = 0 \\
r_2 = x_2 - x_1^2 = 0
\end{cases} \\


Тогда нам нужна мера — насколько хороша наша аппроксимация. Вот она:



F(x) = \frac{1}{2}\sum_i^m r_i^2(x) =  \frac{1}{2} r^T r=   \frac{1}{2} ||r||^2 \tag{5}


Я применил обратную операцию: подстроил вектор-функцию r под целевую F. Но можно и наоборот: если дана вектор-функция r: R^n \to R^m, строим F(x) из (5). Например:



r= \left [ \begin{matrix} x_1^2 \\ x_2^2 \end{matrix} \right ], F(x) = \frac{1}{2} x_1^2 +  \frac{1}{2} x_2^2


Напоследок, один очень важдный момент. Должно выполняться условие m \ge n, иначе методом пользоваться нельзя. В нашем случае условие выполняется



Метод Гаусса-Ньютона



Метод основан на всё той же линейной аппроксимации, только теперь имеем дело с двумя функциями:



r(x+d) \approx l(d) \equiv r(x) + J(x)d \\
F(x+d)  \approx L(d) \equiv \frac{1}{2}l^T(d) l(d)


Далее делаем то же, что и в методе Ньютона — решаем уравнение(только для L(d)):



L^{''}d_{гн} = -L^{'}


Несложно показать, что вблизи \text(\ d \to 0):



L^{''}(d) = J_r^T J_r, \ L^{'}(d) = J_r^T(d) r(d)  \\
(J^T J)d_{гн} = -J^Tr \\
d_{гн} = -(J^T J)^{-1} J^Tr


Код оптимизатора:


class NewtonGaussOptimizer(Optimizer):
    def next_point(self):
        # Solve (J_t * J)d_ng = -J*f
        jacobi = self.jacobi(self.x)
        jacobisLeft = np.dot(jacobi.T, jacobi)
        jacobiLeftInverse = np.linalg.inv(jacobisLeft)
        jjj = np.dot(jacobiLeftInverse, jacobi.T)  # (J_t * J)^-1 * J_t
        nextX = self.x - self.learningRate * np.dot(jjj, self.function_array(self.x)).reshape((-1))
        return self.move_next(nextX)


Результат превысил мои ожидания. Всего за 3 итерации мы пришли в точку (x=1, y=1). Чтобы продемонстрировать траекторию движения я уменьшил learningrate до 0.2




Алгоритм Левенберга — Марквардта



Он основан на одной из версий Методе Гаусса-Ньютона( "damped version" ):



(J^T J + \mu I)d_{лм} = -J^Tr , \mu \ge 0


\mu называется параметром регулизации. Иногда I заменяют на diag(J^T J) для улучшения сходимости.
Диагональные элементы J^T J будут положительны, т.к. элемент a_{ii} матрицы J^T J является скалярным произведением вектора-строки i в J^T на самого себя.



Для больших \mu получается метод наискорейшего спуска, для маленьких — метод Ньютона.
Сам алгоритм в процессе оптимизации подбирает нужный \mu на основе gain ratio, определяющийся как:



g=\frac{F(x) - F(x_{new)}}{L(0) - L(d_{лм})}


Если g>0, то L(d) — хорошая аппроксимация для F(x+d), иначе — нужно увеличить \mu.
Начальное значение \mu задаётся как \tau \cdot  max\{{a_{ij}}\}, где a_{ij} — элементы матрицы J^T J.
\tau рекомендовано назначать за 10^{-3}. Критерием остановки является достижение глобального минимуму, т.е. F^{'}(x^*)=g(x^*)=0




В оптимизаторах я не реализовывал критерий остановки — за это отвечает пользователь. Мне нужно было только движение к следующей точке.


class LevenbergMarquardtOptimizer(Optimizer):
    def __init__(self, function, initialPoint, gradient=None, jacobi=None, hessian=None,
                 interval=None, function_array=None, learningRate=1):
        self.learningRate = learningRate
        functionNew = lambda x: np.array([function(x)])
        super().__init__(functionNew, initialPoint, gradient, jacobi, hessian, interval, function_array=function_array)
        self.v = 2
        self.alpha = 1e-3
        self.m = self.alpha * np.max(self.getA(jacobi(initialPoint)))

   def getA(self, jacobi):
       return np.dot(jacobi.T, jacobi)

   def getF(self, d):
       function = self.function_array(d)
       return 0.5 * np.dot(function.T, function)

   def next_point(self):
       if self.y==0: # finished. Y can't be less than zero
           return self.x, self.y

    jacobi = self.jacobi(self.x)
    A = self.getA(jacobi)
    g = np.dot(jacobi.T, self.function_array(self.x)).reshape((-1, 1))
    leftPartInverse = np.linalg.inv(A + self.m * np.eye(A.shape[0], A.shape[1]))
    d_lm = - np.dot(leftPartInverse, g) # moving direction
    x_new = self.x + self.learningRate * d_lm.reshape((-1)) # line search
    grain_numerator = (self.getF(self.x) - self.getF(x_new))
    gain_divisor = 0.5* np.dot(d_lm.T, self.m*d_lm-g) + 1e-10
    gain = grain_numerator / gain_divisor
    if gain > 0: # it's a good function approximation.
        self.move_next(x_new) # ok, step acceptable
        self.m = self.m * max(1 / 3, 1 - (2 * gain - 1) ** 3)
        self.v = 2
    else:
        self.m *= self.v
        self.v *= 2

    return self.x, self.y


Результат получился тоже хороший:


Итерация X Y Z
0 -2 -2 22.5
4 0.999 0.998 10^{-7}
11 1 1 0

При learningrate =0.2:




Сравнение методов


Название метода Целевая функция Достоинства Недостатки Сходимость
Метод наискорейший спуск дифференцируемая -широкий круг применения
-простая реализация

-низкая цена одной итерации
-глобальный минимум ищется хуже, чем в остальных методах

-низкая скорость сходимости вблизи экстремума
локальная
Метод Нютона дважды дифференцируемая -высокая скорость сходимости вблизи экстремума

-использует информацию о кривизне
-функция должны быть дважды дифференцируема

-вернёт ошибку, если матрица Гессе вырождена ( не имеет обратной)

-есть шанс уйти не туда, если находится далеко от экстремума
локальная
Метод Гаусса-Нютона нелинейный МНК -очень высокая скорость сходимости

-хорошо работает с задачей curve-fitting
-колонки матрицы J должны быть линейно-независимы

-налагает ограничения на вид целевой функции
локальная
Алгоритм Левенберга — Марквардта нелинейный МНК -наибольная устойчивость среди рассмотренных методов

-наибольшие шансы найти глобальный экстремум

-очень высокая скорость сходимости(адаптивная)

-хорошо работает с задачей curve-fitting
-колонки матрицы J должны быть линейно-независимы

-налагает ограничения на вид целевой функции

-сложность в реализации
локальная

Несмотря на хорошие результаты в конкретном примере рассмотренные методы не гарантируют глобальную сходимость(найти которую — крайне трудная задача). Примером из немногих методов, позволяющих всё же достичь этого, является алгоритм basin-hopping



Совмещённый результат(специально понижена скорость последних двух методов):




Исходники можно скачать с github



Источники


  1. K. Madsen, H.B. Nielsen, O. Tingleff(2004): Methods for non-linear least square
  2. Florent Brunet(2011): Basics on Continuous Optimization
  3. Least Squares Problems
Tags:
Hubs:
+76
Comments 28
Comments Comments 28

Articles