Лекція 02 · Лінійна регресія та оптимізація

У першій лекції ми розібрали, з яких задач складається машинне навчання. Сьогодні побудуємо першу повноцінну модель — лінійну регресію — і три різні способи її навчити. Робочий матеріал: 20 640 округів Каліфорнії, для кожного треба передбачити медіанну вартість житла.

Зміст лекції

Лекція проходить один повний цикл: від сирих даних до навченої моделі та розбору її помилок.

1. Дані й задача регресії. Набір California Housing, розподіли ознак і цільової змінної, вади реальних даних.
2. Модель. Лінійна регресія з однією та багатьма ознаками, залишки, матричний запис.
3. Функція втрат. MSE і MAE як два різні погляди на ціну помилки, стійкість до викидів, опуклість.
4. Оптимізація. Перебір по сітці, градієнтний спуск, вибір кроку навчання, роль масштабування ознак.
5. Точний розв'язок. Нормальні рівняння, межі застосовності закритої формули.
6. Стохастичний градієнтний спуск. Як навчатися, коли зразків мільйони.

Порядок не випадковий: спершу з'являється модель, потім міра її помилки, і лише тоді — способи цю помилку мінімізувати.

Чого ви навчитеся

Після цієї лекції ви зможете:

• записати лінійну регресію в матричній формі й пояснити роль стовпця одиниць;
• пояснити, чим квадратична втрата відрізняється від абсолютної та в якій ситуації доречна кожна;
• вивести градієнт середньоквадратичної втрати й реалізувати градієнтний спуск на numpy;
• дібрати крок навчання і розпізнати за кривою втрат, що він завеликий або замалий;
• пояснити, чому стандартизація ознак скорочує кількість ітерацій у десятки разів;
• розв'язати нормальні рівняння через numpy.linalg.solve і назвати випадки, коли єдиного розв'язку немає;
• обґрунтувати перехід до стохастичного градієнтного спуску на великих вибірках;
• прочитати графік залишків і побачити на ньому систематичну ваду моделі.

Регресія: коли відповідь є числом

У задачі регресії (regression) модель має повернути число: ціну квартири, температуру завтра, час доставки, дозу препарату. Це відрізняє її від класифікації, де відповідь — одна з кількох категорій.

Формально нам дано $N$ пар «об'єкт — відповідь», і треба побудувати функцію, яка для нового об'єкта поверне правдоподібне число.

Позначення	Значення
$N$	кількість спостережень у вибірці
$D$	кількість ознак одного спостереження
${\mathbf{x}}_n\in {\mathbb{R}}^D$	вектор ознак $n$-го спостереження
$y_n\in \mathbb{R}$	справжня відповідь для нього
${\hat{y}}_n$	прогноз моделі

Ця нотація витримана в усіх чотирнадцяти лекціях курсу. У багатьох підручниках ваги позначають через $\theta$, а втрати через $J$ — це та сама річ під іншими літерами.

Дані лекції: 20 640 округів Каліфорнії

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

California Housing. 20 640 груп кварталів («переписних округів») за переписом населення США 1990 року, вісім числових ознак на кожну: медіанний дохід мешканців, медіанний вік будинків, середня кількість кімнат і спалень на домогосподарство, населення, середня кількість мешканців на домогосподарство, широта й довгота. Передбачаємо медіанну вартість житла в окрузі, виміряну в сотнях тисяч доларів.

Розподіл цільової змінної скошений праворуч: медіана 1,80, тобто 180 тисяч доларів, а середнє помітно вище — 2,07. Медіанний дохід — найсильніше пов'язана з ціною ознака, кореляція 0,688.

Джерело даних: scikit-learn · California Housing

Матеріали до розділу

• ⌨Géron A. Hands-On Machine Learning — наскрізний проєкт на цьому самому наборікод

  Той самий California Housing розібрано від завантаження до розгортання моделі; ноутбук виконуваний.

• 📊scikit-learn · fetch_california_housingдані

  Опис усіх восьми ознак і застереження про цензуровану цільову змінну.

Стеля 5,0: цільову змінну обрізано зверху

Перш ніж навчати будь-що, варто подивитися на праву частину розподілу. Там видно вертикальну стіну: 992 округи, тобто 4,8 % вибірки, впираються у стелю 5,00.

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

Наслідок буде видно наприкінці лекції: модель систематично занижує ціну там, де вона впирається у стелю, і жодне покращення оптимізації цього не виправить. Медіанний дохід обрізано так само — 51 округ сидить на максимумі 15,0. (У сирому файлі стеля записана двома близькими значеннями, 5,0 і 5,00001, — це артефакт округлення під час підготовки набору, а не два різні рівні цін.)

Викиди: округи, яких майже не буває

Три ознаки набору мають дуже важкі хвости. Середня кількість кімнат на домогосподарство сягає 141,9, а середня кількість мешканців — 1243. Обидва числа не є помилками: перше походить від курортного селища з тридцятьма мешканцями і майже порожнім житловим фондом, друге — від гуртожитку, де сім із половиною тисяч людей припадають на шість домогосподарств.

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

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

Кореляційна матриця: що з чим пов'язано

Наступний обов'язковий крок — подивитися на попарні зв'язки. Матриця кореляцій відповідає одразу на два питання: які ознаки варті уваги і які з них дублюють одна одну.

Ціна найсильніше пов'язана з медіанним доходом, 0,69. Решта ознак поодинці пояснюють мало: вік будинків дає 0,11, кількість кімнат 0,15.

Важливіші за це дві пари всередині самих ознак. Кількість кімнат і кількість спалень корелюють на 0,85, широта й довгота — на −0,92 (Каліфорнія витягнута з північного заходу на південний схід). Такі пари створюють мультиколінеарність (multicollinearity), через яку ваги моделі стають нестійкими; до цього ми повернемось після нормальних рівнянь.

Модель з однією ознакою

Найпростіша модель проводить через хмару точок пряму.

"Linear regression is a model that assumes a linear relationship between inputs and the output." — EPFL CS-433, лекція 01c

$${\hat{y}}_n=w_0+w_1{x}_{n1}$$

Тут $w_1$ — нахил: на скільки змінюється прогноз, коли ознака зростає на одиницю; $w_0$ — зсув (offset, або bias): значення прогнозу при нульовій ознаці. Двох чисел досить, щоб повністю задати модель.

Різницю між справжнім значенням і прогнозом називають залишком (residual): $e_n=y_n-{\hat{y}}_n$. На рисунку залишки показані вертикальними відрізками — саме їх модель намагатиметься зменшити.

Джерело: Dive into Deep Learning · Linear Regression

Пряма найменших квадратів на реальних даних

Проведемо пряму через усі 20 640 округів, узявши єдину ознаку — медіанний дохід. Результат:

$$\hat{y}=\mathrm{0,418}\cdot \text{дохід}+\mathrm{0,451}$$

Читається це так: кожна додаткова одиниця медіанного доходу піднімає прогноз ціни на 0,418 сотні тисяч доларів, тобто приблизно на 42 тисячі. Середня абсолютна помилка такої моделі — 0,626, або 63 тисячі доларів на округ.

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

Модель із кількома ознаками

Одна ознака пояснює далеко не все. Додамо решту сім:

$${\hat{y}}_n=w_0+w_1{x}_{n1}+\cdots +w_D{x}_{nD}={\tilde{\mathbf{x}}}_n^{\mathrm{\top }}\tilde{\mathbf{w}}$$

Компактний запис праворуч отримано простим трюком: до вектора ознак дописують одиницю, ${\tilde{\mathbf{x}}}_n=(1,x_{n1},\dots ,x_{nD}{)}^{\mathrm{\top }}$, а до вектора ваг — зсув, $\tilde{\mathbf{w}}=(w_0,w_1,\dots ,w_D{)}^{\mathrm{\top }}$. Тепер зсув перестає бути окремим випадком і живе разом з іншими вагами. Далі тильду опускаємо, і вся вибірка записується однією матрицею $\mathbf{X}$ розміру $N\times (D+1)$, а всі прогнози — одним добутком $\mathbf{X}\mathbf{w}$.

Схема праворуч показує ту саму модель у мові нейронних мереж: входи $x_1,\dots ,x_d$ сходяться в єдиний вихід. Лінійна регресія — це нейронна мережа з одного нейрона без функції активації, і в лекції 10 ми добудуємо до неї приховані шари.

Джерело: Dive into Deep Learning · Linear Regression

Лінійна за параметрами, а не за ознаками

Слово «лінійна» стосується ваг, а не вхідних змінних. Ніщо не заважає подати в ту саму модель квадрат ознаки, її логарифм чи добуток двох ознак: модель залишиться лінійною за $\mathbf{w}$, а отже всі формули цієї лекції працюватимуть без змін.

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

Плата за гнучкість видна вже тут: криві високого степеня починають повторювати шум замість закономірності. Це перенавчання (overfitting), і воно є темою лекції 03.

Джерело: scikit-learn · Polynomial and Spline interpolation

Функція втрат задає, що означає «добра модель»

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

"A loss function (or energy, cost, training objective) is used to learn parameters that explain the data well. The loss function quantifies how well our model does — or in other words how costly our mistakes are." — EPFL CS-433, лекція 01d

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

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

Середньоквадратична втрата

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

$$\mathcal{L}(\mathbf{w})=\frac{1}{2N}\sum _{n=1}^N(y_n-{\mathbf{x}}_n^{\mathrm{\top }}\mathbf{w}{)}^2=\frac{1}{2N}{\mathbf{e}}^{\mathrm{\top }}\mathbf{e},\mathbf{e}=\mathbf{y}-\mathbf{X}\mathbf{w}$$

Позначення	Значення
$\mathbf{e}\in {\mathbb{R}}^N$	вектор залишків, по одному на спостереження
${\mathbf{e}}^{\mathrm{\top }}\mathbf{e}$	сума квадратів залишків
$1/N$	усереднення, щоб втрата не залежала від розміру вибірки
$1/2$	зручність: двійка скоротиться при диференціюванні

Множник одна друга — суто технічна деталь. У частині джерел середньоквадратичну втрату (mean squared error, MSE) пишуть без нього; на положення мінімуму це не впливає, тому далі скрізь використовуємо запис із $\frac{1}{2N}$.

Числовий приклад на трьох округах із залишками $+\mathrm{1,086}$, $-\mathrm{0,218}$ і $-\mathrm{0,584}$: сума квадратів дорівнює 1,5676, середнє 0,5225, а втрата з множником одна друга — 0,2613.

Форма штрафу вирішує все

Порівняємо три способи перетворити залишок на штраф. Квадратична втрата $e^2$ зростає прискорено: помилка вдвічі більша коштує вчетверо дорожче. Абсолютна втрата $|e|$ зростає рівномірно. Функція Хубера поводиться як квадратична біля нуля і як абсолютна далеко від нього:

$$\mathrm{Huber}(e)=\{\begin{array}{ll}\frac{1}{2}{e}^2,& |e|\le \delta \\ \delta |e|-\frac{1}{2}{\delta }^2,& |e|>\delta \end{array}$$

При $e=3$ квадратична втрата дає 9, абсолютна — 3, Хубера при $\delta =1$ — 2,5. Саме тому квадратична втрата так чутлива до викидів: одне спостереження з помилкою 10 важить у ній стільки ж, скільки сто спостережень з помилкою 1.

Середня абсолютна втрата (mean absolute error, MAE) від цієї вади вільна, але має свою: у нулі вона не диференційовна, тому звичайний градієнт для неї не визначений і доводиться користуватися субградієнтом.

Один викид зсуває мінімум

Найпростіша модель узагалі не має ознак: вона всім передбачає одне й те саме число $w_0$. Візьмімо чотири спостереження зі значеннями 1, 2, 3, 4 і подивимося, де лежить мінімум кожної втрати.

Квадратична втрата мінімальна в точці 2,5 — це середнє арифметичне. Абсолютна втрата має мінімумом весь відрізок від 2 до 3, і будь-яка точка з нього однаково добра; за домовленістю беруть медіану.

Тепер додамо п'яте значення 20. Мінімум квадратичної втрати їде на 6,0 — модель відходить від усіх чотирьох звичайних спостережень, аби трохи наблизитися до одного дикого. Мінімум абсолютної втрати зсувається лише до 3,0 і залишається серед даних.

Куди поїде мінімум, якщо потягнути викид

Чотири спостереження — 1, 2, 3, 4 — і п'яте, яке ви пересуваєте самі. Модель без ознак: усім передбачає одне число, тож уся задача — знайти мінімум втрати по цьому числу.

П'яте значення: 20,00Поріг Хубера δ = 1,00

MSEMAEХубер

6,00мінімум MSE — середнє

3,00мінімум MAE — медіана

3,00мінімум Хубера, δ = 1,00

Тягніть п'яте значення від 4 до 20 і дивіться на перші два числа. Мінімум квадратичної втрати їде слідом за викидом: він завжди дорівнює середньому, а середнє тягнеться за будь-яким одним спостереженням як завгодно далеко. Мінімум абсолютної втрати дорівнює медіані й лишається серед даних — на 3,0, хоч би куди поїхав викид. Хубер стоїть посередині: при малому δ він майже повторює MAE, при великому зливається з MSE. Саме тому вибір функції втрат — це рішення про те, скільки влади ви віддаєте найдальшій точці.

Стійка пряма на даних із викидами

Той самий ефект на двовимірних даних. Зелена пряма — звичайні найменші квадрати — розвертається до чотирьох викидів; прямі за втратою Хубера тримаються тренду. Параметр $\delta$ (у scikit-learn — epsilon) задає, з якої помилки починається аномалія.

Тепер те саме на наших двадцяти тисячах округів — чотири способи поруч:

Спосіб	MAE на тесті	Вага при AveOccup
МНК	0,535	−0,031
Хубер	0,533	−0,761
RANSAC	0,530	−0,370
обрізання 1 %	0,623	−0,271

Точність майже не змінилася, і це головне: стійкі методи купують тут не якість, а правильну інтерпретацію. МНК каже, що кількість мешканців на домогосподарство на ціну майже не впливає — вага −0,031. Але цю вагу він отримав, підлаштовуючись під гуртожиток із 1243 мешканцями; щойно Хубер зменшує внесок таких округів, проступає справжня залежність — чим тісніше живуть, тим дешевше житло.

Останній рядок застерігає: обрізання хвостів прибрало 330 округів і погіршило помилку до 0,623.

Джерело: scikit-learn · HuberRegressor vs Ridge

Опуклість прибирає погані локальні мінімуми

Те, що ми щойно виміняли на стійкість, має точну назву. Функція $h$ опукла (convex), якщо відрізок, проведений між будь-якими двома точками її графіка, лежить не нижче за сам графік:

$$h(\alpha \mathbf{u}+(1-\alpha )\mathbf{v})\le \alpha h(\mathbf{u})+(1-\alpha )h(\mathbf{v}),0\le \alpha \le 1$$

"A strictly convex function has a unique global minimum. For convex functions, every local minimum is a global minimum." — EPFL CS-433, лекція 01d

Для нас це означає таке: у лінійної регресії з квадратичною втратою немає поганих локальних мінімумів — кожен локальний мінімум є глобальним. Сума опуклих функцій опукла, а квадрат лінійної функції опуклий, тому вся втрата опукла за $\mathbf{w}$.

Обережність: опуклість не дає єдиності точки мінімуму — для цього потрібен повний стовпцевий ранг $\mathbf{X}$. При лінійно залежних ознаках мінімальних $\mathbf{w}$ нескінченно багато (див. розділ про межі закритої формули), а сам спуск сходиться лише за достатньо малого кроку.

Абсолютна втрата теж опукла; втрачається лише гладкість.

Матеріали до розділу

• 📄Boyd S., Vandenberghe L. Convex OptimizationPDF

  Стандартний підручник; розділ 3 — означення опуклості та її наслідки для оптимізації. Повний текст на сайті авторів.

Перебір по сітці не масштабується

Найпряміший спосіб знайти ваги — перебрати їх: накриваємо простір параметрів сіткою, рахуємо втрату в кожному вузлі й беремо найкращий. На рисунку так шукають дві ваги моделі з однією ознакою: сітка 7 × 7 дає 49 вузлів і втрату 0,3675 проти точного мінімуму 0,3506.

Проблема в тому, що десять значень на кожен вимір при $D$ параметрах дають ${10}^D$ обчислень: для двох ваг це сто, для наших дев'яти — мільярд, а сучасні моделі мають мільйони параметрів. Перебір придатний лише як ілюстрація і як спосіб побачити ландшафт втрат очима.

Йдеться про перебір самих ваг моделі, а не гіперпараметрів — це різні речі з однаковою назвою.

Градієнт вказує напрям найшвидшого зростання

Замість перебору скористаємося похідними.

"A gradient (at a point) is the slope of the tangent to the function (at that point). It points to the direction of largest increase of the function." — EPFL CS-433, лекція 02a

$$\mathrm{\nabla }\mathcal{L}(\mathbf{w})={(\frac{\partial \mathcal{L}}{\partial w_0},\dots ,\frac{\partial \mathcal{L}}{\partial w_D})}^{\mathrm{\top }}$$

Якщо градієнт показує вгору, то щоб спускатися, треба йти проти нього. Звідси градієнтний спуск (gradient descent):

$${\mathbf{w}}^{(t+1)}:={\mathbf{w}}^{(t)}-\eta \mathrm{\nabla }\mathcal{L}({\mathbf{w}}^{(t)})$$

Тут $t$ — номер ітерації, а $\eta >0$ — крок навчання (step size, або learning rate). Це головне число, яке доводиться підбирати вручну, і воно вирішує, чи навчиться модель узагалі.

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

Градієнт середньоквадратичної втрати

Для лінійної моделі градієнт рахується в один рядок. Продиференціюємо $\mathcal{L}=\frac{1}{2N}{\mathbf{e}}^{\mathrm{\top }}\mathbf{e}$ за $\mathbf{w}$, пам'ятаючи, що $\mathbf{e}=\mathbf{y}-\mathbf{X}\mathbf{w}$:

$$\mathrm{\nabla }\mathcal{L}(\mathbf{w})=-\frac{1}{N}{\mathbf{X}}^{\mathrm{\top }}\mathbf{e}=\frac{1}{N}{\mathbf{X}}^{\mathrm{\top }}(\mathbf{X}\mathbf{w}-\mathbf{y})$$

Покомпонентно це виглядає ще прозоріше: $\partial \mathcal{L}/\partial w_0=-\frac{1}{N}\sum _n{e}_n$, а $\partial \mathcal{L}/\partial w_1=-\frac{1}{N}\sum _n{e}_n{x}_{n1}$. Тобто кожна вага зсувається на середню помилку, зважену відповідною ознакою: якщо модель систематично занижує прогноз там, де ознака велика, відповідна вага зросте.

Вартість одного кроку — $O(ND)$: стільки коштує один прохід по всій матриці ознак. Саме цей множник $N$ стане проблемою наприкінці лекції.

Ця сама структура «прогноз мінус мітка, зважена на ознаки» повториться в логістичній регресії в лекції 05 і в нейронних мережах у лекції 10.

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

Двадцять рядків, у яких немає нічого, крім формул із попереднього розділу.

import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.preprocessing import StandardScaler

data = fetch_california_housing()
X = StandardScaler().fit_transform(data.data)  # 20640 × 8
X = np.c_[X, np.ones(len(X))]  # стовпець одиниць під w0
y = data.target  # ціна, сотні тисяч доларів


def loss(w):  # L = (1/2N)·eᵀe
    e = y - X @ w
    return 0.5 * np.mean(e**2)


def gradient(w):  # ∇L = −(1/N)·Xᵀe
    return -X.T @ (y - X @ w) / len(y)


w = np.zeros(X.shape[1])
for step in range(201):
    if step % 50 == 0:
        print(f"крок {step:3d}   втрата {loss(w):.4f}")
    w -= 0.5 * gradient(w)

Вивід код запущено, вивід не редагувався

крок   0   втрата 2.8052
крок  50   втрата 0.2637
крок 100   втрата 0.2622
крок 150   втрата 0.2622
крок 200   втрата 0.2622

На що дивитися: втрата падає з 2,8052 до 0,2637 уже за 50 кроків, а далі змінюється лише четвертий знак — спуск дійшов до мінімуму.

Початкове значення 2,8052 — це втрата моделі з нульовими вагами, яка всім передбачає нуль.

Спуск у просторі ваг і в просторі даних

Один і той самий процес корисно бачити з двох боків одразу.

Ліворуч — простір параметрів. Кожна точка площини є парою $(w_0,w_1)$, кольором показано втрату, а червона ламана — сорок кроків спуску з початкової точки. Видно характерну поведінку: перші кроки великі, бо градієнт крутий, далі траєкторія входить у пласку долину і сповільнюється.

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

Втрата за ці сорок кроків впала з 0,828 до 0,3506 — того самого значення, яке дає точна формула.

Завеликий крок розхитує спуск

Крок навчання не можна взяти навмання. Класична демонстрація на функції $f(x)=x^2$ зі стартом у точці $x=10$: при $\eta =0,2$ десять ітерацій майже приводять у нуль, при $\eta =0,05$ за ті самі десять ітерацій вдається дійти лише до $x\approx 3,49$.

На рисунку показано третій випадок, $\eta =1,1$. Кожен крок перестрибує мінімум і закидає точку далі, ніж вона була: після десяти ітерацій $x\approx 61,9$ замість нуля. Спуск не просто повільний — він розбігається.

Межа стійкості залежить від кривини функції втрат, а та — від даних. Універсального «правильного» значення $\eta$ немає, тому його підбирають експериментально, зазвичай пробуючи значення, що відрізняються в десять разів.

Джерело: Dive into Deep Learning · Gradient Descent

Темп навчання на власних очах

Та сама демонстрація, що на рисунку: парабола f(x) = x², старт у точці x = 10. Крок спуску — x ← x − η·f′(x) = x(1 − 2η).

Темп навчання η = 0,20Ітерацій: 10

0,06x після 10 ітерацій

0,60множник |1 − 2η| за крок

збігаєтьсяумова: |1 − 2η| < 1, тобто η < 1

Спуск збігається, поки множник за крок менший від одиниці. Для параболи це дає точну межу η < 1; за η = 0,5 мінімум береться за один крок, а далі кожен крок перестрибує його все далі. Для реальної функції втрат кривина залежить від даних — перемкніть режим і подивіться, що з цього виходить.

Крива втрат показує, чи вдалий крок

Підбирати $\eta$ наосліп не доводиться: усе видно на графіку втрати за ітераціями.

Замалий крок дає повільне, майже пряме спадання, яке не встигає дійти до мінімуму. Вдалий крок дає швидке падіння з виходом на пласке плато. Завеликий крок спричиняє передчасну зупинку на високому рівні: спуск скаче через мінімум і не може в нього потрапити. Зовсім великий крок веде до зростання втрати аж до переповнення.

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

Джерело: Stanford CS231n · Neural Networks 3

Масштаб ознак вирішує, скільки кроків потрібно

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

Ліворуч показано лінії рівня втрати для сирих ознак. Вони витягнуті у вузьку долину: число обумовленості (condition number) дорівнює 44, тобто вздовж однієї осі поверхня в сорок чотири рази крутіша, ніж уздовж іншої. Градієнт майже не має складової вздовж плаского напрямку, і спуск повзе по дну — 305 ітерацій.

Праворуч ті самі дані після стандартизації (standardization): кожну ознаку поділено на її стандартне відхилення, тож обидві мають однаковий розмах. Середнє віднято в обох випадках, тому різниця між панелями стосується саме масштабу. Лінії рівня стали майже колами, число обумовленості впало до 1,27, і того самого результату спуск досягає за 5 ітерацій.

Стандартизацію завжди навчають на тренувальній частині даних і лише потім застосовують до тестової.

Нормальні рівняння: точна відповідь без ітерацій

Лінійна регресія з квадратичною втратою — рідкісний випадок, коли мінімум виписується формулою. Прирівняємо градієнт до нуля:

$${\mathbf{X}}^{\mathrm{\top }}(\mathbf{y}-\mathbf{X}\mathbf{w})=\mathbf{0}⟺{\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\mathbf{w}={\mathbf{X}}^{\mathrm{\top }}\mathbf{y}$$

Ці рівності називають нормальними рівняннями (normal equations), а метод, який на них спирається, — методом найменших квадратів (least squares). Матрицю ${\mathbf{X}}^{\mathrm{\top }}\mathbf{X}$ розміру $(D+1)\times (D+1)$ називають матрицею Грама (Gram matrix).

Якщо матриця Грама оборотна, розв'язок єдиний:

$${\mathbf{w}}^{\star }=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathbf{y}$$

Метод має довгу історію. Першим його опублікував Лежандр 1805 року. Гаусс 1809-го застосував його до обчислення орбіт небесних тіл і додав імовірнісне обґрунтування, заявивши при цьому, що користувався методом іще з 1795 року, — звідси й давня суперечка про пріоритет.

Матеріали до розділу

• ▶MIT 18.06 · Linear Algebra, лекція про проєкції та найменші квадративідео

  Гілберт Стренг виводить нормальні рівняння геометрично — через проєкцію на простір стовпців.

• ⚙NumPy · linalg.lstsqдокументація

  Чому в реальному коді беруть саме lstsq, а не обертання матриці.

Межі застосовності закритої формули

Матриця Грама оборотна тоді й лише тоді, коли стовпці $\mathbf{X}$ лінійно незалежні. На практиці це порушується у двох ситуаціях.

Ознак більше, ніж спостережень. При $D>N$ ранг матриці менший за кількість стовпців завжди, і розв'язків стає нескінченно багато. Крайній випадок: через одну точку можна провести безліч прямих. Ліки — регуляризація, тема лекції 03.

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

У нашому наборі така пара теж є: широта й довгота з кореляцією −0,92.

Джерело: scikit-learn · Common pitfalls in the interpretation of coefficients

Нормальні рівняння на numpy

Головне правило реалізації: ніколи не обертати матрицю явно. Лінійний солвер розв'язує систему через розклад матриці, він точніший і щонайменше втричі дешевший за обчислення оберненої.

import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.preprocessing import StandardScaler

data = fetch_california_housing()
X = np.c_[StandardScaler().fit_transform(data.data), np.ones(20640)]
y = data.target

G = X.T @ X  # матриця Грама, 9 × 9
w = np.linalg.solve(G, X.T @ y)  # а не np.linalg.inv(G) @ X.T @ y

print("ранг X:", np.linalg.matrix_rank(X), "з", X.shape[1], "стовпців")
print("втрата в оптимумі:", round(0.5 * np.mean((y - X @ w) ** 2), 4))
for name, v in zip(list(data.feature_names) + ["зсув w0"], w, strict=True):
    print(f"  {name:<11} {v:+.3f}")

Вивід код запущено, вивід не редагувався

ранг X: 9 з 9 стовпців
втрата в оптимумі: 0.2622
  MedInc      +0.830
  HouseAge    +0.119
  AveRooms    -0.266
  AveBedrms   +0.306
  Population  -0.005
  AveOccup    -0.039
  Latitude    -0.900
  Longitude   -0.871
  зсув w0     +2.069

На що дивитися: втрата 0,2622 збігається з тією, до якої дійшов спуск, а найбільші за модулем ваги мають дохід і координати округу.

Та сама модель через scikit-learn

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

from sklearn.datasets import fetch_california_housing
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error, mean_squared_error
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X, y = fetch_california_housing(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)

model = make_pipeline(StandardScaler(), LinearRegression())
model.fit(X_train, y_train)
pred = model.predict(X_test)

print("MSE:", round(mean_squared_error(y_test, pred), 3))
print("MAE:", round(mean_absolute_error(y_test, pred), 3))

Вивід код запущено, вивід не редагувався

MSE: 0.529
MAE: 0.535

На що дивитися: MAE 0,535 означає, що прогноз ціни помиляється в середньому на 53,5 тисячі доларів; MSE 0,529 — середній квадрат тієї самої помилки.

Для порівняння: прогноз-константа, що всім видає середню ціну, дає MAE 0,907.

make_pipeline тут не косметика: стандартизація навчається лише на тренувальній частині, інакше відомості про тестову вибірку просочилися б у модель.

Джерело: scikit-learn · LinearRegression

Прогноз проти факту на відкладених округах

Кожна точка — один з 4 128 округів, яких модель не бачила; пунктир — лінія ідеального прогнозу. Хмара витягнута вздовж нього, отже головний зв'язок уловлено. Але 19 округів дістали від'ємну ціну, окремі прогнози сягають 7,16, а вздовж верхнього краю тягнеться смуга з 185 обрізаних округів.

Замовникові зазвичай потрібне не одне число, а межі, у яких воно тримається. Їх дає квантильна регресія (quantile regression): замість квадрата помилки мінімізується асиметричний модуль, тож модель наближає свій квантиль ціни, а не середнє. Три моделі — на 10-му, 50-му і 90-му квантилях — дають смугу, у яку ціна має потрапляти у 80 % випадків.

# продовження: спирається на попередні блоки лекції
from sklearn.linear_model import QuantileRegressor

band = {}
for q in (0.1, 0.5, 0.9):  # три окремі моделі
    band[q] = (
        make_pipeline(
            StandardScaler(), QuantileRegressor(quantile=q, alpha=0, solver="highs")
        )
        .fit(X_train, y_train)
        .predict(X_test)
    )

lo, mid, hi = band[0.1], band[0.5], band[0.9]
print("покриття смуги 10-90 :", round(((y_test >= lo) & (y_test <= hi)).mean(), 3))
print("середня ширина смуги :", round(float((hi - lo).mean()), 3))
print("MAE медіанної моделі :", round(mean_absolute_error(y_test, mid), 3))

Вивід код запущено, вивід не редагувався

покриття смуги 10-90 : 0.793
середня ширина смуги : 1.647
MAE медіанної моделі : 0.535

На що дивитися: смуга між 10-м і 90-м квантилями накрила 79,3 % тестових округів при обіцяних 80 %, а медіанна модель помиляється так само, як МНК.

Покриття 0,793 при обіцяних 0,80 — і жодного припущення про форму розподілу помилок для цього не знадобилося. MAE медіанної моделі 0,535, стільки ж, скільки в МНК: смуга дістається не ціною точності, а ціною втричі довшого навчання. Ширина смуги росте з прогнозом — 1,33 у найдешевшій чверті округів проти 2,06 у найдорожчій.

Графік залишків виявляє системну ваду

Той самий результат у координатах «прогноз проти залишку». Ідеальна картина — безструктурна хмара навколо нуля однакової товщини по всій ширині.

Ми бачимо інше. Ліворуч, де прогнози малі, залишки переважно додатні; праворуч вони йдуть у мінус — ознака того, що лінійної залежності недостатньо. Червоні точки на чіткій спадній прямій — обрізані округи: у них справжнє значення завжди 5,0, тому залишок дорівнює $5,0-\hat{y}$. Середній залишок на них становить $+1,20$ проти $-0,07$ на решті.

Розподіл помилок праворуч має середнє $-\mathrm{0,010}$ і стандартне відхилення 0,727, але симетричним його не назвеш: правий хвіст помітно довший — знову наслідок стелі.

Головне тут те, що товщина хмари неоднакова: біля дешевих округів вона вужча, біля дорогих ширша. Це гетероскедастичність (heteroscedasticity), і саме через неї одне число «промахується в середньому на 53,5 тисячі» описує ситуацію погано — для дешевого округу похибка менша, для дорогого більша, а середнє не стосується жодного з них. Смуга з попереднього розділу відповідає на це прямо.

Матеріали до розділу

• ⚙Common pitfalls in the interpretation of coefficients of linear modelsдокументація

  Розгорнутий приклад: чому коефіцієнти не можна читати як «важливість ознаки».

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

Кожен крок повного градієнта коштує проходу по всій вибірці. На двадцяти тисячах округів це непомітно, на мільярді записів — неприйнятно: перш ніж зробити один-єдиний крок, треба переглянути всі дані.

Вихід підказує сама структура втрати: вона є середнім по спостереженнях, $\mathcal{L}=\frac{1}{N}\sum _n{\mathcal{L}}_n$. Отже, градієнт, порахований на одному випадково вибраному спостереженні, є незміщеною оцінкою повного градієнта: ${\mathbb{E}}_n[\mathrm{\nabla }{\mathcal{L}}_n]=\mathrm{\nabla }\mathcal{L}$. Крок стає в $N$ разів дешевшим:

$${\mathbf{w}}^{(t+1)}:={\mathbf{w}}^{(t)}-\eta \mathrm{\nabla }{\mathcal{L}}_n({\mathbf{w}}^{(t)}),n\sim \text{рівномірно з }\{1,\dots ,N\}$$

Ціна — шум. На рисунку показано квадратичну функцію $f(x_1,x_2)=x_1^2+2x_2^2$ і п'ятдесят кроків стохастичного спуску з $\eta =0,1$. Повний градієнт на цій самій поверхні з тим самим кроком дає гладку криву; тут траєкторія смикається і зрештою блукає навколо мінімуму, не зупиняючись у ньому.

Джерело: Dive into Deep Learning · Stochastic Gradient Descent

Матеріали до розділу

• 📄Ruder S. An overview of gradient descent optimization algorithmsстаття

  Оглядова стаття про всі варіанти спуску — від пакетного до Adam, з анімаціями.

• ▶Distill · Why Momentum Really Worksдемо

  Інтерактивне пояснення, чому імпульс прискорює спуск у витягнутих ярах.

Крок мусить спадати, але повільно

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

Саме це показано на рисунку: при експоненційному спаданні ${\eta }_t=e^{-0,1t}$ навіть за тисячу кроків точка зупиняється в $(-0,76;-0,12)$ замість початку координат.

Правильний баланс задає умова Роббінса — Монро: сума кроків має розбігатися, а сума їхніх квадратів — збігатися.

$$\sum _t{\eta }^{(t)}=\mathrm{\infty },\sum _t({\eta }^{(t)}{)}^2<\mathrm{\infty }$$

Перша умова гарантує, що запасу кроків вистачить дійти куди завгодно; друга — що шум зрештою згасне. Обом задовольняє степеневе спадання ${\eta }^{(t)}=(t+1{)}^{-r}$ при $r\in (0,5;1]$. Межу зліва не можна брати: при рівно $r=0,5$ сума квадратів перетворюється на гармонічний ряд і розбігається, тому «чистий» розклад $1/\sqrt{t}$ умові не задовольняє. На практиці беруть $r$ трохи більше за 0,5.

Джерело: Dive into Deep Learning · Stochastic Gradient Descent

Міні-батч: компроміс, який лягає на відеокарту

Між одним спостереженням і всією вибіркою є природна середина. Міні-батч (mini-batch) — випадкова підмножина $B$ фіксованого розміру, градієнт усереднюють по ній:

$$\mathbf{g}=\frac{1}{|B|}\sum _{n\in B}\mathrm{\nabla }{\mathcal{L}}_n({\mathbf{w}}^{(t)}),{\mathbf{w}}^{(t+1)}:={\mathbf{w}}^{(t)}-\eta \mathbf{g}$$

Обидві крайності — окремі випадки: $|B|=1$ дає стохастичний спуск, $|B|=N$ — повний градієнт.

Головна перевага не математична, а апаратна: $|B|$ градієнтів рахуються незалежно, тому їх обчислюють одночасно. Саме так глибоке навчання використовує графічні процесори, запускаючи обчислення в $|B|$ паралельних потоках.

Виграш від збільшення батча спадний: перехід від 1 до 2 дає багато, від 100 до 101 — майже нічого. Типові розміри — від кількох десятків до кількох сотень, і обмежує їх здебільшого обсяг пам'яті пристрою.

Повний і стохастичний градієнт на California Housing

Порівняємо обидва методи на нашій задачі чесно — за обсягом обчислень, а не за кількістю ітерацій.

Ліворуч по горизонталі відкладено проходи по всій вибірці. За цією міркою методи майже нерозрізненні: після дванадцяти проходів повний градієнт дає втрату 0,2779, стохастичний — 0,2625 при точному мінімумі 0,2617.

Праворуч по горизонталі відкладено кількість оброблених окремих зразків у логарифмічній шкалі, і картина міняється. Стохастичний спуск встигає збити втрату до 0,768 за перші тисячу зразків — на цей момент повний градієнт ще не зробив жодного кроку, бо йому потрібні всі 16 512 записів навчальної частини.

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

Чотири способи знайти ваги

Метод	Коли застосовний	Вартість	Слабке місце
Перебір по сітці	будь-яка втрата, один-два параметри	${10}^D$ обчислень	не масштабується, відповідь наближена
Нормальні рівняння	лише лінійна модель із квадратичною втратою	$O(N{D}^2+D^3)$	вимагає повного стовпцевого рангу, важкий при великому $D$
Градієнтний спуск	будь-яка диференційовна втрата, будь-який $D$	$O(NDt)$ за $t$ ітерацій	потрібно дібрати $\eta$ і кількість ітерацій
Стохастичний спуск	втрата у вигляді середнього по зразках	$O(D)$ на крок	шум, потрібен розклад кроку

Правило вибору просте. Кілька десятків ознак і кілька тисяч спостережень — беріть закриту формулу. Багато ознак або багато даних — градієнтний спуск. Дані не поміщаються в пам'ять або надходять потоком — стохастичний спуск.

Реальні ландшафти втрат не опуклі

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

На рисунку градієнтний спуск працює на функції $f(x)=x\cos (cx)$, у якої мінімумів кілька. Стартувавши з $x=10$, процес зупиняється в точці $x\approx -1,53$, яка є локальним, але не глобальним мінімумом. Інша початкова точка дала б іншу відповідь.

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

Джерело: Dive into Deep Learning · Gradient Descent

Обмеження лінійної регресії

Модель бачить лише лінійні зв'язки. Криволінійну залежність вона опише тільки якщо нелінійність внести руками, через перетворені ознаки.

Квадратична втрата нестійка до викидів. Один аномальний округ здатен помітно зрушити пряму; ліки — втрата Хубера або відсів аномалій до навчання.

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

Коефіцієнти важко тлумачити при корельованих ознаках. Ваги при широті й довготі великі за модулем, $-0,90$ і $-0,87$, але через кореляцію $-0,92$ між самими ознаками жодну з них окремо не можна читати як «вплив широти на ціну».

Модель не виправляє вад даних. Стеля 5,0 залишиться стелею за будь-якого оптимізатора: якщо справжня відповідь у даних відсутня, вивчити її неможливо.

Підсумок лекції

Лінійна модель — це $\mathbf{X}\mathbf{w}$, і нічого більше. Стовпець одиниць перетворює зсув на звичайну вагу, а нелінійні ознаки дають криві, не змінюючи жодної формули.

Функція втрат — це рішення, а не даність. Квадратична зручна для обчислень, абсолютна стійка до викидів, Хубера поєднує обидві. Один викид зсунув мінімум квадратичної втрати з 2,5 на 6,0 — це й є ціна зручності.

Опуклість — головна причина, чому лінійна регресія навчається надійно. Поганих локальних мінімумів немає, тож початкова точка не важить; єдиність мінімуму дає повний стовпцевий ранг $\mathbf{X}$.

Крок навчання і масштаб ознак важать не менше за сам алгоритм. Стандартизація скоротила шлях з 305 ітерацій до 5, не змінивши в задачі нічого, крім одиниць виміру.

Стохастичний градієнт міняє точність кроку на його ціну — і саме цей обмін робить можливим навчання на великих даних.

Гарна оптимізація не рятує від поганих даних. Наша модель промахується в середньому на 53,5 тисячі доларів, і найбільша частина цієї помилки походить не від алгоритму, а від стелі 5,0.

Джерела й додаткові матеріали

Курси, на яких побудовано лекцію

• EPFL CS-433 — лекції 01c, 01d, 02a, 03a
• UBC CPSC 340 — лекції 11, 13, 22
• Stanford CS229 — розділ 1

Підручники

• Dive into Deep Learning — розділ 3
• Bishop. Pattern Recognition and ML — розділ 3
• Boyd, Vandenberghe. Convex Optimization — розділи 3, 9

Документація й дані

• scikit-learn · лінійні моделі
• numpy.linalg.solve
• California Housing — Pace, Barry, 1997