Розділ 6
Програмуємо МНК
Формулу $\hat{\boldsymbol\theta} = (A^\mathsf{T}A)^{-1}A^\mathsf{T}\mathbf y$ можна записати в один рядок коду. І саме так робити не варто. Комп'ютер рахує з обмеженою кількістю цифр, і вибір алгоритму вирішує, скільки з них доживе до відповіді.
6.1. Формула проти алгоритму
float64 приблизно з 16 значущими цифрами, у float32 (стандарт для відеокарт і нейромереж) —
лише з 7. Кожна операція трохи округлює результат. Чи всі способи обчислити ту саму формулу однаково добрі?
Ні. Математично рівні вирази можуть давати зовсім різні числа в комп'ютері. Мірило «підступності» задачі — число обумовленості $\kappa(A)$. Воно показує, у скільки разів відносна похибка вхідних даних може зрости у відповіді. Грубе правило:
Звідки береться велике $\kappa$, ми вже бачили в розділах 3–4: наприклад, коли всі $x$ далеко від нуля (роки, координати в метрах). Стовпці $\mathbf 1$ і $\mathbf x$ стають майже паралельними, а чаша — вузькою долиною.
6.2. Скільки цифр губиться
Візьмемо 10 точок на прямій з невеликим шумом, але відсунемо всі $x$ на $X_0$ від нуля: $x = X_0,\ X_0 + 0{,}37,\ \dots,\ X_0 + 3{,}33$. Порахуємо нахил трьома способами, округлюючи кожну арифметичну операцію до обраної точності. Так справді рахує процесор.
Що видно на графіку:
-
Нормальні рівняння ламаються першими. Кожен порядок $X_0$ коштує їм приблизно двох цифр.
У
float32за $X_0 = 1000$ лишається одна правильна цифра, за $10\,000$ — жодної (у таблиці «—»: знаменник формули став нулем). Уfloat64запас більший, але за $X_0 = 10^7$ теж лишаються лише дві цифри. - QR тримається набагато довше. Той самий комп'ютер, та сама точність, а відповідь точна до останньої цифри, доки $x$ взагалі вдається записати. Секрет у 6.3: для прямої перший крок QR — це, по суті, центрування.
- Центрування прибирає причину. Якщо спершу відняти середнє $x$, задача стає добре обумовленою, і навіть проста формула дає всі цифри. Найкращий алгоритм — той, що не створює поганої обумовленості.
6.3. QR: спершу зробити стовпці перпендикулярними
Підказку дав розділ 5. Проєктувати на площину легко, якщо в неї є перпендикулярні напрямні вектори одиничної довжини: тоді кожна координата проєкції — просто скалярний добуток. QR-розклад робить саме це:
$$A = Q\,R,$$де стовпці $Q$ ортонормовані (взаємно перпендикулярні, довжини 1), а $R$ — верхньотрикутна матриця $k\times k$. Підставляємо в нормальні рівняння: $R^\mathsf{T}Q^\mathsf{T}QR\,\boldsymbol\theta = R^\mathsf{T}Q^\mathsf{T}\mathbf y$. Оскільки $Q^\mathsf{T}Q = I$, лишається
$$R\,\hat{\boldsymbol\theta} = Q^\mathsf{T}\mathbf y.$$Трикутна система розв'язується знизу вгору, по одному невідомому. І головне: $A^\mathsf{T}A$ ми взагалі не утворюємо, тож обумовленість не підноситься до квадрата.
Приклад: три точки з розділу 1
Ортогоналізуємо стовпці $\mathbf 1 = (1,1,1)$ і $\mathbf x = (1,2,3)$ (процес Грама — Шмідта):
$$\mathbf q_1 = \frac{\mathbf 1}{\sqrt3},\qquad \mathbf x - (\mathbf q_1\!\cdot\mathbf x)\,\mathbf q_1 = \begin{pmatrix}1\\2\\3\end{pmatrix} - 2\begin{pmatrix}1\\1\\1\end{pmatrix} = \begin{pmatrix}-1\\0\\1\end{pmatrix},\qquad \mathbf q_2 = \frac{1}{\sqrt2}\begin{pmatrix}-1\\0\\1\end{pmatrix}.$$Другий вектор — це просто центровані $x$. Центрування, яке рятувало нас у розділах 3–4, виявляється першим кроком QR. Далі:
$$R = \begin{pmatrix}\sqrt3 & 2\sqrt3\\ 0 & \sqrt2\end{pmatrix},\qquad Q^\mathsf{T}\mathbf y = \begin{pmatrix}5/\sqrt3\\ 1/\sqrt2\end{pmatrix}.$$Зворотний хід — спершу останнє рівняння, потім перше:
Знайомі $a = \tfrac23$, $b = \tfrac12$.
6.4. SVD: найнадійніший інструмент
Сингулярний розклад (singular value decomposition) іде ще далі:
$$A = U\,\Sigma\,V^\mathsf{T},$$де $U$ і $V$ мають ортонормовані стовпці, а $\Sigma$ — діагональна матриця сингулярних чисел $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_k \ge 0$. Геометрично: будь-яка матриця — це поворот, розтягування вздовж осей і ще один поворот. Розв'язок МНК:
$$\hat{\boldsymbol\theta} = \sum_{i=1}^{k} \frac{\mathbf u_i^\mathsf{T}\mathbf y}{\sigma_i}\,\mathbf v_i,\qquad \kappa(A) = \frac{\sigma_1}{\sigma_k}.$$Формула показує, де ховається небезпека: ділення на малі $\sigma_i$. Якщо стовпці майже залежні, якесь $\sigma_i$ крихітне, і шум у $\mathbf u_i^\mathsf{T}\mathbf y$ підсилюється в $1/\sigma_i$ разів. SVD дає ліки: доданки з надто малими $\sigma_i$ можна просто відкинути.
Коли $\delta$ мале, повний розв'язок «розхитується»: за $\delta = 0{,}001$ виходить $b_1 \approx -31$, $b_2 \approx 32$, за
$\delta = 0{,}0001$ — уже $\mp 312$. Модель підлаштовує різницю двох майже однакових ознак під шум. Сума квадратів при цьому майже така сама, як у розумного розв'язку $b_1 \approx b_2 \approx 0{,}5$.
Обрізання малих сингулярних чисел обирає з усіх майже однаково добрих розв'язків той, що має найменшу
довжину $|\boldsymbol\theta|$. Саме так працює псевдообернена матриця $A^{+}$ і параметр rcond у
numpy.linalg.lstsq.
6.5. Що обрати на практиці
| Спосіб | Вартість (n ≫ k) | Втрата цифр | Коли доречний |
|---|---|---|---|
| Нормальні рівняння + розклад Холецького | $\approx n k^2$ — найдешевший | $\approx 2\log_{10}\kappa(A)$ | добре обумовлені задачі; дані надходять потоком; дуже великі розріджені системи |
| QR (відбиття Хаусхолдера) | $\approx 2 n k^2$ | $\approx \log_{10}\kappa(A)$ | стандартний вибір «за замовчуванням» |
| SVD | у кілька разів дорожче за QR | $\approx \log_{10}\kappa(A)$ | майже залежні стовпці, потрібна діагностика або розв'язок з найменшою нормою |
| Ітераційні (CG, LSQR, градієнтні) | залежить від кількості ітерацій | залежить від $\kappa$ через швидкість збіжності | гігантські розріджені задачі; нейромережі |
І завжди, незалежно від алгоритму: явно обертати $A^\mathsf{T}A$ не треба. Обернена матриця дорожча й менш точна, ніж розв'язання системи, і нічого не дає, крім красивого запису.
6.6. Геодезія: мережі на сотні тисяч невідомих
Геодезія була першою галуззю, де МНК став обчислювальною проблемою. Державні мережі мають тисячі й сотні тисяч пунктів, кожен пункт дає два-три невідомих. Спасає те, що кожен вимір пов'язує лише кілька сусідніх пунктів: матриця $A$ і $A^\mathsf{T}A$ майже повністю складаються з нулів (розріджені).
- Ще у 1880-х роках Гельмерт запропонував розбивати мережу на блоки, зрівнювати їх окремо й потім «зшивати». Цей підхід (Helmert blocking) використовували для зрівнювання континентальних мереж аж до кінця XX століття, зокрема для північноамериканського датуму NAD 83 із сотнями тисяч пунктів.
- Сьогодні такі системи розв'язують розрідженим розкладом Холецького з вдалим упорядкуванням невідомих, щоб нулі лишалися нулями. Тут нормальні рівняння виграють за швидкістю й пам'яттю, а обумовленість забезпечують центруванням координат і розумним вибором ваг.
- Та сама структура — у фотограмметрії й комп'ютерному зорі: bundle adjustment одночасно уточнює положення тисяч камер і мільйонів точок. Про це — в підсумковому розділі 11.
6.7. Що всередині ML-бібліотек
-
numpy.linalg.lstsqрозв'язує МНК через SVD (процедура LAPACKgelsd).sklearn.linear_model.LinearRegressionдля щільних даних викликаєscipy.linalg.lstsq, теж на основі SVD.statsmodelsза замовчуванням використовує псевдообернену матрицю, аlmу R і операторA\bу MATLAB — QR-розклад. Ніхто з них не обертає $X^\mathsf{T}X$. - Масштабування ознак (розділ 3) — це не лише про швидкість градієнтного спуску, а й про обумовленість: воно зменшує $\kappa$, а отже й втрату цифр у будь-якому алгоритмі.
-
Точність має значення в глибокому навчанні. Відеокарти рахують у
float32, а часто й уfloat16/bfloat16(3–4 значущі цифри). Погана обумовленість, якуfloat64пробачив би, там стає реальною проблемою. Звідси нормалізація шарів і обережність з формулами на кшталт наших «шкільних» сум. - Мультиколінеарність з розділу 5 — це малі сингулярні числа. Обрізання SVD — один спосіб з нею впоратися, гребенева регресія (ridge, розділ 10) — інший, м'якший: вона не відкидає малі $\sigma_i$, а «пригальмовує» їх.
- Розв'язок з найменшою нормою виникає і там, де його не просили: градієнтний спуск для лінійної регресії, стартуючи з нуля, сходиться саме до нього, коли розв'язків багато. Це один із прикладів неявної регуляризації, про яку багато говорять у теорії глибокого навчання.
- Великі дані: $X^\mathsf{T}X$ і $X^\mathsf{T}\mathbf y$ можна накопичувати частинами (розділ 4), якщо $k$ невелике, а дані добре обумовлені. Інакше — ітераційні методи й стохастичний градієнтний спуск.
🔍 Хочете більше подробиць? Число обумовленості, Хаусхолдер проти Грама — Шмідта, псевдообернена матриця
Число обумовленості
Для прямокутної матриці з лінійно незалежними стовпцями $\kappa(A) = \sigma_{\max}/\sigma_{\min}$. Оскільки $A^\mathsf{T}A = V\Sigma^2V^\mathsf{T}$, сингулярні числа $A^\mathsf{T}A$ — це $\sigma_i^2$, звідки $\kappa(A^\mathsf{T}A) = \kappa(A)^2$. Для прямої з $x$ поблизу великого $X_0$ число $\kappa(A)$ росте приблизно як $X_0^2$. Правило «губиться $\log_{10}\kappa$ цифр» — це оцінка для найгіршого випадку. У нашому інтерактиві нормальні рівняння на практиці губили близько двох цифр на порядок $X_0$, тобто менше за межу, але все одно катастрофічно порівняно з QR.
Точніше, для QR похибка має вигляд $\varepsilon\,\bigl(\kappa(A) + \kappa(A)^2\,\frac{|\mathbf r|}{|A||\boldsymbol\theta|}\bigr)$. Квадрат з'являється і тут, але помножений на відносну величину нев'язки. Якщо модель добре описує дані, він малий. Нормальні рівняння платять $\kappa^2$ завжди.
Як будують Q
- Класичний Грам — Шмідт (як у прикладі вище) наочний, але на погано обумовлених матрицях стовпці $Q$ поступово перестають бути перпендикулярними.
- Модифікований Грам — Шмідт віднімає проєкції по одній, і втрата ортогональності значно менша. Саме його використовує інтерактив 6.2.
- Відбиття Хаусхолдера послідовно «дзеркалять» стовпці, обнуляючи все під діагоналлю. Це найстійкіший варіант і стандарт у LAPACK. Його реалізовано в JavaScript-коді нижче.
- QR з вибором головного стовпця (column pivoting) переставляє стовпці так, щоб найменш залежні йшли першими. Діагональ $R$ тоді показує ефективний ранг, майже як SVD, але дешевше.
Псевдообернена матриця
$A^{+} = V\Sigma^{+}U^\mathsf{T}$, де $\Sigma^{+}$ містить $1/\sigma_i$ для ненульових (або не надто малих) $\sigma_i$ і нулі для решти. Розв'язок $\hat{\boldsymbol\theta} = A^{+}\mathbf y$ має дві властивості: він мінімізує $|\mathbf y - A\boldsymbol\theta|$, а серед усіх таких — ще й $|\boldsymbol\theta|$. Для матриці з незалежними стовпцями $A^{+} = (A^\mathsf{T}A)^{-1}A^\mathsf{T}$, тобто це узагальнення нашої формули на вироджені випадки.
Уточнення розв'язку
Якщо розв'язок уже є, але неточний, можна порахувати нев'язку $\mathbf r = \mathbf y - A\boldsymbol\theta$ з підвищеною точністю, розв'язати ту саму задачу для $\mathbf r$ замість $\mathbf y$ і додати поправку. Це ітераційне уточнення. Ним користуються, коли обчислення йдуть у низькій точності (зокрема на відеокартах), а відповідь потрібна у високій.
6.8. Код
Python: п'ять способів розв'язати ту саму задачу з $x \approx 10^6$, від явного оберненого до центрування, і похибка кожного.
JavaScript: власні реалізації нормальних рівнянь з розкладом Холецького та QR на відбиттях Хаусхолдера, перевірені на
$x \approx 10^7$ у звичайній для JS точності float64.
6.9. Перевірте себе
$\kappa(A) = 1000$. Чому дорівнює $\kappa(A^\mathsf{T}A)$?
Нормальні рівняння розв'язують у float64 (≈16 цифр), $\kappa(A^\mathsf{T}A) = 10^{12}$. Скільки правильних цифр лишиться?
Чому дорівнює $R_{11}$ у QR-розкладі матриці з трьома рядками $(1, x_i)$? (десятковим дробом)
Скільки ненульових сингулярних чисел має матриця зі стовпцями $\mathbf 1$, $\mathbf x$ і $2\mathbf x$ (за різних $x_i$)?
Тепер у нас є надійний «двигун»: будь-яка задача вигляду $A\boldsymbol\theta \approx \mathbf y$ розв'язується одним викликом. У розділі 7 ми побачимо, що під цей вигляд підходить набагато більше, ніж прямі: параболи й многочлени, коло, експонента, сезонні коливання. Головне — лінійність за параметрами, а не за $x$. Там же зустрінемо й головну небезпеку гнучких моделей — перенавчання.