ML: Градієнтний метод

Вступ

У цьому документі продовжується обговорення багатовимірної лінійної моделі. На її прикладі ми розглянемо метод стохастичного градієнтного спуску і вибір оптимальних гіперпараметрів. Наведений код на Python з використанням бібліотек Numpy, PyTorch і Keras, можна знайти у файлі: ML_Line_Model.ipynb.

Градієнтний метод

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

Позначимо набір параметрів вектором $\mathbf{w}$. Помилка є їх функцією $L=L(\mathbf{w})$. У просторі параметрів існують поверхні постійного значення $L(\mathbf{w})=\mathrm{const}$. Зсув $d\mathbf{w}$ уздовж таких поверхонь не змінює $L$ і, отже, її диференціал дорівнює нулю:

$$ dL = \frac{\partial L}{\partial \mathbf{w}}\cdot d\mathbf{w} = \mathbf{g}\cdot d\mathbf{w} = 0, $$ $$ ~ $$ $$ \mathbf{g} = \frac{\partial L}{\partial \mathbf{w}} = \Bigr\{\frac{\partial L}{\partial w_1},..., \frac{\partial L}{\partial w_m} \Bigr\} $$

Таким чином, вектор градієнта $\mathbf{g}$ перпендикулярний поверхням $L=\mathrm{const}$ і спрямований у бік збільшення $L$ (як і будь-яка похідна). При наближенні до мінімуму довжина градієнта, зазвичай зменшується, а в точці мінімуму (екстремум) вона дорівнює нулю (вище другий рисунок).

Щоб знайти мінімум $L$, необхідно рухатися у зворотному до градієнта напрямку (уздовж антиградієнта), з кроком пропорційним деякому числу $\lambda$.
Цей гіперпараметр називається швидкістю навчання:

$$ \mathbf{w}^{(t+1)} = \mathbf{w}^{(t)} - \lambda\,\mathbf{g}^{(t)},~~~~~~~~~~~~~\mathbf{g}^{(t)} = \frac{\partial L}{\partial\mathbf{w}^{(t)}}. $$

Чим більша швидкість навчання $\lambda$, тим швидше параметри моделі наближаються до оптимальних значень. Однак при великих $\lambda$ існує ризик проскочити мінімум (незважаючи на зменшення довжини градієнта в його околі). Проілюструємо це прикладами.

✍ В одновимірному випадку рівняння параболи з мінімумом у точці $w=0$ і її "градієнт" $g$ мають вигляд: $$ L(w)=\frac{1}{2}\,\Gamma\,w^2,~~~~~~~~~~~~~g=\frac{dL}{dw}=\Gamma\,w, $$ де $\Gamma$ - константа. Спуск зі швидкістю $\lambda$ з точки $w^{(0)}$, на $t$-тій ітерації призводить у точку: $$w^{(t)}=(1-\lambda\,\Gamma)^t\,w^{(0)}.$$ При $\lambda > 1/\Gamma$ метод стає нестійким: $w^{(t)} \to \pm \infty$. Однак, чим ближче $\lambda$ знизу до критичного значення $\lambda_c = 1/\Gamma$, тим швидше досягається мінімум.

✍ Аналогічно поводиться градієнтний спуск для "багатовимірної параболи" з різними коефіцієнтами $\Gamma_i$ по кожній координаті: $$ L = \frac{1}{2}\sum_i \Gamma_i\,w_i ^2,~~~~~~~~~g_i = \Gamma_i\, w_i,~~~~~~w^{(t)}_i=(1-\lambda\Gamma_i)^t\,w^{(0)}_i. $$ Якщо ми знаходимося в точці $\mathbf{w}$, то найкоротший шлях до мінімуму $\mathbf{w}_{\text{min}}=\mathbf{0}$ буде в напрямку вектора $-\mathbf{w}$. Однак градієнтний метод буде змінювати параметри у відмінному напрямку $-\mathbf{g}$ і поведе до мінімуму "окружним шляхом". Справа наведено траєкторії руху до мінімуму $\{0,0\}$ у двовимірному випадку з різних початкових положень c $\Gamma=\{5, 1\}$ і $\lambda=0.1$. Звернімо увагу не тільки на кривизну траєкторій, але і на помітне зменшення кроку в міру наближення до мінімуму. Критичне значення $\lambda_c$ визначається найбільшим коефіцієнтом параболи: $\lambda_c = 1/\max \Gamma_i = 0.2$.

Stochastic Gradient Descent

SGD (Stochastic Gradient Descent) - кілька більш просунутий градієнтний метод оптимізації. У ньому помилка $L$ обчислюється не по всіх даних, а по вибірці (пачка, batch) розміру batch_size. Це прискорює збіжність, т.к. для обчислення чергового кроку потрібні не всі дані, а лише невелика їх підмножина. Після кожної підправки параметрів, вибірка змінюється (тому метод називається "stochastic"). При проходженні через усі N приклади (одна епоха), відбувається int(N/btach_size) ітерацій (кроків у просторі ознак).

Оскільки даних для обчислення $\mathbf{g}$ небагато (зазвичай batch_size = 10-300), градієнт від ітерації до ітерації може "виляти" в різні боки. Щоб зменшити цей ефект, вектор $\mathbf{g}$ усереднюється $\langle\mathbf{g}\rangle$ за допомогою ковзного середнього. Це ж усереднення допомагає переміщатися по покритій "дрібною брижами" поверхні. Ступінь згладжування (усереднення) регулюється параметром $\beta=[0...1)$. Ітерації $t=0,1,2,..$ методу SGD при $\langle\mathbf{g}\rangle_{(-1)}=0$ мають вигляд:

$$ \mathbf{w}^{(t+1)} = \mathbf{w}^{(t)} - \lambda \,\langle\mathbf{g}\rangle^{(t)},~~~~~~~~~~~~ \langle\mathbf{g}\rangle^{(t)} = \beta\,\langle\mathbf{g}\rangle^{(t-1)} + (1-\beta)\,\mathbf{g}^{(t)}. $$

Величини batch_size і $\beta$ скорельовані. Чим вони більші, тим сильніше згладжений градієнт. Збільшення batch_size уповільнює навчання (за епоху робиться менше ітерацій). Однак на графічних карточках, завдяки швидкому перемноженню великих матриць, збільшення batch_size може зменшувати час проходу по всіх даних, навіть у розрахунку на одну ітерацію.

☝ Ковзне середнє градієнта має такий зміст. З вагами $\beta$ і $1-\beta$ (сума яких дорівнює 1) складаються середнє значення $\langle\mathbf{g}\rangle^{(t-1)}$ (отримане на попередній ітерації) і поточний градієнт $\mathbf{g}^{(t)}$. Чим ближче $\beta$ до 1, тим сильніше згладжування $\langle\mathbf{g}\rangle^{(t)}\approx \langle\mathbf{g}\rangle^{(t-1)}$. При $\beta=0$ усереднення не відбувається і $\langle\mathbf{g}\rangle^{(t)} = \mathbf{g}^{(t)}$.

Нижче наведено приклад згладжування даних за допомогою ковзного середнього. Синя лінія - це одновимірне стохастичне випадкове блукання (подібне ціні на акцію). Гладкі жовта і зелена лінії - це усереднення з різним параметром $\beta$.

Зрозуміло, що чим більше усереднення (ближче $\beta$ до одиниці), тим сильніше відстає середнє від вихідних даних.

Якщо усереднювана величина дорівнює константі $\mathbf{g}^{(t)}=\mathbf{g}=\text{const}$, то її середнє, через достатньо велике число ітерацій, також вийде на цю константу (сума геометричної прогресії): $\langle\mathbf{g}\rangle^{(1)} = \beta\,\mathbf{0} + (1-\beta)\,\mathbf{g}$, $$ \langle\mathbf{g}\rangle^{(2)} =\beta\,\langle\mathbf{g}\rangle^{(1)}+ (1-\beta)\,\mathbf{g}=(1+\beta)(1-\beta)\mathbf{g},....~~~~~~~~ \langle\mathbf{g}\rangle^{(t)} = (1+\beta+...+\beta^{t-1})(1-\beta)\,\mathbf{g} = (1-\beta^t)\,\mathbf{g}. $$

☝ У бібліотеках PyTorch і Keras параметр $\beta$ називається momentum, а швидкість навчання lr - це $\lambda\,(1-\beta)$ і $\langle\mathbf{g}\rangle^{(t)} \mapsto (1-\beta)\langle\mathbf{g}\rangle^{(t)}$. Крім цього, метод SGD "includes support for momentum, learning rate decay, and Nesterov momentum". Розпад означає зменшення lr з кожною епохою iterations відповідно до параметра розпаду decay (код з Keras): 

lr_t = lr * (1. / (1. + decay * iterations ))
Далі відбуваються обчислення: 
v = momentum * v - lr_t * grad
if nesterov:  params = params + momentum * v - lr_t * grad
       else:  params = params + v

За замовчуванням прийняті такі параметри: SGD(lr=0.01, momentum=0.0, decay=0.0, nesterov=False).

Обчислення градієнта

Знайдемо градієнт від середньоквадратичної помилки (MSE) у лінійній моделі. Для цього скористаємося правилом диференціювання складеної функції. У даному випадку $L$ залежить від $y_{i\gamma}$, які вже безпосередньо залежать від параметрів: $y_{i\gamma}=\sum_\mu\,x_{i\mu}w_{\mu\gamma}+b_{\gamma}$. Тому для середньоквадратичної помилки $L= \langle (\mathbf{y}-\hat{\mathbf{y}})^2\rangle$, усередненої по $N$ прикладах (перший індекс) і $m$ виходах моделі (другий індекс) маємо:

$$ \frac{\partial L}{\partial w_{\alpha\beta}} = \sum_{i, \gamma} \frac{\partial L}{\partial y_{i\gamma}} \, \frac{\partial y_{i\gamma}}{\partial w_{\alpha\beta}} = \frac{2}{N\,m} \sum_{i, \gamma} (y_{i\gamma}-\hat{y}_{i\gamma})\,\sum_\mu\,x_{i\mu}\,\delta_{\mu\alpha}\,\delta_{\gamma\beta} = \frac{2}{N\,m} \sum_{i} x_{i\alpha}\, (y_{i\beta}-\hat{y}_{i\beta}). $$

Аналогічно для градієнта за $\mathbf{b}$:

$$ \frac{\partial L}{\partial b_{\alpha}} = \frac{2}{N\,m} \sum_{i} (y_{i\alpha}-\hat{y}_{i\alpha}). $$

У матричному вигляді це можна записати так:

$$ \frac{\partial L}{\partial \mathbf{w} } = \frac{2}{N\,m}~\mathbf{x}^\top\cdot(\mathbf{y}-\hat{\mathbf{y}}),~~~~~~~~~~~~~ \frac{\partial L}{\partial \mathbf{b} } = \frac{2}{N\,m}~\mathrm{sum}(\mathbf{y}-\hat{\mathbf{y}},~ \mathrm{axis}=0). $$

Добуток $N\,m$ дорівнює числу елементів масиву $\mathbf{y}$, які в numpy знаходяться в атрибуті size.

Чисельний пошук мінімуму

Наведемо алгоритм пошуку оптимальних параметрів методом SGD за допомогою бібліотеки numpy. Помилку будемо обчислювати по batch_size прикладах, проходячи по всіх навчальних прикладах epochs разів: 

def Loss(X,Y, W, B):                                     # функція помилки
    return np.mean((X @ W + B - Y) ** 2)

def My_SGD(X, Y, lr, mo, batch_size, epochs):
    
    W = np.random.random ((nX, nY))                      # випадкові початкові значення
    B = np.random.random ((nY, )  )                      # параметрів моделі W,B

    agW, agB, iters = 0, 0, int( len(X)/batch_size )     # середні градієнти, число батчів
        
    for epoch in range(epochs):                          # епоха - прохід по всіх прикладах
        idx = np.random.permutation( len(X) )            # перемішаний список індексів
        X, Y = X[idx], Y[idx]                            # перемішуємо дані 
        
        for i in range(0, iters*batch_size, batch_size): # приклади розбиті на пачки 
            xb = X[i: i+batch_size]                      # входи пачки
            yb = Y[i: i+batch_size]                      # виходи пачки (цільові значення)
        
            y = xb @ W + B                               # виходи моделі
        
            gW = (2 / y.size) * xb.T @ (y - yb)          # градієнт за w
            gB = (2 / y.size) * np.sum(y - yb, axis=0)   # градієнт за b
        
            agW = mo*agW + (1-mo)*gW;   W -= lr * agW    # згладжуємо градієнти і
            agB = mo*agB + (1-mo)*gB;   B -= lr * agB    # підправляємо параметри
        
    return W, B, Loss(X,Y, W, B)                         # параметри і помилка

Звернімо увагу, що навчальні приклади перед початком кожної епохи випадково перемішуються. Для цього, за допомогою функції np.random.permutation створюється перемішаний список цілих чисел від 0 до len(X)-1. Наступний рядок власне виробляє перемішування. Перемішувати дані перед кожною епохою в принципі не обов'язково, але якщо даних небагато це може поліпшити збіжність до мінімуму. Крім цього, якщо число прикладів N націло не ділиться на batch_size, не всі дані потраплять в обчислення градієнта. Перемішування ослаблює цю проблему.

Справа наведено графіки помилки при різних швидкостях навчання і mo=0. Критичне значення lr знаходиться в районі 1.75. При наближенні до цього значення швидкість навчання помітно зростає.

Оптимальні гіперпараметри

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

Намалюємо карту висот помилки Loss лінійної моделі, як функцію гіперпараметрів $\lambda$ (=lr) і $\beta$ (=mo).
Синій колір - мінімум помилки, коричневий - максимум. Кількість епох epochs = 10 і розмір пачки batch_size = 10 фіксовані (точок N=100). Блакитна лінія з точками означає значення $\beta$ відповідне мінімальній помилці при даному $\lambda$. Під картою висот наведено графік помилки як функції $\lambda$ (при оптимальному $\beta$ для даного $\lambda$). Коричневий колір на карті висот відповідає області нестійкості градієнтного методу:

Як видно, для такої простої моделі, при $\beta \sim 0.8$, значення $\lambda$ може змінюватися в широких межах. Оптимальні гіперпараметри: mo=0.63, lr=5.31, loss=0.00013. Зазначимо, що підходяща параметризація гіперпараметрів важлива.
Наприклад, якщо вибрати як у Keras lr=$\lambda\,(1-\beta)$, то вийде менш стійка картина:

Оптимальні гіперпараметри в цьому випадку: mo=0.58, lr=2.25, loss=0.00013.

Градієнт на графі в PyTorch

Бібліотека PyTorch від Facebook є популярним інструментом при розробці складних архітектур нейронних мереж. PyTorch створює динамічний обчислювальний граф (будується в процесі обчислень і не вимагає компіляції).

Відтворимо на PyTorch простий градієнтний метод. Параметри моделі W,B є термінальними вузлами обчислювального графа, для яких буде обчислюватися градієнт (requires_grad=True):

Метод randn повертає тензор нормально розподілених випадкових чисел. За замовчуванням вони мають тип float32. Щоб він збігався з типом навчальних даних X,Y, відбувається явне задання типу аргументом dtype.

Потім йде основний цикл обчислень: 

import torch

def grad_torch(X, Y, lr=1, batch_size=10, epochs=10):
    W = torch.randn(nX, nY, dtype=torch.float64,  requires_grad=True)
    B = torch.randn(nY,     dtype=torch.float64,  requires_grad=True)

    for epoch in range(epochs):                          # епоха - прохід по всіх прикладах
        
        idx = torch.randperm( len(X) )                   # перемішаний список індексів
        X, Y = X[idx], Y[idx]                            # мішаємо дані
        
        sumL, iters = 0, int( len(X)/batch_size)         # сумарна помилка і число батчів 
        
        for i in range(0, iters*batch_size, batch_size): # приклади розбиті на пачки 
         
            xb = torch.from_numpy(X[i: i+batch_size])
            yb = torch.from_numpy(Y[i: i+batch_size])
            
            y = xb.mm(W).add(B)                          # модель y = bx @ W + B
            loss = ((y-yb)**2).mean()                    # mse помилка по батчу
            sumL += loss.data.item()

            loss.backward()                              # обчислення градієнтів
                                                
            with torch.no_grad():   
                W.add_(- lr *  W.grad)                   # без перебудови графа
                B.add_(- lr *  B.grad)            
            
            W.grad.zero_()                               # обнуляємо градієнти
            B.grad.zero_()                               
            
        print(f"{epoch:3d}({iters})  loss:{sumL/iters:0.5f}")

    return W, B

grad_torch(X, Y)

Спочатку формуються torch-тензори пачок з numpy-тензорів (нова пам'ять при цьому не виділяється). Потім будується обчислювальний граф, коренем якого буде помилка моделі loss. Попутно відбуваються обчислення на графі. Потім для кореня запускається метод backward, який повертається по гілках і обчислює градієнти.

Будь-який вираз з тензорами, що мають властивість requires_grad=True призводить до динамічної перебудови обчислювального графа, що веде до його кореня (помилки моделі loss). Щоб цього не сталося, після backward встановлюється оточення no_grad(), яке блокує зміну графа При цьому, змінюючи параметри, ми не створюємо нових тензорів (функція add_ це "інкрементація"). Після зміни W,B необхідно скинути в нуль градієнти (для наступної ітерації).

Нейронна мережа на PyTorch

Параметри лінійної моделі можна також знайти за допомогою нейронної мережі з одного шару без активаційної функції: 

import torch.nn as nn

model = nn.Sequential( nn.Linear(nX, nY) )

Створимо SGD оптимізатор (який буде підправляти параметри) і mse функцію помилки: 

optimizer = torch.optim.SGD(model.parameters(), lr=1, momentum=0.9)
criterion = nn.MSELoss()

Цикл навчання виглядає таким чином. На кожній ітерації на мережу передають пачку прикладів y=model(xb), у результаті чого відбувається пряме поширення (forward). Отриманий вихід мережі y, разом з "істинними" значеннями yb передають функції помилки. Для неї викликають зворотне поширення backward, у результаті якого будуть обчислені градієнти. Метод оптимізатора step, за допомогою цих градієнтів, підправляє параметри моделі. Потім оптимізатор обнуляє градієнти (zero_grad): 

iters = int( len(X)/batch_size )
for epoch in range(epochs):                           # епоха - прохід по всіх прикладах
    for i in range(0, iters*batch_size, batch_size):  # приклади розбиті на пачки 
        bx = torch.Tensor( X[i: i+batch_size] )       # по numpy масивах створюємо
        by = torch.Tensor( Y[i: i+batch_size] )       # тензори torch типу float32
    
        y = model(bx)                                 # пряме поширення
        
        loss = criterion(y, by)                       # обчислюємо помилку
        loss.backward()                               # обчислюємо градієнти
         
        optimizer.step()                              # підправляємо параметри
        optimizer.zero_grad()                         # обнуляємо градієнти
            
    print('epoch: %d  Loss: %.6f'  % (epoch,  loss) )

Після навчання можна вивести параметри моделі: 

print(model[0].weight, model[0].bias)
Матриця ваг weight, як і в sklearn, зберігається в транспонованому вигляді.

Нейронна мережа на Keras

Бібліотека Keras є простим у використанні інструментом для проектування нейронних мереж. У даний час вона є складовою частиною фреймворка tensorflow від Google.
Вступ до Keras можна знайти в документі NN_Base_Keras.

Лінійна модель (один шар без активаційної функції) в Keras реалізована в повнозв'язному шарі Dense. Створення нейронної мережі з одним таким шаром робиться таким чином: 

from   tensorflow       import keras              # keras з tensorflow
from   keras.models     import Sequential         # спосіб формування шарів (стопка)
from   keras.layers     import Dense              # повнозв'язний шар
from   keras.optimizers import SGD                # метод оптимізації

model = Sequential()                              # лінійна стопка шарів
model.add(Dense(units=nY, input_dim=nX))

Після створення моделі, вона компілюється (compile) і запускається її навчанням (fit). При навчанні будемо використовувати SGD-оптимізатор і mse-помилку: 

model.compile(optimizer = SGD(lr=2, momentum=0.5), loss = 'mse') 

res = model.fit(X, Y, batch_size=batch_size, epochs=10,  verbose=0 ) )
Після навчання можна вивести графік помилок, як функцію числа епох: 
plt.plot(res.history['loss'], marker=".")
plt.legend(["loss: %.5f"  % ( res.history['loss'][-1] )])
plt.show()
Параметри кожного шару (у нашому випадку шар один) знаходяться в списку, отримуваному методом шару get_weights(). Тому ваги і зміщення лінійної моделі (або в загальному випадку багатошарової повнозв'язної моделі) виводяться таким чином: 
for lr in model.layers:                           # по всіх шарах
    w = lr.get_weights()[0]                       # ваги синапсів нейронів
    b = lr.get_weights()[1]                       # зміщення нейронів
    ...