Rozwiązujemy równanie prostej regresji liniowej

Artykuł omawia kilka sposobów wyznaczania matematycznego równania linii prostej (parowej) regresji.

Wszystkie omawiane tutaj sposoby rozwiązania równania opierają się na metodzie najmniejszych kwadratów. Oznaczymy metody w następujący sposób:

• Rozwiązanie analityczne
• Spadek gradientu
• Stochastyczny spadek gradientu

Dla każdego ze sposobów rozwiązania równania linii, w artykule przedstawiono różne funkcje, które głównie dzielą się na te, które są napisane bez użycia biblioteki NumPy oraz te, które do obliczeń stosują NumPy. Uważa się, że umiejętne wykorzystanie NumPy pozwoli na zmniejszenie kosztów obliczeń.

Cały kod przedstawiony w artykule został napisany w języku python 2.7 z użyciem Jupyter Notebook. Źródłowy kod i plik z danymi próbki są dostępne na GitHubie

Artykuł jest w większym stopniu skierowany zarówno do początkujących, jak i do tych, którzy już zaczęli oswajać się z bardzo obszernym działem sztucznej inteligencji - uczenia maszynowego.

Aby zilustrować materiał, użyjemy bardzo prostego przykładu.

Warunki przykładu

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

Tabela nr 1 "Warunki przykładu"

Zakładamy, że wartości to miesiąc roku, a to przychody w tym miesiącu. Innymi słowy, przychody zależą od miesiąca roku, a to jedyny czynnik, od którego zależą przychody.

Przykład jest dość uproszczony, 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. Niemniej jednak takie uproszczenie pozwoli, że jak to się mówi, w zrozumiały sposób wyjaśni materiał, który nowicjusze przyswajają często z trudem. A także prostota cyfr pozwoli bez dużych nakładów czasu, chętnym, rozwiązać przykład „na papierze”.

Załóżmy, że przedstawiona w przykładzie zależność może być dostatecznie dobrze aproksymowana matematycznym równaniem linii prostej (parowej) regresji w postaci:

gdzie - to miesiąc, w którym uzyskano przychody, - przychody odpowiadające miesiącowi, i - współczynniki regresji oszacowanej linii.

Zauważmy, że współczynnik często nazywany jest współczynnikiem kierunkowym lub gradientem oszacowanej linii; reprezentuje to wartość, o którą zmieni się przy zmianie .

Oczywiście, nasza zadanie w przykładzie polega na dobraniu w równaniu takich współczynników i , przy których odchylenia naszych obliczonych wartości przychodów w miesiącach od prawdziwych odpowiedzi, tzn. wartości przedstawionych w próbie, będą minimalne.

Metoda najmniejszych kwadratów

Zgodnie z metodą najmniejszych kwadratów, odchylenie należy obliczać, podnosząc je do kwadratu. Taki zabieg pozwala uniknąć wzajemnego wygaszania odchyleń w przypadku, gdy mają one przeciwne znaki. Na przykład, jeśli w jednym przypadku odchylenie wynosi +5 (plus pięć), a w innym -5 (minus pięć), to suma odchyleń wzajemnie się wygasi i wyniesie 0 (zero). Można nie podnosić odchylenia do kwadratu, a skorzystać z własności modułu, a wtedy wszystkie nasze odchylenia będą dodatnie i będą się kumulować. Nie będziemy się szczegółowo nad tym zatrzymywać, a jedynie zaznaczymy, że dla wygody obliczeń przyjęto podnosić odchylenie do kwadratu.

Oto jak wygląda wzór, za pomocą którego określimy najmniejszą sumę kwadratów odchyleń (błędu):

gdzie — to funkcja aproksymacji prawdziwych odpowiedzi (czyli obliczona przez nas przychód),

— to prawdziwe odpowiedzi (przekazane w próbie przychody),

— to indeks próby (numer miesiąca, w którym zachodzi określenie odchylenia)

Zróżniczkujemy funkcję, określimy równania pochodnych cząstkowych i będziemy gotowi przejść do analitycznego rozwiązania. Ale najpierw przeprowadzimy mały wstęp na temat tego, czym jest różniczkowanie, i przypomnimy sobie geometryczne znaczenie pochodnej.

Różniczkowanie

Różniczkowanie to operacja polegająca na znajdowaniu pochodnej funkcji.

Po co potrzebna jest pochodna? Pochodna funkcji określa szybkość zmiany funkcji i wskazuje nam jej kierunek. Jeśli pochodna w danym punkcie jest dodatnia, to funkcja rośnie, w przeciwnym przypadku — funkcja maleje. Im wyższa wartość pochodnej w module, tym wyższa szybkość zmiany wartości funkcji oraz bardziej stromy kąt nachylenia wykresu funkcji.

Na przykład, w warunkach kartezjańskiego układu współrzędnych, wartość pochodnej w punkcie M(0,0) równa +25 oznacza, że w danym punkcie, przy przesunięciu wartości w prawo o jednostkową wartość, wartość wzrasta o 25 jednostek warunkowych. Na wykresie wygląda to jak wystarczająco stromy kąt wzrostu wartości. od danego punktu.

Inny przykład. Wartość pochodnej równa -0,1 oznacza, że przy przesunięciu o jedną jednostkę warunkową, wartość spada jedynie o 0,1 jednostki warunkowej. Przy tym, na wykresie funkcji, możemy obserwować ledwo zauważalny spadek. Porównując to do góry, jakby bardzo powoli schodziliśmy po łagodnym zboczu, w przeciwieństwie do poprzedniego przykładu, gdzie musieliśmy zdobywać bardzo strome szczyty :)

Tak więc, przeprowadzając różniczkowanie funkcji według współczynników i , określimy równania pierwszej pochodnej. Po określeniu równań, otrzymamy system dwóch równań, rozwiązując który będziemy mogli dobrać takie wartości współczynników i , przy których wartości odpowiadających 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 pierwszych pochodnych w tych punktach będą równe zeru.

Zatem, według zasad różniczkowania równanie pierwszej pochodnej według współczynnika przyjmie postać:

również równanie pierwszej pochodnej według przyjmie postać:

W efekcie otrzymaliśmy system równań, który ma dość proste rozwiązanie analityczne:

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*}

Zanim rozwiążemy równanie, najpierw załadujemy, sprawdzimy poprawność załadunku i sformatujemy dane.

Ładowanie i formatowanie danych

Należy zauważyć, że w związku z tym, że dla rozwiązania analitycznego, a następnie dla gradientowego i stochastycznego spadku gradientu, zastosujemy kod w dwóch wariantach: z wykorzystaniem biblioteki NumPy i bez niej, 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 pierwsze, że załadowaliśmy dane, po drugie, sprawdziliśmy poprawność załadunku i wreszcie sformatowaliśmy dane, przeprowadzimy pierwszą wizualizację. Często do tego celu używa się metody pairplot biblioteki Seaborn. W naszym przykładzie, z uwagi na ograniczoną ilość cyfr, nie ma sensu stosować biblioteki. Seaborn. Skorzystamy z tradycyjnej biblioteki. Matplotlib i skupimy się tylko na wykresie punktowym.

Kod wykresu punktowego.

print 'Wykres nr 1 "Zależność przychodu od miesiąca roku"'

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ść przychodu od miesiąca roku”.

Rozwiązanie analityczne

Skorzystamy z najbardziej podstawowych narzędzi w python i rozwiążemy system równań:

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*}

Zgodnie z regułą Cramera, znajdziemy ogólny wyznacznik, a także wyznaczniki według i według , po czym, dzieląc wyznacznik według przez ogólny wyznacznik — znajdziemy współczynnik , podobnie znajdziemy 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 udało nam się osiągnąć:

Tak, wartości współczynników zostały znalezione, suma kwadratów odchyleń została ustalona. Narysujemy na histogramie rozrzutu linię zgodnie z 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 rzeczywiste i obliczone”.

Można zobaczyć wykres odchyleń za każdy miesiąc. W naszym przypadku nie wyciągniemy żadnych znaczących praktycznych wniosków, ale zaspokoimy ciekawość, jak dobrze równanie prostej regresji liniowej opisuje zależność przychodu 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 nasze zadanie zostało wykonane.

Napiszemy funkcję, która do określenia współczynników i używa biblioteki NumPy, a dokładniej — napiszemy dwie funkcje: jedną z wykorzystaniem pseudoodwrotnej macierzy (niezalecana w praktyce, ponieważ proces jest obliczeniowo skomplikowany i niestabilny), drugą z wykorzystaniem równania macierzowego.

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ównamy czas, który został poświęcony na określenie współczynników i , zgodnie z trzema przedstawionymi metodami.

Kod do obliczania czasu obliczeń.

print ' 33[1m' + ' 33[4m' + "Czas wykonania obliczeń współczynników bez użycia biblioteki NumPy:" + ' 33[0m'
% timeit ab_us = Kramer_method(x_us,y_us)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Czas wykonania obliczeń współczynników z użyciem pseudoodwrotnej macierzy:" + ' 33[0m'
%timeit ab_np = pseudoinverse_matrix(x_np, y_np)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Czas wykonania obliczeń współczynników z użyciem równania macierzowego:" + ' 33[0m'
%timeit ab_np = matrix_equation(x_np, y_np)

Na niewielkiej ilości danych na pierwszy plan wychodzi "napisana ręcznie" funkcja, która znajduje współczynniki metodą Cramera.

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

Spadek gradientu

Na początek zdefiniujmy, czym jest gradient. Prosto mówiąc, gradient to odcinek, który wskazuje kierunek maksymalnego wzrostu funkcji. W odniesieniu do wspinaczki w górę, kierunek, w którym spogląda gradient, to najstromszy wzrost w kierunku szczytu góry. Rozwijając przykład z górą, warto pamiętać, że tak naprawdę potrzebujemy najstromniejszego zejścia, aby jak najszybciej osiągnąć dolinę, czyli minimum – miejsce, w którym funkcja ani nie rośnie, ani nie maleje. W tym miejscu pochodna będzie równa zeru. W związku z tym potrzebujemy nie gradientu, a antygradientu. Aby znaleźć antygradient, wystarczy pomnożyć gradient przez -1 (minus jeden).

Zwróćmy uwagę, że funkcja może mieć kilka minimów, a po zejściu w jedno z nich według opisanego dalej algorytmu, nie będziemy w stanie znaleźć innego minimum, które może znajdować się poniżej znalezionego. Nie martwmy się, to nam nie grozi! W naszym przypadku mamy do czynienia z jedynym minimum, ponieważ nasza funkcja na wykresie przedstawia zwykłą parabola. Jak wszyscy powinniśmy wiedzieć z szkolnego kursu matematyki – parabola ma tylko jedno minimum.

Po tym, jak ustaliliśmy, do czego potrzebny nam był gradient oraz to, że gradient to odcinek, a więc wektor z określonymi współrzędnymi, które są dokładnie tymi współczynnikami i możemy zrealizować spadek gradientowy.

Przed uruchomieniem proponuję przeczytać dosłownie kilka zdań na temat algorytmu spadku:

• Określamy współrzędne współczynników w sposób pseudo-losowy. i . W naszym przykładzie będziemy ustalać współczynniki w pobliżu zera. Jest to powszechna praktyka, jednak w każdym przypadku może być przewidziana inna praktyka.
• Od współrzędnej odejmujemy wartość pierwszej pochodnej w punkcie. . Więc jeśli pochodna będzie dodatnia, funkcja rośnie. W związku z tym, odejmując wartość pochodnej, będziemy poruszać się w przeciwną stronę wzrostu, czyli w kierunku spadku. Jeśli pochodna jest ujemna, to znaczy, że funkcja w tym punkcie maleje, a odejmując wartość pochodnej, poruszamy się w kierunku spadku.
• Przeprowadzamy analogiczną operację z współrzędną : odejmujemy wartość częściowej pochodnej w punkcie .
• Aby nie przegapić minimum i nie odlecieć w daleki kosmos, należy ustalić rozmiar kroku w kierunku spadku. Ogólnie rzecz biorąc, można napisać cały artykuł o tym, jak poprawnie ustalić krok i jak go zmieniać w trakcie spadku, aby zredukować koszty obliczeń. Jednak teraz mamy nieco inne zadanie, i metodą naukowego "tykania" lub jak mówią potocznie, empiryczną drogą, ustalimy rozmiar kroku.
• Po tym, jak z zadanych współrzędnych i odejmujemy wartości pochodnych, uzyskujemy nowe współrzędne i . Robimy następny krok (odejmowanie), już z obliczonych współrzędnych. I tak cykl rozpoczyna się na nowo, aż osiągniemy wymaganą zbieżność.

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

Kod dla gradientowego spadku

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

Zanurzyliśmy się na samo dno Rowu Mariańskiego i tam odkryliśmy te same wartości współczynników i , co właściwie było do przewidzenia.

Wykonamy jeszcze jedno zanurzenie, tylko tym razem technologia naszego głębokowodnego aparatu będzie inna, a mianowicie biblioteka NumPy.

Kod dla gradientowego spadku (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 i pozostają niezmienne.

Zobaczymy, jak zmieniała się błąd podczas gradientowego spadku, to znaczy jak zmieniała się suma kwadratów odchyleń z każdym krokiem.

Kod dla wykresu sum kwadratów odchyleń

print 'Wykres nr 4 "Suma kwadratów odchyleń krok po kroku"'
plt.plot(range(len(list_parametres_gradient_descence[1])), list_parametres_gradient_descence[1], color='red', lw=3)
plt.xlabel('Kroki (Iteracja)', size=16)
plt.ylabel('Suma kwadratów odchyleń', size=16)
plt.show()

Wykres nr 4 „Suma kwadratów odchyleń podczas gradientowego spadku”

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

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

Kod do określenia czasu obliczeń dla gradientowego spadku

print ' 33[1m' + ' 33[4m' + "Czas wykonania gradientowego spadku bez użycia biblioteki 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' + "Czas wykonania gradientowego spadku z użyciem biblioteki NumPy:" + ' 33[0m'
%timeit list_parametres_gradient_descence = gradient_descent_numpy(x_np,y_np,l=0.1,tolerance=0.000000000001)

Możliwe, że robimy coś nie tak, ale znowu prosta „samodzielnie napisana” funkcja, która nie używa biblioteki NumPy przekracza czasowo funkcję wykorzystującą bibliotekę NumPy.

Ale nie stoimy w miejscu, lecz zmierzamy w stronę poznania kolejnego fascynującego sposobu rozwiązania równania prostej regresji liniowej. Oto!

Stochastyczny spadek gradientu

Aby szybciej zrozumieć zasadę działania stochastycznego gradientowego spadku, lepiej jest określić jego różnice w porównaniu do zwykłego gradientowego spadku. W przypadku gradientowego spadku w równaniach pochodnych i używaliśmy sum wartości wszystkich cech i rzeczywistych odpowiedzi dostępnych w próbie (to znaczy sum wszystkich i ). W stochastycznym gradientowym spadku nie będziemy korzystać ze wszystkich wartości dostępnych w próbie, lecz zamiast tego losowo wybierzemy tzw. indeks próby i wykorzystamy jego wartości.

Na przykład, jeśli indeks został ustalony jako numer 3 (trzy), to bierzemy wartości i , następnie podstawiamy wartości do równań pochodnych i określamy nowe współrzędne. Następnie, określając współrzędne, znowu losowo wybieramy indeks próby, podstawiamy wartości odpowiadające indeksowi do równań pochodnych cząstkowych, a następnie na nowo określamy współrzędne i i tak dalej, aż do osiągnięcia zbieżności. Na pierwszy rzut oka może się wydawać, że to w ogóle może działać, ale działa. Warto zaznaczyć, że nie z każdym krokiem błąd się zmniejsza, ale tendencja zdecydowanie istnieje.

Jakie są zalety stochastycznego spadku gradientowego w porównaniu do zwykłego? Jeśli mamy bardzo dużą próbkę, mierzona dziesiątkami tysięcy wartości, znacznie łatwiej jest przetworzyć na przykład losową tysiąc z nich, zamiast całej próbki. W takim przypadku uruchamiany jest stochastyczny spadek gradientowy. W naszym przypadku oczywiście nie zauważymy dużej różnicy.

Patrzymy na kod.

Kod dla stochastycznego spadku gradientowego

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

Dokładnie przyglądamy się współczynnikom i łapiemy się na pytaniu „Jak to możliwe?”. Otrzymaliśmy inne wartości współczynników. i . Może stochastyczny spadek gradientowy znalazł 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 spieszmy się z rozczarowaniem. Zbudujemy wykres zmiany błędu.

Kod dla wykresu sumy kwadratów odchyleń przy stochastycznym spadku gradientowym

print 'Wykres nr 5 "Suma kwadratów odchyleń krok po kroku"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1])), list_parametres_stoch_gradient_descence[1], color='red', lw=2)
plt.xlabel('Kroki (Iteracja)', size=16)
plt.ylabel('Suma kwadratów odchyleń', size=16)
plt.show()

Wykres nr 5 „Suma kwadratów odchyleń przy stochastycznym spadku gradientowym”

Patrząc na wykres, wszystko staje się jasne i teraz wszystko poprawimy.

Cóż, co się wydarzyło? Wydarzyło się to, że gdy wybieramy losowo miesiąc, to właśnie dla wybranego miesiąca nasz algorytm stara się zmniejszyć błąd w obliczeniach przychodu. Następnie wybieramy inny miesiąc i powtarzamy obliczenia, ale błąd zmniejszamy już dla drugiego wybranego miesiąca. A teraz przypomnijmy sobie, że nasze pierwsze dwa miesiące znacznie odbiegają od linii równania prostej regresji liniowej. To oznacza, że gdy zostaje wybrany którykolwiek z tych dwóch miesięcy, to zmniejszając błąd każdego z nich, nasz algorytm poważnie zwiększa błąd w całej próbce. Co więc robić? Odpowiedź jest prosta: trzeba zmniejszyć krok spadku. W końcu zmniejszając krok spadku, błąd również przestanie „skakać” to w górę, to w dół. Dokładniej, błąd może nadal „skakać”, ale nie tak żwawo :) Sprawdźmy.

Kod dla uruchomienia SGD z mniejszym krokiem

# запустим функцию, уменьшив шаг в 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ń przy stochastycznym spadku gradientowym (80 tys. kroków)”

Wartości współczynników poprawiły się, ale wciąż nie są idealne. Teoretycznie można to naprawić w ten sposób. Wybieramy, na przykład, z ostatnich 1000 iteracji wartości współczynników, dla których popełniono minimalny błąd. Jednak do tego musielibyśmy także zanotować same wartości współczynników. Nie zrobimy tego, a raczej zwrócimy uwagę na wykres. Wygląda gładko, a błąd zdaje się zmniejszać równomiernie. W rzeczywistości nie jest jednak tak. Spójrzmy na pierwsze 1000 iteracji i porównajmy je z ostatnimi.

Kod do wykresu SGD (pierwsze 1000 kroków)

print 'Wykres nr 7 "Suma kwadratów odchyleń krok po kroku. Pierwsze 1000 iteracji"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][:1000])), 
         list_parametres_stoch_gradient_descence[1][:1000], color='red', lw=2)
plt.xlabel('Kroki (Iteracja)', size=16)
plt.ylabel('Suma kwadratów odchyleń', size=16)
plt.show()

print 'Wykres nr 7 "Suma kwadratów odchyleń krok po kroku. Ostatnie 1000 iteracji"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][-1000:])), 
         list_parametres_stoch_gradient_descence[1][-1000:], color='red', lw=2)
plt.xlabel('Kroki (Iteracja)', size=16)
plt.ylabel('Suma kwadratów odchyleń', 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)»

Na początku spadku obserwujemy dość równomierne i strome zmniejszenie błędu. W ostatnich iteracjach widzimy, że błąd oscyluje w okolicach wartości 1,475 i w niektórych momentach nawet osiąga tę optymalną wartość, ale potem znowu idzie w górę… Powtórzę, można zapisywać wartości współczynników, i , a następnie wybrać te, dla których błąd jest minimalny. Jednak pojawił się poważniejszy problem: musieliśmy wykonać 80 tysięcy kroków (patrz kod), aby uzyskać wartości zbliżone do optymalnych. A to już stoi w sprzeczności z ideą oszczędności czasu obliczeń przy stochastycznym spadku gradientu w porównaniu do klasycznego. Co można naprawić i poprawić? Nie trudno zauważyć, że w pierwszych iteracjach pewnie idziemy w dół i w związku z tym warto pozostawić duży krok na początku, a w miarę postępu zmniejszać krok. Nie zrobimy tego w tym artykule — już się przeciąga. Chętni mogą sami pomyśleć, jak to zrobić, to nie jest trudne 🙂

Teraz przeprowadzimy stochastyczny spadek gradientu, wykorzystując bibliotekę NumPy (i nie będziemy potykać się o kamienie, które wcześniej zidentyfikowaliśmy)

Kod dla stochastycznego spadku 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

Wyniki są prawie takie same, jak przy spadku bez użycia NumPy. W każdym razie to logiczne.

Dowiedzmy się, ile czasu zajęły nam stochastyczne spadki gradientu.

Kod do określenia czasu obliczeń SGD (80 tys. kroków)

print ' 33[1m' + ' 33[4m' +
"Czas wykonania stochastycznego spadku gradientu bez użycia biblioteki 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' +
"Czas wykonania stochastycznego spadku gradientu przy użyciu biblioteki 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: znów 'autorska' formuła pokazuje lepszy wynik. Wszystko to skłania do myślenia, że istnieją jeszcze bardziej zaawansowane sposoby użycia biblioteki NumPy, które rzeczywiście przyspieszają obliczenia. O tym w tym artykule się nie dowiemy. Będzie o czym pomyśleć w wolnym czasie :)

Podsumujmy

Zanim podsumujemy, chciałbym odpowiedzieć na pytanie, które z pewnością pojawiło się u naszego drogiego czytelnika. Po co takie 'męki' ze spadkami, dlaczego musimy chodzić po górach w górę i w dół (głównie w dół), aby znaleźć upragnioną dolinę, jeśli mamy potężne i proste narzędzie w postaci rozwiązania analitycznego, które natychmiast teleportuje nas w potrzebne miejsce?

Odpowiedź na to pytanie jest oczywista. Teraz omówiliśmy bardzo prosty przykład, w którym prawdziwa odpowiedź zależy od jednego atrybutu . W życiu rzadko się to zdarza, dlatego wyobraźmy sobie, że mamy 2, 30, 50 lub więcej atrybutó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ć próby i zawieść. Z kolei spadek gradientu i jego warianty będą powoli, ale pewnie przybliżać nas do celu — minimum funkcji. A jeśli chodzi o prędkość, nie martwcie się — z pewnością jeszcze raz omówimy sposoby, które pozwolą nam ustalać i regulować długość kroku (czyli prędkość).

A teraz, w skrócie, krótkie podsumowanie.

Po pierwsze, mam nadzieję, że materiały przedstawione w artykule pomogą początkującym "data scientistom" w zrozumieniu, jak rozwiązywać równania prostej (i nie tylko) regresji liniowej.

Po drugie, omówiliśmy kilka sposobów rozwiązania równania. Teraz, w zależności od sytuacji, możemy wybrać ten, który najlepiej pasuje do rozwiązania postawionego zadania.

Po trzecie, zobaczyliśmy moc dodatkowych ustawień, a mianowicie długości kroku spadku gradientu. Tego parametru nie można lekceważyć. Jak wspomniano wcześniej, aby ograniczyć koszty obliczeń, warto dostosować długość kroku w trakcie spadku.

Po czwarte, w naszym przypadku "autorskie" funkcje wykazały lepszy czas obliczeń. Prawdopodobnie wynika to z nieco nieprofesjonalnego zastosowania możliwości biblioteki. NumPy. Niemniej jednak, kolejny wniosek jest następujący. Z jednej strony czasami warto poddać w wątpliwość ugruntowane opinie, z drugiej zaś - nie zawsze warto wszystko komplikować - przeciwnie, czasami prostsze podejście jest bardziej efektywne. A ponieważ naszym celem było omówienie trzech podejść do rozwiązania równania prostej regresji liniowej, użycie "autorskich" funkcji w zupełności wystarczyło.

Literatura (lub coś w tym stylu)

1. Regresja liniowa

http://statistica.ru/theory/osnovy-lineynoy-regressii/

2. Metoda najmniejszych kwadratów

mathprofi.ru/metod_naimenshih_kvadratov.html

3. Pochodna

www.mathprofi.ru/chastnye_proizvodnye_primery.html

4. Gradient

mathprofi.ru/proizvodnaja_po_napravleniju_i_gradient.html

5. Spadek gradientu

habr.com/ru/post/471458

habr.com/ru/post/307312

artemarakcheev.com/2017-12-31/linear_regression

6. Biblioteka NumPy

docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.linalg.solve.html

docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.linalg.pinv.html

pythonworld.ru/numpy/2.html

Źródło: habr.com