import numpy as np
%matplotlib inline
import pylab as plt
x = np.array([[1, 2, 3, 4],
[1, 3, 3, 5]])
w = np.array([0.1, 0.2, 0.3, 0.4])
y = np.array([0, 1])
Как посчитать квадратичный функционал если X — матрица объектов, y — ответы, а w — некоторый вектор весов:
q = 0
for i in range(x.shape[0]):
q += (np.sum(w * x[i]) - y[i]) ** 2
print q
15.76
q = np.sum((np.sum(w * x, axis=1) - y) ** 2)
print q
15.76
Градиентный спуск позволяет оптимизировать (обычно — минимизировать) некоторый функционал.
В начале просто попробуем найти минимум некоторой функции.
Предположим что у нас есть зависимость $y = 2 * x_1^2 + 3 * x_2^2$. Давайте найдем минимум этой функции (он просто вычисляется аналитически, однако попробуем сделать это с помощью данного метода). То есть для заданной функции мы хотим найти точку $(x_1^*, x_2^*)$, в которой достигается минимум.
Для начала опишем функции, которые нам понадобятся:
def function(x1, x2):
return 2 * x1 ** 2 + 3 * x2 ** 2
def grad(x1, x2):
return np.array([4 * x1, 6 * x2])
Тогда градиентный спуск будет выглядеть как функция (уже из питона, которая принимает на вход):
Будем считать что шаг — это константа и не зависит от номера итерации.
И которая будет возвращать:
Так как мы решаем задачу минимизации функции, то хотим, чтобы с номером итерации значение функции уменьшалось.
def grad_descent(x1, x2, step, max_iters):
points = [(x1, x2)]
values = [function(x1, x2)]
for _ in range(max_iters):
x1_new, x2_new = np.array([x1, x2]) - step * grad(x1, x2)
x1, x2 = x1_new, x2_new
points.append((x1, x2))
values.append(function(x1, x2))
return points, values
Запустим наш градиентный спуск из случайной точки и построим график:
points, values = grad_descent(10, 20, 0.01, 100)
print points[-1], values[-1]
plt.plot(values)
(0.16870319358849653, 0.041097495410471981) 0.061988547441
[<matplotlib.lines.Line2D at 0x10c912b50>]
Заметим, что значение градиентный спуск не успевает найти минимум. Возможно стоит увелчить число итераций:
points, values = grad_descent(10, 20, 0.01, 1000)
print points[-1], values[-1]
plt.plot(values)
(1.8673814466701882e-17, 2.6846247849867386e-26) 6.97422693474e-34
[<matplotlib.lines.Line2D at 0x10ca6cd10>]
Или можно было увеличить шаг:
points, values = grad_descent(10, 20, 0.1, 100)
print points[-1], values[-1]
plt.plot(values)
(6.5331862350006843e-22, 3.2138760885179519e-39) 8.53650447624e-43
[<matplotlib.lines.Line2D at 0x10cad1210>]
Если посмотреть на графики, то видно что начиная с некоторой итерациизначение функции не меняется. То есть скорее всего мы уже нашли оптимальную точку. Так что давайте добавим следующий критерий остановки: если предыдущая и следующая точки отличаются не сильно, то градиентный спуск останавливается. Для этого нужно всего лишь дописать пару строчек:
def grad_descent(x1, x2, step, max_iters):
points = [(x1, x2)]
values = [function(x1, x2)]
for iters in range(max_iters):
x1_new, x2_new = np.array([x1, x2]) - step * grad(x1, x2)
if np.linalg.norm(np.array([x1_new, x2_new]) - np.array([x1, x2])) < 1e-5:
break
x1, x2 = x1_new, x2_new
points.append((x1, x2))
values.append(function(x1, x2))
print iters
return points, values
points, values = grad_descent(10, 20, 0.1, 100)
print points[-1], values[-1]
plt.plot(values)
26 (1.7058172817957805e-05, 9.0071992547409624e-10) 5.81962522209e-10
[<matplotlib.lines.Line2D at 0x10cc5d610>]
Теперь градиентный спуск проделал всего 55 итераций.
По факту нахождение оптимального вектора весов (например, в случае линейных алгоритмов), это тоже оптимизационная задача, но уже мы пытаемся оптимизировать функционал качества.
В качестве функционала можно взять упомянутый ранее функционал квадратичных потерь. Пусть мы считаем, что зависимость описывается уравнением $ y = a + b * x_1^2$ и нам необходимо найти коэффициенты a и b. Для этого введем фиктивный признак, везде равный 1. Тогда можно переписать уравнение следующим образом: $ y = a * x_1 + b * x_2^2$.
Так как мы оптимизируем функционал и пытаемся найти коэффициенты a и b, необходимо записать градиент Q по этим переменным.
def function(x1, x2, a, b):
return a * x1 + b * x2 ** 2
def Q(x1, x2, a, b, y):
return np.sum((function(x1, x2, a, b) - y) ** 2)
def grad_Q(x1, x2, a, b, y):
return np.array([np.sum(2 * (function(x1, x2, a, b) - y) * x1),
np.sum(2 * (function(x1, x2, a, b) - y) * x2 ** 2)])
Отлично, теперь нам будет передаваться в градиентный спуск некоторая выборка, а мы будем по ней подбирать коэффициенты a и b. Небольшой модификацией предыдущего кода, получаем код, решающий нашу задачу:
def grad_descent(x1, x2, y, a0, b0, step, max_iters):
points = [(a0, b0)]
values = [Q(x1, x2, a0, b0, y)]
a, b = a0, b0
for iters in range(max_iters):
a_new, b_new = np.array([a, b]) - step * grad_Q(x1, x2, a, b, y)
if np.linalg.norm(np.array([a_new, b_new]) - np.array([a, b])) < 1e-5:
break
a, b = a_new, b_new
points.append((a, b))
values.append(Q(x1, x2, a, b, y))
print iters
return points, values
Только нужно сгенерировать выборку. Предполагаем, что истинная зависимость $ y = 2 + 3 * x_1^2$ (то есь настоящие a и b равны 2 и 3 соответственно, но об этом никто не знает). И добавим немного шума в ответы.
n = 20
x1 = np.ones((1, 20))
x2 = np.linspace(0, 10, 20)
y = function(x1, x2, 2, 3) + np.random.normal(0, 1)
Запустим наш спуск из случайной начальной точки 0, 0.
points, values = grad_descent(x1, x2, y, 0, 0, 0.01, 100)
print points[-1], values[-1]
plt.plot(values)
99 (-1.8552849864902168e+292, -1.1704828397605956e+294) inf
/usr/local/lib/python2.7/site-packages/ipykernel/__main__.py:5: RuntimeWarning: overflow encountered in square
[<matplotlib.lines.Line2D at 0x10725d810>]
Получилось что-то явно не то (значение функционала равно +inf). Сократим шаг.
points, values = grad_descent(x1, x2, y, 0, 0, 0.00001, 100)
print points[-1], values[-1]
plt.plot(values)
99 (0.0799334813069459, 3.0273948768068668) 27.3532007379
[<matplotlib.lines.Line2D at 0x10721ba90>]
Уже лучше, однако кажется что не хватило итераций:
points, values = grad_descent(x1, x2, y, 0, 0, 0.00001, 1000)
print points[-1], values[-1]
plt.plot(values)
999 (0.34253214066478249, 3.023232531716197) 19.6726346491
[<matplotlib.lines.Line2D at 0x107373490>]
points, values = grad_descent(x1, x2, y, 0, 0, 0.00001, 10000)
print points[-1], values[-1]
plt.plot(values)
9999 (1.5262040841834328, 3.0044706267045096) 0.728460931851
[<matplotlib.lines.Line2D at 0x107610090>]
График получился не очень репрезентативен. Давайте посмотрим в последнем случае на каждое сотое значение функционала.
print values[::100]
[323297.57383659371, 27.461800541681825, 26.474322124700553, 25.522351708096949, 24.604612486151403, 23.71987356484934, 22.866948310976838, 22.044692760579082, 21.252004084648654, 20.487819109983999, 19.751112893237053, 19.04089734623286, 18.356219910723159, 17.69616228079153, 17.059839171199489, 16.446397130020777, 15.85501339397095, 15.28489478489689, 14.735276645947607, 14.205421815998115, 13.694619640950711, 13.202185020588491, 12.727457489702996, 12.269800332261717, 11.82859972742928, 11.403263926296411, 10.993222458211823, 10.597925365653026, 10.216842466610231, 9.8494626434933341, 9.4952931576087316, 9.1538589882863466, 8.8247021957700458, 8.5073813070182869, 8.2014707235895727, 7.9065601508194003, 7.6222540475233425, 7.3481710954877109, 7.083943688036296, 6.8292174369869976, 6.5836506973385776, 6.3469141090476526, 6.1186901552829696, 5.8986727365647234, 5.686566760215797, 5.4820877445768792, 5.2849614374524894, 5.0949234482777008, 4.9117188935104528, 4.7351020547762888, 4.5648360493044935, 4.4006925122154428, 4.2424512902319682, 4.089900146404374, 3.9428344754534685, 3.8010570293485793, 3.6643776527542444, 3.5326130279886137, 3.4055864291539821, 3.2831274851073577, 3.1650719509547276, 3.051261487761797, 2.9415434501860527, 2.8357706817449468, 2.7338013174462925, 2.63549859351546, 2.5407306639643266, 2.4493704237565042, 2.3612953383308062, 2.2763872792550481, 2.1945323657892679, 2.1156208121460778, 2.0395467802435103, 1.9662082377521981, 1.8955068212473267, 1.8273477042811495, 1.7616394701994864, 1.6982939895314009, 1.6372263017879596, 1.578354501511243, 1.5215996284204532, 1.4668855615088179, 1.4141389169480298, 1.3632889496642051, 1.3142674584528251, 1.2670086945056902, 1.2214492732269306, 1.1775280892201385, 1.1351862343323449, 1.0943669186448488, 1.055015394305455, 1.017078882099038, 0.98050650065901157, 0.94524919822390074, 0.91125968684798597, 0.87849237897747279, 0.84690332630753606, 0.81645016083761146, 0.78709203804647498, 0.75878958211068481, 0.73150483309262149]
Как видно, фукционал действительно уменьшается, то есть все верно (но видим, что даже 10000 итераций оказалось недостаточно).