Введение в машинное обучение.
Поиск "оптимального" решения путем минимизации функционала.
Это предварительная версия! Любые замечания приветствуются.
import numpy as np
import matplotlib.pyplot as plt
import numpy.polynomial.polynomial as poly
%matplotlib inline
Пусть ${\{ (x_i,y_i) \}}_{i=1}^{i=N}$, где N количество данных. При решении задачи регресии (например, линейной) ищется такая функция $f_{a,b}(x)=ax+b$ дабы минимизировать среднеквадротичное отклонение от исходного значения, а именно -- решается задача: $$ \min_{a,b}\sum_{i=1}^{i=N}(f_{a,b}(x_i) - y_i )^2. $$ Точнее ищется именно параметр ($a,b$) для которого данный минимум достигается: $$ a^*, b^* = \underset{a,b}{\operatorname{argmin}}\sum_{i=1}^{i=N}(f_{a,b}(x_i) - y_i )^2. $$ Последнее показывает важность решения задачи минимизации.
Следует отметить, что в простых случаях, например, как ранее задачу можно решить явно. Но бывают случаи, когда этого сделать нельзя.
Поэтому, ставится задача о приближенном поиске минимума. Подразумевается итеративное решение, т.е. когда приближение к искомому минимумому достигается за некотрое (может и сотни, тысячи) количество шагов. Так, решение начинается с некоторго $\theta_0$. Далее, $$ \theta_{i+1} = G(X, \theta_i). $$ Обновление $\theta_i$ происходит пока не будет достигнута нужное приближение к минимуму.
import scipy.optimize as opt
Парабола
Парабола нам хорошо известна. Например, то что экстремум (минимум/максимум) достигается в середине между корнями. Поэтому начнем с неё.
def fun(x):
return (x-5)*(x+11) # x*x + 6x -55
# Корни: 5 и -11. Значит миниум в -3.
def fun(x):
print(x)#
return (x-5)*(x+11) # x*x + 6x -55
# Корни: 5 и -11. Значит миниум в -3.
def dfun(x):
return (x+11) + (x-5)
sol = opt.minimize( fun, 0, method="CG", jac = dfun )
sol
[0.] [-1.01] [-3.]
fun: -64.0 jac: array([-1.59872116e-13]) message: 'Optimization terminated successfully.' nfev: 3 nit: 1 njev: 3 status: 0 success: True x: array([-3.])
sol = opt.minimize( fun, 0, method="CG" )
sol # solution, решение
[0.] [1.49011612e-08] [-1.01] [-1.00999999] [-3.0000002] [-3.00000018]
fun: -63.999999999999964 jac: array([0.]) message: 'Optimization terminated successfully.' nfev: 6 nit: 1 njev: 3 status: 0 success: True x: array([-3.0000002])
# Значение аргумента при котором "достигается" (локальный) минимум.
sol.x # Приближенное решение.
array([-3.0000002])
x = np.arange(-13, 7, 0.1)
plt.plot(x, fun(x), sol.x, fun(sol.x), 'r*')
# Отрисуем параболу с найденным минимумом.
[-1.30000000e+01 -1.29000000e+01 -1.28000000e+01 -1.27000000e+01 -1.26000000e+01 -1.25000000e+01 -1.24000000e+01 -1.23000000e+01 -1.22000000e+01 -1.21000000e+01 -1.20000000e+01 -1.19000000e+01 -1.18000000e+01 -1.17000000e+01 -1.16000000e+01 -1.15000000e+01 -1.14000000e+01 -1.13000000e+01 -1.12000000e+01 -1.11000000e+01 -1.10000000e+01 -1.09000000e+01 -1.08000000e+01 -1.07000000e+01 -1.06000000e+01 -1.05000000e+01 -1.04000000e+01 -1.03000000e+01 -1.02000000e+01 -1.01000000e+01 -1.00000000e+01 -9.90000000e+00 -9.80000000e+00 -9.70000000e+00 -9.60000000e+00 -9.50000000e+00 -9.40000000e+00 -9.30000000e+00 -9.20000000e+00 -9.10000000e+00 -9.00000000e+00 -8.90000000e+00 -8.80000000e+00 -8.70000000e+00 -8.60000000e+00 -8.50000000e+00 -8.40000000e+00 -8.30000000e+00 -8.20000000e+00 -8.10000000e+00 -8.00000000e+00 -7.90000000e+00 -7.80000000e+00 -7.70000000e+00 -7.60000000e+00 -7.50000000e+00 -7.40000000e+00 -7.30000000e+00 -7.20000000e+00 -7.10000000e+00 -7.00000000e+00 -6.90000000e+00 -6.80000000e+00 -6.70000000e+00 -6.60000000e+00 -6.50000000e+00 -6.40000000e+00 -6.30000000e+00 -6.20000000e+00 -6.10000000e+00 -6.00000000e+00 -5.90000000e+00 -5.80000000e+00 -5.70000000e+00 -5.60000000e+00 -5.50000000e+00 -5.40000000e+00 -5.30000000e+00 -5.20000000e+00 -5.10000000e+00 -5.00000000e+00 -4.90000000e+00 -4.80000000e+00 -4.70000000e+00 -4.60000000e+00 -4.50000000e+00 -4.40000000e+00 -4.30000000e+00 -4.20000000e+00 -4.10000000e+00 -4.00000000e+00 -3.90000000e+00 -3.80000000e+00 -3.70000000e+00 -3.60000000e+00 -3.50000000e+00 -3.40000000e+00 -3.30000000e+00 -3.20000000e+00 -3.10000000e+00 -3.00000000e+00 -2.90000000e+00 -2.80000000e+00 -2.70000000e+00 -2.60000000e+00 -2.50000000e+00 -2.40000000e+00 -2.30000000e+00 -2.20000000e+00 -2.10000000e+00 -2.00000000e+00 -1.90000000e+00 -1.80000000e+00 -1.70000000e+00 -1.60000000e+00 -1.50000000e+00 -1.40000000e+00 -1.30000000e+00 -1.20000000e+00 -1.10000000e+00 -1.00000000e+00 -9.00000000e-01 -8.00000000e-01 -7.00000000e-01 -6.00000000e-01 -5.00000000e-01 -4.00000000e-01 -3.00000000e-01 -2.00000000e-01 -1.00000000e-01 -4.61852778e-14 1.00000000e-01 2.00000000e-01 3.00000000e-01 4.00000000e-01 5.00000000e-01 6.00000000e-01 7.00000000e-01 8.00000000e-01 9.00000000e-01 1.00000000e+00 1.10000000e+00 1.20000000e+00 1.30000000e+00 1.40000000e+00 1.50000000e+00 1.60000000e+00 1.70000000e+00 1.80000000e+00 1.90000000e+00 2.00000000e+00 2.10000000e+00 2.20000000e+00 2.30000000e+00 2.40000000e+00 2.50000000e+00 2.60000000e+00 2.70000000e+00 2.80000000e+00 2.90000000e+00 3.00000000e+00 3.10000000e+00 3.20000000e+00 3.30000000e+00 3.40000000e+00 3.50000000e+00 3.60000000e+00 3.70000000e+00 3.80000000e+00 3.90000000e+00 4.00000000e+00 4.10000000e+00 4.20000000e+00 4.30000000e+00 4.40000000e+00 4.50000000e+00 4.60000000e+00 4.70000000e+00 4.80000000e+00 4.90000000e+00 5.00000000e+00 5.10000000e+00 5.20000000e+00 5.30000000e+00 5.40000000e+00 5.50000000e+00 5.60000000e+00 5.70000000e+00 5.80000000e+00 5.90000000e+00 6.00000000e+00 6.10000000e+00 6.20000000e+00 6.30000000e+00 6.40000000e+00 6.50000000e+00 6.60000000e+00 6.70000000e+00 6.80000000e+00 6.90000000e+00] [-3.0000002]
[<matplotlib.lines.Line2D at 0x7f45ae742a90>, <matplotlib.lines.Line2D at 0x7f45ae742ac8>]
# Для отображения информации связанной со сходимостью нужно указать дополнительный параметр
sol = opt.minimize( fun, 0, method="CG", options = {'disp':True} )
# Будет напечатана дополнительная информация.
[0.] [1.49011612e-08] [-1.01] [-1.00999999] [-3.0000002] [-3.00000018] Optimization terminated successfully. Current function value: -64.000000 Iterations: 1 Function evaluations: 6 Gradient evaluations: 3
Iterations - Сколько в итоге было итераций. Current function value - Значение функции в последней точке. С точки зрения производительности важны: Function evaluations - Сколько раз было вычисленно значение функции. Gradient evaluations - Сколько раз был вычислен градиаент.
Наша функция очень простая, поэтому и минимум был найден за один шаг. Усложним функцию.
Полином четвертой степени
def fun2(x):
return (x-5)*(x+11)*(x+5)*(x+3) # x*x + 6x -55
x = np.arange(-13, 7, 0.1)
plt.plot(x, fun2(x))
[<matplotlib.lines.Line2D at 0x7f45a6636f98>]
sol = opt.minimize( fun2, 0, method="CG", options = {'disp':True} )
Optimization terminated successfully. Current function value: -1392.564272 Iterations: 5 Function evaluations: 26 Gradient evaluations: 13
x = np.arange(-13, 7, 0.1)
plt.plot(x, fun2(x))
plt.plot(sol.x, fun2(sol.x), 'r*')
[<matplotlib.lines.Line2D at 0x7f45a6603358>]
Теперь мы видим, что было 5 шагов. И вычисления подросли.
Как вывести промежуточные шаги. Это может быть интересно с точки зрения понимания что происходит при поиске минимума. И когда минимум будет плохо искаться увидеть наглядно в чем затык.
def printx(x):
print(x)
opt.minimize( fun2, 0, method="CG", callback = printx );
[2.23596463] [2.43404495] [2.44736661] [2.44954084] [2.44955921]
steps = []
def savex(x):
steps.append(x[0])
sol = opt.minimize( fun2, 0, method="CG", callback = savex )
steps = np.array(steps)
labels = list(map(str,range(1,len(steps)+1)))
labels
['1', '2', '3', '4', '5']
pnts = list(zip(steps,fun2(steps)))
pnts # Вычислим значение в точках итерации.
[(2.23596463167933, -1386.0928830542748), (2.434044953172843, -1392.5290073297645), (2.4473666094713526, -1392.563566390497), (2.449540835296399, -1392.5642722397167), (2.449559214909005, -1392.5642722893724)]
lpnts = list(zip(labels, pnts))
lpnts # Объединим вместе.
[('1', (2.23596463167933, -1386.0928830542748)), ('2', (2.434044953172843, -1392.5290073297645)), ('3', (2.4473666094713526, -1392.563566390497)), ('4', (2.449540835296399, -1392.5642722397167)), ('5', (2.449559214909005, -1392.5642722893724))]
dpnts = np.array( pnts )
dpnts.T
array([[ 2.23596463, 2.43404495, 2.44736661, 2.44954084, 2.44955921], [-1386.09288305, -1392.52900733, -1392.56356639, -1392.56427224, -1392.56427229]])
# Два аргумента.
np.power(2,3)
8
arg = [2,3]
np.power( *arg ) # Распаковка аргумента.
8
x = np.arange(2.2, 2.6, 0.02)
plt.plot(x, fun2(x))
plt.plot( *dpnts.T, '*' ) # Распакова отдельных точек.
#plt.plot( dpnts[:,0], dpnts[:,1], '*' ) # Так без распаковки.
p=list(map( lambda d: plt.annotate(d[0], d[1]), lpnts[0:3] ));
Не только полнином!
def fun3(x):
return (x-2)**4 + (x+2)**2 + x + np.sin(x*3)*3 # x*x + 6x -55
x = np.arange(0.5, 2, 0.05)
plt.plot(x, fun3(x))
[<matplotlib.lines.Line2D at 0x7f45a6508550>]
steps=[]
sol = opt.minimize( fun3, 0, method="CG", callback = savex )
steps = np.array(steps)
labels = list(map(str,range(1,len(steps)+1)))
pnts = list(zip(steps,fun3(steps))) # fun3!
lpnts = list(zip(labels, pnts))
dpnts = np.array( pnts )
x = np.arange(0.98, 1.45, 0.02)
plt.plot(x, fun3(x))
plt.plot( *dpnts.T, '*' ) # Распакова отдельных точек.
#plt.plot( dpnts[:,0], dpnts[:,1], '*' )
p=list(map( lambda d: plt.annotate(d[0], d[1]), lpnts[0:3] ));
Увеличим масштаб.
x = np.arange(1.29, 1.32, 0.002)
plt.plot(x, fun3(x))
plt.plot( *dpnts[1:].T, '*' ) # Распакова отдельных точек.
#plt.plot( dpnts[:,0], dpnts[:,1], '*' )
p=list(map( lambda d: plt.annotate(d[0], d[1]), lpnts[1:] ));
График похож, но на самом деле это увеличение возле экстремума. Оно должно быть "впукло", т.е. яма. А ямы похожи друг на друга. Смотрим на ост X.
def dfun(x):
return 2*x+6
Поиск нулей
Будем искать минимум там, где производная равна нулю.
# Там где 0
opt.newton(fun, 0, dfun)
0.0 9.166666666666666 5.713470319634704 5.029209940375455 5.000053132289669 5.000000000176438
5.0
По аналогии с одномерным случаем покажем как отрисовывать функции от двух переменных, т.е. 3-хмерные графики.
# Нужно только для того, чтобы 3d графика заработала.
from mpl_toolkits.mplot3d import Axes3D
# т.е. да, сам Axes3D не будет использован.
Кривая
Полным аналогом plot в 3трехмерном пространстве является тоже функция plot. Она отвечает за отрисовку одномерной кривой. Так, по аналогии с двумерным графиком функции подаются точки в прострастве которые нужно соеденить подряд. Токчи задаются тремя списками.
plt.gca(projection='3d') # Вклюаем 3хмерный режим
theta = np.linspace(-8 * np.pi, 8 * np.pi, 200)
z = np.linspace(-2, 2, 200)
r = z**2 + 1
x = r * np.sin(theta)
y = r * np.cos(theta)
plt.plot(x, y, z, label='parametric curve')
plt.legend();
Трехмарность обеспечена вызовом plt.gca(projection='3d'). Функции plot были переданы три списка: x, y, z.
plt.gca(projection='3d') # Вклюаем 3хмерный режим
plt.plot(x, y, z, '*--',label='parametric curve')
plt.legend();# Стиль отрисовки кривой задается таким же образом как и для двумерных графиков.
x1 = z*z/2 + 2 *z - 5
y1 = z - 1
# Можно и две кривые на одном графике.
plt.gca(projection='3d') # Вклюаем 3хмерный режим
plt.plot(x, y, z, '*--',label='parametric curve')
plt.plot(x1, y1, z, '*--',label='parametric curve2')
plt.legend();# Стиль отрисовки кривой задается таким же образом как и для двумерных графиков.
Поверхность
Для создание поверхности нам нужно задать сетку и значение функции в каждом узле данной сетки. Сетка задается прямоугольником: отсечками вдоль одной и дрйгой оси (X и Y).
# Сетки принято делать с равномерным шагом и совпадающим шагов вдоль каждой из осей.
X = [5, 7, 9, 11] # Отсечки вдоль оси X.
Y = [3, 5, 7] # Отсечки вдоль оси y.
Сетку можно визуализировать так:
list( map( lambda y:plt.plot( X, [y]*len(X) ), Y ));
list( map( lambda x:plt.plot( [x]*len(Y), Y ), X ));
plt.xticks(X)
plt.yticks(Y)
plt.axis('equal')
(4.7, 11.3, 2.8, 7.2)
Теперь необходимо получить координаты пересечения.
В качестве упражнения получити их самостоятельно через map или цикл for.
# Есть готовая функция:
grid_x, grid_y = np.meshgrid(X, Y)
grid_x # Значения x-ой координаты точек сетки.
# Видно, что размер матрицы соответствует ранее нарисованной сетки.
array([[ 5, 7, 9, 11], [ 5, 7, 9, 11], [ 5, 7, 9, 11]])
grid_y # Аналогично с координатой y.
array([[3, 3, 3, 3], [5, 5, 5, 5], [7, 7, 7, 7]])
i, j = 1, 2 # Индексы как в матрице.
# Иначе говоря координаты узла с индексом i,j задается как
grid_x[i][j], grid_y[i][j]
# Отмечу, что индекс i является строчкой матрицы,
# а индекс j -- столбцом матрцы.
# Координаты же задаются как обычно: вдоль осей x и y.
(9, 5)
list( map( lambda y:plt.plot( X, [y]*len(X) ), Y ));
list( map( lambda x:plt.plot( [x]*len(Y), Y ), X ));
plt.xticks(X)
plt.yticks(Y)
plt.plot( grid_x[i][j], grid_y[i][j], 'o' ) # Отметим данный узел.
plt.axis('equal')
(4.7, 11.3, 2.8, 7.2)
Z = np.abs(grid_x-8) + np.abs(grid_y-5)
Z
array([[5, 3, 3, 5], [3, 1, 1, 3], [5, 3, 3, 5]])
Z.shape # По размеру матрицы зачений видно что она соответствует сетки.
(3, 4)
ax = plt.gca(projection='3d') # Вклюаем 3хмерный режим
# Отрисовываем поверхность.
surf = ax.plot_surface(grid_x, grid_y, Z ) #, cmap=cm.coolwarm, linewidth=0, antialiased=False)
from matplotlib import cm
ax = plt.gca(projection='3d') # Вклюаем 3хмерный режим
X = np.arange(-5, 5, 0.25)
Y = np.arange(-5, 5, 0.25)
X, Y = np.meshgrid(X, Y)
R = np.sin(X/4)*np.cos(Y)
Z = np.sin(R)
# Отрисовываем поверхность.
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
ax = plt.gca(projection='3d')
X = np.arange(-5, 5, 0.25)
Y = np.arange(-5, 5, 0.25)
X, Y = np.meshgrid(X, Y)
R = np.sqrt(X**2 + Y**2)
Z = np.sin(R)
# Отрисовываем поверхность.
surf = ax.plot_surface(X, Y, Z)
def fun2(x):
return (x[0]-5)**2*(x[1]-7)**2
#def funff(x):
# return ( (x[0]*x[2]-3)**2 + (x[0]*x[3]-4)**2 + (x[1]*x[2]-6)**2 + (x[1]*x[3] - 8)**2 )
grid = np.stack( [X, Y])
grid.shape
(2, 40, 40)
Z = fun2( grid )
Z.shape
(40, 40)
Z[3,3]
10829.00390625
X[3,3], Y[3,3]
(-4.25, -4.25)
fun2( [4.75, 6.75] )
0.00390625
ax = plt.gca(projection='3d')
X = np.arange(4, 6, 0.25)
Y = np.arange(6, 8, 0.25)
X, Y = np.meshgrid(X, Y)
grid = np.stack( [X, Y])
Z = fun2( grid )
surf = ax.plot_surface(X, Y, Z)