Розділ 4

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

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

4.1. Від «чаші» до рівнянь

Проблема Ми знаємо, що $S(a, b)$ — чаша з одним дном. Кнопка «На дно» в розділі 3 якось його знаходила. Як саме обчислити координати дна точно, без перебору й без поступового спуску?

Згадайте перерізи з розділу 3: на дні обидві параболи мають вершину, тобто нахил уздовж $a$ і нахил уздовж $b$ дорівнюють нулю. Будь-де в іншому місці хоча б один нахил ненульовий, і можна спуститися нижче. Мовою похідних:

$$\frac{\partial S}{\partial a} = 0, \qquad \frac{\partial S}{\partial b} = 0.$$

Два рівняння, два невідомих. А оскільки $S$ — многочлен другого степеня, його похідні мають перший степінь. Тож рівняння будуть лінійними, а лінійні системи ми розв'язувати вміємо. У цьому й полягає головна перевага суми квадратів, яку ми обіцяли в розділі 2.

4.2. Дві умови рівноваги

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

$$\sum_{i=1}^n r_i = 0, \qquad \sum_{i=1}^n x_i\, r_i = 0, \qquad\text{де } r_i = y_i - a - b\,x_i.$$

$\sum r_i$

$\sum x_i\,r_i$

Для точок, які ви щойно рухали, ось «ручна» таблиця сум, з якої складаються нормальні рівняння. Так МНК рахували геодезисти й астрономи до появи комп'ютерів: стовпчики чисел, підсумки внизу, дві формули.

4.3. Вивід крок за кроком

Тепер власне вивід. Кожен крок показано двічі: загальною формулою і з числами для точок з інтерактиву вище. Якщо змінити точки, числа тут теж зміняться.

Система з кроку 4 і є нормальні рівняння. Для прямої їх два, бо параметрів два. Для параболи $y = a + b\,x + c\,x^2$ їх буде три, для $k$ параметрів — $k$. Щоразу це лінійна система, яку розв'язує стандартний алгоритм (метод Гаусса), навіть якщо рівнянь сотні.

4.4. Що означають ці умови

$\sum r_i = 0$: пряма проходить через центр ваги

Перше нормальне рівняння, поділене на $n$, дає $\bar y = a + b\,\bar x$. Отже, пряма МНК завжди проходить через точку $(\bar x;\,\bar y)$ — середнє всіх точок. Це прямий аналог розділу 2, де для одного невідомого оцінкою МНК було середнє арифметичне. Кнопка «Лише Σr = 0» в інтерактиві зсуває пряму до центру ваги, не змінюючи нахилу.

$\sum x_i r_i = 0$: у нев'язках не лишилося тренду

Друга умова каже, що нев'язки не «корелюють» з $x$. Якби праворуч нев'язки були переважно додатні, а ліворуч від'ємні, пряму варто було б трохи повернути. Тоді $S$ зменшилася б, і це не мінімум. Графік нев'язок у інтерактиві показує саме це: пунктир — пряма МНК через точки $(x_i, r_i)$. Для оптимальної прямої він лежить точно на нулі. Уся лінійна залежність від $x$ вичерпана, у нев'язках лишився тільки «шум».

Чому «нормальні»?

Обидві умови мають однаковий вигляд: сума нев'язок, помножених на відповідний стовпчик матриці $A$ з розділу 1 (стовпчик одиниць і стовпчик $x_i$), дорівнює нулю. Скалярний добуток, рівний нулю, означає перпендикулярність. Отже, вектор нев'язок перпендикулярний (нормальний) до стовпців матриці. Звідси найпоширеніше пояснення назви. Сам термін походить від Гаусса, і історики сперечаються, чи мав він на увазі саме перпендикулярність, чи просто «стандартну» форму рівнянь. Але геометричне тлумачення прижилося, і в розділі 5 ми його розгорнемо.

4.5. Геодезія: зрівнювання нівелірної мережі

Нормальні рівняння потрібні не лише для прямих. Ось класична геодезична задача. Від репера Rp з висотою 100,000 м виміряли перевищення до пунктів B і C, а також між ними. Невідомих два ($H_B$ і $H_C$), вимірів три. Кожен вимір дає рівняння:

$$H_B - 100 = h_1, \qquad H_C - H_B = h_2, \qquad H_C - 100 = h_3.$$

Це та сама система $A\boldsymbol\theta = \mathbf l$, лише матриця $A$ складена з коефіцієнтів $\pm1$ і $0$, а не з $1$ і $x_i$. Нормальні рівняння мають той самий вигляд: $A^\mathsf{T}A\,\boldsymbol\theta = A^\mathsf{T}\mathbf l$. Спробуйте змінити виміряні перевищення.

У реальних мережах пунктів сотні, вимірів тисячі, до них додаються ваги (розділ 8). Але схема та сама: скласти $A$, обчислити $A^\mathsf{T}A$ і $A^\mathsf{T}\mathbf l$, розв'язати систему. Метод виключення, який ми звемо методом Гаусса, Гаусс описав саме для розв'язання нормальних рівнянь у 1810 році, коли обчислював орбіту астероїда Паллада. А розклад Холецького, стандартний спосіб розв'язувати такі системи й сьогодні, придумав французький військовий топограф Андре-Луї Холецький для зрівнювання геодезичних мереж.

4.6. Normal equation у машинному навчанні

Зв'язок з ML

У курсах машинного навчання ця формула так і називається: normal equation. Її зазвичай пишуть у матричному вигляді, з матрицею ознак $X$ замість $A$:

$$\boldsymbol\theta = \bigl(X^\mathsf{T}X\bigr)^{-1} X^\mathsf{T}\mathbf y.$$

Два способи дістатися дна. Градієнтний спуск (кулька з розділу 3) іде до мінімуму кроками, доки градієнт не стане нулем. Нормальні рівняння ставлять умову «градієнт дорівнює нулю» одразу й розв'язують її за один раз. Градієнт тут — це ті самі дві суми з індикаторів: $\nabla S = -2A^\mathsf{T}\mathbf r$.
Коли що обирати. Щоб скласти $X^\mathsf{T}X$ для $n$ прикладів і $k$ ознак, потрібно порядку $n k^2$ операцій, а щоб розв'язати систему — порядку $k^3$. Для десятків чи тисяч ознак це миттєво. Для мільйонів параметрів (нейромережі) неможливо, тож там лише градієнтні методи.
Дані частинами. $X^\mathsf{T}X$ і $X^\mathsf{T}\mathbf y$ — суми внесків окремих прикладів, як наша таблиця сум. Їх можна накопичувати потоком або рахувати паралельно на багатьох машинах і потім додати.
Графік нев'язок (residual plot) — стандартна діагностика регресії. Якщо в ньому видно криву чи «віяло», модель щось пропустила. Для прямої МНК лінійного тренду там немає за побудовою, тож помітна закономірність означає, що залежність нелінійна (розділ 7).
Застереження: формулу з оберненою матрицею зручно писати, але обчислювати $(X^\mathsf{T}X)^{-1}$ явно не варто. Бібліотеки (зокрема sklearn і numpy.linalg.lstsq) роблять інакше. Чому — у розділі 6.

🔍 Хочете більше подробиць? Похідні, матричний вивід, гарантія мінімуму, чисельні пастки

Похідні детально

Позначимо $r_i = y_i - a - b\,x_i$. Тоді $\partial r_i/\partial a = -1$, $\partial r_i/\partial b = -x_i$, і за правилом диференціювання складеної функції

$$\frac{\partial S}{\partial a} = \sum 2 r_i \cdot (-1) = -2\sum r_i, \qquad \frac{\partial S}{\partial b} = \sum 2 r_i \cdot (-x_i) = -2\sum x_i r_i.$$

Підставляючи $r_i$ і розкриваючи суми, отримуємо $-2\bigl(\sum y_i - n a - b\sum x_i\bigr) = 0$ і $-2\bigl(\sum x_i y_i - a\sum x_i - b\sum x_i^2\bigr) = 0$. Це нормальні рівняння.

Матричний вивід для будь-якої кількості параметрів

Для моделі $A\boldsymbol\theta \approx \mathbf y$ з $k$ параметрами $S(\boldsymbol\theta) = \|\mathbf y - A\boldsymbol\theta\|^2 = \mathbf y^\mathsf{T}\mathbf y - 2\boldsymbol\theta^\mathsf{T}A^\mathsf{T}\mathbf y + \boldsymbol\theta^\mathsf{T}A^\mathsf{T}A\boldsymbol\theta$. Градієнт:

$$\nabla S = -2A^\mathsf{T}\mathbf y + 2A^\mathsf{T}A\boldsymbol\theta = -2A^\mathsf{T}\bigl(\mathbf y - A\boldsymbol\theta\bigr) = -2A^\mathsf{T}\mathbf r.$$

Прирівнюючи до нуля:

$$A^\mathsf{T}A\,\boldsymbol\theta = A^\mathsf{T}\mathbf y \qquad\Longleftrightarrow\qquad A^\mathsf{T}\mathbf r = \mathbf 0.$$

Для прямої $A^\mathsf{T}A$ — та сама матриця $M$ з розділу 3, а два рядки рівності $A^\mathsf{T}\mathbf r = \mathbf 0$ — наші умови $\sum r_i = 0$ і $\sum x_i r_i = 0$.

Чому знайдена точка — справді мінімум

Нульовий градієнт буває і в максимумі, і в сідлі. Але матриця других похідних (гессіан) тут дорівнює $2A^\mathsf{T}A$. У розділі 3 ми бачили, що вона додатно визначена, якщо стовпці $A$ лінійно незалежні (для прямої — якщо не всі $x_i$ однакові). Тоді функція опукла, і єдина точка з нульовим градієнтом — глобальний мінімум. Якщо ж стовпці залежні, $A^\mathsf{T}A$ вироджена, і розв'язків нескінченно багато (жолоб).

Зв'язок зі статистикою

Формулу з кроку 6 можна записати через вибіркові коваріацію й дисперсію, а також через коефіцієнт кореляції $\rho_{xy}$:

$$b = \frac{\operatorname{cov}(x, y)}{\operatorname{var}(x)} = \rho_{xy}\,\frac{s_y}{s_x}.$$

Звідси, зокрема, видно, що нахил регресії $y$ на $x$ і нахил регресії $x$ на $y$ — різні прямі (якщо $|\rho_{xy}| < 1$). Це наслідок того, що ми мінімізуємо вертикальні, а не перпендикулярні відстані.

Чисельна пастка «шкільної» формули

Формула з кроку 5 має в знаменнику різницю $n\sum x_i^2 - (\sum x_i)^2$. Якщо $x_i$ великі й близькі одне до одного (скажімо, роки 2001…2010 або координати в метрах у державній системі, порядку $5\,000\,000$), обидва доданки величезні й майже рівні. Їхня різниця втрачає більшість значущих цифр (катастрофічне скорочення). Формула з центрованими даними (крок 6) рахує те саме без цієї проблеми. Це та сама «вузька долина» з розділу 3, тільки в арифметиці. Звідси практичне правило: спершу відняти середнє, потім рахувати.

4.7. Код

Python: таблиця сум, нормальні рівняння, центровані формули, перевірка обох умов і матричний запис $A^\mathsf{T}A\boldsymbol\theta = A^\mathsf{T}\mathbf y$ (порівняно з np.polyfit). JavaScript: універсальна функція leastSquares(A, y) з методом Гаусса, застосована і до прямої, і до нівелірної мережі.

4.8. Перевірте себе

Точки: $(0;\,1)$, $(1;\,3)$, $(2;\,4)$.

Знайдіть нахил $b$ прямої МНК.

Знайдіть зсув $a$ (можна десятковим дробом).

Чому дорівнює сума нев'язок прямої МНК для будь-яких даних?

Скільки нормальних рівнянь буде для параболи $y = a + b\,x + c\,x^2$?

Далі: геометрія

Ми з'ясували, що в розв'язку вектор нев'язок перпендикулярний до стовпців матриці $A$: $A^\mathsf{T}\mathbf r = \mathbf 0$. Це підказка, що за нормальними рівняннями ховається проста геометрична картина: найкоротший перпендикуляр з точки на площину. У розділі 5 ми намалюємо її у тривимірному просторі: вектор $\mathbf y$, площина стовпців $A$ і його проєкція.