🥇Rozwiązywanie równania prostej regresji liniowej

W artykule omówiono kilka sposobów wyznaczania równania matematycznego prostej (sparowanej) linii regresji.

Wszystkie metody rozwiązywania omawianego równania opierają się na metodzie najmniejszych kwadratów. Oznaczmy metody następująco:

Rozwiązanie analityczne
Zejście gradientowe
Stochastyczne zejście gradientowe

Dla każdej metody rozwiązywania równania prostej w artykule podano różne funkcje, które dzielą się głównie na te, które są pisane bez użycia biblioteki numpy i te, które wykorzystują do obliczeń numpy. Uważa się, że umiejętne wykorzystanie numpy obniży koszty obliczeń.

Cały kod podany w artykule jest napisany w tym języku python 2.7 za pomocą Notebook Jupyter. Kod źródłowy i plik z przykładowymi danymi są publikowane na stronie GitHub

Artykuł jest bardziej skierowany zarówno do początkujących, jak i tych, którzy już stopniowo zaczęli opanowywać naukę w bardzo szerokim dziale sztucznej inteligencji - uczeniu maszynowym.

Aby zilustrować materiał, posłużmy się bardzo prostym przykładem.

Przykładowe warunki

Mamy pięć wartości, które charakteryzują zależność Y od X (Tabela nr 1):

Tabela nr 1 „Przykładowe warunki”

Zakładamy, że wartości jest miesiącem roku i — przychody w tym miesiącu. Innymi słowy, przychody zależą od miesiąca roku i - jedyny znak, od którego zależy dochód.

Przykład jest taki sobie, zarówno z punktu widzenia warunkowej zależności przychodów od miesiąca roku, jak i z punktu widzenia liczby wartości – jest ich bardzo mało. Jednak takie uproszczenie pozwoli, jak mówią, wyjaśnić, nie zawsze łatwo, materiał, który przyswajają początkujący. A także prostota liczb pozwoli tym, którzy chcą rozwiązać przykład na papierze bez znacznych kosztów pracy.

Załóżmy, że zależność podaną w przykładzie można dość dobrze przybliżyć równaniem matematycznym prostej (sparowanej) linii regresji o postaci:

gdzie to miesiąc, w którym uzyskano przychód, — przychód odpowiadający miesiącowi, и są współczynnikami regresji oszacowanej linii.

Należy pamiętać, że współczynnik często nazywany nachyleniem lub nachyleniem szacowanej linii; oznacza kwotę, o jaką kiedy to się zmienia .

Oczywiście naszym zadaniem w przykładzie jest wybranie takich współczynników w równaniu и , przy czym odchylenia naszych obliczonych wartości przychodów w poszczególnych miesiącach od prawdziwych odpowiedzi, tj. wartości prezentowane w próbce będą minimalne.

Metoda najmniejszych kwadratów

Zgodnie z metodą najmniejszych kwadratów odchylenie należy obliczyć podnosząc je do kwadratu. Ta technika pozwala uniknąć wzajemnego anulowania odchyleń, jeśli mają one przeciwne znaki. Na przykład, jeśli w jednym przypadku odchylenie wynosi +5 (plus pięć) i w drugim -5 (minus pięć), to suma odchyleń zniesie się i wyniesie 0 (zero). Można nie podnosić odchylenia do kwadratu, ale skorzystać z właściwości modułu, a wtedy wszystkie odchylenia będą dodatnie i będą się kumulować. Nie będziemy szczegółowo rozwodzić się nad tym punktem, ale po prostu wskażemy, że dla wygody obliczeń zwyczajowo podnosi się odchylenie do kwadratu.

Tak wygląda wzór, za pomocą którego wyznaczymy najmniejszą sumę kwadratów odchyleń (błędów):

gdzie jest funkcją przybliżenia prawdziwych odpowiedzi (czyli obliczonego przez nas przychodu),

są prawdziwymi odpowiedziami (przychody podane w próbie),

to indeks próbki (numer miesiąca, w którym określono odchylenie)

Zróżniczkujmy funkcję, zdefiniujmy równania różniczkowe cząstkowe i przejdźmy do rozwiązania analitycznego. Ale najpierw zróbmy krótką wycieczkę po tym, czym jest różniczkowanie i pamiętajmy o geometrycznym znaczeniu pochodnej.

Różnicowanie

Różniczkowanie to operacja znajdowania pochodnej funkcji.

Do czego służy pochodna? Pochodna funkcji charakteryzuje szybkość zmian funkcji i mówi nam o jej kierunku. Jeśli pochodna w danym punkcie jest dodatnia, to funkcja rośnie, w przeciwnym razie funkcja maleje. Im większa wartość pochodnej bezwzględnej, tym większa szybkość zmian wartości funkcji i tym bardziej strome nachylenie wykresu funkcji.

Przykładowo w warunkach kartezjańskiego układu współrzędnych wartość pochodnej w punkcie M(0,0) jest równa + 25 oznacza, że w danym momencie, gdy wartość zostanie przesunięta po prawej stronie konwencjonalną jednostką – wartość wzrasta o 25 jednostek konwencjonalnych. Na wykresie wygląda to na dość gwałtowny wzrost wartości z danego punktu.

Inny przykład. Wartość pochodnej jest równa -0,1 oznacza to, że po przemieszczeniu na jedną jednostkę konwencjonalną, wartość zmniejsza się jedynie o 0,1 jednostki konwencjonalnej. Jednocześnie na wykresie funkcji możemy zaobserwować ledwo zauważalne nachylenie w dół. Rysując analogię z górą, to tak, jakbyśmy bardzo powoli schodzili z góry łagodnym zboczem, w przeciwieństwie do poprzedniego przykładu, gdzie musieliśmy wspinać się na bardzo strome szczyty :)

Zatem po zróżnicowaniu funkcji według szans и , definiujemy równania różniczkowe cząstkowe pierwszego rzędu. Po ustaleniu równań otrzymamy układ dwóch równań, rozwiązując które będziemy mogli wybrać takie wartości współczynników и , dla których wartości odpowiednich pochodnych w danych punktach zmieniają się o bardzo, bardzo małą wielkość, a w przypadku rozwiązania analitycznego nie zmieniają się wcale. Innymi słowy, funkcja błędu przy znalezionych współczynnikach osiągnie minimum, ponieważ wartości pochodnych cząstkowych w tych punktach będą równe zero.

Zatem zgodnie z zasadami różniczkowania równanie pochodnej cząstkowej pierwszego rzędu względem współczynnika przyjmie postać:

Równanie pochodnej cząstkowej pierwszego rzędu w odniesieniu do przyjmie postać:

W rezultacie otrzymaliśmy układ równań, który ma dość proste rozwiązanie analityczne:

rozpocznij{równanie*}
rozpocząć{przypadki}
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
koniec{przypadków}
koniec{równania*}

Przed rozwiązaniem równania załadujmy wstępnie, sprawdźmy czy ładowanie jest prawidłowe i sformatujmy dane.

Ładowanie i formatowanie danych

Należy zaznaczyć, że ze względu na to, że dla rozwiązania analitycznego, a następnie dla gradientowego i stochastycznego opadania gradientu, kod wykorzystamy w dwóch wariantach: z wykorzystaniem biblioteki numpy i bez jego użycia będziemy potrzebować odpowiedniego formatowania danych (patrz kod).

Kod ładowania i przetwarzania danych

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

Wizualizacja

Teraz, po tym jak po pierwsze wczytaliśmy dane, po drugie sprawdziliśmy poprawność wczytania i na koniec sformatowaliśmy dane, przeprowadzimy pierwszą wizualizację. Często stosowaną w tym celu metodą jest wykres par Biblioteka Dno morskie. W naszym przykładzie ze względu na ograniczoną liczbę nie ma sensu korzystać z biblioteki Dno morskie. Będziemy korzystać ze zwykłej biblioteki Biblioteki Matplotu i po prostu spójrz na wykres rozrzutu.

Kod wykresu punktowego

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

Wykres nr 1 „Zależność przychodów od miesiąca roku”

Rozwiązanie analityczne

Użyjmy najpopularniejszych narzędzi w pyton i rozwiąż układ równań:

rozpocznij{równanie*}
rozpocząć{przypadki}
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
koniec{przypadków}
koniec{równania*}

Zgodnie z regułą Cramera znajdziemy wyznacznik ogólny, a także wyznaczniki według i , po czym dzieląc wyznacznik przez do ogólnego wyznacznika - znajdź współczynnik , podobnie znajdujemy współczynnik .

Kod rozwiązania analitycznego

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

Oto co otrzymaliśmy:

Znaleziono więc wartości współczynników, ustalono sumę kwadratów odchyleń. Narysujmy linię prostą na histogramie rozproszenia zgodnie ze znalezionymi współczynnikami.

Kod linii regresji

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

Wykres nr 2 „Odpowiedzi prawidłowe i obliczone”

Możesz spojrzeć na wykres odchyleń dla każdego miesiąca. W naszym przypadku nie wyciągniemy z tego żadnej istotnej wartości praktycznej, ale zaspokoimy naszą ciekawość, jak dobrze proste równanie regresji liniowej charakteryzuje zależność przychodów od miesiąca roku.

Kod wykresu odchyleń

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

Wykres nr 3 „Odchylenia, %”

Nie idealnie, ale wykonaliśmy swoje zadanie.

Napiszmy funkcję, która w celu wyznaczenia współczynników и korzysta z biblioteki numpya dokładniej napiszemy dwie funkcje: jedną wykorzystując macierz pseudoodwrotną (w praktyce nie jest to zalecane, gdyż proces jest obliczeniowo skomplikowany i niestabilny), drugą wykorzystując równanie macierzowe.

Kod rozwiązania analitycznego (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

Porównajmy czas spędzony na wyznaczaniu współczynników и zgodnie z 3 przedstawionymi metodami.

Kod do obliczania czasu obliczeń

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)

Przy małej ilości danych wychodzi naprzeciw „samodzielnie napisanej” funkcji, która znajduje współczynniki metodą Cramera.

Teraz możesz przejść do innych sposobów znajdowania współczynników и .

Zejście gradientowe

Najpierw zdefiniujmy, czym jest gradient. Mówiąc najprościej, gradient to odcinek wskazujący kierunek maksymalnego wzrostu funkcji. Przez analogię do wspinaczki górskiej, gdzie pochyłe ściany są najbardziej stromym podejściem na szczyt góry. Rozwijając przykład z górą pamiętamy, że tak naprawdę potrzebujemy jak najbardziej stromego zjazdu, aby jak najszybciej dotrzeć na nizinę, czyli do minimum – miejsca, w którym funkcja nie rośnie ani nie maleje. W tym momencie pochodna będzie równa zeru. Dlatego nie potrzebujemy gradientu, ale antygradientu. Aby znaleźć antygradient, wystarczy pomnożyć gradient przez -1 (minus jeden).

Zwróćmy uwagę na to, że funkcja może mieć kilka minimów i po zejściu do jednego z nich za pomocą zaproponowanego poniżej algorytmu nie uda nam się znaleźć kolejnego minimum, które może być mniejsze od znalezionego. Spokojnie, nam to nie grozi! W naszym przypadku mamy do czynienia z jednym minimum, ponieważ jest to nasza funkcja na wykresie jest parabola foremna. A jak wszyscy powinniśmy dobrze wiedzieć z naszego szkolnego kursu matematyki, parabola ma tylko jedno minimum.

Po tym jak dowiedzieliśmy się po co nam gradient i że gradient jest odcinkiem, czyli wektorem o podanych współrzędnych, które mają dokładnie takie same współczynniki и możemy wdrożyć opadanie gradientowe.

Zanim zacznę, sugeruję przeczytanie kilku zdań na temat algorytmu opadania:

Współrzędne współczynników wyznaczamy w sposób pseudolosowy и . W naszym przykładzie zdefiniujemy współczynniki bliskie zeru. Jest to powszechna praktyka, ale każdy przypadek może mieć swoją własną praktykę.
Ze współrzędnych odejmij w punkcie wartość pochodnej cząstkowej pierwszego rzędu . Jeśli więc pochodna jest dodatnia, to funkcja rośnie. Odejmując zatem wartość pochodnej, będziemy poruszać się w przeciwnym kierunku wzrostu, czyli w kierunku opadania. Jeśli pochodna jest ujemna, to funkcja w tym miejscu maleje i odejmując wartość pochodnej przesuwamy się w kierunku malejącym.
Podobną operację wykonujemy ze współrzędnymi : odejmij wartość pochodnej cząstkowej w punkcie .
Aby nie przeskoczyć minimum i nie polecieć w głęboką przestrzeń kosmiczną, konieczne jest ustawienie wielkości kroku w kierunku opadania. W sumie można by napisać cały artykuł o tym jak prawidłowo ustawić stopień i jak go zmienić w trakcie schodzenia, aby obniżyć koszty obliczeniowe. Ale teraz przed nami nieco inne zadanie, a wielkość kroku ustalimy za pomocą naukowej metody „szturchania” lub, jak to się mówi w potocznym języku, empirycznie.
Gdy już znajdziemy się pod podanymi współrzędnymi и odejmij wartości pochodnych, otrzymamy nowe współrzędne и . Wykonujemy kolejny krok (odejmowanie) już od obliczonych współrzędnych. I tak cykl zaczyna się od nowa, aż do osiągnięcia wymaganej zbieżności.

Wszystko! Teraz jesteśmy gotowi wyruszyć na poszukiwania najgłębszego wąwozu Rowu Mariańskiego. Zacznijmy.

Kod opadania gradientowego

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

Zanurkowaliśmy na samo dno rowu Mariana i tam znaleźliśmy te same wartości współczynników и , czyli dokładnie to, czego można było się spodziewać.

Zanurkujmy jeszcze raz, tylko tym razem nasz pojazd głębinowy będzie wypełniony innymi technologiami, a mianowicie biblioteką numpy.

Kod opadania gradientu (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

Wartości współczynników и niezmienny.

Przyjrzyjmy się, jak zmieniał się błąd podczas opadania gradientu, czyli jak zmieniała się suma kwadratów odchyleń z każdym krokiem.

Kod do wykreślania sum kwadratów odchyleń

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

Wykres nr 4 „Suma kwadratów odchyleń podczas opadania gradientu”

Na wykresie widzimy, że z każdym krokiem błąd maleje, a po określonej liczbie iteracji obserwujemy niemal poziomą linię.

Na koniec oszacujmy różnicę w czasie wykonania kodu:

Kod określający czas obliczania opadania gradientowego

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)

Być może robimy coś źle, ale znowu jest to prosta, „domowo napisana” funkcja, która nie korzysta z biblioteki numpy przewyższa czas obliczeń funkcji przy użyciu biblioteki numpy.

Ale nie stoimy w miejscu, ale zmierzamy w kierunku zbadania innego ekscytującego sposobu rozwiązania prostego równania regresji liniowej. Poznać!

Stochastyczne zejście gradientowe

Aby szybko zrozumieć zasadę działania stochastycznego opadania gradientowego, lepiej określić jego różnice w stosunku do zwykłego opadania gradientowego. My, w przypadku gradientu, w równaniach pochodnych и wykorzystano sumy wartości wszystkich cech i prawdziwych odpowiedzi dostępnych w próbie (czyli sumy wszystkich и ). W stochastycznym spadku gradientowym nie będziemy wykorzystywać wszystkich wartości występujących w próbce, lecz pseudolosowo wybieramy tzw. indeks próbki i wykorzystujemy jego wartości.

Na przykład, jeśli indeks zostanie określony jako numer 3 (trzy), wówczas przyjmujemy wartości и , następnie podstawiamy wartości do równań pochodnych i wyznaczamy nowe współrzędne. Następnie, po ustaleniu współrzędnych, ponownie pseudolosowo wyznaczamy indeks próbki, podstawiamy wartości odpowiadające indeksowi do równań różniczkowych cząstkowych i wyznaczamy współrzędne w nowy sposób и itp. aż zbieżność zmieni kolor na zielony. Na pierwszy rzut oka może się wydawać, że to w ogóle nie działa, a jednak tak jest. Co prawda warto zauważyć, że błąd nie maleje z każdym krokiem, ale na pewno jest tendencja.

Jakie są zalety stochastycznego opadania gradientowego w porównaniu z konwencjonalnym? Jeśli wielkość naszej próby jest bardzo duża i mierzona w dziesiątkach tysięcy wartości, wówczas znacznie łatwiej jest przetworzyć, powiedzmy, losowy tysiąc z nich niż całą próbę. W tym miejscu wchodzi w grę stochastyczne zejście gradientowe. W naszym przypadku oczywiście dużej różnicy nie zauważymy.

Spójrzmy na kod.

Kod stochastycznego spadku gradientu

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

Przyglądamy się uważnie współczynnikom i przyłapujemy się na zadawaniu pytania: „Jak to możliwe?” Otrzymaliśmy inne wartości współczynników и . Może stochastyczne zejście gradientowe znalazło bardziej optymalne parametry równania? Niestety nie. Wystarczy spojrzeć na sumę kwadratów odchyleń i zobaczyć, że przy nowych wartościach współczynników błąd jest większy. Nie spieszymy się do rozpaczy. Zbudujmy wykres zmiany błędu.

Kod do wykreślania sumy kwadratów odchyleń w stochastycznym spadku gradientu

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

Wykres nr 5 „Suma kwadratów odchyleń podczas opadania w gradiencie stochastycznym”

Patrząc na harmonogram, wszystko się układa i teraz wszystko naprawimy.

Więc co się stało? Wydarzyło się co następuje. Kiedy losowo wybieramy miesiąc, to właśnie dla wybranego miesiąca nasz algorytm stara się zredukować błąd w obliczaniu przychodów. Następnie wybieramy kolejny miesiąc i powtarzamy obliczenia, ale zmniejszamy błąd dla drugiego wybranego miesiąca. Pamiętajmy teraz, że pierwsze dwa miesiące znacznie odbiegają od linii prostego równania regresji liniowej. Oznacza to, że w przypadku wybrania któregokolwiek z tych dwóch miesięcy, zmniejszając błąd każdego z nich, nasz algorytm poważnie zwiększa błąd dla całej próby. Co więc zrobić? Odpowiedź jest prosta: musisz zmniejszyć stopień zejścia. Przecież zmniejszając stopień zejścia, błąd przestanie także „skakać” w górę i w dół. A raczej błąd „skakania” nie ustanie, ale nie zrobi tego tak szybko :) Sprawdźmy.

Kod uruchamiający SGD z mniejszymi przyrostami

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

Wykres nr 6 „Suma kwadratów odchyleń podczas opadania w gradiencie stochastycznym (80 tys. kroków)”

Współczynniki poprawiły się, ale nadal nie są idealne. Hipotetycznie można to skorygować w ten sposób. Wybieramy np. w ostatnich 1000 iteracjach wartości współczynników, przy których popełniono błąd minimalny. To prawda, że w tym celu będziemy musieli również zapisać wartości samych współczynników. Nie zrobimy tego, ale raczej zwrócimy uwagę na harmonogram. Wygląda gładko, a błąd wydaje się równomiernie zmniejszać. W rzeczywistości nie jest to prawdą. Przyjrzyjmy się pierwszym 1000 iteracjom i porównajmy je z ostatnimi.

Kod wykresu SGD (pierwsze 1000 kroków)

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

Wykres nr 7 „Suma kwadratów odchyleń SGD (pierwsze 1000 kroków)”

Wykres nr 8 „Suma kwadratów odchyleń SGD (ostatnie 1000 kroków)”

Już na samym początku opadania obserwujemy dość równomierny i stromy spadek błędu. W ostatnich iteracjach widzimy, że błąd oscyluje wokół wartości 1,475 i w niektórych momentach jest nawet równy tej wartości optymalnej, ale potem wciąż rośnie... Powtarzam, możesz zapisać wartości współczynniki и , a następnie wybierz te, dla których błąd jest minimalny. Mieliśmy jednak poważniejszy problem: musieliśmy wykonać 80 tysięcy kroków (patrz kod), aby uzyskać wartości zbliżone do optymalnych. A to już jest sprzeczne z ideą oszczędzania czasu obliczeń dzięki stochastycznemu zejściu gradientu w stosunku do zejścia gradientu. Co można poprawić i ulepszyć? Nietrudno zauważyć, że w pierwszych iteracjach pewnie schodzimy w dół, dlatego w pierwszych iteracjach powinniśmy zostawić duży krok i w miarę postępów go zmniejszać. Nie będziemy tego robić w tym artykule – jest już za długi. Kto chce, może sam pomyśleć, jak to zrobić, nie jest to trudne :)

Teraz przeprowadźmy stochastyczne opadanie w gradiencie, korzystając z biblioteki numpy (i nie potykajmy się o kamienie, które zidentyfikowaliśmy wcześniej)

Kod stochastycznego opadania gradientu (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

Wartości okazały się prawie takie same jak przy zejściu bez użycia numpy. Jest to jednak logiczne.

Dowiedzmy się, ile czasu zajęło nam zejście w gradiencie stochastycznym.

Kod określający czas obliczenia SGD (80 tys. kroków)

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 dalej w las, tym ciemniejsze chmury: ponownie „samodzielnie napisany” wzór daje najlepszy wynik. Wszystko to sugeruje, że muszą istnieć jeszcze bardziej subtelne sposoby korzystania z biblioteki numpy, które naprawdę przyspieszają operacje obliczeniowe. W tym artykule nie dowiemy się o nich. Będzie o czym myśleć w wolnym czasie :)

Podsumowujemy

Zanim zacznę podsumowywać, chciałbym odpowiedzieć na pytanie, które najprawdopodobniej padło od naszego drogiego czytelnika. Po co w ogóle takie „męki” ze zjazdami, dlaczego musimy chodzić po górach (głównie w dół), aby znaleźć cenną nizinę, skoro mamy w rękach tak potężne i proste urządzenie, w formę rozwiązania analitycznego, które błyskawicznie teleportuje nas we właściwe miejsce?

Odpowiedź na to pytanie leży na powierzchni. Teraz przyjrzeliśmy się bardzo prostemu przykładowi, w którym prawdziwa odpowiedź brzmi zależy od jednego znaku . Nie widujesz tego często w życiu, więc wyobraźmy sobie, że mamy 2, 30, 50 lub więcej znaków. Dodajmy do tego tysiące, a nawet dziesiątki tysięcy wartości dla każdego atrybutu. W takim przypadku rozwiązanie analityczne może nie wytrzymać testu i zakończyć się niepowodzeniem. Z kolei opadanie gradientowe i jego zmiany będą powoli, ale pewnie przybliżać nas do celu - minimum funkcji. I nie martw się o prędkość – prawdopodobnie przyjrzymy się sposobom, które pozwolą nam ustawić i regulować długość kroku (czyli prędkość).

A teraz właściwe krótkie podsumowanie.

Po pierwsze, mam nadzieję, że materiał przedstawiony w artykule pomoże początkującym „data science” w zrozumieniu sposobu rozwiązywania prostych (i nie tylko) równań regresji liniowej.

Po drugie, przyjrzeliśmy się kilku sposobom rozwiązania równania. Teraz w zależności od sytuacji możemy wybrać ten, który najlepiej rozwiąże problem.

Po trzecie, zobaczyliśmy siłę dodatkowych ustawień, a mianowicie długość kroku opadania gradientu. Nie można pominąć tego parametru. Jak wspomniano powyżej, aby obniżyć koszty obliczeń, należy zmieniać długość kroku podczas schodzenia.

Po czwarte, w naszym przypadku najlepsze wyniki czasowe obliczeń wykazały funkcje „samodzielnie napisane”. Prawdopodobnie wynika to z niezbyt profesjonalnego wykorzystania możliwości biblioteki numpy. Tak czy inaczej, następujący wniosek nasuwa się sam. Z jednej strony czasami warto kwestionować utarte opinie, z drugiej nie zawsze warto wszystko komplikować – wręcz przeciwnie, czasami prostszy sposób rozwiązania problemu jest skuteczniejszy. A ponieważ naszym celem była analiza trzech podejść do rozwiązania prostego równania regresji liniowej, wystarczyło nam użycie „samodzielnie napisanych” funkcji.