Решаваме уравнението на проста линейна регресия

В статията се разглеждат няколко начина за определяне на математическото уравнение на правата линия на проста (двойна) регресия.

Всички разглеждани тук методи за решаване на уравнението се основават на метода на най-малките квадрати. Нека обозначим методите по следния начин:

• Аналитично решение
• Градиентен спуск
• Стохастичен градиентен спуск

За всеки от методите за решаване на уравнението на правата, в статията са представени различни функции, които основно се разделят на тези, написани без използване на библиотека NumPy и тези, които за извършване на изчисления използват NumPy. Смята се, че умелото използване на NumPy ще доведе до намаляване на разходите за изчисления.

Целият код, представен в статията, е написан на езика python 2.7 с използването на Jupyter Notebook. Изходният код и файлът с данни от извадката са публикувани на Гитхаб.

Статията е ориентирована както към начинаещи, така и към тези, които вече са започнали да изучават доста обширната област в изкуствения интелект — машинното обучение.

За илюстрация на материала ще използваме много прост пример.

Условия на примера

Имаме пет стойности, които характеризират зависимостта 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 библиотеки Seaborn. В нашия пример, поради ограничението на цифрите, няма смисъл да се прилага библиотеката Seaborn. Ще използваме стандартната библиотека Matplotlib и ще разгледаме само диаграмата на разсейване.

Код на диаграмата на разсейване

print 'Графика №1 "Зависимост на приходите от месец на годината"'

plt.plot(x_us,y_us,'o',color='green',markersize=16)
plt.xlabel('$Месеци$', size=16)
plt.ylabel('$Приходи$', size=16)
plt.show()

Графика №1 «Зависимост на приходите от месец на годината»

Аналитично решение

Ще използваме най-обикновените инструменти в python и ще решим системата от уравнения:

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

Да сравним времето, което е било похарчено за определяне на коефициентите и , в съответствие с 3-те представени метода.

Код за изчисляване на времето на изчисленията

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 (минус едно).

Обърнете внимание, че функцията може да има няколко минимума, а следвайки предложената по-долу алгоритъм, не можем да намерим друг минимум, който евентуално да е по-нисък от намерения. Нека се успокоим, това не ни заплашва! В нашия случай имаме работа с единствен минимум, тъй като нашата функция на графика представлява обикновена парабола. Както всички трябва да знаем от училищния курс по математика - параболата има само един минимум.

След като установихме защо ни е нужен градиентът, както и че градиентът е отсечка, тоест вектор с зададени координати, които точно представляват същите тези коефициенти и можем да реализираме градиентния спуск.

Преди да стартираме, предлагам да прочетем буквално няколко изречения за алгоритъма на спускането:

• Определяме псевдослучайно координатите на коефициентите и . В нашия пример ще определяме коефициентите в близост до нула. Това е разпространена практика, но за всеки случай може да има своя собствена практика.
• От координатата изваждаме стойността на първичната производна в точка . Така, ако производната е положителна, функцията расте. Следователно, изваждайки стойността на производната, ще се движим в обратна посока на растежа, тоест в посока на спускането. Ако производната е отрицателна, значи функцията в тази точка намалява и изваждайки стойността на производната, се движим в посока на спускането.
• Извършваме аналогична операция с координатата : извличаме стойността на частната производна в точката .
• За да не прескочим минимума и да не излетим в далечния космос, трябва да зададем размера на стъпката в посока на спускане. По принцип, може да се напише цяла статия за това как правилно да установим стъпката и как да я променяме в процеса на спускане, за да намалим разходите за изчисления. Но в момента пред нас е малко по-различна задача, и с научния метод на "тикане", или както казват в обикновеното народ, емпиричен път, ще установим размера на стъпката.
• След като от зададените координати и извлякохме стойностите на производните, получаваме нови координати и . Правим следващата стъпка (изваждане), вече от изчислените координати. И така цикълът се стартира отново и отново, докато не бъде постигната необходимата конвергенция.

Всичко! Сега сме готови да се отправим на търсене на най-дълбокото дефиле на Марианската впадина. Започваме.

Код за градиентен спуск

# напишем функцию градиентного спуска без использования библиотеки 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('Стъпки (Итерация)', size=16)
plt.ylabel('Сума на квадратите на отклоненията', 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('Стъпки (Итерация)', size=16)
plt.ylabel('Сума на квадратите на отклоненията', 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('Стъпки (Итерация)', size=16)
plt.ylabel('Сума на квадратите на отклоненията', 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('Стъпки (Итерация)', size=16)
plt.ylabel('Сума на квадратите на отклоненията', 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Но каквото и да е, извода е следният. От една страна, понякога е необходимо да се поставят под съмнение утвърдените мнения, а от друга - не винаги е нужно да се усложняват нещата — понякога по-простият метод за решаване на задача се оказва по-ефективен. И тъй като нашата цел беше да разгледаме три подхода за решаване на уравнението на проста линейна регресия, то използването на „самописни“ функции беше напълно достатъчно.

Литература (или нещо подобно)

1. Линейна регресия

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

2. Метод на най-малките квадрати

mathprofi.ru/metod_naimenshih_kvadratov.html

3. Производна

www.mathprofi.ru/chastnye_proizvodnye_primery.html

4. Градиент

mathprofi.ru/proizvodnaja_po_napravleniju_i_gradient.html

5. Градиентен спуск

habr.com/ru/post/471458

habr.com/ru/post/307312

artemarakcheev.com/2017-12-31/linear_regression

6. Библиотека 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

Източник: habr.com