Članek obravnava več načinov za določitev matematične enačbe preproste (parne) regresijske črte.
Vse tukaj obravnavane metode reševanja enačbe temeljijo na metodi najmanjših kvadratov. Označimo metode na naslednji način:
- Analitična rešitev
- Gradientni spust
- Stohastični gradientni spust
Za vsako metodo reševanja enačbe premice so v članku podane različne funkcije, ki so v glavnem razdeljene na tiste, ki so napisane brez uporabe knjižnice numpy in tiste, ki se uporabljajo za izračune numpy. Menijo, da je spretna uporaba numpy zmanjša stroške računalništva.
Vsa koda v članku je napisana v jeziku pyton 2.7 z uporabo Jupyter Notebook. Izvorna koda in datoteka z vzorčnimi podatki sta objavljeni na
Članek je bolj namenjen tako začetnikom kot tistim, ki so že postopoma začeli obvladovati študij zelo širokega dela umetne inteligence - strojnega učenja.
Za ponazoritev gradiva uporabljamo zelo preprost primer.
Primeri pogojev
Imamo pet vrednosti, ki označujejo odvisnost Y od X (Tabela št. 1):
Tabela št. 1 "Primeri pogojev"
Predpostavili bomo, da so vrednosti je mesec v letu in — prihodek ta mesec. Z drugimi besedami, prihodki so odvisni od meseca v letu in - edini znak, od katerega je odvisen prihodek.
Primer je tako tako, tako z vidika pogojne odvisnosti prihodkov od meseca v letu kot z vidika števila vrednosti - zelo malo jih je. Vendar pa bo takšna poenostavitev omogočila, kot pravijo, razložiti, ne vedno z lahkoto, gradivo, ki ga začetniki usvojijo. In tudi preprostost številk bo tistim, ki želijo rešiti primer na papirju, omogočila brez večjih stroškov dela.
Predpostavimo, da je odvisnost, podano v primeru, mogoče precej dobro aproksimirati z matematično enačbo preproste (parne) regresijske črte oblike:
če je mesec, v katerem je bil prejet prihodek, — prihodek, ki ustreza mesecu, и so regresijski koeficienti ocenjene črte.
Upoštevajte, da je koeficient pogosto imenovan naklon ali gradient ocenjene črte; predstavlja znesek, za katerega je ko se spremeni .
Očitno je naša naloga v primeru izbrati takšne koeficiente v enačbi и , pri kateri odstopanja naših izračunanih vrednosti prihodkov po mesecih od pravih odgovorov, t.j. vrednosti, predstavljene v vzorcu, bodo minimalne.
Metoda najmanjših kvadratov
Po metodi najmanjših kvadratov naj bi odstopanje izračunali tako, da ga kvadriramo. Ta tehnika vam omogoča, da se izognete medsebojnemu preklicu odstopanj, če imajo nasprotne znake. Na primer, če je v enem primeru odstopanje +5 (plus pet), v drugi pa -5 (minus pet), potem se bo vsota odstopanj med seboj izničila in znašala 0 (nič). Možno je, da odstopanja ne kvadriramo, ampak uporabimo lastnost modula in potem bodo vsa odstopanja pozitivna in se bodo kopičila. Na tej točki se ne bomo podrobneje ukvarjali, ampak preprosto navedli, da je za udobje izračunov običajno kvadrat odstopanja.
Takole izgleda formula, s katero bomo določili najmanjšo vsoto kvadratov odstopanj (napak):
če je funkcija približka resničnih odgovorov (tj. prihodka, ki smo ga izračunali),
so pravi odgovori (prihodki navedeni v vzorcu),
je vzorčni indeks (številka meseca, v katerem je ugotovljeno odstopanje)
Diferencirajmo funkcijo, definirajmo parcialne diferencialne enačbe in bodimo pripravljeni, da nadaljujemo z analitično rešitvijo. Najprej pa naredimo kratek izlet o tem, kaj je diferenciacija, in se spomnimo geometrijskega pomena odvoda.
Diferenciacija
Diferenciranje je operacija iskanja odvoda funkcije.
Za kaj se uporablja derivat? Odvod funkcije označuje hitrost spreminjanja funkcije in nam pove njeno smer. Če je odvod v dani točki pozitiven, potem funkcija narašča, sicer pa pada. In večja kot je vrednost absolutnega odvoda, višja je hitrost spreminjanja vrednosti funkcije, pa tudi bolj strm je naklon grafa funkcije.
Na primer, v pogojih kartezičnega koordinatnega sistema je vrednost odvoda v točki M(0,0) enaka +25 pomeni, da na dani točki, ko je vrednost premaknjena na desno s konvencionalno enoto, vrednost poveča za 25 konvencionalnih enot. Na grafu je videti kot dokaj strm porast vrednosti od dane točke.
Še en primer. Izpeljanka je enaka -0,1 pomeni, da ko je premaknjen na eno konvencionalno enoto, vrednost zmanjša le za 0,1 konvencionalne enote. Hkrati lahko na grafu funkcije opazimo komaj opazen naklon navzdol. Če potegnemo analogijo z goro, je videti, kot da se zelo počasi spuščamo po blagem pobočju z gore, za razliko od prejšnjega primera, kjer smo se morali povzpeti na zelo strme vrhove :)
Tako po diferenciaciji funkcije po kvotah и , definiramo parcialne diferencialne enačbe 1. reda. Po določitvi enačb bomo prejeli sistem dveh enačb, z reševanjem katerih bomo lahko izbrali takšne vrednosti koeficientov и , pri kateri se vrednosti ustreznih derivatov v danih točkah spremenijo za zelo, zelo majhno količino, v primeru analitične rešitve pa se sploh ne spremenijo. Z drugimi besedami, funkcija napake pri najdenih koeficientih bo dosegla minimum, saj bodo vrednosti delnih derivatov na teh točkah enake nič.
Torej, v skladu s pravili diferenciacije, enačba delnega odvoda 1. reda glede na koeficient bo imel obliko:
Enačba parcialnega odvoda 1. reda glede na bo imel obliko:
Kot rezultat smo dobili sistem enačb, ki ima dokaj preprosto analitično rešitev:
začetek{enačba*}
začetek{primeri}
na + bsumlimits_{i=1}^nx_i — sumlimits_{i=1}^ny_i = 0
sumlimits_{i=1}^nx_i(a +bsumlimits_{i=1}^nx_i — sumlimits_{i=1}^ny_i) = 0
konec{primeri}
konec{enačba*}
Preden rešimo enačbo, prednaložimo, preverimo pravilnost nalaganja in oblikujmo podatke.
Nalaganje in oblikovanje podatkov
Opozoriti je treba, da bomo zaradi dejstva, da bomo za analitično rešitev in nato za gradientni in stohastični gradientni spust uporabili kodo v dveh različicah: z uporabo knjižnice numpy in brez uporabe, potem bomo potrebovali ustrezno oblikovanje podatkov (glejte kodo).
Koda za nalaganje in obdelavo podatkov
# импортируем все нужные нам библиотеки
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import math
import pylab as pl
import random
# графики отобразим в Jupyter
%matplotlib inline
# укажем размер графиков
from pylab import rcParams
rcParams['figure.figsize'] = 12, 6
# отключим предупреждения Anaconda
import warnings
warnings.simplefilter('ignore')
# загрузим значения
table_zero = pd.read_csv('data_example.txt', header=0, sep='t')
# посмотрим информацию о таблице и на саму таблицу
print table_zero.info()
print '********************************************'
print table_zero
print '********************************************'
# подготовим данные без использования NumPy
x_us = []
[x_us.append(float(i)) for i in table_zero['x']]
print x_us
print type(x_us)
print '********************************************'
y_us = []
[y_us.append(float(i)) for i in table_zero['y']]
print y_us
print type(y_us)
print '********************************************'
# подготовим данные с использованием NumPy
x_np = table_zero[['x']].values
print x_np
print type(x_np)
print x_np.shape
print '********************************************'
y_np = table_zero[['y']].values
print y_np
print type(y_np)
print y_np.shape
print '********************************************'Vizualizacija
Zdaj, ko smo prvič naložili podatke, drugič preverili pravilnost nalaganja in končno formatirali podatke, bomo izvedli prvo vizualizacijo. Metoda, ki se pogosto uporablja za to, je pairplot knjižnice rojen na morju. V našem primeru uporaba knjižnice zaradi omejenega števila nima smisla rojen na morju. Uporabljali bomo običajno knjižnico matplotlib in samo poglejte razpršeni grafikon.
Koda razpršitve
print 'График №1 "Зависимость выручки от месяца года"'
plt.plot(x_us,y_us,'o',color='green',markersize=16)
plt.xlabel('$Months$', size=16)
plt.ylabel('$Sales$', size=16)
plt.show()Grafikon št. 1 "Odvisnost prihodkov od meseca v letu"
Analitična rešitev
Uporabimo najpogostejša orodja v python in reši sistem enačb:
začetek{enačba*}
začetek{primeri}
na + bsumlimits_{i=1}^nx_i — sumlimits_{i=1}^ny_i = 0
sumlimits_{i=1}^nx_i(a +bsumlimits_{i=1}^nx_i — sumlimits_{i=1}^ny_i) = 0
konec{primeri}
konec{enačba*}
Po Cramerjevem pravilu bomo našli splošno determinanto, pa tudi determinante po in , nato pa deljenje determinante z na splošno determinanto - poiščite koeficient , podobno najdemo koeficient .
Koda analitične rešitve
# определим функцию для расчета коэффициентов a и b по правилу Крамера
def Kramer_method (x,y):
# сумма значений (все месяца)
sx = sum(x)
# сумма истинных ответов (выручка за весь период)
sy = sum(y)
# сумма произведения значений на истинные ответы
list_xy = []
[list_xy.append(x[i]*y[i]) for i in range(len(x))]
sxy = sum(list_xy)
# сумма квадратов значений
list_x_sq = []
[list_x_sq.append(x[i]**2) for i in range(len(x))]
sx_sq = sum(list_x_sq)
# количество значений
n = len(x)
# общий определитель
det = sx_sq*n - sx*sx
# определитель по a
det_a = sx_sq*sy - sx*sxy
# искомый параметр a
a = (det_a / det)
# определитель по b
det_b = sxy*n - sy*sx
# искомый параметр b
b = (det_b / det)
# контрольные значения (прооверка)
check1 = (n*b + a*sx - sy)
check2 = (b*sx + a*sx_sq - sxy)
return [round(a,4), round(b,4)]
# запустим функцию и запишем правильные ответы
ab_us = Kramer_method(x_us,y_us)
a_us = ab_us[0]
b_us = ab_us[1]
print ' 33[1m' + ' 33[4m' + "Оптимальные значения коэффициентов a и b:" + ' 33[0m'
print 'a =', a_us
print 'b =', b_us
print
# определим функцию для подсчета суммы квадратов ошибок
def errors_sq_Kramer_method(answers,x,y):
list_errors_sq = []
for i in range(len(x)):
err = (answers[0] + answers[1]*x[i] - y[i])**2
list_errors_sq.append(err)
return sum(list_errors_sq)
# запустим функцию и запишем значение ошибки
error_sq = errors_sq_Kramer_method(ab_us,x_us,y_us)
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений" + ' 33[0m'
print error_sq
print
# замерим время расчета
# print ' 33[1m' + ' 33[4m' + "Время выполнения расчета суммы квадратов отклонений:" + ' 33[0m'
# % timeit error_sq = errors_sq_Kramer_method(ab,x_us,y_us)Evo, kaj imamo:
Torej so bile ugotovljene vrednosti koeficientov, ugotovljena je bila vsota kvadratov odstopanj. Na histogramu sipanja narišimo ravno črto v skladu z najdenimi koeficienti.
Koda regresijske črte
# определим функцию для формирования массива рассчетных значений выручки
def sales_count(ab,x,y):
line_answers = []
[line_answers.append(ab[0]+ab[1]*x[i]) for i in range(len(x))]
return line_answers
# построим графики
print 'Грфик№2 "Правильные и расчетные ответы"'
plt.plot(x_us,y_us,'o',color='green',markersize=16, label = '$True$ $answers$')
plt.plot(x_us, sales_count(ab_us,x_us,y_us), color='red',lw=4,
label='$Function: a + bx,$ $where$ $a='+str(round(ab_us[0],2))+',$ $b='+str(round(ab_us[1],2))+'$')
plt.xlabel('$Months$', size=16)
plt.ylabel('$Sales$', size=16)
plt.legend(loc=1, prop={'size': 16})
plt.show()Tabela št. 2 "Pravilni in izračunani odgovori"
Za vsak mesec si lahko ogledate graf odstopanj. V našem primeru iz tega ne bomo izpeljali nobene pomembne praktične vrednosti, bomo pa zadovoljili našo radovednost o tem, kako dobro enačba preproste linearne regresije označuje odvisnost prihodkov od meseca v letu.
Koda grafikona odstopanj
# определим функцию для формирования массива отклонений в процентах
def error_per_month(ab,x,y):
sales_c = sales_count(ab,x,y)
errors_percent = []
for i in range(len(x)):
errors_percent.append(100*(sales_c[i]-y[i])/y[i])
return errors_percent
# построим график
print 'График№3 "Отклонения по-месячно, %"'
plt.gca().bar(x_us, error_per_month(ab_us,x_us,y_us), color='brown')
plt.xlabel('Months', size=16)
plt.ylabel('Calculation error, %', size=16)
plt.show()Grafikon št. 3 “Odstopanja, %”
Ni popolno, vendar smo svojo nalogo opravili.
Napišimo funkcijo, ki določa koeficiente и uporablja knjižnico numpy, natančneje, zapisali bomo dve funkciji: eno z uporabo psevdoinverzne matrike (ni priporočljivo v praksi, ker je proces računsko zapleten in nestabilen), drugo z uporabo matrične enačbe.
Koda analitične rešitve (NumPy)
# для начала добавим столбец с не изменяющимся значением в 1.
# Данный столбец нужен для того, чтобы не обрабатывать отдельно коэффицент a
vector_1 = np.ones((x_np.shape[0],1))
x_np = table_zero[['x']].values # на всякий случай приведем в первичный формат вектор x_np
x_np = np.hstack((vector_1,x_np))
# проверим то, что все сделали правильно
print vector_1[0:3]
print x_np[0:3]
print '***************************************'
print
# напишем функцию, которая определяет значения коэффициентов a и b с использованием псевдообратной матрицы
def pseudoinverse_matrix(X, y):
# задаем явный формат матрицы признаков
X = np.matrix(X)
# определяем транспонированную матрицу
XT = X.T
# определяем квадратную матрицу
XTX = XT*X
# определяем псевдообратную матрицу
inv = np.linalg.pinv(XTX)
# задаем явный формат матрицы ответов
y = np.matrix(y)
# находим вектор весов
return (inv*XT)*y
# запустим функцию
ab_np = pseudoinverse_matrix(x_np, y_np)
print ab_np
print '***************************************'
print
# напишем функцию, которая использует для решения матричное уравнение
def matrix_equation(X,y):
a = np.dot(X.T, X)
b = np.dot(X.T, y)
return np.linalg.solve(a, b)
# запустим функцию
ab_np = matrix_equation(x_np,y_np)
print ab_npPrimerjajmo čas, porabljen za določanje koeficientov и , v skladu s 3 predstavljenimi metodami.
Koda za izračun časa izračuna
print ' 33[1m' + ' 33[4m' + "Время выполнения расчета коэффициентов без использования библиотеки NumPy:" + ' 33[0m'
% timeit ab_us = Kramer_method(x_us,y_us)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Время выполнения расчета коэффициентов с использованием псевдообратной матрицы:" + ' 33[0m'
%timeit ab_np = pseudoinverse_matrix(x_np, y_np)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Время выполнения расчета коэффициентов с использованием матричного уравнения:" + ' 33[0m'
%timeit ab_np = matrix_equation(x_np, y_np)
Z majhno količino podatkov se pojavi "samonapisana" funkcija, ki poišče koeficiente po Cramerjevi metodi.
Zdaj lahko nadaljujete z drugimi načini iskanja koeficientov и .
Gradientni spust
Najprej opredelimo, kaj je gradient. Preprosto povedano, gradient je segment, ki označuje smer največje rasti funkcije. Po analogiji s plezanjem na goro je tam, kjer je strm vzpon, tisti, kjer je najstrmejši vzpon na vrh gore. Če razvijemo primer z goro, se spomnimo, da pravzaprav potrebujemo najstrmejši spust, da čim hitreje dosežemo nižino, torej minimum - mesto, kjer se funkcija ne poveča ali zmanjša. Na tej točki bo odvod enak nič. Zato ne potrebujemo gradienta, ampak antigradient. Če želite najti antigradient, morate samo pomnožiti gradient s -1 (minus ena).
Bodimo pozorni na dejstvo, da ima lahko funkcija več minimumov in ko se spustimo v enega od njih z uporabo spodaj predlaganega algoritma, ne bomo mogli najti drugega minimuma, ki je lahko nižji od najdenega. Sprostimo se, to nas ne ogroža! V našem primeru imamo opravka z enim samim minimumom, saj je naša funkcija na grafu je pravilna parabola. In kot bi morali vsi zelo dobro vedeti iz šolskega tečaja matematike, ima parabola samo en minimum.
Ko smo ugotovili, zakaj potrebujemo gradient in tudi, da je gradient segment, to je vektor z danimi koordinatami, ki so popolnoma enaki koeficienti и lahko izvajamo gradientni spust.
Preden začnete, predlagam, da preberete le nekaj stavkov o algoritmu spuščanja:
- Na psevdonaključni način določimo koordinate koeficientov и . V našem primeru bomo definirali koeficiente blizu ničle. To je običajna praksa, vendar ima lahko vsak primer svojo prakso.
- Iz koordinat odštejte vrednost delnega odvoda 1. reda v točki . Torej, če je odvod pozitiven, potem funkcija narašča. Zato se bomo z odštevanjem vrednosti derivata premaknili v nasprotni smeri rasti, torej v smeri padanja. Če je odvod negativen, potem funkcija na tej točki pada in se z odštevanjem vrednosti odvoda premaknemo v smeri padanja.
- Podobno operacijo izvedemo s koordinato : odštejte vrednost delnega odvoda v točki .
- Da ne bi preskočili minimuma in poleteli v globoko vesolje, je treba nastaviti velikost koraka v smeri spuščanja. Na splošno bi lahko napisali cel članek o tem, kako pravilno nastaviti korak in kako ga spremeniti med postopkom spuščanja, da zmanjšate računske stroške. Sedaj pa nas čaka nekoliko drugačna naloga in velikost koraka bomo ugotavljali z znanstveno metodo »poke« oziroma po domače empirično.
- Ko smo iz podanih koordinat и odštejemo vrednosti derivatov, dobimo nove koordinate и . Naredimo naslednji korak (odštevanje), že od izračunanih koordinat. In tako se cikel začne znova in znova, dokler ni dosežena zahtevana konvergenca.
Vse! Zdaj smo pripravljeni na iskanje najgloblje soteske Marianskega jarka. Začnimo.
Koda za gradientni spust
# напишем функцию градиентного спуска без использования библиотеки NumPy.
# Функция на вход принимает диапазоны значений x,y, длину шага (по умолчанию=0,1), допустимую погрешность(tolerance)
def gradient_descent_usual(x_us,y_us,l=0.1,tolerance=0.000000000001):
# сумма значений (все месяца)
sx = sum(x_us)
# сумма истинных ответов (выручка за весь период)
sy = sum(y_us)
# сумма произведения значений на истинные ответы
list_xy = []
[list_xy.append(x_us[i]*y_us[i]) for i in range(len(x_us))]
sxy = sum(list_xy)
# сумма квадратов значений
list_x_sq = []
[list_x_sq.append(x_us[i]**2) for i in range(len(x_us))]
sx_sq = sum(list_x_sq)
# количество значений
num = len(x_us)
# начальные значения коэффициентов, определенные псевдослучайным образом
a = float(random.uniform(-0.5, 0.5))
b = float(random.uniform(-0.5, 0.5))
# создаем массив с ошибками, для старта используем значения 1 и 0
# после завершения спуска стартовые значения удалим
errors = [1,0]
# запускаем цикл спуска
# цикл работает до тех пор, пока отклонение последней ошибки суммы квадратов от предыдущей, не будет меньше tolerance
while abs(errors[-1]-errors[-2]) > tolerance:
a_step = a - l*(num*a + b*sx - sy)/num
b_step = b - l*(a*sx + b*sx_sq - sxy)/num
a = a_step
b = b_step
ab = [a,b]
errors.append(errors_sq_Kramer_method(ab,x_us,y_us))
return (ab),(errors[2:])
# запишем массив значений
list_parametres_gradient_descence = gradient_descent_usual(x_us,y_us,l=0.1,tolerance=0.000000000001)
print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_gradient_descence[0][0],3)
print 'b =', round(list_parametres_gradient_descence[0][1],3)
print
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_gradient_descence[1][-1],3)
print
print ' 33[1m' + ' 33[4m' + "Количество итераций в градиентном спуске:" + ' 33[0m'
print len(list_parametres_gradient_descence[1])
print
Potopili smo se do samega dna Marianskega jarka in tam našli vse enake vrednosti koeficientov и , kar je bilo točno to, kar je bilo pričakovati.
Potopimo se še enkrat, le da bo tokrat naše globokomorsko vozilo polno drugih tehnologij, in sicer knjižnice numpy.
Koda za gradientni spust (NumPy)
# перед тем определить функцию для градиентного спуска с использованием библиотеки NumPy,
# напишем функцию определения суммы квадратов отклонений также с использованием NumPy
def error_square_numpy(ab,x_np,y_np):
y_pred = np.dot(x_np,ab)
error = y_pred - y_np
return sum((error)**2)
# напишем функцию градиентного спуска с использованием библиотеки NumPy.
# Функция на вход принимает диапазоны значений x,y, длину шага (по умолчанию=0,1), допустимую погрешность(tolerance)
def gradient_descent_numpy(x_np,y_np,l=0.1,tolerance=0.000000000001):
# сумма значений (все месяца)
sx = float(sum(x_np[:,1]))
# сумма истинных ответов (выручка за весь период)
sy = float(sum(y_np))
# сумма произведения значений на истинные ответы
sxy = x_np*y_np
sxy = float(sum(sxy[:,1]))
# сумма квадратов значений
sx_sq = float(sum(x_np[:,1]**2))
# количество значений
num = float(x_np.shape[0])
# начальные значения коэффициентов, определенные псевдослучайным образом
a = float(random.uniform(-0.5, 0.5))
b = float(random.uniform(-0.5, 0.5))
# создаем массив с ошибками, для старта используем значения 1 и 0
# после завершения спуска стартовые значения удалим
errors = [1,0]
# запускаем цикл спуска
# цикл работает до тех пор, пока отклонение последней ошибки суммы квадратов от предыдущей, не будет меньше tolerance
while abs(errors[-1]-errors[-2]) > tolerance:
a_step = a - l*(num*a + b*sx - sy)/num
b_step = b - l*(a*sx + b*sx_sq - sxy)/num
a = a_step
b = b_step
ab = np.array([[a],[b]])
errors.append(error_square_numpy(ab,x_np,y_np))
return (ab),(errors[2:])
# запишем массив значений
list_parametres_gradient_descence = gradient_descent_numpy(x_np,y_np,l=0.1,tolerance=0.000000000001)
print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_gradient_descence[0][0],3)
print 'b =', round(list_parametres_gradient_descence[0][1],3)
print
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_gradient_descence[1][-1],3)
print
print ' 33[1m' + ' 33[4m' + "Количество итераций в градиентном спуске:" + ' 33[0m'
print len(list_parametres_gradient_descence[1])
print
Vrednosti koeficientov и nespremenljivo.
Poglejmo, kako se je spreminjala napaka med gradientnim spuščanjem, torej kako se je spreminjala vsota kvadratov odstopanj z vsakim korakom.
Koda za risanje vsot kvadratov odklonov
print 'График№4 "Сумма квадратов отклонений по-шагово"'
plt.plot(range(len(list_parametres_gradient_descence[1])), list_parametres_gradient_descence[1], color='red', lw=3)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()Graf št. 4 “Vsota kvadratov odstopanj med gradientnim spustom”
Na grafu vidimo, da se z vsakim korakom napaka zmanjšuje, po določenem številu iteracij pa opazimo skoraj vodoravno črto.
Nazadnje ocenimo razliko v času izvajanja kode:
Koda za določanje časa izračuna gradientnega spusta
print ' 33[1m' + ' 33[4m' + "Время выполнения градиентного спуска без использования библиотеки NumPy:" + ' 33[0m'
%timeit list_parametres_gradient_descence = gradient_descent_usual(x_us,y_us,l=0.1,tolerance=0.000000000001)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Время выполнения градиентного спуска с использованием библиотеки NumPy:" + ' 33[0m'
%timeit list_parametres_gradient_descence = gradient_descent_numpy(x_np,y_np,l=0.1,tolerance=0.000000000001)
Morda delamo kaj narobe, a spet gre za preprosto "doma napisano" funkcijo, ki ne uporablja knjižnice numpy prekaša čas izračuna funkcije, ki uporablja knjižnico numpy.
Vendar ne stojimo na mestu, ampak gremo k preučevanju še enega vznemirljivega načina za reševanje preproste enačbe linearne regresije. Srečati!
Stohastični gradientni spust
Da bi hitro razumeli princip delovanja stohastičnega gradientnega spuščanja, je bolje ugotoviti njegove razlike od običajnega gradientnega spuščanja. Mi, v primeru gradientnega spusta, v enačbah odvodov od и uporabil vsote vrednosti vseh lastnosti in resničnih odgovorov, ki so na voljo v vzorcu (to je vsoto vseh и ). Pri stohastičnem gradientnem spuščanju ne bomo uporabili vseh vrednosti, ki so prisotne v vzorcu, ampak namesto tega psevdo-naključno izberemo tako imenovani vzorčni indeks in uporabimo njegove vrednosti.
Na primer, če je indeks določen kot številka 3 (tri), potem vzamemo vrednosti и , nato vrednosti nadomestimo v izvedene enačbe in določimo nove koordinate. Potem, ko določimo koordinate, ponovno psevdo-naključno določimo vzorčni indeks, nadomestimo vrednosti, ki ustrezajo indeksu, v parcialne diferencialne enačbe in določimo koordinate na nov način и itd. dokler konvergenca ne postane zelena. Na prvi pogled se morda zdi, da to sploh ne bi moglo delovati, a je. Res je, da velja poudariti, da se napaka ne zmanjšuje z vsakim korakom, a tendenca vsekakor obstaja.
Kakšne so prednosti stohastičnega gradientnega spuščanja pred običajnim? Če je velikost našega vzorca zelo velika in se meri v več deset tisoč vrednostih, potem je veliko lažje obdelati, recimo, naključnih tisoč izmed njih, kot celotnega vzorca. Tu pride v poštev stohastični gradientni spust. V našem primeru seveda ne bomo opazili velike razlike.
Poglejmo kodo.
Koda za stohastični gradientni spust
# определим функцию стох.град.шага
def stoch_grad_step_usual(vector_init, x_us, ind, y_us, l):
# выбираем значение икс, которое соответствует случайному значению параметра ind
# (см.ф-цию stoch_grad_descent_usual)
x = x_us[ind]
# рассчитывыаем значение y (выручку), которая соответствует выбранному значению x
y_pred = vector_init[0] + vector_init[1]*x_us[ind]
# вычисляем ошибку расчетной выручки относительно представленной в выборке
error = y_pred - y_us[ind]
# определяем первую координату градиента ab
grad_a = error
# определяем вторую координату ab
grad_b = x_us[ind]*error
# вычисляем новый вектор коэффициентов
vector_new = [vector_init[0]-l*grad_a, vector_init[1]-l*grad_b]
return vector_new
# определим функцию стох.град.спуска
def stoch_grad_descent_usual(x_us, y_us, l=0.1, steps = 800):
# для самого начала работы функции зададим начальные значения коэффициентов
vector_init = [float(random.uniform(-0.5, 0.5)), float(random.uniform(-0.5, 0.5))]
errors = []
# запустим цикл спуска
# цикл расчитан на определенное количество шагов (steps)
for i in range(steps):
ind = random.choice(range(len(x_us)))
new_vector = stoch_grad_step_usual(vector_init, x_us, ind, y_us, l)
vector_init = new_vector
errors.append(errors_sq_Kramer_method(vector_init,x_us,y_us))
return (vector_init),(errors)
# запишем массив значений
list_parametres_stoch_gradient_descence = stoch_grad_descent_usual(x_us, y_us, l=0.1, steps = 800)
print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_stoch_gradient_descence[0][0],3)
print 'b =', round(list_parametres_stoch_gradient_descence[0][1],3)
print
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_stoch_gradient_descence[1][-1],3)
print
print ' 33[1m' + ' 33[4m' + "Количество итераций в стохастическом градиентном спуске:" + ' 33[0m'
print len(list_parametres_stoch_gradient_descence[1])
Pozorno pogledamo koeficiente in se ujamemo, da postavljamo vprašanje "Kako je to mogoče?" Imamo druge vrednosti koeficientov и . Mogoče je stohastični gradientni spust našel bolj optimalne parametre za enačbo? Žal ne. Dovolj je, da pogledamo vsoto kvadratov odstopanj in vidimo, da je z novimi vrednostmi koeficientov napaka večja. Ne mudi se nam v obup. Zgradimo graf spremembe napake.
Koda za risanje vsote kvadratov odstopanj pri stohastičnem gradientnem spuščanju
print 'График №5 "Сумма квадратов отклонений по-шагово"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1])), list_parametres_stoch_gradient_descence[1], color='red', lw=2)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()Graf št. 5 “Vsota kvadratov odklonov med spuščanjem stohastičnega gradienta”
Če pogledamo urnik, se vse postavi na svoje mesto in zdaj bomo vse popravili.
Torej kaj se je zgodilo? Zgodilo se je naslednje. Ko naključno izberemo mesec, potem naš algoritem za izbrani mesec skuša zmanjšati napako pri izračunu prihodkov. Nato izberemo drug mesec in ponovimo izračun, vendar zmanjšamo napako za drugi izbrani mesec. Ne pozabite, da prva dva meseca močno odstopata od linije preproste linearne regresijske enačbe. To pomeni, da ko je izbran kateri koli od teh dveh mesecev, naš algoritem z zmanjšanjem napake vsakega od njih resno poveča napako za celoten vzorec. Torej, kaj narediti? Odgovor je preprost: zmanjšati morate korak spusta. Navsezadnje bo z zmanjšanjem koraka spusta tudi napaka prenehala "skakati" gor in dol. Oziroma napaka "skakanje" se ne bo ustavila, vendar ne bo tako hitro :) Preverimo.
Koda za izvajanje SGD z manjšimi koraki
# запустим функцию, уменьшив шаг в 100 раз и увеличив количество шагов соответсвующе
list_parametres_stoch_gradient_descence = stoch_grad_descent_usual(x_us, y_us, l=0.001, steps = 80000)
print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_stoch_gradient_descence[0][0],3)
print 'b =', round(list_parametres_stoch_gradient_descence[0][1],3)
print
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_stoch_gradient_descence[1][-1],3)
print
print ' 33[1m' + ' 33[4m' + "Количество итераций в стохастическом градиентном спуске:" + ' 33[0m'
print len(list_parametres_stoch_gradient_descence[1])
print 'График №6 "Сумма квадратов отклонений по-шагово"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1])), list_parametres_stoch_gradient_descence[1], color='red', lw=2)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()
Graf št. 6 “Vsota kvadratov odstopanj med stohastičnim gradientnim spuščanjem (80 tisoč korakov)”
Koeficienti so se izboljšali, vendar še vedno niso idealni. Hipotetično je to mogoče popraviti na ta način. Izberemo, na primer, v zadnjih 1000 ponovitvah vrednosti koeficientov, s katerimi je bila narejena najmanjša napaka. Res je, za to bomo morali zapisati tudi vrednosti samih koeficientov. Tega ne bomo storili, ampak bomo raje pozorni na urnik. Videti je gladko in zdi se, da se napaka enakomerno zmanjšuje. Pravzaprav to ni res. Poglejmo prvih 1000 ponovitev in jih primerjajmo z zadnjimi.
Koda za grafikon SGD (prvih 1000 korakov)
print 'График №7 "Сумма квадратов отклонений по-шагово. Первые 1000 итераций"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][:1000])),
list_parametres_stoch_gradient_descence[1][:1000], color='red', lw=2)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()
print 'График №7 "Сумма квадратов отклонений по-шагово. Последние 1000 итераций"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][-1000:])),
list_parametres_stoch_gradient_descence[1][-1000:], color='red', lw=2)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()Graf št. 7 “Vsota kvadratov odstopanj SGD (prvih 1000 korakov)”
Graf št. 8 “Vsota kvadratov odstopanj SGD (zadnjih 1000 korakov)”
Na samem začetku spusta opazimo dokaj enakomerno in strmo zmanjšanje napake. V zadnjih iteracijah vidimo, da gre napaka okrog in okoli vrednosti 1,475 in je v nekaterih trenutkih celo enaka tej optimalni vrednosti, potem pa gre še vedno navzgor ... Ponavljam, lahko si zapišete vrednosti koeficientov и , nato pa izberite tiste, pri katerih je napaka minimalna. Vendar smo imeli resnejšo težavo: narediti smo morali 80 tisoč korakov (glej kodo), da smo dobili vrednosti blizu optimalnih. In to je že v nasprotju z idejo o prihranku računalnega časa s stohastičnim gradientnim spuščanjem glede na gradientni spust. Kaj je mogoče popraviti in izboljšati? Ni težko opaziti, da se v prvih iteracijah samozavestno spuščamo navzdol, zato bi morali v prvih iteracijah pustiti velik korak in zmanjševati korak naprej. V tem članku tega ne bomo storili - že predolg je. Tisti, ki želijo, se lahko sami zamislijo, kako to narediti, ni težko :)
Zdaj pa izvedimo stohastični gradientni spust z uporabo knjižnice numpy (in ne spotikajmo se ob kamne, ki smo jih identificirali prej)
Koda za stohastični gradientni spust (NumPy)
# для начала напишем функцию градиентного шага
def stoch_grad_step_numpy(vector_init, X, ind, y, l):
x = X[ind]
y_pred = np.dot(x,vector_init)
err = y_pred - y[ind]
grad_a = err
grad_b = x[1]*err
return vector_init - l*np.array([grad_a, grad_b])
# определим функцию стохастического градиентного спуска
def stoch_grad_descent_numpy(X, y, l=0.1, steps = 800):
vector_init = np.array([[np.random.randint(X.shape[0])], [np.random.randint(X.shape[0])]])
errors = []
for i in range(steps):
ind = np.random.randint(X.shape[0])
new_vector = stoch_grad_step_numpy(vector_init, X, ind, y, l)
vector_init = new_vector
errors.append(error_square_numpy(vector_init,X,y))
return (vector_init), (errors)
# запишем массив значений
list_parametres_stoch_gradient_descence = stoch_grad_descent_numpy(x_np, y_np, l=0.001, steps = 80000)
print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_stoch_gradient_descence[0][0],3)
print 'b =', round(list_parametres_stoch_gradient_descence[0][1],3)
print
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_stoch_gradient_descence[1][-1],3)
print
print ' 33[1m' + ' 33[4m' + "Количество итераций в стохастическом градиентном спуске:" + ' 33[0m'
print len(list_parametres_stoch_gradient_descence[1])
print
Vrednosti so se izkazale za skoraj enake kot pri spustu brez uporabe numpy. Vendar je to logično.
Ugotovimo, koliko časa so nam vzeli stohastični gradientni spusti.
Koda za določanje časa izračuna SGD (80 tisoč korakov)
print ' 33[1m' + ' 33[4m' +
"Время выполнения стохастического градиентного спуска без использования библиотеки NumPy:"
+ ' 33[0m'
%timeit list_parametres_stoch_gradient_descence = stoch_grad_descent_usual(x_us, y_us, l=0.001, steps = 80000)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' +
"Время выполнения стохастического градиентного спуска с использованием библиотеки NumPy:"
+ ' 33[0m'
%timeit list_parametres_stoch_gradient_descence = stoch_grad_descent_numpy(x_np, y_np, l=0.001, steps = 80000)
Čim globlje v gozd, temnejši so oblaki: ponovno »samonapisana« formula kaže najboljši rezultat. Vse to nakazuje, da morajo obstajati še bolj subtilni načini uporabe knjižnice numpy, ki resnično pospešijo računske operacije. V tem članku se ne bomo učili o njih. V prostem času bo nekaj za razmišljati :)
Povzemamo
Preden povzamem, bi rad odgovoril na vprašanje, ki se je najverjetneje pojavilo pri naši dragi bralki. Zakaj pravzaprav takšno »mučenje« s spusti, zakaj moramo hoditi gor in dol po gori (večinoma navzdol), da bi našli dragoceno nižino, če imamo v rokah tako zmogljivo in preprosto napravo, v obliki analitične rešitve, ki nas v trenutku teleportira na Pravo mesto?
Odgovor na to vprašanje leži na površini. Zdaj smo pogledali zelo preprost primer, v katerem je pravi odgovor odvisno od enega znaka . Tega v življenju ne vidite pogosto, zato si predstavljajmo, da imamo 2, 30, 50 ali več znamenj. Temu dodamo na tisoče ali celo na desettisoče vrednosti za vsak atribut. V tem primeru analitična raztopina morda ne bo zdržala preskusa in bo padla. Gradientni spust in njegove različice pa nas bodo počasi, a zanesljivo pripeljali bližje cilju - minimumu funkcije. In ne skrbite za hitrost - verjetno bomo pogledali načine, ki nam bodo omogočili nastavitev in uravnavanje dolžine koraka (torej hitrosti).
In zdaj pravi kratek povzetek.
Prvič, upam, da bo gradivo, predstavljeno v članku, pomagalo začetnikom »podatkovnim znanstvenikom« pri razumevanju reševanja preprostih (in ne samo) enačb linearne regresije.
Drugič, pogledali smo več načinov za rešitev enačbe. Zdaj lahko glede na situacijo izberemo tisto, ki je najprimernejša za rešitev težave.
Tretjič, videli smo moč dodatnih nastavitev, in sicer dolžino koraka gradientnega spuščanja. Tega parametra ni mogoče zanemariti. Kot je navedeno zgoraj, je treba za zmanjšanje stroškov izračunov spremeniti dolžino koraka med spustom.
Četrtič, v našem primeru so »doma napisane« funkcije pokazale najboljše časovne rezultate za izračune. To je verjetno posledica ne najbolj profesionalne uporabe zmogljivosti knjižnice numpy. Kakor koli že, naslednji sklep se nakazuje sam od sebe. Po eni strani je včasih vredno dvomiti o ustaljenih mnenjih, po drugi strani pa se ne splača vedno vsega komplicirati – nasprotno, včasih je preprostejši način reševanja problema učinkovitejši. In ker je bil naš cilj analizirati tri pristope k reševanju enostavne linearne regresijske enačbe, nam je bila uporaba »samonapisanih« funkcij povsem dovolj.
Literatura (ali kaj podobnega)
1. Linearna regresija
2. Metoda najmanjših kvadratov
3. Izpeljanka
4. Gradient
5. Gradientni spust
6. Knjižnica NumPy
Vir: www.habr.com