🥇Вырашаем раўнанне простай лінейнай рэгрэсіі

У артыкуле разглядаецца некалькі спосабаў вызначэння матэматычнага ўраўнення лініі простай (парнай) рэгрэсіі.

Усе разгляданыя тут спосабы рашэння раўнання заснаваны на метадзе найменшых квадратаў. Абазначым спосабы наступным чынам:

Аналітычнае рашэнне
Градыентны спуск
Выпадковы градыентны спуск

Для кожнага са спосабаў рашэння ўраўненні прамой, у артыкуле прыведзены розныя функцыі, якія ў асноўным дзеляцца на тыя, якія напісаны без выкарыстання бібліятэкі. numpy і тыя, якія для правядзення разлікаў ужываюць numpy. Лічыцца, што ўмелае выкарыстанне numpy дазволіць скараціць выдаткі на вылічэнні.

Увесь код, прыведзены ў артыкуле, напісаны на мове python 2.7 з выкарыстаннем Jupyter наўтбукаў. Зыходны код і файл з дадзенымі выбаркі выкладзены на Гітхабе

Артыкул у большай ступені арыентаваны як на пачаткоўцаў, так і на тых, хто ўжо патроху пачаў асвойваць вывучэнне вельмі шырокага падзелу ў штучным інтэлекце - машыннага навучання.

Для ілюстрацыі матэрыялу выкарыстоўваем вельмі просты прыклад.

Умовы прыкладу

У нас ёсць пяць значэнняў, якія характарызуюць залежнасць Y ад X (Табліца №1):

Табліца №1 "Умовы прыкладу"

Будзем лічыць, што значэння - гэта месяц года, а - Выручка ў гэтым месяцы. Іншымі словамі, выручка залежыць ад месяца года, а - Адзіная прыкмета, ад якога залежыць выручка.

Прыклад так сабе, як з пункту гледжання ўмоўнай залежнасці выручкі ад месяца года, так і з пункту гледжання колькасці значэнняў іх вельмі мала. Аднак такое спрашчэнне дазволіць, што завецца на пальцах, растлумачыць, не заўсёды з лёгкасцю, які засвойваецца пачаткоўцамі матэрыял. А таксама прастата лікаў дазволіць без важкіх працавыдаткаў, жадаючым, вырашыць прыклад на «паперы».

Выкажам здагадку, што прыведзеная ў прыкладзе залежнасць, можа быць досыць добра апраксімаваная матэматычным раўнаннем лініі простай (парнай) рэгрэсіі выгляду:

дзе - Гэта месяц, у якім была атрымана выручка, - Выручка, якая адпавядае месяцу, и - Каэфіцыенты рэгрэсіі ацэненай лініі.

Адзначым, што каэфіцыент часта называюць кутнім каэфіцыентам або градыентам ацэненай лініі; уяўляе сабой велічыню, на якую зменіцца пры змене .

Відавочна, што наша задача ў прыкладзе - падабраць ва ўраўненні такія каэфіцыенты и , пры якіх адхіленні нашых разліковых значэнняў выручкі па месяцах ад сапраўдных адказаў, г.зн. значэнняў, прадстаўленых у выбарцы, будуць мінімальныя.

Метад найменшых квадратаў

У адпаведнасці з метадам найменшых квадратаў, адхіленне варта разлічваць, узводзячы яго ў квадрат. Падобны прыём дазваляе пазбегнуць узаемнага пагашэння адхіленняў, у тым выпадку, калі яны маюць супрацьлеглыя знакі. Напрыклад, калі ў адным выпадку, адхіленне складае +5 (плюс пяць), а ў іншым -5 (мінус пяць), то сума адхіленняў узаемна пагасіцца і складзе 0 (нуль). Можна і не ўзводзіць адхіленне ў квадрат, а скарыстацца ўласцівасцю модуля і тады ў нас усе адхіленні будуць дадатнымі і будуць назапашвацца. Мы не будзем спыняцца на гэтым моманце падрабязна, а проста абазначым, што для зручнасці разлікаў, прынята ўзводзіць адхіленне ў квадрат.

Вось так выглядае формула, з дапамогай якой мы вызначым найменшую суму квадратаў адхіленняў (памылкі):

дзе - Гэта функцыя апраксімацыі сапраўдных адказаў (гэта значыць палічаная намі выручка),

- Гэта сапраўдныя адказы (прадстаўленая ў выбарцы выручка),

- Гэта індэкс выбаркі (нумар месяца, у якім адбываецца вызначэнне адхіленні)

Прадыферэнцыруем функцыю, вызначым раўнанні прыватных вытворных і будзем гатовыя перайсці да аналітычнага рашэння. Але для пачатку правядзем невялікі экскурс аб тым, што такое дыферэнцыяванне і ўспомнім геаметрычны сэнс вытворнай.

Дыферэнцыяванне

Дыферэнцыяваннем называецца аперацыя па знаходжанні вытворнай функцыі.

Навошта патрэбна вытворная? Вытворная функцыі характарызуе хуткасць змены функцыі і паказвае нам яе напрамак. Калі вытворная ў зададзенай кропцы дадатная, то функцыя ўзрастае, у адваротным выпадку - функцыя спадае. І чым больш значэнне вытворнай па модулі, тым вышэй хуткасць змены значэнняў функцыі, а таксама страмчэй кут нахілу графіка функцыі.

Напрыклад, ва ўмовах дэкартавай сістэмы каардынат, значэнне вытворнай у кропцы M(0,0) роўнае +25 азначае, што ў зададзенай кропцы, пры зняцці значэння направа на ўмоўную адзінку, значэнне узрастае на 25 умоўных адзінак. На графіцы гэта выглядае, як дастаткова круты кут ўздыму значэнняў з зададзенай кропкі.

Іншы прыклад. Значэнне вытворнай роўнае 0,1- азначае, што пры зняцці на адну ўмоўную адзінку, значэнне меншае ўсяго толькі на 0,1 умоўную адзінку. Пры гэтым, на графіцы функцыі, мы можам назіраць ледзь прыкметны нахіл уніз. Праводзячы аналогію з гарой, то мы як быццам вельмі павольна спускаемся па спадзістым схіле з гары, у адрозненне ад папярэдняга прыкладу, дзе нам даводзілася браць вельмі стромкія вяршыні:)

Такім чынам, правёўшы дыферэнцыяванне функцыі па каэфіцыентах и , вызначым раўнанні прыватных вытворных 1-га парадку. Пасля вызначэння раўнанняў, мы атрымаем сістэму з двух раўнанняў, вырашыўшы якую мы зможам падабраць такія значэнні каэфіцыентаў и , Пры якіх значэння адпаведных вытворных у зададзеных кропках змяняюцца на вельмі і вельмі малую велічыню, а ў выпадку з аналітычным рашэннем не змяняюцца зусім. Іншымі словамі, функцыя памылкі пры знойдзеных каэфіцыентах дасягне мінімуму, бо значэнні прыватных вытворных у гэтых кропках будуць роўныя нулю.

Такім чынам, па правілах дыферэнцыявання раўнанне прыватнай вытворнай 1-га парадку па каэфіцыенце прыме выгляд:

раўнанне прыватнай вытворнай 1-га парадку па прыме выгляд:

У выніку мы атрымалі сістэму раўнанняў, якая мае дастаткова простае аналітычнае рашэнне:

begin{equation*}
begin{cases}
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
end{cases}
end{equation*}

Перш чым вырашаць раўнанне, папярэдне загрузім, праверым правільнасць загрузкі і адфарматуем дадзеныя.

Загрузка і фарматаванне даных

Неабходна адзначыць, што ў сувязі з тым, што для аналітычнага рашэння, а ў далейшым для градыентнага і стахастычнага градыентнага спуску, мы будзем прымяняць код у двух варыяцыях: з выкарыстаннем бібліятэкі numpy і без яе выкарыстання, то нам запатрабуецца якое адпавядае фарматаванне дадзеных (гл. код).

Код загрузкі і апрацоўкі дадзеных

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

Візуалізацыя

Цяпер, пасля таго, як мы, па-першае, загрузілі дадзеныя, па-другое, праверылі правільнасць загрузкі і нарэшце адфарматавалі дадзеныя, правядзём першую візуалізацыю. Часта для гэтага выкарыстоўваюць метад pairplot бібліятэкі Марскі народ. У нашым прыкладзе, з прычыны абмежаванасці лічбаў няма сэнсу прымяняць бібліятэку Марскі народ. Мы скарыстаемся звычайнай бібліятэкай matplotlib і паглядзім толькі на дыяграму рассейвання.

Код дыяграмы рассейвання

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

Графік №1 "Залежнасць выручкі ад месяца года"

Аналітычнае рашэнне

Скарыстаемся самымі звычайнымі прыладамі ў пітон і вырашым сістэму раўнанняў:

begin{equation*}
begin{cases}
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
end{cases}
end{equation*}

Па правіле Крамера знойдзем агульны вызначальнік, а таксама вызначальнікі па і па , пасля чаго, падзяліўшы вызначальнік па на агульны вызначальнік - знойдзем каэфіцыент , аналагічна знойдзем каэфіцыент .

Код аналітычнага рашэння

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

Вось, што ў нас атрымалася:

Такім чынам, значэнні каэфіцыентаў знойдзены, сума квадратаў адхіленняў устаноўлена. Намалюем на гістаграме рассейвання прамую лінію ў адпаведнасці са знойдзенымі каэфіцыентамі.

Код лініі рэгрэсіі

# определим функцию для формирования массива рассчетных значений выручки
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()

Графік №2 "Правільныя і разліковыя адказы"

Можна паглядзець на графік адхіленняў за кожны месяц. У нашым выпадку, які-небудзь значнай практычнай каштоўнасці мы з яго не вынесем, але задаволім цікаўнасць у тым, наколькі добра, раўнанне простай лінейнай рэгрэсіі характарызуе залежнасць выручкі ад месяца года.

Код графіка адхіленняў

# определим функцию для формирования массива отклонений в процентах
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()

Графік №3 "Адхіленні,%"

Не ідэальна, але нашу задачу мы выканалі.

Напішам функцыю, якая для вызначэння каэфіцыентаў и выкарыстоўвае бібліятэку numpy, Дакладней - напішам дзве функцыі: адну з выкарыстаннем псеўдазваротнай матрыцы (не рэкамендуецца на практыцы, бо працэс вылічальна складаны і нестабільны), іншую з выкарыстаннем матрычнага раўнання.

Код аналітычнага рашэння (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

Параўнальны час, які быў затрачаны на вызначэнне каэфіцыентаў и , у адпаведнасці з трыма прадстаўленымі спосабамі.

Код для вылічэння часу разлікаў

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)

На невялікай колькасці дадзеных, наперад выходзіць "самапісная" функцыя, якая знаходзіць каэфіцыенты метадам Крамера.

Цяпер можна перайсці да іншых спосабаў знаходжання каэфіцыентаў и .

Градыентны спуск

Для пачатку вызначым, што такое градыент. Па-простаму, градыент - гэта адрэзак, які паказвае кірунак максімальнага росту функцыі. Па аналогіі з уздымам у гару, то куды глядзіць градыент, тамака і ёсць самы круты ўздым да вяршыні гары. Развіваючы прыклад з гарой, успамінаем, што насамрэч нам патрэбен самы круты спуск, каб як мага хутчэй дасягнуць нізіны, гэта значыць мінімуму месцы дзе функцыя не ўзрастае і не меншае. У гэтым месцы вытворная будзе роўная нулю. Такім чынам, нам патрэбен не градыент, а антыградыент. Для знаходжання антыградыенту трэба ўсяго толькі памножыць градыент на -1 (мінус адзін).

Зварачальны ўвага на тое, што функцыя можа мець некалькі мінімумаў, і апусціўшыся ў адзін з іх па прапанаваным далей алгарытме, мы не зможам знайсці іншы мінімум, які магчыма знаходзіцца ніжэй знойдзенага. Паслабімся, нам гэта не пагражае! У нашым выпадку мы маем справу з адзіным мінімумам, бо наша функцыя на графіцы ўяўляе сабой звычайную парабалу. А як мы ўсе павінны выдатна ведаць са школьнага курсу матэматыкі - у парабалы існуе толькі адзін мінімум.

Пасля таго, як мы высветлілі для чаго нам запатрабаваўся градыент, а таксама тое, што градыент - гэта адрэзак, гэта значыць вектар з зададзенымі каардынатамі, якія як раз з'яўляюцца тымі самымі каэфіцыентамі и мы можам рэалізаваць градыентны спуск.

Перад запускам, прапаную прачытаць літаральна некалькі прапаноў аб алгарытме спуску:

Вызначаем псеўдавыпадковым чынам каардынаты каэфіцыентаў и . У нашым прыкладзе, мы будзем вызначаць каэфіцыенты паблізу нуля. Гэта з'яўляецца распаўсюджанай практыкай, аднак для кожнага выпадку можа быць прадугледжана свая практыка.
Ад каардынаты адымаем значэнне прыватнай вытворнай 1-га парадку ў кропцы . Так, калі вытворная будзе дадатная, то функцыя ўзрастае. Такім чынам, адымаючы значэнне вытворнай, мы будзем рухацца ў зваротны бок росту, гэта значыць у бок спуску. Калі вытворная адмоўная, значыць функцыя ў гэтай кропцы меншае і адымаючы значэнне вытворнай мы рухаемся ў бок спуску.
Праводзім аналагічную аперацыю з каардынатай : адымаем значэнне прыватнай вытворнай у кропцы .
Для таго, каб не пераскочыць мінімум і не паляцець у далёкі космас, неабходна ўстанавіць памер кроку ў бок спуску. Увогуле і цэлым, можна напісаць цэлы артыкул аб тым, як правільней усталяваць крок і як яго змяняць падчас спуску, каб зменшыць выдаткі на вылічэнні. Але цяпер перад намі некалькі іншая задача, і мы навуковым метадам "тыка" ці як кажуць у просты люд, эмпірычным шляхам, усталюем памер кроку.
Пасля таго, як мы з зададзеных каардынат и аднялі значэнні вытворных, атрымліваем новыя каардынаты и . Які робіцца наступны крок (адніманне), ужо з разлічаных каардынат. І так цыкл запускаецца ізноў і ізноў, датуль, пакуль не будзе дасягнута патрабаваная збежнасць.

Усё! Цяпер мы гатовыя адправіцца на пошукі самай глыбокай цясніны Марыянскай упадзіны. Прыступаем.

Код для градыентнага спуску

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

Мы пагрузіліся на самае дно Марыянскай западзіны і там выявілі ўсё тыя ж значэнні каэфіцыентаў и , што ўласна і трэба было чакаць.

Здзейснім яшчэ адно апусканне, толькі на гэты раз, начыннем нашага глыбакаводнага апарата будуць іншыя тэхналогіі, а менавіта бібліятэка numpy.

Код для градыентнага спуску (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

Значэнні каэфіцыентаў и нязменныя.

Паглядзім на тое, як змянялася памылка пры градыентнай спуску, гэта значыць як змянялася сума квадратаў адхіленняў з кожным крокам.

Код для графіка сум квадратаў адхіленняў

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

Графік №4 "Сума квадратаў адхіленняў пры градыентным спуску"

На графіцы мы бачым, што з кожным крокам памылка памяншаецца, а праз нейкую колькасць ітэрацый назіраем практычна гарызантальную лінію.

Напрыканцы ацэнім розніцу ў часе выканання кода:

Код для вызначэння часу вылічэнні градыентнага спуску

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)

Магчыма мы робім нешта не тое, але зноў простая “самапісная” функцыя, якая не выкарыстоўвае бібліятэку numpy апярэджвае па часе выканання разлікаў функцыю, якая выкарыстоўвае бібліятэку numpy.

Але мы не стаім на месцы, а рухаемся ў бок вывучэння яшчэ аднаго займальнага спосабу рашэння раўнання простай лінейнай рэгрэсіі. Сустракайце!

Выпадковы градыентны спуск

Для таго, каб хутчэй зразумець прынцып працы стахастычнага градыентнага спуску, лепш вызначыць яго адрозненні ад звычайнага градыентнага спуску. Мы, у выпадку з градыентным спускам, ва ўраўненнях вытворных ад и выкарыстоўвалі сумы значэнняў усіх прыкмет і сапраўдных адказаў, якія ёсць у выбарцы (гэта значыць сумы ўсіх и ). У стахастычным градыентным спуску мы не будзем выкарыстоўваць усе значэнні, наяўныя ў выбарцы, а замест гэтага, псеўдавыпадковай выявай абярэм так званы азначнік выбаркі і выкарыстаны яго значэнні.

Напрыклад, калі індэкс вызначыўся за нумарам 3 (тры), то мы бярэм значэнні и , Далей падстаўляем значэння ва ўраўненні вытворных і вызначаем новыя каардынаты. Затым, вызначыўшы каардынаты, мы зноў псеўдавыпадковым чынам вызначаем азначнік выбаркі, падстаўляемы значэнні, якія адпавядаюць азначніку ў раўнанні дзеляў вытворных, па новым вызначаем каардынаты и і г.д. да пазелянення збежнасці. На першы погляд, можа падацца, як гэта ўвогуле можа працаваць, аднак працуе. Праўда варта адзначыць, што не з кожным крокам змяншаецца памылка, але тэндэнцыя безумоўна ёсць.

Якія перавагі стахастычнага градыентнага спуску перад звычайным? У выпадку, калі ў нас памер выбаркі вельмі вялікі і вымяраецца дзясяткамі тысяч значэнняў, то значна прасцей апрацаваць, дапусцім выпадковую тысячу з іх, чым усю выбарку. Вось у гэтым выпадку і запускаецца стахастычны градыентны спуск. У нашым выпадку мы вядома ж вялікай розніцы не заўважым.

Глядзім код.

Код для стахастычнага градыентнага спуску

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

Глядзім уважліва на каэфіцыенты і ловім сябе на пытанні "Як жа так?". У нас атрымаліся іншыя значэнні каэфіцыентаў и . Можа быць стахастычны градыентны спуск знайшоў больш аптымальныя параметры раўнання? Нажаль, не. Дастаткова паглядзець на суму квадратаў адхіленняў і ўбачыць, што пры новых значэннях каэфіцыентаў, памылка большая. Не спяшаемся адчайвацца. Пабудуем графік змены памылкі.

Код для графіка сумы квадратаў адхіленняў пры стахастычным градыентным спуску

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

Графік №5 «Сума квадратаў адхіленняў пры стахастычным градыентным спуску»

Паглядзеўшы на графік, усё становіцца на свае месцы і зараз мы ўсё выправім.

Што ж адбылося? Адбылося наступнае. Калі мы выбіраем выпадковым чынам месяц, то менавіта для абранага месяца наш алгарытм імкнецца паменшыць памылку ў разліку выручкі. Затым выбіраемы іншы месяц і паўтараем разлік, але памылку памяншаем ужо для другога абранага месяца. А зараз успомнім, што ў нас першыя два месяцы істотна адхіляюцца ад лініі раўнання простай лінейнай рэгрэсіі. Гэта значыць, што калі выбіраецца любы з гэтых двух месяцаў, тое памяншаючы памылку кожнага з іх, наш алгарытм сур'ёзна павялічвае памылку па ўсёй выбарцы. Дык што ж рабіць? Адказ просты: трэба зменшыць крок спуску. Бо паменшыўшы крок спуску, памылка гэтак жа перастане "скакаць" то ўгару, то ўніз. Дакладней, памылка "скакаць" не перастане, але будзе гэта рабіць не так жвава:) Праверым.

Код для запуску SGD з меншым крокам

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

Графік №6 "Сума квадратаў адхіленняў пры стахастычным градыентным спуску (80 тыс. крокаў)"

Значэнні каэфіцыентаў палепшыліся, але ўсё роўна не ідэальныя. Гіпатэтычна гэта можна паправіць такім чынам. Выбіраемы, напрыклад, на апошніх 1000 ітэрацыях значэння каэфіцыентаў, з якімі была дапушчаная мінімальная памылка. Праўда нам для гэтага давядзецца запісваць яшчэ і самі значэнні каэфіцыентаў. Мы не будзем гэтага рабіць, а лепш звернем увагу на графік. Ён выглядае гладкім, і памылка нібыта памяншаецца раўнамерна. Насамрэч гэта не так. Паглядзім на першыя 1000 ітэрацый і параўнаем іх з апошнімі.

Код для графіка SGD (першыя 1000 крокаў)

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

Графік №7 "Сума квадратаў адхіленняў SGD (першыя 1000 крокаў)"

Графік №8 "Сума квадратаў адхіленняў SGD (апошнія 1000 крокаў)"

У самым пачатку спуску мы назіраем дастаткова раўнамернае і крутое памяншэнне памылкі. На апошніх ітэрацыях мы бачым, што памылка ходзіць вакол ды каля значэння ў 1,475 і ў некаторыя моманты нават ураўноўваецца гэтаму аптымальнаму значэнню, але потым усё роўна сыходзіць увысь… Паўтаруся, можна запісваць значэнні каэфіцыентаў и , а потым абраць тыя, пры якіх памылка мінімальная. Аднак у нас узнікла праблема больш сур'ёзны: нам прыйшлося зрабіць 80 тыс. крокаў (гл. код), каб атрымаць значэння, блізкія да аптымальных. А гэта, ужо супярэчыць ідэі аб эканоміі часу вылічэнняў пры стахастычным градыентным спуску адносна градыентнага. Што можна паправіць і палепшыць? Не цяжка заўважыць, што на першых ітэрацыях мы ўпэўнена ідзём уніз і, такім чынам, нам варта пакінуць вялікі крок на першых ітэрацыях і па меры пасоўвання наперад крок памяншаць. Мы не будзем гэтага рабіць у гэтым артыкуле - ён і так ужо зацягнуўся. Ахвочыя могуць і самі падумаць, як гэта зрабіць, гэта не складана 🙂

Цяпер выканаем стахастычны градыентны спуск, выкарыстоўваючы бібліятэку numpy (і не будзем спатыкацца аб камяні, якія мы выявілі ранняе)

Код для стахастычнага градыентнага спуску (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

Значэнні атрымаліся амаль такімі ж, як і пры спуску без выкарыстання numpy. Зрэшты, гэта лагічна.

Даведаемся колькі ж часу займалі ў нас стахастычныя градыентныя спускі.

Код для вызначэння часу вылічэнні SGD (80 тыс. крокаў)

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)

Чым далей у лес, тым цямней хмары: зноў "самапісная" формула паказвае лепшы вынік. Усё гэта наводзіць на думкі аб тым, што павінны існаваць яшчэ больш тонкія спосабы выкарыстання бібліятэкі. numpy, якія сапраўды паскараюць аперацыі вылічэнняў. У гэтым артыкуле мы пра іх ужо не даведаемся. Будзе пра што падумаць у вольны час:)

рэзюмуецца

Перад тым як рэзюмаваць, хацелася б адказаць на пытанне, якое хутчэй за ўсё, узнікла ў нашага дарагога чытача. Для чаго, уласна, такія "пакуты" са спускамі, навошта нам хадзіць па гары ўверх і ўніз (пераважна ўніз), каб знайсці запаветную нізіну, калі ў нашых руках такі магутны і просты прыбор, у выглядзе аналітычнага рашэння, які імгненна тэлепартуе нас у патрэбнае месца?

Адказ на гэтае пытанне ляжыць на паверхні. Цяпер мы разбіралі вельмі просты прыклад, у якім праўдзівы адказ залежыць ад адной прыкметы . У жыцці такое сустрэнеш не часта, таму ўявім, што ў нас прыкмет 2, 30, 50 ці больш. Дадамо да гэтага тысячы, а то і дзясяткі тысяч значэнняў для кожнай прыкметы. У гэтым выпадку аналітычнае рашэнне можа не вытрымаць выпрабаванні і даць збой. У сваю чаргу градыентны спуск і яго варыяцыі будуць павольна, але дакладна набліжаць нас да мэты - мінімуму функцыі. А на рахунак хуткасці не хвалюйцеся - мы напэўна яшчэ разбяром спосабы, якія дазволяць нам задаваць і рэгуляваць даўжыню кроку (гэта значыць хуткасць).

А зараз уласна кароткае рэзюмэ.

Па-першае, спадзяюся, што выкладзены ў артыкуле матэрыял, дапаможа пачаткоўцам «дата сайнтыстам» у разуменні таго, як вырашаць раўнанні простай (і не толькі) лінейнай рэгрэсіі.

Па-другое, мы разгледзелі некалькі спосабаў рашэння раўнання. Цяпер, у залежнасці ад сітуацыі, мы можам выбраць той, які лепш за ўсё падыходзіць для вырашэння пастаўленай задачы.

Па-трэцяе, мы ўбачылі сілу дадатковых налад, а менавіта даўжыні кроку градыентнага спуску. Гэтым параметрам нельга грэбаваць. Як было прыкмечана вышэй, з мэтай скарачэння затрат на правядзенне вылічэнняў, даўжыню кроку варта змяняць па ходзе спуску.

Па-чацвёртае, у нашым выпадку, самапісныя функцыі паказалі лепшы часавы вынік вылічэнняў. Верагодна, гэта звязана з не самым прафесійным прымяненнем магчымасцей бібліятэкі. numpy. Але як бы там ні было, выснова напрошваецца наступны. З аднаго боку, часам варта падвяргаць сумневу ўстояныя меркаванні, а з другога - не заўсёды варта ўсё ўскладняць - наадварот часам больш эфектыўна аказваецца больш просты спосаб рашэння задачы. А бо мэта ў нас была разабраць тры падыходу ў рашэнні раўнання простай лінейнай рэгрэсіі, то выкарыстанне "самапісных" функцый нам суцэль хапіла.