Заметка 2. Оптимизационная задача. курса Введение в машинное обучение. Шокуров Антон В. shokurov.anton.v@yandex.ru http://машинноезрение.рф Версия 0.11

Аннотация

Поиск "оптимального" решения путем минимизации функционала.

Это предварительная версия! Любые замечания приветствуются.

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)