Розділ 11

Підсумок

Десять розділів в одному: шпаргалка, яку можна тримати поруч під час роботи, інтерактивний вибір методу, наскрізні задачі з перевіркою і два напрями, куди МНК веде далі, — фільтр Калмана і bundle adjustment.

11.1. Шпаргалка: МНК за сім кроків

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

1Модель

Кожен вимір записати як рівняння: $y_i \approx a_{i1}\theta_1 + \dots + a_{ik}\theta_k$. Модель має бути лінійною за параметрами, не обов'язково за $x$ (розділ 7). Нелінійну лінеаризуємо (розділ 8).
2Матриці $A$, $\mathbf y$, $P$

Рядок $A$ — один вимір, стовпець — один параметр. Ваги $p_i = \sigma_0^2/\sigma_i^2$ (розділ 8).
3Надлишковість

$r = n - k > 0$ і стовпці $A$ незалежні. Без надлишкових вимірів немає ні вибору, ні контролю (розділ 1).
4Розв'язок

$\hat{\boldsymbol\theta}$ з нормальних рівнянь $A^\mathsf{T}PA\,\boldsymbol\theta = A^\mathsf{T}P\mathbf y$ (розділи 4–5), а в програмі — через QR або SVD, np.linalg.lstsq (розділ 6).
5Поправки й $\sigma_0$

$\mathbf v = A\hat{\boldsymbol\theta} - \mathbf y$, $\ \hat\sigma_0 = \sqrt{\mathbf v^\mathsf{T}P\mathbf v / r}$. Це оцінка точності одного виміру з одиничною вагою, отримана з самих даних (розділ 9).
6Точність параметрів

$\mathrm{Cov}(\hat{\boldsymbol\theta}) = \hat\sigma_0^2\,Q$, $\ Q = (A^\mathsf{T}PA)^{-1}$; $\ \sigma_{\theta_j} = \hat\sigma_0\sqrt{Q_{jj}}$. Еліпс похибок — з власних чисел блоку $2\times 2$ (розділ 9).
7Діагностика

Тест Баарди $w_i = v_i / (\sigma\sqrt{q_{v_iv_i}})$, $|w_i| > 3{,}29$ — підозра на промах (розділи 9–10). $R^2$ і перевірка на нових даних — чи не перенавчилися (розділи 5, 7).

Головні формули

Що	Формула	Розділ
критерій	$S(\boldsymbol\theta) = \sum p_i r_i^2 = (\mathbf y - A\boldsymbol\theta)^\mathsf{T}P(\mathbf y - A\boldsymbol\theta)$	2
пряма $y = a + bx$	$b = \dfrac{\sum (x_i - \bar x)(y_i - \bar y)}{\sum (x_i - \bar x)^2},\ a = \bar y - b\bar x$	4
нормальні рівняння	$A^\mathsf{T}PA\,\hat{\boldsymbol\theta} = A^\mathsf{T}P\mathbf y$	4, 5, 8
проєкція	$\hat{\mathbf y} = H\mathbf y,\ H = A(A^\mathsf{T}A)^{-1}A^\mathsf{T},\ A^\mathsf{T}\mathbf v = 0$	5
частка поясненого	$R^2 = 1 - \mathrm{SS_{res}}/\mathrm{SS_{tot}}$	5
обумовленість	$\kappa(A^\mathsf{T}A) = \kappa(A)^2$	6
Гаусс — Ньютон	$J^\mathsf{T}PJ\,\Delta\boldsymbol\theta = J^\mathsf{T}P\,(\mathbf y - f(\boldsymbol\theta))$	8
середнє $n$ вимірів	$\sigma_{\bar y} = \sigma/\sqrt n$	9
точність	$\hat\sigma_0^2 = \mathbf v^\mathsf{T}P\mathbf v/(n-k),\ \mathrm{Cov} = \hat\sigma_0^2(A^\mathsf{T}PA)^{-1}$	9
градієнтний спуск	$\boldsymbol\theta \leftarrow \boldsymbol\theta - \eta\nabla\mathrm{MSE},\ \eta < 2/L$	10
ridge	$(A^\mathsf{T}A + \lambda I)\,\hat{\boldsymbol\theta} = A^\mathsf{T}\mathbf y$	10
Губер (IRLS)	$p_i = 1$ при $\|r_i\| \le \delta$, інакше $\delta/\|r_i\|$	10

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

11.2. Виміри надходять потоком: фільтр Калмана

Проблема GNSS-приймач на опорі мосту видає нову висоту щосекунди. Перераховувати МНК по всіх вимірах з початку спостережень щоразу марнотратно. І гірше: якщо опора раптом просіла, середнє за всі дні ще довго «пам'ятатиме» стару висоту.

Перша проблема розв'язується просто: суми $A^\mathsf{T}A$ і $A^\mathsf{T}\mathbf y$ можна накопичувати, додаючи по одному доданку на кожен вимір. Для однієї величини це відоме з розділу 9 середнє, що оновлюється:

$$\hat x_t = \hat x_{t-1} + K_t\,(z_t - \hat x_{t-1}),\qquad K_t = \frac1t .$$

Нова оцінка — стара плюс частка «здивування» $z_t - \hat x_{t-1}$. Це рекурсивний МНК: результат точно той самий, що й пакетний, але кожен крок коштує кілька множень. Проте $K_t = 1/t$ прямує до нуля: після сотні вимірів новий вимір майже нічого не змінює. Для незмінної величини це правильно, а для опори, що просідає, — ні.

Фільтр Калмана (1960) додає одну ідею: між вимірами величина могла змінитися, тому невизначеність оцінки перед кожним новим виміром трохи зростає на $q$, дисперсію «шуму процесу»:

$$P \leftarrow P + q,\qquad K = \frac{P}{P + R},\qquad \hat x \leftarrow \hat x + K\,(z - \hat x),\qquad P \leftarrow (1 - K)\,P .$$

Тут $R$ — дисперсія шуму виміру, $P$ — дисперсія оцінки. Коефіцієнт $K$ — це вага нового виміру порівняно з усім, що ми вже знаємо, тобто те саме зважування з розділу 8. При $q = 0$ формули дають рівно рекурсивний МНК.

$q = 0$ — чистий МНК: на спокійній ділянці найточніше, але після стрибка оцінка повзе до нового значення десятками кроків.
Великий $q$ — фільтр майже не усереднює: швидко реагує, але тремтить разом із шумом.
Компроміс ($q \approx 0{,}01$–$0{,}1$) — та сама боротьба зсуву й розкиду, що в розділі 7 між простою й гнучкою моделлю. Для сценарію «Стала» найкраще $q = 0$: якщо величина справді не змінюється, МНК неперевершений.

У реальних задачах стан — вектор (координати, швидкості, орієнтація), а прогноз задає рівняння руху. Саме так працюють навігація в смартфоні, автопілоти й GNSS-приймачі, а одним з перших застосувань фільтра була навігація «Аполлона».

11.3. Мільйони невідомих: bundle adjustment і SLAM

Інший напрям — масштаб. У фотограмметрії та комп'ютерному зорі за тисячами фотографій відновлюють 3D-модель місцевості чи будівлі. Невідомі — положення й орієнтація кожної камери та координати мільйонів точок. Виміри — пікселі, де кожна точка видна на кожному знімку. Модель проєкції нелінійна, тож це МНК з Гауссом — Ньютоном з розділу 8, лише величезний. Його називають bundle adjustment, зрівнюванням «в'язок» променів.

Роботи й дрони розв'язують ту саму задачу на ходу (SLAM, одночасна локалізація й побудова карти), а геодезисти так само зрівнюють загальнодержавні мережі. Усе це можливо, бо матриця $A$ розріджена: кожен вимір стосується лише однієї камери й однієї точки, тож майже всі елементи рядка нульові. Спеціальні розв'язувачі (розріджений Холецький, доповнення Шура, метод спряжених градієнтів) використовують цю структуру й не зберігають нулі.

Куди МНК веде в машинному навчанні

Логістична регресія (класифікація «так/ні») навчається методом Ньютона, який на кожному кроці розв'язує зважений МНК (IRLS з розділу 10). Так само влаштовані всі узагальнені лінійні моделі (GLM).
Гауссові процеси й ядрові методи — ridge-регресія в нескінченновимірному просторі ознак; нормальні рівняння в них розв'язуються для матриці «схожостей» між прикладами.
Нейромережа з MSE-втратою — нелінійний МНК з мільйонами параметрів. Навчають її SGD і Adam, а не Гауссом — Ньютоном, бо матриця $J^\mathsf{T}J$ завелика. Але останній лінійний шар часто можна дорахувати звичайним МНК, а ідеї обумовленості, нормалізації ознак, регуляризації й перевірки на нових даних перейшли без змін.
Фільтр Калмана — окремий випадок байєсівського оновлення: апріорні знання плюс новий вимір дають апостеріорні. Це основа онлайн-навчання та ймовірнісного ML.

11.4. Наскрізні задачі

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

1. Середнє та його точність. Довжину лінії виміряли п'ять разів: 25,31; 25,35; 25,28; 25,34; 25,32 м. Найкраща оцінка — середнє 25,32 м (розділи 1, 4). Яка його стандартна похибка в міліметрах?

Підказка

Відхилення від середнього: −1, 3, −4, 2, 0 см. $s = \sqrt{\sum v_i^2 / (n-1)}$, бо одну невідому вже оцінено; далі $\sigma_{\bar y} = s/\sqrt n$ (розділ 9).

2. Пряма. Точки $(0;\,1{,}1)$, $(1;\,2{,}9)$, $(2;\,5{,}2)$, $(3;\,6{,}8)$, $(4;\,9{,}0)$. Знайдіть нахил $b$ прямої $y = a + bx$ за МНК.

Підказка

$\bar x = 2$, $\bar y = 5$. $b = \sum (x_i - \bar x)(y_i - \bar y) / \sum (x_i - \bar x)^2$ (розділ 4), знаменник дорівнює 10.

3. Ваги. Один і той самий зсув виміряли двома приладами: 100,0 мм ($\sigma = 1$ мм) і 103,0 мм ($\sigma = 2$ мм). Яка оцінка за зваженим МНК (мм)?

Підказка

Ваги $p_i = 1/\sigma_i^2$, тобто 1 і 0,25. Зважене середнє $\sum p_i y_i / \sum p_i$ (розділ 8).

4. Нівелірний хід. Замкнений хід з трьох секцій довжиною 2, 1 і 1 км. Сума виміряних перевищень дорівнює −12 мм, хоча має бути 0. Яку поправку (мм) отримає перша секція, якщо дисперсія пропорційна довжині?

Підказка

Вага секції $1/L_i$. МНК розподіляє нев'язку з протилежним знаком пропорційно $L_i$: $v_i = -f\,L_i / \sum L$ (розділ 8).

5. Якість моделі. Сума квадратів відхилень $y$ від середнього — 40, сума квадратів нев'язок після МНК — 4. Чому дорівнює $R^2$?

Підказка

$R^2 = 1 - \mathrm{SS_{res}}/\mathrm{SS_{tot}}$ (розділ 5).

6. Точність параметра. Після зрівнювання $\hat\sigma_0 = 2$ мм, а діагональний елемент матриці $Q = (A^\mathsf{T}PA)^{-1}$ для висоти репера $Q_{jj} = 0{,}25$. Яка стандартна похибка висоти (мм)?

Підказка

$\sigma_{\theta_j} = \hat\sigma_0\sqrt{Q_{jj}}$ (розділ 9).

7. Калман. Поточна оцінка висоти 10 см з дисперсією $P = 1$. Надходить вимір 12 см з дисперсією $R = 3$. Яка нова оцінка (см)? Вважайте $q = 0$.

Підказка

$K = P/(P + R)$, $\hat x \leftarrow \hat x + K(z - \hat x)$ (розділ 11.2). Те саме дає зважене середнє з вагами $1/P$ і $1/R$.

11.5. Що читати і чим рахувати

Книги

G. Strang, «Introduction to Linear Algebra» — проєкції, $A^\mathsf{T}A$, SVD з геометричною інтуїцією (розділи 5–6).
S. Boyd, L. Vandenberghe, «Introduction to Applied Linear Algebra: Vectors, Matrices, and Least Squares» — вільно доступна, багато прикладів з даних.
C. Ghilani, «Adjustment Computations: Spatial Data Analysis» — зрівнювання геодезичних мереж від простого до складного (розділ 8).
P. Teunissen, «Adjustment Theory» і «Testing Theory» — теорія зрівнювання й тестування в традиції Баарди (розділи 9–10).
G. Golub, C. Van Loan, «Matrix Computations»; Å. Björck, «Numerical Methods for Least Squares Problems» — чисельні методи всерйоз (розділ 6).
T. Hastie, R. Tibshirani, J. Friedman, «The Elements of Statistical Learning» — МНК, ridge, LASSO і далі в ML (розділ 10).
D. Simon, «Optimal State Estimation» — фільтр Калмана та його нелінійні варіанти (розділ 11.2).

Інструменти

Python: numpy.linalg.lstsq, scipy.linalg.lstsq, scipy.optimize.least_squares (нелінійний МНК, стійкі втрати), scipy.sparse.linalg.lsqr (великі розріджені задачі).
Статистика: statsmodels (OLS, WLS з повним звітом: похибки, довірчі інтервали, $R^2$).
ML: scikit-learn (LinearRegression, Ridge, HuberRegressor, RANSACRegressor).
Великі нелінійні задачі: Ceres Solver, g2o, GTSAM (bundle adjustment, SLAM).
JavaScript: інтерактиви цього підручника обходяться без математичних бібліотек, МНК там — Хаусхолдерів QR на кількох десятках рядків.

🔍 Хочете більше подробиць? Калман як рекурсивний МНК, формула Шермана — Моррісона, згладжування, нелінійні фільтри

Рекурсивний МНК у матричному вигляді

Нехай уже оброблено виміри з матрицею $A$ і відомо $\hat{\boldsymbol\theta}$ та $P = (A^\mathsf{T}A)^{-1}$. Надходить ще один вимір $z = \mathbf a^\mathsf{T}\boldsymbol\theta + \varepsilon$. Нова матриця нормальних рівнянь $A^\mathsf{T}A + \mathbf a\mathbf a^\mathsf{T}$ відрізняється на матрицю рангу 1, і її обернену можна отримати без обернення (формула Шермана — Моррісона):

$$\mathbf K = \frac{P\mathbf a}{R + \mathbf a^\mathsf{T}P\mathbf a},\qquad \hat{\boldsymbol\theta} \leftarrow \hat{\boldsymbol\theta} + \mathbf K\,(z - \mathbf a^\mathsf{T}\hat{\boldsymbol\theta}),\qquad P \leftarrow (I - \mathbf K\mathbf a^\mathsf{T})\,P .$$

Кожен крок коштує $O(k^2)$ операцій замість $O(k^3)$, а результат збігається з пакетним МНК. Це і є фільтр Калмана без прогнозу. Повний фільтр додає крок прогнозу $\boldsymbol\theta \leftarrow F\boldsymbol\theta$, $P \leftarrow FPF^\mathsf{T} + Q$, де $F$ описує, як стан змінюється між вимірами (наприклад, координата + швидкість × час).

Інформаційна форма

Замість $P$ можна зберігати обернену матрицю $P^{-1}$, «інформацію», тобто просто накопичені $A^\mathsf{T}PA$ і $A^\mathsf{T}P\mathbf y$. Оновлення тоді зводиться до додавання, як у коді цього розділу. Обидві форми еквівалентні, а вибір між ними залежить від того, чого більше: вимірів чи параметрів.

Фільтр і згладжувач

Фільтр дає найкращу оцінку на момент $t$ за вимірами до $t$ включно. Якщо обробка не термінова, можна використати і пізніші виміри: згладжувач Рауха — Тунга — Штрібеля проходить дані назад і отримує рівно той самий результат, що й пакетний МНК по всіх вимірах одразу з урахуванням моделі руху. Тобто фільтр Калмана — це МНК, розкладений у часі.

Нелінійні моделі

Якщо рівняння вимірів чи руху нелінійні, їх лінеаризують у поточній оцінці (розширений фільтр Калмана, EKF), тобто роблять один крок Гаусса — Ньютона на кожен вимір. Ітераційний EKF робить кілька кроків, а сучасні методи факторних графів (GTSAM) перезрівнюють усе вікно останніх станів, стираючи межу між фільтром і bundle adjustment.

11.6. Код

Увесь конвеєр з шпаргалки в одній функції adjust: розв'язок, поправки, $\hat\sigma_0$, точність параметрів і тест Баарди. Потім рекурсивний МНК, що дає той самий результат, і одновимірний фільтр Калмана з різними $q$.

Кінець підручника

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