Розділ 5

Матричний погляд і геометрія

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

5.1. Третій простір: простір спостережень

Проблема Нормальні рівняння $A^\mathsf{T}A\boldsymbol\theta = A^\mathsf{T}\mathbf y$ ми отримали через похідні. Вони працюють, але чому саме така формула? І як зрозуміти МНК, коли параметрів не два, а десять, і чашу вже не намалювати?

Досі ми бачили два простори. У просторі даних (осі $x$, $y$) кожна точка — один вимір. У просторі параметрів (осі $a$, $b$) кожна точка — ціла пряма. Тепер додамо третій.

Запишемо всі виміряні $y_i$ в один стовпчик $\mathbf y = (y_1, y_2, y_3)$. Для трьох точок це вектор у тривимірному просторі, і весь набір даних стає однією точкою. Це простір спостережень. Будь-яка пряма з параметрами $a$, $b$ «передбачає» у точках $x_i$ значення $a + b\,x_i$, і це теж вектор:

$$A\boldsymbol\theta = a\begin{pmatrix}1\\1\\1\end{pmatrix} + b\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}.$$

Це «погляд на систему по стовпчиках», який ми бачили в подробицях розділу 1. Усі можливі передбачення — це всі комбінації двох векторів-стовпців, тобто площина через початок координат. Якщо $\mathbf y$ лежить у ній, система сумісна. Шум виштовхує $\mathbf y$ з площини, і задача МНК стає геометричною:

МНК одним реченням Знайти в площині стовпців $A$ точку, найближчу до $\mathbf y$. Відстань між ними — це $\sqrt{S}$, бо $|\mathbf y - A\boldsymbol\theta|^2 = \sum (y_i - a - b\,x_i)^2 = S$.

5.2. Проєкція в тривимірній сцені

Наблизити до $\mathbf y$

Погравшись, зверніть увагу на кілька речей.

Червоний відрізок — нев'язка $\mathbf r = \mathbf y - A\boldsymbol\theta$. Його довжина дорівнює $\sqrt S$, тож мінімізувати суму квадратів означає вкоротити відрізок.
Найкоротший відрізок перпендикулярний площині. Це шкільна геометрія: перпендикуляр коротший за будь-яку похилу. Натисніть «Спроєктувати», і біля основи з'явиться прямий кут.
Перпендикулярний площині означає перпендикулярний обом стовпцям. Скалярні добутки $\mathbf r\cdot\mathbf 1$ і $\mathbf r\cdot\mathbf x$ дорівнюють нулю. Але це ж $\sum r_i$ і $\sum x_i r_i$, дві умови рівноваги з розділу 4. Нормальні рівняння — це умова перпендикулярності, і похідні для неї не потрібні.
Поставте точки ліворуч на одну пряму: вектор $\mathbf y$ ляже в площину, і нев'язка зникне. Поставте всі три точки на одну вертикаль: вектор $\mathbf x$ стане паралельним вектору $\mathbf 1$, і площина сплющиться в пряму. Це той самий «жолоб» з розділу 3.

Запишемо умову перпендикулярності формулою. Нев'язка ортогональна кожному стовпцю $A$:

$$A^\mathsf{T}\bigl(\mathbf y - A\hat{\boldsymbol\theta}\bigr) = \mathbf 0 \quad\Longleftrightarrow\quad A^\mathsf{T}A\,\hat{\boldsymbol\theta} = A^\mathsf{T}\mathbf y.$$

Це нормальні рівняння, отримані в один рядок. Той самий результат, що й у розділі 4, але тепер зрозуміло, чому він такий.

5.3. Піфагор і коефіцієнт детермінації

Прямий кут дає теорему Піфагора. Вектори $\hat{\mathbf y}$ і $\mathbf r$ перпендикулярні, тому

$$|\mathbf y|^2 = |\hat{\mathbf y}|^2 + |\mathbf r|^2.$$

Для наших трьох точок $9 = 8\tfrac56 + \tfrac16$. Цікавіша версія тієї самої теореми виходить, якщо спершу відняти від усіх $y_i$ середнє $\bar y$. Оскільки модель має вільний член, вектор $\mathbf 1$ лежить у площині, і розклад зберігається:

$$\underbrace{\sum (y_i - \bar y)^2}_{\text{увесь розкид}} = \underbrace{\sum (\hat y_i - \bar y)^2}_{\text{пояснений моделлю}} + \underbrace{\sum r_i^2}_{\text{залишковий}}.$$

Частку поясненого розкиду називають коефіцієнтом детермінації:

$$R^2 = 1 - \frac{\sum r_i^2}{\sum (y_i - \bar y)^2}.$$

Геометрично $R^2 = \cos^2\varphi$, де $\varphi$ — кут між центрованим вектором вимірів і центрованим вектором передбачень. Для наших трьох точок $R^2 = 1 - \tfrac{1/6}{2/3} = 0{,}75$: пряма пояснює три чверті розкиду $y$. Для прямої $R^2$ збігається з квадратом звичайного коефіцієнта кореляції між $x$ і $y$. Кореляція — це косинус кута між центрованими векторами.

5.4. Матриця проєкції

Проєкція — лінійна операція: $\hat{\mathbf y}$ виходить з $\mathbf y$ множенням на матрицю. Підставимо розв'язок нормальних рівнянь:

$$\hat{\mathbf y} = A\hat{\boldsymbol\theta} = \underbrace{A\bigl(A^\mathsf{T}A\bigr)^{-1}A^\mathsf{T}}_{P}\,\mathbf y.$$

Матриця $P$ залежить лише від $x_i$, не від $y_i$. У статистиці її називають hat matrix, бо вона «надягає капелюшок» на $\mathbf y$. Нижче вона порахована для точок з інтерактиву вище. Посуньте їх і подивіться, як вона зміниться.

$P^2 = P$. Спроєктувати вдруге нічого не змінює: точка вже в площині.
$P$ симетрична, а її слід (сума діагоналі) дорівнює кількості параметрів, тут 2. Це розмірність площини.
Діагональ — важелі точок. $\hat y_i$ залежить від власного $y_i$ з вагою $P_{ii}$. Крайні точки мають більший важіль (для $x = 1, 2, 3$ це $\tfrac56$, $\tfrac13$, $\tfrac56$): зсув крайньої точки сильно «тягне» пряму, середньої — слабше. Відтягніть одну точку далеко вбік, і її важіль наблизиться до 1. У розділі 10 ми побачимо, чому такі точки особливо небезпечні, коли в них промах.

5.5. Будь-яка кількість вимірів і параметрів

Для трьох точок простір спостережень тривимірний, і його можна намалювати. Для ста точок він стовимірний, і намалювати його неможливо. Але нічого в міркуваннях не змінюється:

	3 точки, пряма	$n$ вимірів, $k$ параметрів
вектор вимірів $\mathbf y$	точка в $\mathbb R^3$	точка в $\mathbb R^n$
усі передбачення $A\boldsymbol\theta$	площина (2 виміри)	підпростір розмірності $k$
нев'язка $\mathbf r$	пряма, перпендикулярна площині	живе в $(n - k)$-вимірному доповненні
умова мінімуму	$\mathbf r \perp \mathbf 1,\ \mathbf r \perp \mathbf x$	$\mathbf r \perp$ кожного стовпця: $A^\mathsf{T}\mathbf r = \mathbf 0$

Число $n - k$ ми вже зустрічали: у розділі 1 це була кількість надлишкових вимірів $r$. Тепер у нього є геометричний сенс. Це розмірність простору, в якому «живуть» нев'язки, тобто кількість незалежних напрямків, у яких шум може відхилити виміри від моделі.

5.6. Геодезія: два боки однієї проєкції

Повернемося до нівелірної мережі Rp–B–C з розділу 4. Три виміряні перевищення — вектор $\mathbf h = (h_1, h_2, h_3)$ у тривимірному просторі. Безпомилкові перевищення мусять задовольняти умову $h_1 + h_2 - h_3 = 0$ (обійшли трикутник, повернулися в ту саму точку). Усі такі вектори утворюють площину.

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

Параметричний спосіб (розділ 4): невідомі висоти $H_B$, $H_C$, рівняння поправок, нормальні рівняння $A^\mathsf{T}A\boldsymbol\theta = A^\mathsf{T}\mathbf l$. Геометрично це проєкція на площину стовпців $A$.
Корелатний (умовний) спосіб: невідомих не вводимо. Записуємо лише умову $\mathbf c^\mathsf{T}\mathbf h = 0$ з $\mathbf c = (1, 1, -1)$. Вектор $\mathbf c$ перпендикулярний площині правильних вимірів, тож поправки — це проєкція $\mathbf h$ на напрямок $\mathbf c$ зі зворотним знаком:

$$w = \mathbf c^\mathsf{T}\mathbf h = h_1 + h_2 - h_3 = -6\ \text{мм},\qquad \mathbf v = -\frac{w}{\mathbf c^\mathsf{T}\mathbf c}\,\mathbf c = -\frac{-6}{3}\begin{pmatrix}1\\1\\-1\end{pmatrix} = \begin{pmatrix}+2\\+2\\-2\end{pmatrix}\ \text{мм}.$$

Ті самі поправки, що й у розділі 4, але отримані без жодних невідомих висот. Параметричний спосіб проєктує на площину (невідомих $k$), корелатний проєктує на перпендикуляр до неї (умов $n - k$). Геодезисти обирали той, у якому менше рівнянь. Коли пунктів багато, а надлишкових вимірів мало, умовний спосіб вимагав розв'язувати значно меншу систему, а в епоху ручних обчислень це вирішувало все.

5.7. Геометрія лінійної регресії

Зв'язок з ML

Стовпці матриці ознак $X$ задають, що модель здатна виразити. Будь-яке передбачення лінійної моделі — це комбінація ознак, тобто точка в підпросторі стовпців. Навчена модель дає проєкцію вектора міток на цей підпростір.
Нев'язки некорельовані з кожною ознакою (якщо є вільний член). Це не властивість даних, а наслідок побудови: $X^\mathsf{T}\mathbf r = \mathbf 0$. Тож перевіряти «чи корелюють залишки з ознакою» для навчальних даних безглуздо, відповідь завжди «ні».
$R^2$ — стандартна метрика регресії, sklearn.metrics.r2_score і метод .score() у LinearRegression. Геометрично це квадрат косинуса кута між центрованими векторами $\mathbf y$ і $\hat{\mathbf y}$.
Нова ознака розширює підпростір, тому відстань до $\mathbf y$ може лише зменшитися. Помилка на навчальних даних ніколи не зростає від додавання ознак, навіть випадкових. Звідси спокуса додавати ознаки й небезпека перенавчання (розділ 7). Тому для порівняння моделей використовують скоригований $R^2$ або перевірку на окремих даних.
Мультиколінеарність: якщо стовпці майже паралельні, площина «погано натягнута»: маленька зміна $\mathbf y$ сильно змінює $\boldsymbol\theta$, хоча $\hat{\mathbf y}$ майже не змінюється. Від цього лікують регуляризацією (ridge, розділ 10) і стійкими алгоритмами (розділ 6).
Leverage з діагоналі hat matrix — стандартний інструмент діагностики впливових точок (відстань Кука та інші міри).

🔍 Хочете більше подробиць? Доведення без похідних, властивості P, чотири підпростори, ступені свободи

Чому основа перпендикуляра — найближча точка

Нехай $\hat{\mathbf y} = A\hat{\boldsymbol\theta}$ і $\mathbf r = \mathbf y - \hat{\mathbf y} \perp$ площини. Для будь-якої іншої точки площини $A\boldsymbol\theta$ вектор $\hat{\mathbf y} - A\boldsymbol\theta$ лежить у площині, отже, перпендикулярний $\mathbf r$. За Піфагором

$$|\mathbf y - A\boldsymbol\theta|^2 = |\mathbf r|^2 + |\hat{\mathbf y} - A\boldsymbol\theta|^2 \;\ge\; |\mathbf r|^2,$$

і рівність лише при $A\boldsymbol\theta = \hat{\mathbf y}$. Це повне доведення мінімальності без жодної похідної. Воно ж пояснює формулу з розділу 3: $S(\boldsymbol\theta) = S_{\min} + (\boldsymbol\theta - \hat{\boldsymbol\theta})^\mathsf{T}A^\mathsf{T}A(\boldsymbol\theta - \hat{\boldsymbol\theta})$, бо $|A(\boldsymbol\theta - \hat{\boldsymbol\theta})|^2$ — це другий доданок. Цю тотожність показує й інтерактив («Піфагор» у таблиці).

Властивості матриці проєкції

$P^2 = P$ (ідемпотентність) і $P^\mathsf{T} = P$ (симетричність). Ці дві властивості повністю характеризують ортогональні проєктори.
$I - P$ — теж проєктор, на ортогональне доповнення: $\mathbf r = (I - P)\mathbf y$.
Власні числа $P$ дорівнюють лише 0 або 1: одиниці відповідають напрямкам у площині, нулі — перпендикулярним. Тому $\operatorname{tr} P = k$.
$0 \le P_{ii} \le 1$, і сума важелів дорівнює $k$. Середній важіль дорівнює $k/n$. Точки з важелем, значно більшим за середній, варто перевіряти окремо.

Чотири підпростори

Для матриці $A$ розміром $n\times k$ простір $\mathbb R^n$ розпадається на два перпендикулярні підпростори: простір стовпців (розмірність $k$, тут живуть $\hat{\mathbf y}$) і лівий нуль-простір $\{\mathbf z : A^\mathsf{T}\mathbf z = \mathbf 0\}$ (розмірність $n - k$, тут живуть нев'язки). Базис лівого нуль-простору — якраз коефіцієнти умовних рівнянь геодезії. Для нівелірної мережі це один вектор $\mathbf c = (1, 1, -1)$: перевірте, що $A^\mathsf{T}\mathbf c = \mathbf 0$ для $A$ з розділу 4.

Ступені свободи й оцінка шуму

Якщо кожен вимір має випадкову похибку з дисперсією $\sigma^2$, то нев'язки «бачать» лише $n - k$ вимірів шуму: решта $k$ напрямків поглинула модель. Звідси $\mathbb E\,|\mathbf r|^2 = (n - k)\,\sigma^2$, а незміщена оцінка дисперсії —

$$\hat\sigma^2 = \frac{|\mathbf r|^2}{n - k}.$$

Саме тому в геодезії середню квадратичну похибку одиниці ваги рахують діленням на кількість надлишкових вимірів, а не на $n$. Докладно — у розділі 9.

Якщо стовпці залежні

Проєкція $\hat{\mathbf y}$ єдина завжди: найближча точка площини одна, хоч як площину натягнуто. Але якщо стовпці $A$ лінійно залежні, то одну й ту саму точку площини дають різні $\boldsymbol\theta$, і параметри визначені неоднозначно. Передбачення модель зробить, а параметрам вірити не можна. Серед усіх розв'язків зазвичай обирають той, що має найменшу довжину; його дає псевдообернена матриця (розділ 6).

5.8. Код

Python: проєкція, перевірка $A^\mathsf{T}\mathbf r = \mathbf 0$ і теореми Піфагора, матриця $P$ з її властивостями, важелі, $R^2$. JavaScript: проєкція на площину двох векторів «вручну» і нівелірна мережа, розв'язана корелатним способом.

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

Для точок $(1;1)$, $(2;2)$, $(3;2)$: $|\mathbf y|^2 = 9$, $|\hat{\mathbf y}|^2 = 8\tfrac56$. Чому дорівнює $|\mathbf r|^2$?

Пряму проводять через 10 точок. Яка розмірність підпростору, в якому лежать нев'язки?

Чому дорівнює слід матриці $P$ для параболи $y = a + b\,x + c\,x^2$ (за будь-якою кількістю точок, більшою за 3)?

У мережі Rp–B–C нев'язка $h_1 + h_2 - h_3 = +9$ мм. Яка поправка до $h_1$ (у мм)?

Далі: як це рахує комп'ютер

Формулу $\hat{\boldsymbol\theta} = (A^\mathsf{T}A)^{-1}A^\mathsf{T}\mathbf y$ легко написати, але в програмах її майже ніколи не обчислюють буквально. Геометрія підказує кращий шлях: якщо спершу зробити стовпці $A$ взаємно перпендикулярними, проєкція рахується тривіально. Це QR-розклад. А ще загальніший інструмент — SVD. У розділі 6 ми запрограмуємо МНК «як слід» і побачимо, чому $A^\mathsf{T}A$ може підвести.