🥇Rješavanje jednadžbe jednostavne linearne regresije

U članku se razmatra nekoliko načina za određivanje matematičke jednadžbe jednostavne (uparene) regresijske linije.

Sve metode rješavanja jednadžbe o kojima se ovdje raspravlja temelje se na metodi najmanjih kvadrata. Označimo metode na sljedeći način:

Analitičko rješenje
Gradient Descent
Stohastički gradijentni pad

Za svaku od metoda za rješavanje jednačine prave, u članku su navedene različite funkcije, koje se uglavnom dijele na one koje se pišu bez korištenja biblioteke numpy i one koje koriste za proračune numpy. Vjeruje se da je to vješto korištenje numpy smanjiće troškove računara.

Sav kod dat u članku je napisan na jeziku python 2.7 koristeći Jupyter Notebook. Izvorni kod i datoteka sa uzorcima podataka su objavljeni Github

Članak je više namijenjen i početnicima i onima koji su već postupno počeli savladavati proučavanje vrlo širokog dijela umjetne inteligencije – strojnog učenja.

Za ilustraciju materijala koristimo vrlo jednostavan primjer.

Primer uslova

Imamo pet vrijednosti koje karakteriziraju ovisnost Y iz X (Tabela br. 1):

Tabela br. 1 “Primjer uvjeta”

Pretpostavit ćemo da su vrijednosti je mjesec u godini, i — prihod ovog mjeseca. Drugim riječima, prihod zavisi od mjeseca u godini, i - jedini znak od kojeg zavisi prihod.

Primjer je tako-tako, kako sa stanovišta uslovne zavisnosti prihoda od mjeseca u godini, tako i sa stanovišta broja vrijednosti - vrlo ih je malo. Međutim, takvo pojednostavljenje će omogućiti, kako kažu, da se objasni, ne uvijek s lakoćom, materijal koji početnici usvajaju. A i jednostavnost brojeva omogućit će onima koji žele riješiti primjer na papiru bez značajnih troškova rada.

Pretpostavimo da se zavisnost data u primjeru može prilično dobro aproksimirati matematičkom jednadžbom jednostavne (uparene) regresijske linije oblika:

gdje je mjesec u kojem je prihod primljen, — prihod koji odgovara mjesecu, и su koeficijenti regresije procijenjene linije.

Imajte na umu da je koeficijent često se naziva nagib ili gradijent procijenjene linije; predstavlja iznos za koji je kad se promijeni .

Očigledno, naš zadatak u primjeru je da izaberemo takve koeficijente u jednadžbi и , pri čemu su odstupanja naših izračunatih vrijednosti prihoda po mjesecima od tačnih odgovora, tj. vrijednosti prikazane u uzorku će biti minimalne.

Metoda najmanjeg kvadrata

Prema metodi najmanjih kvadrata, odstupanje treba izračunati kvadriranjem. Ova tehnika vam omogućava da izbjegnete međusobno otkazivanje odstupanja ako imaju suprotne znakove. Na primjer, ako je u jednom slučaju, odstupanje je +5 (plus pet), au drugom -5 (minus pet), tada će se zbir odstupanja međusobno poništavati i iznositi 0 (nula). Moguće je ne kvadrirati odstupanje, već koristiti svojstvo modula i tada će sva odstupanja biti pozitivna i akumulirati. Nećemo se detaljnije zadržavati na ovoj točki, već jednostavno naznačimo da je za praktičnost proračuna uobičajeno kvadrirati odstupanje.

Ovako izgleda formula pomoću koje ćemo odrediti najmanji zbir kvadrata odstupanja (greške):

gdje je funkcija aproksimacije tačnih odgovora (tj. prihoda koji smo izračunali),

su tačni odgovori (prihod naveden u uzorku),

je indeks uzorka (broj mjeseca u kojem je utvrđeno odstupanje)

Hajde da izdiferenciramo funkciju, definišemo parcijalne diferencijalne jednadžbe i budimo spremni da pređemo na analitičko rešenje. Ali prvo, hajde da napravimo kratak izlet o tome šta je diferencijacija i prisjetimo se geometrijskog značenja derivacije.

Diferencijacija

Diferencijacija je operacija pronalaženja derivacije funkcije.

Za šta se koristi derivat? Derivat funkcije karakterizira brzinu promjene funkcije i govori nam njen smjer. Ako je derivacija u datoj tački pozitivna, tada se funkcija povećava; u suprotnom se funkcija smanjuje. I što je veća vrijednost apsolutne derivacije, to je veća brzina promjene vrijednosti funkcije, kao i strmiji nagib grafa funkcije.

Na primjer, pod uslovima kartezijanskog koordinatnog sistema, vrijednost derivacije u tački M(0,0) jednaka je + 25 znači da u datoj tački, kada se vrijednost pomjeri desno za konvencionalnu jedinicu, vrijednost povećava za 25 konvencionalnih jedinica. Na grafikonu to izgleda kao prilično strm porast vrijednosti iz date tačke.

Još jedan primjer. Vrijednost derivata je jednaka -0,1 znači da kada se raseljava po jednoj konvencionalnoj jedinici, vrijednost smanjuje se za samo 0,1 konvencionalnu jedinicu. Istovremeno, na grafu funkcije možemo uočiti jedva primjetan silazni nagib. Povodeći analogiju sa planinom, kao da se vrlo polako spuštamo niz blagu padinu sa planine, za razliku od prethodnog primera gde smo morali da se penjemo na veoma strme vrhove :)

Dakle, nakon diferenciranja funkcije po izgledima и , definiramo parcijalne diferencijalne jednadžbe 1. reda. Nakon određivanja jednadžbi, dobićemo sistem od dvije jednadžbe, rješavanjem kojih ćemo moći odabrati takve vrijednosti koeficijenata и , za koje se vrijednosti odgovarajućih derivacija u datim tačkama mijenjaju za vrlo, vrlo mali iznos, a u slučaju analitičkog rješenja uopće se ne mijenjaju. Drugim riječima, funkcija greške pri pronađenim koeficijentima će dostići minimum, jer će vrijednosti parcijalnih izvoda u tim točkama biti jednake nuli.

Dakle, prema pravilima diferencijacije, jednadžba parcijalnog derivata 1. reda u odnosu na koeficijent će poprimiti oblik:

Jednačina parcijalnog izvoda 1. reda u odnosu na će poprimiti oblik:

Kao rezultat, dobili smo sistem jednačina koji ima prilično jednostavno analitičko rješenje:

početi{jednačina*}
započeti{slučajevi}
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
kraj{case}
kraj{jednačina*}

Prije rješavanja jednadžbe izvršimo prethodno učitavanje, provjerimo da li je učitavanje ispravno i formatirajmo podatke.

Učitavanje i formatiranje podataka

Treba napomenuti da ćemo zbog činjenice da ćemo za analitičko rješenje, a potom i za gradijentni i stohastički gradijentni spust, koristiti kod u dvije varijacije: korištenjem biblioteke numpy i bez upotrebe, tada će nam trebati odgovarajuće formatiranje podataka (vidi kod).

Kod za učitavanje i obradu podataka

# импортируем все нужные нам библиотеки
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

Sada, nakon što smo prvo učitali podatke, drugo, provjerili ispravnost učitavanja i konačno formatirali podatke, izvršit ćemo prvu vizualizaciju. Metoda koja se često koristi za ovo je parplot biblioteke rođen na moru. U našem primjeru, zbog ograničenog broja, nema smisla koristiti biblioteku rođen na moru. Koristićemo redovnu biblioteku matplotlib i samo pogledajte dijagram raspršenosti.

Scatterplot code

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 br. 1 “Zavisnost prihoda od mjeseca u godini”

Analitičko rješenje

Koristimo najčešće alate u python i riješi sistem jednačina:

početi{jednačina*}
započeti{slučajevi}
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
kraj{case}
kraj{jednačina*}

Prema Cramerovom pravilu naći ćemo opštu odrednicu, kao i determinante po i po , nakon čega, dijeljenje determinante sa na opštu odrednicu - nađi koeficijent , na sličan način nalazimo koeficijent .

Kod analitičkog rješenja

# определим функцию для расчета коэффициентов 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 šta smo dobili:

Dakle, pronađene su vrijednosti koeficijenata, utvrđen je zbir kvadrata odstupanja. Nacrtajmo pravu liniju na histogramu raspršenja u skladu sa pronađenim koeficijentima.

Kod regresijske linije

# определим функцию для формирования массива рассчетных значений выручки
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 br. 2 “Tačni i proračunati odgovori”

Možete pogledati grafikon odstupanja za svaki mjesec. U našem slučaju, nećemo izvući nikakvu značajnu praktičnu vrijednost iz toga, ali ćemo zadovoljiti našu znatiželju o tome koliko dobro prosta jednačina linearne regresije karakterizira ovisnost prihoda od mjeseca u godini.

Kod grafikona odstupanja

# определим функцию для формирования массива отклонений в процентах
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 br. 3 “Odstupanja, %”

Nije savršeno, ali smo završili svoj zadatak.

Napišimo funkciju koja će odrediti koeficijente и koristi biblioteku numpy, preciznije, napisaćemo dvije funkcije: jednu koristeći pseudoinverznu matricu (ne preporučuje se u praksi, budući da je proces računski složen i nestabilan), drugu koristeći matričnu jednadžbu.

Kod analitičkog rješenja (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_np

Uporedimo vrijeme utrošeno na određivanje koeficijenata и , u skladu sa 3 predstavljene metode.

Šifra za izračunavanje vremena računanja

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)

Sa malom količinom podataka, ispred se pojavljuje „samonapisana“ funkcija koja pronalazi koeficijente koristeći Cramerovu metodu.

Sada možete preći na druge načine pronalaženja koeficijenata и .

Gradient Descent

Prvo, hajde da definišemo šta je gradijent. Jednostavno rečeno, gradijent je segment koji označava smjer maksimalnog rasta funkcije. Po analogiji sa penjanjem na planinu, gdje je uspon tamo gdje je najstrmiji uspon na vrh planine. Razvijajući primjer s planinom, sjećamo se da nam je zapravo potreban najstrmiji spust da bismo što brže stigli do nizije, odnosno minimuma - mjesta gdje se funkcija ne povećava niti smanjuje. U ovom trenutku derivacija će biti jednaka nuli. Dakle, ne treba nam gradijent, već antigradijent. Da biste pronašli antigradijent, samo trebate pomnožiti gradijent sa -1 (minus jedan).

Obratimo pažnju na činjenicu da funkcija može imati nekoliko minimuma, a spuštajući se u jedan od njih koristeći algoritam predložen u nastavku, nećemo moći pronaći drugi minimum, koji može biti niži od pronađenog. Opustimo se, ovo nam nije prijetnja! U našem slučaju imamo posla sa jednim minimumom, budući da je naša funkcija na grafu je pravilna parabola. I kao što svi dobro znamo iz našeg školskog kursa matematike, parabola ima samo jedan minimum.

Nakon što smo saznali zašto nam je potreban gradijent, a takođe i da je gradijent segment, odnosno vektor sa datim koordinatama, koje su potpuno isti koeficijenti и možemo implementirati gradijentni spust.

Prije početka, predlažem da pročitate samo nekoliko rečenica o algoritmu spuštanja:

Određujemo koordinate koeficijenata na pseudo-slučajan način и . U našem primjeru ćemo definirati koeficijente blizu nule. Ovo je uobičajena praksa, ali svaki slučaj može imati svoju praksu.
Od koordinata oduzmi vrijednost parcijalnog izvoda 1. reda u tački . Dakle, ako je izvod pozitivan, tada se funkcija povećava. Dakle, oduzimanjem vrijednosti derivacije kretat ćemo se u suprotnom smjeru rasta, odnosno u smjeru pada. Ako je izvod negativan, tada funkcija u ovoj tački opada i oduzimanjem vrijednosti izvoda krećemo se u smjeru pada.
Izvodimo sličnu operaciju s koordinatom : oduzmi vrijednost parcijalnog izvoda u tački .
Kako ne biste preskočili minimum i odletjeli u duboki svemir, potrebno je podesiti veličinu koraka u smjeru spuštanja. Općenito, možete napisati cijeli članak o tome kako ispravno postaviti korak i kako ga promijeniti tokom procesa spuštanja kako biste smanjili troškove računanja. Ali sada je pred nama nešto drugačiji zadatak, a veličinu koraka ćemo utvrditi naučnom metodom “bockanja” ili, kako se u narodu kaže, empirijski.
Jednom smo iz datih koordinata и oduzimamo vrijednosti izvoda, dobivamo nove koordinate и . Poduzimamo sljedeći korak (oduzimanje), već od izračunatih koordinata. I tako ciklus počinje iznova i iznova, sve dok se ne postigne potrebna konvergencija.

Sve! Sada smo spremni da krenemo u potragu za najdubljom klisurom Marijanskog rova. Hajde da počnemo.

Kod za gradijentni pad

# напишем функцию градиентного спуска без использования библиотеки 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

Zaronili smo do samog dna Marijanskog rova i tamo smo našli sve iste vrijednosti koeficijenata и , što je upravo ono što se moglo očekivati.

Hajdemo još jednom zaroniti, samo što će ovaj put naše dubokomorsko vozilo biti ispunjeno drugim tehnologijama, odnosno bibliotekom numpy.

Kod za gradijentni pad (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

Vrijednosti koeficijenata и nepromjenjivo.

Pogledajmo kako se greška promijenila tokom gradijenta spuštanja, odnosno kako se mijenja zbir kvadrata odstupanja sa svakim korakom.

Kod za crtanje zbira kvadrata odstupanja

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()

Grafikon br. 4 „Zbir kvadrata odstupanja tokom gradijenta spuštanja“

Na grafikonu vidimo da se sa svakim korakom greška smanjuje, a nakon određenog broja iteracija uočavamo gotovo horizontalnu liniju.

Na kraju, procijenimo razliku u vremenu izvršenja koda:

Kod za određivanje vremena izračunavanja gradijenta spuštanja

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)

Možda radimo nešto pogrešno, ali opet to je jednostavna "kućno napisana" funkcija koja ne koristi biblioteku numpy nadmašuje vrijeme izračunavanja funkcije koja koristi biblioteku numpy.

Ali mi ne stojimo mirno, već se krećemo ka proučavanju još jednog uzbudljivog načina rješavanja jednostavne linearne regresijske jednadžbe. Upoznajte nas!

Stohastički gradijentni pad

Da bi se brzo razumio princip rada stohastičkog gradijentnog spuštanja, bolje je utvrditi njegove razlike od običnog gradijentnog spuštanja. Mi, u slučaju gradijentnog spuštanja, u jednadžbama izvoda od и koristio zbroje vrijednosti svih karakteristika i istinitih odgovora dostupnih u uzorku (tj. zbira svih и ). U stohastičkom gradijentnom spuštanju, nećemo koristiti sve vrijednosti prisutne u uzorku, već ćemo umjesto toga pseudo-slučajno odabrati takozvani indeks uzorka i koristiti njegove vrijednosti.

Na primjer, ako se odredi da je indeks broj 3 (tri), onda uzimamo vrijednosti и , zatim zamjenjujemo vrijednosti u jednadžbe derivata i određujemo nove koordinate. Zatim, nakon što smo odredili koordinate, ponovo pseudo-slučajno određujemo indeks uzorka, zamjenjujemo vrijednosti koje odgovaraju indeksu u parcijalne diferencijalne jednadžbe i određujemo koordinate na novi način и itd. sve dok konvergencija ne postane zelena. Na prvi pogled možda izgleda da ovo uopšte ne bi moglo da funkcioniše, ali jeste. Istina je da je vrijedno napomenuti da se greška ne smanjuje sa svakim korakom, ali tendencija svakako postoji.

Koje su prednosti stohastičkog gradijentnog spuštanja u odnosu na konvencionalni? Ako je veličina našeg uzorka vrlo velika i mjeri se u desetinama hiljada vrijednosti, onda je mnogo lakše obraditi, recimo, nasumično hiljadu njih, a ne cijeli uzorak. Ovdje dolazi u obzir stohastički gradijentni pad. U našem slučaju, naravno, nećemo primijetiti veliku razliku.

Pogledajmo kod.

Kod za stohastički gradijentni pad

# определим функцию стох.град.шага
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])

Pažljivo gledamo koeficijente i uhvatimo se da postavljamo pitanje "Kako je to moguće?" Dobili smo druge vrijednosti koeficijenta и . Možda je stohastički gradijentni pad pronašao optimalnije parametre za jednadžbu? Nažalost nema. Dovoljno je pogledati zbir kvadrata odstupanja i vidjeti da je s novim vrijednostima koeficijenata greška veća. Ne žurimo u očajanje. Napravimo grafikon promjene greške.

Kod za crtanje zbira kvadrata odstupanja u stohastičkom gradijentu spuštanja

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()

Grafikon br. 5 “Zbir kvadrata odstupanja tokom stohastičkog gradijenta spuštanja”

Gledajući raspored, sve dolazi na svoje mjesto i sad ćemo sve popraviti.

Šta se desilo? Desilo se sljedeće. Kada nasumično odaberemo mjesec, tada naš algoritam za odabrani mjesec nastoji smanjiti grešku u izračunavanju prihoda. Zatim biramo drugi mjesec i ponavljamo obračun, ali smanjujemo grešku za drugi odabrani mjesec. Sada zapamtite da prva dva mjeseca značajno odstupaju od linije jednostavne linearne regresijske jednačine. To znači da kada se odabere bilo koji od ova dva mjeseca, smanjenjem greške svakog od njih, naš algoritam ozbiljno povećava grešku za cijeli uzorak. Pa šta da radimo? Odgovor je jednostavan: morate smanjiti korak spuštanja. Uostalom, smanjenjem koraka spuštanja, greška će također prestati "skakati" gore-dolje. Tačnije, greška "skakanja" neće prestati, ali neće to učiniti tako brzo :) Hajde da provjerimo.

Kod za pokretanje SGD sa manjim koracima

# запустим функцию, уменьшив шаг в 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()

Grafikon br. 6 “Zbir kvadrata odstupanja tokom stohastičkog gradijenta spuštanja (80 hiljada koraka)”

Koeficijenti su poboljšani, ali još uvijek nisu idealni. Hipotetički, ovo se može ispraviti na ovaj način. Odabiremo, na primjer, u zadnjih 1000 iteracija vrijednosti koeficijenata s kojima je napravljena minimalna greška. Istina, za to ćemo morati zapisati i vrijednosti samih koeficijenata. Nećemo to raditi, već ćemo obratiti pažnju na raspored. Izgleda glatko i čini se da se greška ravnomjerno smanjuje. Zapravo to nije istina. Pogledajmo prvih 1000 iteracija i uporedimo ih s posljednjim.

Kod za SGD grafikon (prvih 1000 koraka)

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()

Grafikon br. 7 “Zbir kvadratnih devijacija SGD (prvih 1000 koraka)”

Grafikon br. 8 “Zbir kvadratnih devijacija SGD (zadnjih 1000 koraka)”

Na samom početku spuštanja uočavamo prilično ujednačeno i strmo smanjenje greške. U zadnjim iteracijama vidimo da se greška vrti okolo i oko vrijednosti 1,475 iu nekim trenucima se čak i izjednačava sa ovoj optimalnoj vrijednosti, ali onda se i dalje povećava... Ponavljam, možete zapisati vrijednosti koeficijenti и , a zatim odaberite one za koje je greška minimalna. Međutim, imali smo ozbiljniji problem: morali smo napraviti 80 hiljada koraka (vidi kod) da bismo vrijednosti približili optimalnim. A to je već u suprotnosti s idejom uštede vremena računanja sa stohastičkim gradijentom spuštanja u odnosu na gradijentni pad. Šta se može ispraviti i poboljšati? Nije teško primijetiti da u prvim iteracijama samopouzdano idemo dolje i stoga bismo trebali ostaviti veliki korak u prvim iteracijama i smanjivati korak kako idemo naprijed. Nećemo to raditi u ovom članku - već je predugačak. Oni koji žele mogu sami da smisle kako da ovo urade, nije teško :)

Sada izvršimo stohastički gradijentni spust koristeći biblioteku numpy (i nemojmo se spotaknuti o kamenje koje smo ranije identifikovali)

Kod za stohastički gradijentni pad (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

Ispostavilo se da su vrijednosti gotovo iste kao kod spuštanja bez korištenja numpy. Međutim, ovo je logično.

Hajde da saznamo koliko su nam dugo trajali stohastički gradijenti.

Kod za određivanje vremena izračunavanja SGD (80 hiljada koraka)

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)

Što dalje u šumu, to su oblaci tamniji: opet, „samonapisana“ formula pokazuje najbolji rezultat. Sve ovo sugerira da moraju postojati još suptilniji načini korištenja biblioteke numpy, što zaista ubrzava računske operacije. U ovom članku nećemo učiti o njima. Imat ćete o čemu razmišljati u slobodno vrijeme :)

Mi sumiramo

Prije rezimiranja, želio bih odgovoriti na pitanje koje je najvjerovatnije poteklo od našeg dragog čitatelja. Čemu, zapravo, takvo „mučenje“ sa spuštanjima, zašto moramo hodati gore-dolje po planini (uglavnom dolje) da bismo pronašli dragocenu niziju, ako u rukama imamo tako moćnu i jednostavnu spravu, u oblik analitičkog rješenja, koje nas momentalno teleportuje na pravo mjesto?

Odgovor na ovo pitanje leži na površini. Sada smo pogledali vrlo jednostavan primjer u kojem je pravi odgovor zavisi od jednog znaka . Ovo se ne viđa često u životu, pa zamislimo da imamo 2, 30, 50 ili više znakova. Dodajmo ovome hiljade, ili čak desetine hiljada vrijednosti za svaki atribut. U tom slučaju, analitičko rješenje možda neće izdržati test i neće uspjeti. Zauzvrat, gradijentni spust i njegove varijacije će nas polako ali sigurno dovesti bliže cilju – minimumu funkcije. I ne brinite o brzini - vjerovatno ćemo pogledati načine koji će nam omogućiti da podesimo i regulišemo dužinu koraka (tj. brzinu).

A sada pravi kratak sažetak.

Prvo, nadam se da će materijal predstavljen u članku pomoći početnicima „naučnika podataka“ u razumijevanju kako riješiti jednostavne (i ne samo) jednačine linearne regresije.

Drugo, pogledali smo nekoliko načina za rješavanje jednačine. Sada, ovisno o situaciji, možemo odabrati onu koja je najprikladnija za rješavanje problema.

Treće, vidjeli smo snagu dodatnih postavki, odnosno gradijentnu dužinu koraka spuštanja. Ovaj parametar se ne može zanemariti. Kao što je gore navedeno, kako bi se smanjili troškovi proračuna, dužinu koraka treba promijeniti tokom spuštanja.

Četvrto, u našem slučaju, “kućne” funkcije su pokazale najbolje vremenske rezultate za proračune. To je vjerovatno zbog neprofesionalne upotrebe bibliotečkih mogućnosti numpy. Ali kako god bilo, sljedeći zaključak se nameće sam od sebe. S jedne strane, ponekad vrijedi preispitivati ustaljena mišljenja, a s druge strane, ne vrijedi uvijek sve komplikovati – naprotiv, ponekad je jednostavniji način rješavanja problema efikasniji. A budući da nam je cilj bio analizirati tri pristupa rješavanju jednostavne linearne regresijske jednadžbe, upotreba “samopisnih” funkcija nam je bila sasvim dovoljna.