Заявление об ограничении ответственности: это очень объемное сообщение в блоге, включающее математические доказательства и реализации на Python для различных алгоритмов оптимизации.
Оптимизация, одна из самых интересных тем в области машинного обучения. Большинство проблем, с которыми мы сталкиваемся в повседневной жизни, решаются с помощью методов численной оптимизации. В этом блоге давайте подробно рассмотрим некоторые базовые алгоритмы численной оптимизации для поиска локальных оптимумов любой заданной функции (которая лучше всего работает с выпуклыми функциями). Давайте начнем с простых выпуклых функций, в которых локальный и глобальный минимумы совпадают, и перейдем к сильно нелинейным функциям с множественными локальными
и глобальными минимумами. Вся оптимизация вращается вокруг основных понятий линейной алгебры и дифференциального исчисления. Недавние обновления в области глубокого обучения вызвали огромный интерес к алгоритмам численной и стохастической оптимизации, чтобы теоретически обосновать потрясающие качественные результаты, демонстрируемые сетями глубокого обучения. В таких алгоритмах обучения не существует какой-либо явно известной функции для оптимизации, но у нас есть доступ только к оракулам 0-го и 1-го порядка. Оракулы - это черные ящики, которые возвращают значение функции (0-й порядок), градиент (1-й порядок) или гессиан (2-й порядок) функции в любой заданной точке. Этот блог предоставляет базовое теоретическое и численное понимание функций неограниченной и ограниченной оптимизации, а также включает их реализацию на языке Python.
Необходимые и достаточные условия для того, чтобы точка была локальным минимумом:
Пусть f (.) - непрерывная дважды дифференцируемая функция. Чтобы любая точка была минимумом, она должна удовлетворять следующим условиям:
- Необходимое условие первого порядка:
- Необходимое условие второго порядка:
- Условие достаточности второго порядка:
Градиентный спуск:
Градиентный спуск является основой всех достижений в области алгоритмов обучения (машинное обучение, глубокое обучение или глубокое обучение с подкреплением). В этом разделе мы увидим различные модификации градиентного спуска для более быстрой и лучшей сходимости. Рассмотрим случай линейной регрессии, где оценим коэффициенты уравнения:
Предполагая, что функция представляет собой линейную комбинацию всех функций. Оптимальный набор коэффициентов определяется путем минимизации функции потерь:
что оказывается оценкой максимального правдоподобия для задачи линейной регрессии. Минимизация OLS (обыкновенного наименьшего квадрата) включает в себя нахождение обратной матрицы признаков, заданной следующим образом:
Для практической задачи размерность данных огромна, что легко подорвет вычислительную мощность. Например, давайте рассмотрим проблему анализа признаков изображения: обычно изображения имеют размер 1024x1024, что означает, что количество признаков порядка 10⁶. Из-за большого количества функций такие задачи оптимизации решаются только итеративным способом, и это приводит нас к очень широко известному алгоритму, называемому градиентным спуском и методу Ньютона-Рафсона.
Алгоритм градиентного спуска:
Алгоритм градиентного спуска обновляет итерацию (X) в направлении отрицательного градиента (следовательно, в направлении наискорейшего спуска) с ранее указанной скоростью обучения (eta). Скорость обучения используется для предотвращения превышения
локальных минимумов на любой данной итерации.
На рисунке ниже показана сходимость алгоритма градиентного спуска для функции f (x) = x² с eta = 0,25
Поиск оптимального eta - это непосредственная задача, которая требует предварительных знаний о понимании функций и области действия.
# python implementation of vanilla gradient descent update ruledef def gradient_descent_update (x, eta):
"""
get_gradient is 1st order oracle
"""
return x - eta*get_gradient(x)
Градиентный спуск с условием Армийо Гольдштейна:
Чтобы сократить объем ручной настройки, применяются условия Армиджо Гольдштейна (AG), чтобы найти (eta) для следующей итерации. Формальный вывод условия AG требует знания линейного приближения, условий Липшица и базового исчисления, о которых я не буду рассказывать в этом блоге.
Давайте определим две функции f1 (x) и f2 (x) как линейную аппроксимацию f (x) с двумя разными коэффициентами альфа и бета, задаваемыми следующим образом:
В каждой итерации условия AG конкретная эта, удовлетворяющая соотношению:
найден, и текущая итерация обновляется соответствующим образом.
На рисунке ниже показан диапазон следующей итерации для сходимости функции f (x) = x² с альфа = 0,25 и бета = 0,5:
Красная, синяя и зеленая линии на рисунке выше соответствуют
зеленому диапазону следующей возможной итерации.
# python implementation of gradient descent with AG condition update rule
def gradient_descent_update_AG(x, alpha =0.5, beta =0.25):
eta =0.5
max_eta =np.inf
min_eta =0.
value = get_value(x)
grad = get_gradient(x)
while True :
x_cand = x - (eta)*grad
f = get_value(x_cand)
f1 = value - eta *a lpha*np.sum(np.abs(grad)* *2 )
f2 = value - eta *be ta *np.sum(np.abs(grad)* *2 )
if f<=f2 and f>=f1:
return x_cand
if f <= f1:
if eta == max_eta:
eta = np.min([2.*eta, eta + max_eta/2.])
else:
eta = np.min([2.*eta, (eta + max_eta)/2.])
if f>=f2:
max_eta = eta
eta = (eta+min_eta)/2.0
Градиентный спуск с условием полной релаксации:
В случае условия полной релаксации новая функция g (eta) минимизируется для получения eta, которая впоследствии используется для поиска следующей итерации.
Этот метод включает решение задачи оптимизации для поиска каждой следующей итерации, которая выполняется следующим образом:
# The python implementation for FR is shown below:
def gradient_descent_update_FR (x):
eta = 0.5
thresh = 1e-6
max_eta = np.inf
min_eta = 0.
while True :
x_cand = x - eta*get_gradient(x)
g_diff = -1. *get_gradient(x)*get_gradient(x_cand)
if np.sum(np.abs(g_diff)) < thresh and eta > 0 :
return x_cand
if g_diff > 0:
if eta == max_eta:
eta = np.min([2.*eta, eta + max_eta/2.])
else:
eta = np.min([2.*eta, (eta + max_eta)/2.])
else:
max_eta = eta
eta = (min_eta + eta)/2.0
Стохастический градиентный спуск
Как мы знаем, в реальной установке размер данных будет довольно большим, что делает дальнейшее вычисление градиента для всех функций дорогостоящим. В таких случаях выбирается случайный набор точек (объектов) и вычисляется ожидаемый градиент. Вся эта установка сводится только к математическому ожиданию.
Математически это будет означать:
Случайным образом выбрать точку (p) для оценки градиентов.
В приведенной выше итерации wt можно рассматривать как шум. Эта итерация приводит к локальным оптимумам только тогда, когда E (wt) стремится к 0.
Точно так же можно показать, что дисперсия wt также конечна. Приведенное выше доказательство гарантирует сходимость SGD.
# SGD implementation in python
def SGD (self, X, Y, batch_size, thresh= 1 ):
loss = 100
step = 0
if self.plot: losses = []
while loss >= thresh:
# mini_batch formation
index = np.random.randint( 0 , len(X), size = batch_size)
trainX, trainY = np.array(X)[index], np.array(Y)[index]
self.forward(trainX)
loss = self.loss(trainY)
self.gradients(trainY)
# update
self.w0 -= np.squeeze(self.alpha*self.grad_w0)
self.weights -= np.squeeze(self.alpha*self.grad_w)
if self.plot: losses.append(loss)
if step % 1000 == 999 :
print "Batch number: {}" .format(step)+ " current loss: {}" .format(loss)
step += 1
if self.plot : self.visualization(X, Y, losses)
pass
AdaGrad
Adagrad - это оптимизатор, который помогает в автоматической адаптации скорости обучения для каждой функции, задействованной в задаче оптимизации. Это достигается за счет отслеживания истории всех градиентов. Этот метод также сходится только в смысле ожидания.
def AdaGrad (self, X, Y, batch_size,thresh=0.5 ,epsilon=1e-6 ):
loss = 100
step = 0
if self.plot: losses = []
G = np.zeros((X.shape[ 1 ], X.shape[ 1 ]))
G0 = 0
while loss >= thresh:
# mini_batch formation
index = np.random.randint( 0 , len(X), size = batch_size)
trainX, trainY = np.array(X)[index], np.array(Y)[index]
self.forward(trainX)
loss = self.loss(trainY)
self.gradients(trainY)
G += self.grad_w.T*self.grad_w
G0 += self.grad_w0** 2
den = np.sqrt(np.diag(G)+epsilon)
delta_w = self.alpha*self.grad_w / den
delta_w0 = self.alpha*self.grad_w0 / np.sqrt(G0 + epsilon)
# update parameters
self.w0 -= np.squeeze(delta_w0)
self.weights -= np.squeeze(delta_w)
if self.plot: losses.append(loss)
if step % 500 == 0 :
print "Batch number: {}".format (step)+ " current loss: {}".format(loss)
step += 1
if self.plot : self.visualization(X, Y, losses)
pass
Давайте перейдем к методам на основе Гессе, Ньютона и квазиньютона:
Методы на основе Гессе - это методы на основе градиента второго порядка для оптимизации, геометрически они включают информацию о градиенте и кривизне для обновления весов. По этой причине сходимость происходит намного быстрее по сравнению с методами, основанными только на градиенте. Правило обновления для метода Ньютона определяется как :
Сходимость этого алгоритма намного быстрее, чем у градиентных методов. Математически скорость сходимости градиентного спуска пропорциональна O (1 / t), тогда как для метода Ньютона она пропорциональна O (1 / t²). Но в случае данных более высокой размерности оценка оракула 2-го порядка для каждой итерации требует больших вычислительных ресурсов, что приводит к моделированию оракула 2-го порядка с использованием оракула 1-го порядка. Это дает нам квазиньютоновский класс алгоритмов. Чаще всего в этом классе квазиньютоновских методов используются алгоритмы BFGS и LMFGS. В этом разделе мы обсуждаем только алгоритм BFGS, который включает обновление матрицы ранга один гессиана функции. Общая идея метода состоит в том, чтобы случайным образом инициализировать Гессен и продолжать обновлять Гессе с каждой итерацией, используя правило обновления первого ранга. Математически это показано как:
# python implementation for BFGS
def BFGS_update (H, s, y):
smooth = 1e-3
s = np.expand_dims(s, axis= -1 )
y = np.expand_dims(y, axis= -1 )
rho = 1. /(s.T.dot(y) + smooth)
p = (np.eye(H.shape[ 0 ]) - rho*s.dot(y.T))
return p.dot(H.dot(p.T)) + rho*s.dot(s.T)
def quasi_Newton (x0, H0, num_iter= 100 , eta= 1 ):
xo, H = x0, H0
for _ in range(num_iter):
xn = xo - eta*H.dot(get_gradient(xo))
y = get_gradient(xn) - get_gradient(xo)
s = xn - xo
H = BFGS_update(H, s, y)
xo = xn
return xo
Ограниченная оптимизация
Теперь пришло время обсудить несколько ключевых концепций, связанных с оптимизацией с ограничениями (которая включает в себя формулировку проблемы и стратегии решения). В этом разделе также обсуждаются теория и реализация Python алгоритма, известного как SVM (машина опорных векторов). Когда мы сталкиваемся с реальными проблемами, придумать идеальную функцию оптимизации довольно сложно, а иногда и невозможно, поэтому мы ослабляем функцию оптимизации, налагая дополнительные ограничения на проблему, и оптимизация этой ограниченной настройки обеспечит приближение, которое мы могли бы быть как можно ближе к реальному решению проблемы, но также и выполнимо. Для решения задач оптимизации с ограничениями используются такие методы, как лагранжева формулировка, методы штрафов, прогнозируемый градиентный спуск, внутренние точки и многие другие методы. В этом разделе мы рассмотрим лагранжеву формулировку и прогнозируемый градиентный спуск. В этом разделе также рассказывается о наборе инструментов с открытым исходным кодом, используемом для оптимизации CVXOPT, а также рассматривается реализация SVM с использованием этого набора инструментов.
Общая форма задач оптимизации с ограничениями:
где f (x) - целевая функция, g (x) и h (x) - ограничения неравенства и равенства соответственно. Если f (x) выпукло, а ограничения образуют выпуклое множество (т.е. g (x) выпукло, а h (x) аффинно), оптимизация гарантированно сходится в глобальных минимумах. Для любого другого набора задач он сходится к
локальным минимумам.
Прогнозируемый градиентный спуск
Предварительный (и наиболее очевидный) шаг при решении настройки оптимизации с ограничениями - это проекция итерации по ограниченному набору. Оказывается, это самый мощный алгоритм в решении задачи оптимизации с ограничениями
. Это включает в себя два шага (1) для поиска следующей возможной итерации в направлении минимизации (спуска), (2) поиск проекции итерации на ограниченное множество.
# projected gradient descent implementation
def projection_oracle_l2 (w, l2_norm):
return (l2_norm/np.linalg.norm(w))*w
def projection_oracle_l1 (w, l1_norm):
# first remember signs and store them. Modify w
signs = np.sign(w)
w = w*signs
# project this modified w onto the simplex in first orthant.
d=len(w)
if np.sum(w)<=l1_norm:
return w*signs
for i in range(d):
w_next = w+ 0
w_next[w> 1e-7 ] = w[w> 1e-7 ] - np.min(w[w> 1e-7 ])
if np.sum(w_next)<=l1_norm:
w = ((l1_norm - np.sum(w_next))*w + (np.sum(w) -
l1_norm)*w_next)/(np.sum(w)-np.sum(w_next))
return w*signs
else :
w=w_next
def main ():
eta= 1.0 /smoothness_coeff
for l1_norm in np.arange( 0 , 10 , 0.5 ):
w=np.random.randn(data_dim)
for i in range( 1000 ):
w = w - eta*get_gradient(w)
w = projection_oracle_l1(w, l1_norm)
pass
Понимание лагранжевой формулировки
В большинстве задач оптимизации поиск проекции итерации на ограниченный набор является сложной задачей (особенно в случае сложного ограниченного набора). Это похоже на решение задачи оптимизации на каждой итерации, которая в большинстве случаев была бы невыпуклой. На самом деле люди пытаются избавиться от ограничений, решая двойную проблему, а не основную проблему.
Прежде чем погрузиться в лагранжеву дуальную и первичную формулировки, давайте разберемся с условиями KKT и их значение.
- Чтобы любая точка была локальным / глобальным минимумом с ограничениями равенства:
- Аналогично в случае ограничений неравенства:
Эти два условия могут быть легко соблюдены, если рассматривать единичный круг как ограниченное множество. В первой части рассмотрим только границу, где знак (\ mu) не имеет значения, это результат ограничений равенства. Во втором случае рассмотрим внутренний набор единичного круга, где -ve знак для (\ lambda) означает область допустимого решения.
Условия ККТ (Каруш – Куна – Таккера) рассматриваются как необходимые условия первого порядка, которым должна удовлетворять точка, чтобы считаться стационарной точкой (локальные минимумы, локальные максимумы, седловая точка. При x ^ *
чтобы быть локальными минимумами:
Условие LICQ: все активные ограничения
должны быть линейно независимыми друг от друга.
Лагранжиан функции
Для любой функции f (x) с ограничениями-неравенствами g_i (x) ≤ 0 и ограничениями-равенствами h_i (x) = 0 лагранжиан L (x, \ lambda, \ mu)
Функция оптимизации:
вышеуказанная оптимизация эквивалентна:
Вышеупомянутая формулировка известна как основная проблема (p ^ *) в соответствии с утверждением теории игр (т.е. у второго игрока всегда будет больше шансов на оптимизацию), легко увидеть, что:
Эта формулировка известна как двойная задача (d ^ *).
И прямая, и двойственная формулировки приводят к одному и тому же оптимальному значению тогда и только тогда, когда целевая функция является выпуклой и ограниченное множество также является выпуклым. Доказательство этого утверждения обосновывается ниже:
пусть функция выпуклая, что влечет
SVM (машины опорных векторов)
SVM принадлежит к классу алгоритмов контролируемого обучения, используемых в задачах классификации и регрессии. SVM легко масштабируются и могут решать линейные и нелинейные задачи (используя уловку на основе ядра). В большинстве задач мы не можем провести единую линию, разграничивающую два разных класса данных, из-за чего нам нужно обеспечить небольшой запас при построении границ принятия решений, что легко формулируется с помощью SVM. Идея SVM состоит в том, чтобы создать разделяющую гиперплоскость между двумя разными наборами точек данных. После получения разделяющей гиперплоскости становится относительно проще классифицировать точки данных (контрольный пример) в один из классов. SVM хорошо работает даже с данными более высокой размерности. Преимущества SVM в том, что они эффективны с точки зрения памяти, точны и быстры по сравнению с другими моделями машинного обучения. Давайте посмотрим на математику SVM
Основная проблема SVM
Двойная задача SVM с использованием лагранжиана
Происхождение Dual
CVXOPT
В этом разделе мы обсудим реализацию описанного выше двойного алгоритма SVM на Python с использованием библиотеки CVXOPT. Вы можете узнать больше о библиотеке (http://cvxopt.org/userguide/intro.html)
Общий формат CVXOPT для таких проблем
# SVM using CVXOPT
import numpy as np
from cvxopt import matrix,solvers
def solve_SVM_dual_CVXOPT (x_train, y_train, x_test):
"""
Solves the SVM training optimisation problem (the Arguments:
x_train: A numpy array with shape (n,d), denoting \R^d.
y_train: numpy array with shape (n,) Each element
x_train: A numpy array with shape (m,d), denoting dual) using cvxopt.
n training samples in takes +1 or -1
m test samples in \R^d.
Limits:
n<200, d<100, m<1000
Returns:
y_pred_test : A numpy array with shape (m,). This is the result of running the learned model on the
test instances x_test. Each element is +1 or -1.
"""
n, d = x_train.shape
c = 10 # let max \lambda value be 10
y_train = np.expand_dims(y_train, axis= 1 )* 1.
P = (y_train * x_train).dot((y_train * x_train).T)
q = -1. *np.ones((n, 1 ))
G = np.vstack((np.eye(n)* -1 ,np.eye(n)))
h = np.hstack((np.zeros(n), np.ones(n) * c))
A = y_train.reshape( 1 , -1 )
b = np.array([ 0.0 ])
P = matrix(P); q = matrix(q)
G = matrix(G); h = matrix(h)
A = matrix(A); b = matrix(b)
sol = solvers.qp(P, q, G, h, A, b)
lambdas = np.array(sol[ 'x' ])
w = ((y_train * lambdas).T.dot(x_train)).reshape( -1 , 1 )
b = y_train - np.dot(x_train, w)
prediction = lambda x, w, b: np.sign(np.sum(w.T.dot(x) + b))
y_test = np.array([prediction(x_, w, b) for x_ in x_test])
return y_test
if __name__ == "__main__" :
# Example format of input to the functions
n= 100
m= 100
d= 10
x_train = np.random.randn(n,d)
x_test = np.random.randn(m,d)
w = np.random.randn(d)
y_train = np.sign(np.dot(x_train, w))
y_test = np.sign(np.dot(x_test, w))
y_pred_test = solve_SVM_dual_CVXOPT(x_train, y_train, x_test)
check1 = np.sum(y_pred_test == y_test)
print ( "Score: {}/{}" .format(check1, len(y_Test)))
Для более глубокого понимания см .:
● Численная оптимизация (Хорхе Нокедал и Стивен Дж. Райт)
● Нелинейное программирование (Ю. Нестров)
Полный код, использованный в этом посте, можно найти
(https://github.com/koriavinash1/Numerical-Optimization)
