🥇Вирішуємо рівняння простої лінійної регресії

У статті розглядається кілька способів визначення математичного рівняння лінії простої (парної) регресії.

Всі способи розв'язання рівняння, що розглядаються тут, засновані на методі найменших квадратів. Позначимо способи наступним чином:

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

Для кожного із способів вирішення рівняння прямої, у статті наведено різні функції, які в основному поділяються на ті, що написані без використання бібліотеки numpy та ті, які для проведення розрахунків застосовують numpy. Вважається, що вміле використання numpy дозволить скоротити витрати на обчислення.

Весь код, наведений у статті, написаний мовою пітон 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 '********************************************'

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

Тепер після того, як ми, по-перше, завантажили дані, по-друге, перевірили правильність завантаження і нарешті відформатували дані, проведемо першу візуалізацію. Часто для цього використовують метод парний сюжет бібліотеки Морський народжений. У нашому прикладі, через обмеженість цифр немає сенсу застосовувати бібліотеку Морський народжений. Ми скористаємося звичайною бібліотекою Матплотліб і подивимося лише на діаграму розсіювання.

Код діаграми розсіювання

print 'График №1 "Зависимость выручки от месяца года"'

plt.plot(x_us,y_us,'o',color='green',markersize=16)
plt.xlabel('$Months$', size=16)
plt.ylabel('$Sales$', size=16)
plt.show()

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

Аналітичне рішення

Скористаємося звичайними інструментами в пітон і розв'яжемо систему рівнянь:

begin{equation*}
begin{cases}
na + bsumlimits_{i=1}^nx_i - sumlimits_{i=1}^ny_i = 0

sumlimits_{i=1}^nx_i(a +bsumlimits_{i=1}^nx_i — sumlimits_{i=1}^ny_i) = 0
end{cases}
end{equation*}

За правилом Крамера знайдемо загальний визначник, а також визначники по і , після чого, розділивши визначник по на загальний визначник - знайдемо коефіцієнт , аналогічно знайдемо коефіцієнт .

Код аналітичного рішення

# определим функцию для расчета коэффициентов a и b по правилу Крамера
def Kramer_method (x,y):
        # сумма значений (все месяца)
    sx = sum(x)
        # сумма истинных ответов (выручка за весь период)
    sy = sum(y)
        # сумма произведения значений на истинные ответы
    list_xy = []
    [list_xy.append(x[i]*y[i]) for i in range(len(x))]
    sxy = sum(list_xy)
        # сумма квадратов значений
    list_x_sq = []
    [list_x_sq.append(x[i]**2) for i in range(len(x))]
    sx_sq = sum(list_x_sq)
        # количество значений
    n = len(x)
        # общий определитель
    det = sx_sq*n - sx*sx
        # определитель по a
    det_a = sx_sq*sy - sx*sxy
        # искомый параметр a
    a = (det_a / det)
        # определитель по b
    det_b = sxy*n - sy*sx
        # искомый параметр b
    b = (det_b / det)
        # контрольные значения (прооверка)
    check1 = (n*b + a*sx - sy)
    check2 = (b*sx + a*sx_sq - sxy)
    return [round(a,4), round(b,4)]

# запустим функцию и запишем правильные ответы
ab_us = Kramer_method(x_us,y_us)
a_us = ab_us[0]
b_us = ab_us[1]
print ' 33[1m' + ' 33[4m' + "Оптимальные значения коэффициентов a и b:"  + ' 33[0m' 
print 'a =', a_us
print 'b =', b_us
print

# определим функцию для подсчета суммы квадратов ошибок
def errors_sq_Kramer_method(answers,x,y):
    list_errors_sq = []
    for i in range(len(x)):
        err = (answers[0] + answers[1]*x[i] - y[i])**2
        list_errors_sq.append(err)
    return sum(list_errors_sq)

# запустим функцию и запишем значение ошибки
error_sq = errors_sq_Kramer_method(ab_us,x_us,y_us)
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений" + ' 33[0m'
print error_sq
print

# замерим время расчета
# print ' 33[1m' + ' 33[4m' + "Время выполнения расчета суммы квадратов отклонений:" + ' 33[0m'
# % timeit error_sq = errors_sq_Kramer_method(ab,x_us,y_us)

Ось що у нас вийшло:

Отже, значення коефіцієнтів знайдено, суму квадратів відхилень встановлено. Намалюємо на гістограмі розсіювання пряму лінію відповідно до знайдених коефіцієнтів.

Код лінії регресії

# определим функцию для формирования массива рассчетных значений выручки
def sales_count(ab,x,y):
    line_answers = []
    [line_answers.append(ab[0]+ab[1]*x[i]) for i in range(len(x))]
    return line_answers

# построим графики
print 'Грфик№2 "Правильные и расчетные ответы"'
plt.plot(x_us,y_us,'o',color='green',markersize=16, label = '$True$ $answers$')
plt.plot(x_us, sales_count(ab_us,x_us,y_us), color='red',lw=4,
         label='$Function: a + bx,$ $where$ $a='+str(round(ab_us[0],2))+',$ $b='+str(round(ab_us[1],2))+'$')
plt.xlabel('$Months$', size=16)
plt.ylabel('$Sales$', size=16)
plt.legend(loc=1, prop={'size': 16})
plt.show()

Графік №2 «Правильні та розрахункові відповіді»

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

Код графіка відхилень

# определим функцию для формирования массива отклонений в процентах
def error_per_month(ab,x,y):
    sales_c = sales_count(ab,x,y)
    errors_percent = []
    for i in range(len(x)):
        errors_percent.append(100*(sales_c[i]-y[i])/y[i])
    return errors_percent

# построим график
print 'График№3 "Отклонения по-месячно, %"'
plt.gca().bar(x_us, error_per_month(ab_us,x_us,y_us), color='brown')
plt.xlabel('Months', size=16)
plt.ylabel('Calculation error, %', size=16)
plt.show()

Графік №3 «Відхилення, %»

Не ідеально, але наше завдання ми виконали.

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

Код аналітичного рішення (NumPy)

# для начала добавим столбец с не изменяющимся значением в 1. 
# Данный столбец нужен для того, чтобы не обрабатывать отдельно коэффицент a
vector_1 = np.ones((x_np.shape[0],1))
x_np = table_zero[['x']].values # на всякий случай приведем в первичный формат вектор x_np
x_np = np.hstack((vector_1,x_np))

# проверим то, что все сделали правильно
print vector_1[0:3]
print x_np[0:3]
print '***************************************'
print

# напишем функцию, которая определяет значения коэффициентов a и b с использованием псевдообратной матрицы
def pseudoinverse_matrix(X, y):
    # задаем явный формат матрицы признаков
    X = np.matrix(X)
    # определяем транспонированную матрицу
    XT = X.T
    # определяем квадратную матрицу
    XTX = XT*X
    # определяем псевдообратную матрицу
    inv = np.linalg.pinv(XTX)
    # задаем явный формат матрицы ответов
    y = np.matrix(y)
    # находим вектор весов
    return (inv*XT)*y

# запустим функцию
ab_np = pseudoinverse_matrix(x_np, y_np)
print ab_np
print '***************************************'
print

# напишем функцию, которая использует для решения матричное уравнение
def matrix_equation(X,y):
    a = np.dot(X.T, X)
    b = np.dot(X.T, y)
    return np.linalg.solve(a, b)

# запустим функцию
ab_np = matrix_equation(x_np,y_np)
print ab_np

Порівняємо час, який було витрачено на визначення коефіцієнтів и , Відповідно до трьох представлених способів.

Код для обчислення часу розрахунків

print ' 33[1m' + ' 33[4m' + "Время выполнения расчета коэффициентов без использования библиотеки NumPy:" + ' 33[0m'
% timeit ab_us = Kramer_method(x_us,y_us)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Время выполнения расчета коэффициентов с использованием псевдообратной матрицы:" + ' 33[0m'
%timeit ab_np = pseudoinverse_matrix(x_np, y_np)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Время выполнения расчета коэффициентов с использованием матричного уравнения:" + ' 33[0m'
%timeit ab_np = matrix_equation(x_np, y_np)

На невеликій кількості даних вперед виходить «самописна» функція, яка знаходить коефіцієнти методом Крамера.

Тепер можна перейти до інших способів знаходження коефіцієнтів и .

Градієнтний спуск

Спочатку визначимо, що таке градієнт. По-простому, градієнт - це відрізок, який показує напрямок максимального зростання функції. За аналогією з підйомом у гору, то куди дивиться градієнт, там і є найкрутіший підйом до вершини гори. Розвиваючи приклад із горою, згадуємо, що насправді нам потрібен найкрутіший спуск, щоб якнайшвидше досягти низини, тобто мінімуму — місця, де функція не зростає і не зменшується. У цьому місці похідна дорівнюватиме нулю. Отже, нам потрібний не градієнт, а антиградієнт. Для знаходження антиградієнта потрібно лише помножити градієнт на -1 (мінус один).

Звернемо увагу на те, що функція може мати кілька мінімумів, і опустившись в один з них за запропонованим алгоритмом, ми не зможемо знайти інший мінімум, який можливо знаходиться нижче знайденого. Розслабимося, нам це не загрожує! У нашому випадку ми маємо справу з єдиним мінімумом, оскільки наша функція на графіці є звичайною параболу. А як ми всі маємо чудово знати зі шкільного курсу математики — парабола має лише один мінімум.

Після того, як ми з'ясували для чого нам знадобився градієнт, а також те, що градієнт — це відрізок, тобто вектор із заданими координатами, які якраз є тими самими коефіцієнтами и ми можемо продати градієнтний спуск.

Перед запуском пропоную прочитати буквально кілька пропозицій про алгоритм спуску:

Визначаємо псевдовипадковим чином координати коефіцієнтів и . У нашому прикладі ми визначатимемо коефіцієнти поблизу нуля. Це поширена практика, проте для кожного випадку може бути передбачена своя практика.
Від координати віднімаємо значення приватної похідної 1-го порядку в точці . Тож якщо похідна буде позитивна, то функція зростає. Отже, забираючи значення похідної, ми рухатимемося у зворотний бік зростання, тобто у бік спуску. Якщо похідна негативна, отже функція у цій точці зменшується і забираючи значення похідної ми рухаємося у бік спуску.
Проводимо аналогічну операцію з координатою : віднімаємо значення приватної похідної в точці .
Для того, щоб не перескочити мінімум і не відлетіти в далекий космос, необхідно встановити розмір кроку у бік спуску. Загалом, можна написати цілу статтю про те, як правильніше встановити крок і як його змінювати в процесі спуску, щоб знизити витрати на обчислення. Але зараз перед нами трохи інше завдання, і ми науковим методом «тику» або як кажуть у народі, емпіричним шляхом, встановимо розмір кроку.
Після того, як ми із заданих координат и відняли значення похідних, отримуємо нові координати и . Робимо наступний крок (віднімання), вже з розрахованих координат. І так цикл запускається знову і знову, доки не буде досягнута необхідна збіжність.

Всі! Тепер ми готові вирушити на пошуки найглибшої ущелини Маріанської западини. Приступаємо.

Код для градієнтного спуску

# напишем функцию градиентного спуска без использования библиотеки NumPy. 
# Функция на вход принимает диапазоны значений x,y, длину шага (по умолчанию=0,1), допустимую погрешность(tolerance)
def gradient_descent_usual(x_us,y_us,l=0.1,tolerance=0.000000000001):
    # сумма значений (все месяца)
    sx = sum(x_us)
    # сумма истинных ответов (выручка за весь период)
    sy = sum(y_us)
    # сумма произведения значений на истинные ответы
    list_xy = []
    [list_xy.append(x_us[i]*y_us[i]) for i in range(len(x_us))]
    sxy = sum(list_xy)
    # сумма квадратов значений
    list_x_sq = []
    [list_x_sq.append(x_us[i]**2) for i in range(len(x_us))]
    sx_sq = sum(list_x_sq)
    # количество значений
    num = len(x_us)
    # начальные значения коэффициентов, определенные псевдослучайным образом
    a = float(random.uniform(-0.5, 0.5))
    b = float(random.uniform(-0.5, 0.5))
    # создаем массив с ошибками, для старта используем значения 1 и 0
    # после завершения спуска стартовые значения удалим
    errors = [1,0]
    # запускаем цикл спуска
    # цикл работает до тех пор, пока отклонение последней ошибки суммы квадратов от предыдущей, не будет меньше tolerance
    while abs(errors[-1]-errors[-2]) > tolerance:
        a_step = a - l*(num*a + b*sx - sy)/num
        b_step = b - l*(a*sx + b*sx_sq - sxy)/num
        a = a_step
        b = b_step
        ab = [a,b]
        errors.append(errors_sq_Kramer_method(ab,x_us,y_us))
    return (ab),(errors[2:])

# запишем массив значений 
list_parametres_gradient_descence = gradient_descent_usual(x_us,y_us,l=0.1,tolerance=0.000000000001)


print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_gradient_descence[0][0],3)
print 'b =', round(list_parametres_gradient_descence[0][1],3)
print


print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_gradient_descence[1][-1],3)
print



print ' 33[1m' + ' 33[4m' + "Количество итераций в градиентном спуске:" + ' 33[0m'
print len(list_parametres_gradient_descence[1])
print

Ми поринули на саме дно Маріанської западини і там виявили ті самі значення коефіцієнтів и , Що власне і слід було очікувати.

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

Код для градієнтного спуску (NumPy)

# перед тем определить функцию для градиентного спуска с использованием библиотеки NumPy, 
# напишем функцию определения суммы квадратов отклонений также с использованием NumPy
def error_square_numpy(ab,x_np,y_np):
    y_pred = np.dot(x_np,ab)
    error = y_pred - y_np
    return sum((error)**2)

# напишем функцию градиентного спуска с использованием библиотеки NumPy. 
# Функция на вход принимает диапазоны значений x,y, длину шага (по умолчанию=0,1), допустимую погрешность(tolerance)
def gradient_descent_numpy(x_np,y_np,l=0.1,tolerance=0.000000000001):
    # сумма значений (все месяца)
    sx = float(sum(x_np[:,1]))
    # сумма истинных ответов (выручка за весь период)
    sy = float(sum(y_np))
    # сумма произведения значений на истинные ответы
    sxy = x_np*y_np
    sxy = float(sum(sxy[:,1]))
    # сумма квадратов значений
    sx_sq = float(sum(x_np[:,1]**2))
    # количество значений
    num = float(x_np.shape[0])
    # начальные значения коэффициентов, определенные псевдослучайным образом
    a = float(random.uniform(-0.5, 0.5))
    b = float(random.uniform(-0.5, 0.5))
    # создаем массив с ошибками, для старта используем значения 1 и 0
    # после завершения спуска стартовые значения удалим
    errors = [1,0]
    # запускаем цикл спуска
    # цикл работает до тех пор, пока отклонение последней ошибки суммы квадратов от предыдущей, не будет меньше tolerance
    while abs(errors[-1]-errors[-2]) > tolerance:
        a_step = a - l*(num*a + b*sx - sy)/num
        b_step = b - l*(a*sx + b*sx_sq - sxy)/num
        a = a_step
        b = b_step
        ab = np.array([[a],[b]])
        errors.append(error_square_numpy(ab,x_np,y_np))
    return (ab),(errors[2:])

# запишем массив значений 
list_parametres_gradient_descence = gradient_descent_numpy(x_np,y_np,l=0.1,tolerance=0.000000000001)

print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_gradient_descence[0][0],3)
print 'b =', round(list_parametres_gradient_descence[0][1],3)
print


print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_gradient_descence[1][-1],3)
print

print ' 33[1m' + ' 33[4m' + "Количество итераций в градиентном спуске:" + ' 33[0m'
print len(list_parametres_gradient_descence[1])
print

значення коефіцієнтів и незмінні.

Подивимося те що, як змінювалася помилка при градієнтному спуску, тобто як змінювалася сума квадратів відхилень з кожним кроком.

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

print 'График№4 "Сумма квадратов отклонений по-шагово"'
plt.plot(range(len(list_parametres_gradient_descence[1])), list_parametres_gradient_descence[1], color='red', lw=3)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()

Графік №4 "Сума квадратів відхилень при градієнтному спуску"

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

Насамкінець оцінимо різницю в часі виконання коду:

Код визначення часу обчислення градієнтного спуску

print ' 33[1m' + ' 33[4m' + "Время выполнения градиентного спуска без использования библиотеки NumPy:" + ' 33[0m'
%timeit list_parametres_gradient_descence = gradient_descent_usual(x_us,y_us,l=0.1,tolerance=0.000000000001)
print '***************************************'
print

print ' 33[1m' + ' 33[4m' + "Время выполнения градиентного спуска с использованием библиотеки NumPy:" + ' 33[0m'
%timeit list_parametres_gradient_descence = gradient_descent_numpy(x_np,y_np,l=0.1,tolerance=0.000000000001)

Можливо, ми робимо щось не те, але знову проста «самописна» функція, яка не використовує бібліотеку numpy випереджає за часом виконання розрахунків функцію, яка використовує бібліотеку numpy.

Але ми не стоїмо на місці, а рухаємось у бік вивчення ще одного захоплюючого способу розв'язання рівняння простої лінійної регресії. Зустрічайте!

Стохастичний градієнтний спуск

Щоб швидше зрозуміти принцип роботи стохастичного градієнтного спуску, краще визначити його відмінності від звичайного градієнтного спуску. Ми, у випадку з градієнтним спуском, у рівняннях похідних від и використовували суми значень всіх ознак та справжніх відповідей, що є у вибірці (тобто суми всіх и ). У стохастичному градієнтному спуску ми не будемо використовувати всі значення, що є у вибірці, а натомість, псевдовипадковим чином виберемо так званий індекс вибірки та використовуємо його значення.

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

Які переваги стохастичного градієнтного спуску перед звичайним? Якщо у нас розмір вибірки дуже великий і вимірюється десятками тисяч значень, то значно простіше обробити, припустімо випадкову тисячу з них, ніж усю вибірку. Ось тут і запускається стохастичний градієнтний спуск. У нашому випадку ми звичайно великої різниці не помітимо.

Дивимося код.

Код для стохастичного градієнтного спуску

# определим функцию стох.град.шага
def stoch_grad_step_usual(vector_init, x_us, ind, y_us, l):
#     выбираем значение икс, которое соответствует случайному значению параметра ind 
# (см.ф-цию stoch_grad_descent_usual)
    x = x_us[ind]
#     рассчитывыаем значение y (выручку), которая соответствует выбранному значению x
    y_pred = vector_init[0] + vector_init[1]*x_us[ind]
#     вычисляем ошибку расчетной выручки относительно представленной в выборке
    error = y_pred - y_us[ind]
#     определяем первую координату градиента ab
    grad_a = error
#     определяем вторую координату ab
    grad_b = x_us[ind]*error
#     вычисляем новый вектор коэффициентов
    vector_new = [vector_init[0]-l*grad_a, vector_init[1]-l*grad_b]
    return vector_new


# определим функцию стох.град.спуска
def stoch_grad_descent_usual(x_us, y_us, l=0.1, steps = 800):
#     для самого начала работы функции зададим начальные значения коэффициентов
    vector_init = [float(random.uniform(-0.5, 0.5)), float(random.uniform(-0.5, 0.5))]
    errors = []
#     запустим цикл спуска
# цикл расчитан на определенное количество шагов (steps)
    for i in range(steps):
        ind = random.choice(range(len(x_us)))
        new_vector = stoch_grad_step_usual(vector_init, x_us, ind, y_us, l)
        vector_init = new_vector
        errors.append(errors_sq_Kramer_method(vector_init,x_us,y_us))
    return (vector_init),(errors)


# запишем массив значений 
list_parametres_stoch_gradient_descence = stoch_grad_descent_usual(x_us, y_us, l=0.1, steps = 800)

print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_stoch_gradient_descence[0][0],3)
print 'b =', round(list_parametres_stoch_gradient_descence[0][1],3)
print


print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_stoch_gradient_descence[1][-1],3)
print

print ' 33[1m' + ' 33[4m' + "Количество итераций в стохастическом градиентном спуске:" + ' 33[0m'
print len(list_parametres_stoch_gradient_descence[1])

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

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

print 'График №5 "Сумма квадратов отклонений по-шагово"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1])), list_parametres_stoch_gradient_descence[1], color='red', lw=2)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()

Графік №5 "Сума квадратів відхилень при стохастичному градієнтному спуску"

Подивившись на графік, все стає на свої місця, і зараз ми всі виправимо.

Отже, що сталося? Сталося таке. Коли ми вибираємо випадковим чином місяць, саме для обраного місяця наш алгоритм прагне зменшити помилку в розрахунку виручки. Потім вибираємо інший місяць і повторюємо розрахунок, але помилку зменшуємо вже другого обраного місяця. А тепер згадаємо, що у нас перші два місяці суттєво відхиляються від лінії рівняння простої лінійної регресії. Це означає, що коли вибирається кожен із цих двох місяців, то зменшуючи помилку кожного з них, наш алгоритм серйозно збільшує помилку по всій вибірці. То що робити? Відповідь проста: треба зменшити крок спуску. Адже зменшивши крок спуску, помилка так само перестане «стрибати» то вгору, то вниз. Точніше, помилка «скакати» не перестане, але це робитиме не так швидко:) Перевіримо.

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

# запустим функцию, уменьшив шаг в 100 раз и увеличив количество шагов соответсвующе 
list_parametres_stoch_gradient_descence = stoch_grad_descent_usual(x_us, y_us, l=0.001, steps = 80000)

print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_stoch_gradient_descence[0][0],3)
print 'b =', round(list_parametres_stoch_gradient_descence[0][1],3)
print


print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_stoch_gradient_descence[1][-1],3)
print



print ' 33[1m' + ' 33[4m' + "Количество итераций в стохастическом градиентном спуске:" + ' 33[0m'
print len(list_parametres_stoch_gradient_descence[1])

print 'График №6 "Сумма квадратов отклонений по-шагово"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1])), list_parametres_stoch_gradient_descence[1], color='red', lw=2)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()

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

Значення коефіцієнтів покращали, але все одно не ідеальні. Гіпотетично це можна виправити в такий спосіб. Вибираємо, наприклад, на останніх 1000 ітераціях значення коефіцієнтів, з якими було допущено мінімальну помилку. Щоправда, нам для цього доведеться записувати ще й самі значення коефіцієнтів. Ми не робитимемо цього, а краще звернемо увагу на графік. Він виглядає гладким, і помилка зменшується рівномірно. Насправді, це не так. Подивимося перші 1000 ітерацій і порівняємо їх із останніми.

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

print 'График №7 "Сумма квадратов отклонений по-шагово. Первые 1000 итераций"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][:1000])), 
         list_parametres_stoch_gradient_descence[1][:1000], color='red', lw=2)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()

print 'График №7 "Сумма квадратов отклонений по-шагово. Последние 1000 итераций"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][-1000:])), 
         list_parametres_stoch_gradient_descence[1][-1000:], color='red', lw=2)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()

Графік №7 «Сума квадратів відхилень SGD (перші 1000 кроків)»

Графік №8 "Сума квадратів відхилень SGD (останні 1000 кроків)"

На початку спуску ми спостерігаємо досить рівномірне і круте зменшення помилки. На останніх ітераціях ми бачимо, що помилка ходить навколо та близько значення в 1,475 і в деякі моменти навіть дорівнює цьому оптимальному значенню, але потім все одно йде вгору ... Повторюся, можна записувати значення коефіцієнтів и , а потім вибрати ті, у яких помилка мінімальна. Однак у нас виникла проблема серйозніша: нам довелося зробити 80 тис. кроків, щоб отримати значення, близькі до оптимальних. А це вже суперечить ідеї про економію часу обчислень при стохастичному градієнтному спуску щодо градієнтного. Що можна поправити та покращити? Не важко помітити, що на перших ітераціях ми впевнено йдемо вниз і, отже, варто залишити великий крок на перших ітераціях і в міру просування вперед крок зменшувати. Ми не робитимемо цього в цій статті — вона й так уже затягнулася. Бажаючі можуть самі подумати, як це зробити, це не складно 🙂

Тепер виконаємо стохастичний градієнтний спуск, використовуючи бібліотеку numpy (і не спотикатимемося про каміння, яке ми виявили раніше)

Код для стохастичного градієнтного спуску (NumPy)

# для начала напишем функцию градиентного шага
def stoch_grad_step_numpy(vector_init, X, ind, y, l):
    x = X[ind]
    y_pred = np.dot(x,vector_init)
    err = y_pred - y[ind]
    grad_a = err
    grad_b = x[1]*err
    return vector_init - l*np.array([grad_a, grad_b])

# определим функцию стохастического градиентного спуска
def stoch_grad_descent_numpy(X, y, l=0.1, steps = 800):
    vector_init = np.array([[np.random.randint(X.shape[0])], [np.random.randint(X.shape[0])]])
    errors = []
    for i in range(steps):
        ind = np.random.randint(X.shape[0])
        new_vector = stoch_grad_step_numpy(vector_init, X, ind, y, l)
        vector_init = new_vector
        errors.append(error_square_numpy(vector_init,X,y))
    return (vector_init), (errors)

# запишем массив значений 
list_parametres_stoch_gradient_descence = stoch_grad_descent_numpy(x_np, y_np, l=0.001, steps = 80000)

print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_stoch_gradient_descence[0][0],3)
print 'b =', round(list_parametres_stoch_gradient_descence[0][1],3)
print


print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_stoch_gradient_descence[1][-1],3)
print



print ' 33[1m' + ' 33[4m' + "Количество итераций в стохастическом градиентном спуске:" + ' 33[0m'
print len(list_parametres_stoch_gradient_descence[1])
print

Значення вийшли майже такими, як і при спуску без використання numpy. Втім, це логічно.

Дізнаємось скільки ж часу займали у нас стохастичні градієнтні спуски.

Код визначення часу обчислення SGD (80 тис. кроків)

print ' 33[1m' + ' 33[4m' +
"Время выполнения стохастического градиентного спуска без использования библиотеки NumPy:"
+ ' 33[0m'
%timeit list_parametres_stoch_gradient_descence = stoch_grad_descent_usual(x_us, y_us, l=0.001, steps = 80000)
print '***************************************'
print

print ' 33[1m' + ' 33[4m' +
"Время выполнения стохастического градиентного спуска с использованием библиотеки NumPy:"
+ ' 33[0m'
%timeit list_parametres_stoch_gradient_descence = stoch_grad_descent_numpy(x_np, y_np, l=0.001, steps = 80000)

Чим далі в ліс, тим темніші хмари: знову «самописна» формула показує найкращий результат. Все це наводить на думку про те, що повинні існувати ще тонші способи використання бібліотеки numpy, які справді прискорюють операції обчислень. У цій статті ми про них не дізнаємося. Буде про що подумати на дозвіллі :)

Резюме

Перед тим як резюмувати, хотілося б відповісти на запитання, яке швидше за все виникло у нашого дорогого читача. Для чого, власне, такі «муки» зі спусками, навіщо нам ходити горою вгору і вниз (переважно вниз), щоб знайти заповітну низину, якщо в наших руках такий потужний і простий прилад, у вигляді аналітичного рішення, яке миттєво телепортує нас у потрібне місце?

Відповідь це питання лежить на поверхні. Зараз ми розбирали дуже простий приклад, у якому справжня відповідь залежить від однієї ознаки . У житті таке зустрінеш не часто, тому уявімо, що у нас ознак 2, 30, 50 або більше. Додамо до цього тисячі, а то й десятки тисяч значень для кожної ознаки. У цьому випадку аналітичне рішення може не витримати випробування та дати збій. У свою чергу градієнтний спуск та його варіації будуть повільно, але правильно наближати нас до мети – мінімуму функції. А щодо швидкості не хвилюйтеся — ми напевно ще розберемо способи, які дозволять нам задавати і регулювати довжину кроку (тобто швидкість).

А тепер власне коротке резюме.

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

По-друге, ми розглянули кілька способів розв'язання рівняння. Тепер, залежно від ситуації, ми можемо вибрати той, який найкраще підходить для вирішення поставленого завдання.

По-третє, ми побачили силу додаткових налаштувань, саме довжини кроку градієнтного спуску. Цим параметром не можна нехтувати. Як було зазначено вище, з метою скорочення витрат на проведення обчислень, довжину кроку варто змінювати в процесі спуску.

По-четверте, у нашому випадку, «самописні» функції показали найкращий часовий результат обчислень. Ймовірно, це пов'язано з не професійним застосуванням можливостей бібліотеки. numpy. Але як би там не було, висновок напрошується наступний. З одного боку, іноді варто ставити під сумнів усталені думки, а з іншого — не завжди варто все ускладнювати — навпаки іноді ефективніше виявляється простіший спосіб вирішення завдання. Оскільки мета в нас була розібрати три підходи у вирішенні рівняння простої лінійної регресії, то використання «самописних» функцій нам цілком вистачило.