Заметка 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) 

opt.minimize( fun2, (0,0) )

      fun: 6.6642318134004586e-15
 hess_inv: array([[ 0.00431369, -0.00612637],
       [-0.00612637,  0.01655506]])
      jac: array([-3.06893558e-08, -1.23806936e-15])
  message: 'Optimization terminated successfully.'
     nfev: 33
      nit: 7
     njev: 11
   status: 0
  success: True
        x: array([ 4.99999999, -3.76552244])

#ans = opt.minimize( funff, (2,1,1,-1) )
#ans

#ans.x[2]/ans.x[3]

def fxy(x):
    return (2-x[0])**2 + 100*(x[1]-x[0]**2)**2

ll=[]
for i in range(-10, 10, 1):
    for j in range(-10, 10, 1):
        res = opt.minimize( fxy, (i,j), method="Powell");
        cnt = res.nit
        if cnt == 0:
            cnt = 17
        if cnt < 15:
            ll.append( (cnt, i, j) )

np.min(ll), np.max(ll)

(-9, 14)

ll;

opt.minimize( fxy, (1,4), method="Powell", options = {'disp':True} )

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 2
         Function evaluations: 39

   direc: array([[1., 0.],
       [0., 1.]])
     fun: 0.0
 message: 'Optimization terminated successfully.'
    nfev: 39
     nit: 2
  status: 0
 success: True
       x: array([2., 4.])

opt.minimize( fxy, (1,4) )

      fun: 2.19088020981211e-10
 hess_inv: array([[0.51006097, 2.03921704],
       [2.03921704, 8.15783707]])
      jac: array([ 9.76923749e-06, -2.38829632e-06])
  message: 'Optimization terminated successfully.'
     nfev: 57
      nit: 13
     njev: 19
   status: 0
  success: True
        x: array([1.9999852 , 3.99994078])

Регуляризация

Сформируем данные

x = np.linspace(10, 20, 5)

y = x * 2 - 15

yy = y + np.random.randn(5) # Добавим шум к каждому значению (измерению)

plt.plot( x, yy, '*' ) # Кривая как-то странно выглядит.

[<matplotlib.lines.Line2D at 0x7f45a63aaf28>]

# -5 + 5x
d = np.array( [-5, 5] ) # Параметры задающие линию.

poly.polyval( 0, d ), poly.polyval( 2, d )

(-5.0, 5.0)

# Значение в точки с индексом 1.
poly.polyval( x[1], d )

57.5

# Значение во всех точках.
poly.polyval( x, d )

array([45. , 57.5, 70. , 82.5, 95. ])

Допустим у нас есть кандидат на линейное уравнение. Как вычислить среднеквадратичную ошибку?

# Квадрат отклонения от правильного ответа в каждой точки.
err = (poly.polyval( x, d ) - y)**2

# Среднеквадратичное отклонение:
np.mean( err )

3137.5

# По параметрам прямой (d) и точкам (x,y) вычислим общую ошибку.
# Отмечу, что точки (x,y) в функцию не передаются. Они глобальны по отношению к ней.
def errLine(d):
    err = (poly.polyval( x, d ) - y)**2 # Ошибка на прямой d.
    return np.mean( err )

opt.minimize( errLine, (-5,4), method="CG" )

     fun: 2.1911050158702597e-13
     jac: array([1.57485331e-09, 2.95533678e-08])
 message: 'Optimization terminated successfully.'
    nfev: 24
     nit: 3
    njev: 8
  status: 0
 success: True
       x: array([-14.99999802,   1.99999987])

Теперь можно использовать любую ошибку.

Упр. Сведите все вместе. Применb идеалогию машинного обучения: а) случайное разбиение данных на "обучающие" и "проверочные", б) подбор параметров функции по "обучающим" данным, в) оценка делается по тестовым даным, г) интегральная оценка вычисляется как средняя ошибка по случайным разбиениям. Выбирите несколько функций приближения (например, приближение линейной, квадратичной, кубической, но можно и гораздо изошренее) и посчитайте среднею ошибку на случайных выборках. Та функция, на которой наилучший результат и считается искомой с точки зрнеия машинного обучения.

Функтор

Как мы ранее видели, функция котрую мы минимизировали использовала глобальные переменные, т.е. переменные которые ей не передавались как аргумент. Такой подход негодится в большом объеме кода.

Как сделать иначе? Как сделать так, чтобы к функции были "привязаны" данные?

# Делается это посредстовм класса у которого есть специальные метод __call__.
# Последний и имеет полный аналог функции, но раз это метод,
# то данная "функция" может хранить данные в классе.

class myCallable:
    def __init__(self, x, y):
        self.x = x # Сохраняем для дальнейшего использования.
        self.y = y
        
    def __call__(self, d):
        err = (poly.polyval( self.x, d ) - self.y)**2 # Ошибка на прямой d.
        return np.mean( err )

my = myCallable(x, y)

# Теперь можно к my относится как к функции
my( d )

3137.5

f'На коэффициентах {d} достигается среднеквадратичная ошибка {my( d )}'

'На коэффициентах [-5  5] достигается среднеквадратичная ошибка 3137.5'

opt.minimize( my, (-5,4), method="CG" )

     fun: 2.1911050158702597e-13
     jac: array([1.57485331e-09, 2.95533678e-08])
 message: 'Optimization terminated successfully.'
    nfev: 24
     nit: 3
    njev: 8
  status: 0
 success: True
       x: array([-14.99999802,   1.99999987])

+регулиризация... вывод формул

  File "<ipython-input-79-25784b6616f2>", line 1
    +регулиризация... вывод формул
                    ^
SyntaxError: invalid syntax

Поиск минимума

Линейный способ

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 0x7f45a602d5f8>]

Совсем ручной

y = []
for i in x:
    y.append([i,fun2(i)])

j=0
for ii, i in enumerate(y):
    if y[j][1] > i[1]:
        j = ii

j, y[j]

(154, [2.3999999999999453, -1392.2063999999993])

xfmin = y[j]
xfmin

[2.3999999999999453, -1392.2063999999993]

x = np.arange(-13, 7, 0.1)
plt.plot(x, fun2(x))
plt.plot( xfmin[0], xfmin[1], "r.")

[<matplotlib.lines.Line2D at 0x7f45a5f22c50>]

Ускоряемся

y = fun2(x)
y.shape

(200,)

ymin = np.min(y)
ymin

-1392.2063999999993

y[y<np.min(y)+10]

array([-1383.7824, -1389.3579, -1392.2064, -1392.1875, -1389.1584,
       -1382.9739])

y[y<np.min(y)+5]

array([-1389.3579, -1392.2064, -1392.1875, -1389.1584])

xx = np.mean(x[y<np.min(y)+5])
xx

2.449999999999945

x = np.arange(-13, 7, 0.1)
plt.plot(x, fun2(x))
plt.plot( xx, fun2(xx), "r.")

[<matplotlib.lines.Line2D at 0x7f45a5ef7cf8>]

Градиентный</3>

def fun2(x):
    return (x-5)*(x+11)*(x+5)*(x+3)

def dfun2(x):
    return (x+11)*(x+5)*(x+3) + (x-5)*( (x+5)*(x+3) + (x+11)*((x+3) + (x+5)) )

x = np.arange(-13, 7, 0.1)
plt.plot(x, fun2(x))
plt.plot(x, dfun2(x))

[<matplotlib.lines.Line2D at 0x7f45a5e5ccf8>]

xi = -2.5
j = 0
while j < 100:
    xi -= dfun2(xi)
    j += 1
xi

nan

xi = -2.5
j = 0
while j < 10000:
    xi -= 0.01*dfun2(xi)
    j += 1
xi

0.9291440685528078

xi = 2.
j = 0
t = []
while j < 100:
    xi -= 0.01*dfun2(xi)#0.01 0.001
    j += 1
    t.append(xi)
xi
t = np.array(t)

plt.plot(t[:30])

[<matplotlib.lines.Line2D at 0x7f45a5ddda20>]

plt.plot(t[80:])

[<matplotlib.lines.Line2D at 0x7f45a5dafe80>]

Достигли цикла.

уменьшим шаг.

xi = 2.
j = 0
t = []
while j < 100:
    xi -= 0.001*dfun2(xi)#0.01 0.001
    j += 1
    t.append(xi)
t = np.array(t)
xi

2.4495592217553024

plt.plot(t)

[<matplotlib.lines.Line2D at 0x7f45a63b5438>]

plt.plot(abs(ymin-fun2(t[:20])))

[<matplotlib.lines.Line2D at 0x7f45a5e0a9b0>]

plt.plot(abs(ymin-fun2(t[20:])))

[<matplotlib.lines.Line2D at 0x7f45a6563b00>]

plt.plot(abs(ymin-fun2(t[40:])))

[<matplotlib.lines.Line2D at 0x7f45a64ad588>]

plt.plot(abs(ymin-fun2(t[60:])))

[<matplotlib.lines.Line2D at 0x7f45a6163208>]

plt.plot(np.log(abs(ymin-fun2(t))))

[<matplotlib.lines.Line2D at 0x7f45a6072be0>]

2D

Тестовое множество

Test

x = np.random.uniform( 3, 10, 10 ) # 10 точек на отрезке [3, 10]
y = x * 2 - 15 # Вычисляем значение в каждой из точек. Функция линейная.

xx = np.sort( x ) # Дело в том, что необходимо отсортировать числа по оси x.
yy = xx * 2 - 15 + np.random.randn(10) # Теперь заново вычисляем.

yy.shape, xx.shape # Массивы соответсвуют друг другу по размеру.

((10,), (10,))

%matplotlib inline

plt.plot( x, y, '-' )
plt.plot( xx, yy, '*' )

[<matplotlib.lines.Line2D at 0x7f45a5d701d0>]

Случайная выборка

perm = np.random.permutation( xx.shape[0] ) # Создаем перестановку числе от 0 до xx.shape[0] не включительно.
perm

array([8, 7, 0, 1, 9, 4, 6, 3, 5, 2])

#np.random.randint(0, xx.shape[0], 5)
ii = perm[:5] # Берем первые 5 чисел. # learn
ii

array([8, 7, 0, 1, 9])

нужно доделать

# Вычисляем параметры линейной регрессии для подмножества числе (ii)
d1 = fitLin( xx[ii], yy[ii])
d1

array([[-14.69147132],
       [  2.05556844]])

# Вычисляем ошибку на подмножестве.
ee = errLine(d1, xx[ii], yy[ii])
ee

(array([5.86197757e-15]), array([0.10379148]))

jj=perm[5:] # test
jj

array([2, 0, 5, 6, 3])

errLine(d1, xx[jj], yy[jj]) # Вычисляем ошибку на всем множестве точек.

(array([1.16575575]), array([2.27008794]))

def experimentLin(n, d, p, k ):
    x = np.random.uniform( 3, 10, (n) )
    y = d[1] * x + d[0] + np.random.randn( x.shape[0] )/2
    d0 = fitLin( x, y)
    e0 = errLine(d0, x, y)
    
    ee10 = np.array([0.])
    ee11 = np.array([0.])
    for j in range(k):
        perm = np.random.permutation( x.shape[0] )
        ii = perm[:int(x.shape[0] * p)]#np.random.randint(0, x.shape[0], int(x.shape[0] * p) )
        #print(ii)
        d1 = fitLin( x[ii], y[ii])
        e1 = errLine(d1, x, y)
        ee10 += e1[0]
        ee11 += e1[1]
    print( e0, [ee10 / k, ee11 / k] )

experimentLin(10, np.array([-15, 2.5]), 0.85, 1000)

(array([2.97539771e-15]), array([0.23625543])) [array([0.00348373]), array([0.25332236])]