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



Знаходження екстремуму (мінімуму або максимуму) цільовій функції є важливою задачею в математиці і її додатках (зокрема, в машинному навчанні є завдання curve-fitting). Напевно кожен чув про методі найшвидшого спуску (МПЗ) і методі Ньютона (МН). На жаль, ці методи мають ряд істотних недоліків, зокрема — метод найшвидшого спуску може дуже довго сходитися в кінці оптимізації, а метод Ньютона вимагає обчислення других похідних, для чого потрібно дуже багато обчислень.

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


Алгоритм Левенберга — Марквардта використовується для розв'язання нелінійного методу найменших квадратів. Стаття містить:

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

У коді використані додаткові бібліотеки, такі як numpy, matplotlib. Якщо у вас їх немає — дуже рекомендую встановити їх з пакету Anaconda for Python

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



Тому, спочатку, необхідно вивчити їх. Цим і займемося

Визначення

  • inline_formula— наша цільова функція. Ми будемо мінімізувати inline_formula. У цьому випадку, inline_formulaє функцією втрат
  • inline_formula— градієнт функції inline_formulaв точці inline_formula
  • inline_formulainline_formula, при якому inline_formulaє локальним мінімумом, тобто якщо існує проколотая околиця inline_formula, така що inline_formula
  • inline_formula— глобальний мінімум, якщо inline_formula, тобто inline_formulaне має значень менших, ніж inline_formula
  • inline_formulaматриця Якобі для функції inline_formulaв точці inline_formula. Тобто це таблиця всіх частинних похідних першого порядку. По суті, це аналог градієнта для inline_formula, так як в цьому випадку ми маємо справу з відображенням inline_formula-мірного вектора inline_formula-мірний, тому не можна просто порахувати перші похідні по одному вимірюванню, як це відбувається в градієнті. Подробнее
  • inline_formulaматриця Гессе (матриця других похідних). Необхідна для квадратичної апроксимації
Вибір функції

У математичної оптимізації є функції, на яких тестуються нові методи. Одна з таких функція — Функція Розенброка. У випадку функції двох змінних вона визначається як


Я взяв inline_formula.Додано множник 0.5 перед першою частиною для зручності. Тобто функція має вигляд:


Будемо розглядати поведінку функції на інтервалі inline_formula


Ця функція визначена неотрицательно, має мінімум inline_formulaв точці inline_formula

У коді простіше инкапсулировать всі дані про функції в один клас і брати клас тієї функції, яка буде потрібна. Результат залежить від початкової точки оптимізації. Виберемо її як inline_formula. Як видно з графіка, в цій точці функція приймає найбільше значення на інтервалі.

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)
Сам метод вкрай простий. Приймаємо inline_formula, тобто цільова функція збігається з заданої.

Потрібно знайти inline_formula— напрям найшвидшого спуску функції inline_formulaв точці inline_formula.

inline_formulaможе бути лінійно апроксимована в точці inline_formula:



де inline_formula— кут між вектором inline_formula. inline_formulaпотрібно скалярного добутку

Так як ми мінімізуємо inline_formula, то чим більше різниця в inline_formula, тим краще. При inline_formulaвираз буде максимально( inline_formula, норма вектора завжди невід'ємна), а inline_formulaтільки якщо вектора inline_formulaбудуть протилежні, тому


Напрям у нас вірне, але роблячи крок довжиною inline_formulaможна піти не туди. Робимо крок менше:


Теоретично, чим менше крок, тим краще. Але тоді постраждає швидкість збіжності. Рекомендоване значення inline_formula

У коді це виглядає так: спочатку базовий клас-оптимізатор. Передаємо все, що знадобиться в подальшому (матриці Гессе, Якобі, зараз не потрібні, але знадобляться для інших методів)

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.<strong>class</strong>.<strong>name</strong>.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 inline_formula
Впадає в очі: як швидко йшла оптимізація в 0-25 ітераціях, в 25-75 вже повільніше, а в кінці знадобилося 457 ітерацій, щоб наблизитися до нуля впритул. Така поведінка дуже властиво для МПЗ: дуже хороша швидкість збіжності спочатку, погана в кінці.

Метод Ньютона
Сам Метод Ньютона шукає корінь рівняння, тобто такий inline_formula, inline_formula. Це не зовсім те, що нам потрібно, оскільки функція може мати екстремум не обов'язково в нулі.

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

Розглянемо для inline_formula

Приймаємо inline_formula, тобто цільова функція збігається з заданої.

Розтліваємо inline_formulaв ряд Тейлора, тільки на відміну від МНС нам потрібно квадратичне наближення:


Нескладно показати, що якщо inline_formula, то функція не може мати екстремум в inline_formula. Точка inline_formulaв якій inline_formulaназивається стаціонарної.

Продиференціюємо обидві частини за inline_formula. Наша мета, щоб inline_formula, тому вирішуємо рівняння:


inline_formula— це напрямок екстремуму, але воно може бути як великим, так і мінімумом. Щоб дізнатися — чи є точка inline_formulaмінімумом — потрібно проаналізувати другу похідну. Якщо inline_formulainline_formulaє локальним мінімумом, якщо inline_formula— максимумом.

У багатовимірному випадку перша похідна замінюється на градієнт, друга — на матрицю Гессе. Ділити матриці можна, замість цього множать на зворотний (дотримуючи бік, т. к. комутативність відсутня):


Аналогічно одновимірного випадку — потрібно перевірити, чи правильно ми йдемо? Якщо матриця Гессе позитивно визначена, означає напрямок вірний, інакше використовуємо МНС.

У коді:

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)
# Якщо H - позитивно певна - Ок, інакше ми йдемо не в тому напрямку, використовуємо МНС
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 inline_formula
Порівняйте з МНС. Там був дуже сильний спуск до 25 ітерації (практично впали з гори), але потім збіжність сильно сповільнилася. В МН, навпаки, ми спочатку поволі спускаємося з гори, але потім рухаємося швидше. У МНС пішло з 25 по 532 ітерацію, щоб дійти до нуля з inline_formula. МН ж оптимізував inline_formulaза 154 останніх ітерацій.

Це часте поведінка: МН має квадратичною швидкістю збіжності, якщо починати з точки, близької до локального екстремуму. МНС ж добре працює далеко від екстремуму.

МН використовує інформацію кривизни, що було видно на малюнку вище (плавний спуск з гірки). Ще приклад, що демонструє цю ідею: на малюнку нижче червоний вектор — це напрямок МНС, а зелений — МН



[Нелінійний vs лінійний] метод найменших квадратів
В МНК у нас є модель inline_formula, що має inline_formulaпараметрів, що настроюються так, щоб мінімізувати inline_formula, де inline_formulainline_formula-е спостереження.

лінійному МНК у нас є inline_formulaрівнянь, кожне з яких ми можемо уявити як лінійне рівняння


Для лінійного МНК рішення єдино. Існують потужні методи, такі як QR декомпозиція, SVD декомпозиція, здатні знайти рішення для лінійного МНК за 1 наближене рішення матричного рівняння inline_formula.

нелінійному МНК параметр inline_formulaможе сам бути представлений функцією, наприклад inline_formula. Так само, може бути твір параметрів, наприклад inline_formulaі т. д.

Тут же доводиться знаходити рішення ітеративне, причому рішення залежить від вибору початкової точки.

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


Нічого не нагадує? Це якраз форма МНК! Введемо вектор-функцію inline_formula


і будемо підбирати inline_formulaтак, щоб вирішити систему рівнянь (хоча б наближено):


Тоді нам потрібна міра — наскільки гарна наша апроксимація. Ось вона:


Я застосував зворотну операцію: підстроїв вектор-функцію inline_formulaпід цільову inline_formula. Але можна і навпаки: якщо дано вектор-функція inline_formula, будуємо inline_formulaз (5). Наприклад:


Наостанок, один дуже важдный момент. Повинно виконуватися умова inline_formula, інакше методом користуватися не можна. У нашому випадку умова виконується

Метод Гаусса-Ньютона
Метод заснований на все тій же лінійної апроксимації, тільки тепер маємо справу з двома функціями:


Далі робимо те ж, що і в методі Ньютона — вирішуємо рівняння (тільки для inline_formula):


Нескладно показати, що поблизу inline_formula:


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

class NewtonGaussOptimizer(Optimizer):
def next_point(self):
# Вирішуємо (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 ітерації ми прийшли в точку inline_formula. Щоб продемонструвати траєкторію руху я зменшив learningrate до 0.2



Алгоритм Левенберга — Марквардта
Він заснований на одній з версій Методі Гаусса-Ньютона ("damped version" ):


Для великих inline_formulaвиходить метод найшвидшого спуску, для маленьких — метод Ньютона.
Сам алгоритм у процесі оптимізації підбирає потрібний inline_formulaна основі gain ratio, який визначається як:


Якщо inline_formulainline_formula— апроксимація гарна для inline_formula, інакше — потрібно збільшити inline_formula.

Початкове значення inline_formulaзадається як inline_formula, де inline_formula— елементи матриці inline_formula.

inline_formulaрекомендовано призначати за inline_formula. Критерієм зупинки є досягнення глобального мінімуму, тобто inline_formula


У оптимізаторах я не реалізовував критерій зупинки — за це відповідає користувач. Мені потрібно було тільки рух до наступної точки.

class LevenbergMarquardtOptimizer(Optimizer):
def __init__(self, function, initialPoint, gradient=None, jacobi=None, hesse=None,
interval=None, function_array=None, learningRate=1):
self.learningRate = learningRate
functionNew = lambda x: np.array([function(x)])
super().__init__(functionNew, initialPoint, gradient, jacobi, hesse, 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: # Закінчено. Y не може бути менше нуля
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)) # лінійний пошук
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: # Ok, хороша оптимізація
self.move_next(x_new) # ok, крок прийнятий
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 inline_formula
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

Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.