Метод наименьших квадратов

Заметка о ближайших соседях закончилась моделью, которая отвечает на каждый запрос по горстке близких точек и ничего не узнаёт о зависимости в целом; в пространствах большой размерности близкие точки исчезают. Эта заметка делает противоположную ставку. Она предполагает, что целевая переменная с точностью до шума равна взвешенной сумме признаков, по одному весу на признак, и подбирает веса сразу по всем данным, минимизируя квадраты ошибок. У решения есть замкнутая формула, геометрический смысл и вероятностный смысл, и каждый из трёх объясняет свой способ сломаться: на коллинеарных признаках, на плохо масштабированных числах и на выбросах. Эти поломки задают остальную часть раздела.

Обозначения. Обучающая выборка — $(x_1, y_1), \dots, (x_n, y_n)$, где $x_i \in \mathbb R^d$ и $y_i \in \mathbb R$. Строки матрицы признаков (design matrix) $X \in \mathbb R^{n \times d}$ — векторы $x_i^\top$, а $y \in \mathbb R^n$ собирает целевые значения. Если в модели есть свободный член, первый признак — константа 1, так что $d$ учитывает и свободный член. Вектор весов $w \in \mathbb R^d$ даёт прогноз $\hat y_i = \langle w, x_i\rangle$, для всех объектов сразу — $Xw$. Остатки (residuals) — $r = y - Xw$, сумма квадратов остатков (residual sum of squares) — $\mathrm{RSS}(w) = \lVert y - Xw\rVert^2$, среднеквадратичная ошибка — $\mathrm{MSE} = \mathrm{RSS}/n$. $\hat w$ — точка минимума, $\hat y = X\hat w$.

Линейная модель#

Модель имеет вид

$$a(x) = \langle w, x\rangle = w_1 x^{(1)} + \dots + w_d x^{(d)} .$$

Каждый вес говорит, насколько меняется прогноз, когда один признак растёт на единицу, а остальные не меняются. Модель линейна по $w$, но не обязательно по исходным входам: признаками могут быть степени числа, логарифмы, произведения двух входов или one-hot-столбцы категории из заметки о признаках. Всё, что сказано ниже, верно для любого такого фиксированного набора признаков, потому что использует только матрицу $X$.

Функция потерь — квадрат ошибки, усреднённый по обучающей выборке:

$$Q(w) = \frac1n \sum_{i=1}^n \big(\langle w, x_i\rangle - y_i\big)^2 = \frac1n \lVert Xw - y\rVert^2 .$$

Почему квадраты, а не модули или другой штраф? В этой заметке будет три ответа: у точки минимума есть формула, она является ортогональной проекцией, и она совпадает с оценкой максимального правдоподобия при гауссовском шуме. Последний раздел показывает цену этого выбора.

Нормальные уравнения#

Раскроем сумму квадратов остатков: $\lVert y - Xw\rVert^2 = y^\top y - 2 w^\top X^\top y + w^\top X^\top X w$. Градиент линейного члена равен $-2X^\top y$, а градиент квадратичной формы при симметричной $X^\top X$ равен $2X^\top X w$, поэтому

$$\nabla_w \lVert y - Xw\rVert^2 = 2X^\top (Xw - y) .$$

Гессиан равен $2X^\top X$, и $v^\top X^\top X v = \lVert Xv\rVert^2 \ge 0$ для любого $v$: функция выпукла, и каждая точка, где градиент обращается в ноль, — глобальный минимум. Приравняв градиент к нулю, получаем нормальные уравнения (normal equations)

$$X^\top X\, \hat w = X^\top y .$$

Если столбцы $X$ линейно независимы (для этого нужно $n \ge d$), то $\lVert Xv\rVert^2 > 0$ для любого $v \ne 0$, матрица $X^\top X$ положительно определена и обратима, а точка минимума единственна:

$$\hat w = (X^\top X)^{-1} X^\top y .$$

Если столбцы зависимы, например когда один признак равен сумме двух других или когда $d > n$, у нормальных уравнений бесконечно много решений, и все они дают одни и те же прогнозы; к этому случаю вернётся раздел об устойчивости.

Разобранный пример. Четыре точки: $x = (0, 1, 2, 3)$ и $y = (1, 3, 2, 5)$, модель $\hat y = w_1 + w_2 x$. Столбцы матрицы признаков — $(1, 1, 1, 1)$ и $(0, 1, 2, 3)$, поэтому

$$X^\top X = \begin{pmatrix} 4 & 6 \\ 6 & 14 \end{pmatrix}, \qquad X^\top y = \begin{pmatrix} 11 \\ 22 \end{pmatrix} .$$

Определитель равен $4 \cdot 14 - 6 \cdot 6 = 20$, и $\hat w = \frac{1}{20}\begin{pmatrix} 14 & -6 \\ -6 & 4\end{pmatrix}\begin{pmatrix} 11 \\ 22\end{pmatrix} = (1.1, 1.1)$. Подобранная прямая — $\hat y = 1.1 + 1.1x$, остатки — $(-0.1, 0.8, -1.3, 0.6)$, и $\mathrm{RSS} = 0.01 + 0.64 + 1.69 + 0.36 = 2.7$.

Решение нормальных уравнений стоит $O(nd^2)$ на построение $X^\top X$ и $O(d^3)$ на решение системы размера $d \times d$. Для десяти признаков это ничто; для миллиона разреженных текстовых признаков это невозможно, и с этого начинается следующая заметка.

Геометрия: проекция#

Нормальные уравнения можно переписать как $X^\top (y - X\hat w) = 0$: вектор остатков ортогонален каждому столбцу $X$. Будем считать $n$ целевых значений одной точкой $y \in \mathbb R^n$. Когда $w$ меняется, прогнозы $Xw$ пробегают пространство столбцов (column space) $\operatorname{col}(X)$ — $d$-мерное подпространство $\mathbb R^n$. Метод наименьших квадратов выбирает ближайшую к $y$ точку этого подпространства, а ближайшая точка — основание перпендикуляра.

Рисунок

Метод наименьших квадратов как проекция

Ортогональность даёт доказательство оптимальности без дифференцирования. Для любого $w$ разложим $y - Xw = (y - X\hat w) + X(\hat w - w)$. Первое слагаемое ортогонально пространству столбцов, а второе лежит в нём, поэтому по теореме Пифагора

$$\lVert y - Xw\rVert^2 = \lVert y - X\hat w\rVert^2 + \lVert X(\hat w - w)\rVert^2 \ \ge\ \lVert y - X\hat w\rVert^2 ,$$

причём равенство достигается ровно тогда, когда $Xw = X\hat w$. В примере с тремя точками на рисунке добавочное слагаемое для $w = (1, 1)$ равно $0.5$.

Отображение из $y$ в $\hat y$ линейно: $\hat y = Hy$, где $H = X(X^\top X)^{-1}X^\top$ — матрица-шляпа (hat matrix). Это ортогональный проектор на $\operatorname{col}(X)$, поэтому он симметричен и идемпотентен ($H^2 = H$: повторная проекция ничего не меняет), а его след равен $\operatorname{tr}\big((X^\top X)^{-1}X^\top X\big) = d$, размерности подпространства. Диагональный элемент $h_{ii}$, рычаг (leverage) точки $i$, показывает, насколько сдвигается $\hat y_i$, когда $y_i$ сдвигается на единицу. В примере с четырьмя точками рычаги равны $(0.7, 0.3, 0.3, 0.7)$ и в сумме дают $d = 2$: крайние точки тянут прямую сильнее средних.

Два следствия ортогональности используются постоянно. Если в модели есть свободный член, константный столбец — один из столбцов, поэтому $\sum_i r_i = 0$: остатки в среднем равны нулю, и среднее $\hat y$ равно среднему $y$. В разобранном примере $-0.1 + 0.8 - 1.3 + 0.6 = 0$, а также $0 \cdot (-0.1) + 1 \cdot 0.8 + 2 \cdot (-1.3) + 3 \cdot 0.6 = 0$. Кроме того, $y = \hat y + r$ и $\hat y \perp r$, поэтому общий разброс раскладывается на объяснённую и необъяснённую части — это и измеряет $R^2$ ниже.

Вероятность: гауссовский шум и условное среднее#

Пусть данные порождены как $y_i = \langle w, x_i\rangle + \varepsilon_i$, где $\varepsilon_i$ независимы и имеют распределение $N(0, \sigma^2)$. Отрицательный логарифм правдоподобия равен $\mathrm{RSS}(w)/(2\sigma^2)$ плюс константа, как выведено в заметке о методе максимального правдоподобия, поэтому минимизировать квадратичную ошибку — значит максимизировать правдоподобие. И без всякой модели у квадратичной ошибки есть другое обоснование: её минимум по всем функциям на генеральной совокупности — условное среднее $E[Y \mid X = x]$ (наилучший прогноз). Метод наименьших квадратов, ограниченный линейными функциями, оценивает линейную функцию, ближайшую к $E[Y \mid X]$ в среднеквадратичном смысле, — то есть саму $E[Y \mid X]$, если она линейна.

Распределение коэффициентов. Здесь матрица признаков $X$ считается фиксированной, а случаен только шум: утверждения ниже — о повторных реализациях $\varepsilon$ при тех же входах. Подставим $y = Xw + \varepsilon$ в формулу:

$$\hat w = (X^\top X)^{-1}X^\top (Xw + \varepsilon) = w + (X^\top X)^{-1}X^\top \varepsilon .$$

Оценка — это истинный вес плюс фиксированное линейное отображение шума. Отсюда следуют три факта, каждый при своём предположении:

• Несмещённость, если $E[\varepsilon] = 0$: $E[\hat w] = w$.
• Ковариационная матрица $\sigma^2 (X^\top X)^{-1}$, если шум некоррелирован и имеет одинаковые дисперсии, $\operatorname{Var}\varepsilon = \sigma^2 I$: при $A = (X^\top X)^{-1}X^\top$ получаем $\operatorname{Var}\hat w = A\,\sigma^2 I\,A^\top = \sigma^2 (X^\top X)^{-1}$.
• Нормальность, если шум гауссовский, поскольку аффинное отображение гауссовского вектора гауссовское: $\hat w \sim N\big(w, \sigma^2 (X^\top X)^{-1}\big)$. Для другого независимого шума с конечной дисперсией распределение при больших $n$ приближённо нормально по центральной предельной теореме, если ни одна точка не доминирует в подгонке: наибольший рычаг $\max_i h_{ii}$ (определён выше) должен стремиться к 0.

Дисперсия шума неизвестна и оценивается по остаткам. Поскольку $r = (I - H)y = (I - H)\varepsilon$ (так как $(I - H)X = 0$), получаем $\mathrm{RSS} = \varepsilon^\top (I - H)\varepsilon$ и $E[\mathrm{RSS}] = \sigma^2 \operatorname{tr}(I - H) = \sigma^2 (n - d)$. Поэтому оценка

$$\hat\sigma^2 = \frac{\mathrm{RSS}}{n - d}$$

несмещённая: $d$ степеней свободы тратятся на подгонку, так же как $n - 1$ в выборочной дисперсии отражает одну степень свободы, потраченную на среднее. Квадратные корни из диагонали $\hat\sigma^2 (X^\top X)^{-1}$ — стандартные ошибки (standard errors) коэффициентов.

Разобранный пример. В примере с четырьмя точками $\hat\sigma^2 = 2.7/(4 - 2) = 1.35$ и $(X^\top X)^{-1} = \frac{1}{20}\begin{pmatrix} 14 & -6 \\ -6 & 4\end{pmatrix}$, так что у наклона $1.1$ стандартная ошибка $\sqrt{1.35 \cdot 4/20} = 0.520$, а у свободного члена — $\sqrt{1.35 \cdot 14/20} = 0.972$. При четырёх точках наклон отстоит от нуля едва на две стандартные ошибки. Моделирование подтверждает формулу: 30 равноотстоящих точек на $[0, 1]$, истинный $w = (1, 2)$, $\sigma = 0.5$ и 20 000 реализаций шума дают средние коэффициентов $(1.001, 2.000)$ и стандартные отклонения $(0.177, 0.303)$ против $(0.178, 0.306)$ по формуле $\sigma\sqrt{\operatorname{diag}(X^\top X)^{-1}}$ (раздел «Код»).

Теорема Гаусса–Маркова (Gauss–Markov theorem) добавляет, что при некоррелированном шуме с одинаковыми дисперсиями никакая другая линейная по $y$ несмещённая оценка не имеет меньшей дисперсии ни для какой линейной комбинации коэффициентов (Hastie, Tibshirani and Friedman, 2009, раздел 3.2.2). Оговорка важна: оценки, допускающие небольшое смещение, могут иметь намного меньшую дисперсию и меньшую среднеквадратичную ошибку — это компромисс из заметки об оценках. Такие оценки строит заметка о регуляризации.

Когда решение неустойчиво#

Ковариация $\sigma^2 (X^\top X)^{-1}$ велика в тех направлениях, где у $X^\top X$ малые собственные числа, то есть там, где некоторая комбинация признаков почти постоянна на данных. Именно это дают два признака, которые почти повторяют друг друга.

Коллинеарность (collinearity). Возьмём 100 точек с $x^{(2)} = x^{(1)} + 0.01\,\xi$, где $\xi$ — стандартный нормальный шум, и $y = x^{(1)} + x^{(2)} + \varepsilon$ с $\sigma = 1$. Корреляция двух признаков равна $0.99994$, а число обусловленности $X^\top X$ — $3.5 \cdot 10^4$. Три реализации шума дают веса $(-5.00, 6.88)$, $(11.88, -9.90)$ и $(-8.83, 10.94)$ с суммами $1.88$, $1.98$ и $2.11$. Данные фиксируют сумму $w_1 + w_2$ — направление, в котором признаки меняются, — и почти ничего не говорят о разности — направлении, в котором они не меняются. Стандартная ошибка каждого веса равна $9.74$, против $0.10$ у модели только с одним из двух признаков. Прогнозы на данных, похожих на обучающие, от реализации к реализации почти не меняются; коэффициенты, а с ними и любая их интерпретация, — шум. Типичный симптом — вес $-9.90$ у признака, который помогает предсказывать целевую переменную.

Неполный ранг (rank deficiency). Когда столбцы точно зависимы или $d > n$, матрица $X^\top X$ вырождена, и у нормальных уравнений целое аффинное множество решений. Стандартный выбор — решение с наименьшей нормой, $\hat w = X^+ y$ с псевдообратной матрицей Мура–Пенроуза $X^+$; его и возвращает numpy.linalg.lstsq. Этот выбор — предпочтение, а не то, что решают данные, и заметка о регуляризации делает предпочтение явным.

Численная обусловленность. Даже когда $X^\top X$ обратима в точной арифметике, её построение возводит число обусловленности в квадрат: $\kappa(X^\top X) = \kappa(X)^2$. Двойная точность хранит около 16 значащих цифр, а решение линейной системы теряет примерно $\log_{10}\kappa$ из них. Для полиномиальных признаков $1, t, \dots, t^{10}$ на 50 точках отрезка $[0, 1]$ получаем $\kappa(X) = 2.0 \cdot 10^7$ и $\kappa(X^\top X) = 4.1 \cdot 10^{14}$. Решение нормальных уравнений восстанавливает коэффициенты, которые все равны 1, с ошибкой $4 \cdot 10^{-3}$; решатель на основе SVD, применённый прямо к $X$, как это делает lstsq, — с ошибкой $4 \cdot 10^{-10}$ (Trefethen and Bau, 1997, лекции 11 и 19). Формула $(X^\top X)^{-1}X^\top y$ — для выводов; численный код раскладывает $X$. Центрирование и масштабирование признаков, как в заметке о признаках, тоже уменьшают число обусловленности, а следующая заметка показывает, что то же число управляет скоростью градиентного спуска.

Как измерить качество регрессии#

Среднеквадратичная ошибка измеряется в квадратах единиц целевой переменной; корень из неё, RMSE, — в единицах целевой переменной. Коэффициент детерминации (coefficient of determination) сравнивает модель с лучшей константой — средним $\bar y$:

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

На обучающих данных модели со свободным членом разложение по теореме Пифагора из предыдущих разделов помещает $R^2$ между 0 и 1: это доля разброса $y$, которую объясняет модель. В примере с четырьмя точками $\bar y = 2.75$, полная сумма квадратов равна $8.75$, и $R^2 = 1 - 2.7/8.75 = 0.691$.

На данных о диабете из scikit-learn (442 пациента, 10 центрированных и масштабированных измерений, целевая переменная — мера прогрессирования болезни через год) у метода наименьших квадратов $R^2$ на обучении равен $0.518$, а $R^2$ на 5-кратной кросс-валидации — $0.489$. Ближайшие соседи на стандартизованных признаках достигают в лучшем случае $0.454$ (при $k = 20$; $0.382$, $0.444$ и $0.428$ при $k = 5$, 10 и 40). Одиннадцать коэффициентов, подобранных по всем данным, здесь побеждают локальное голосование, потому что зависимость близка к аддитивной, а 442 точек мало для окрестностей в десятимерном пространстве.

У $R^2$ три ловушки:

• $R^2$ на обучении не убывает с числом признаков. Добавленный столбец может только сохранить или расширить пространство столбцов, поэтому остаток проекции никогда не растёт. При 50 обучающих точках, одном полезном признаке и $k$ столбцах чистого шума $R^2$ на обучении равен $0.572$, $0.649$, $0.747$, $0.921$ при $k = 0, 10, 20, 40$, а на 50 новых точках — $0.325$, $0.154$, $-0.158$, $-0.71$. Отрицательный $R^2$ означает, что модель хуже прогноза средним, и он появляется только на данных, которых модель не видела. Сообщайте $R^2$ на валидационных данных, как заметка о валидации требует для любой метрики.
• $R^2$ зависит от разброса $x$. Для $y = x + \varepsilon$ с $\sigma = 1$ одна и та же модель с тем же шумом имеет $R^2 = (1/12)/(1/12 + 1) = 0.077$, когда $x$ равномерно распределён на $[0, 1]$, и $0.893$, когда на $[0, 10]$; RMSE в обоих случаях равен 1. $R^2$ сравнивает модель со средним, а среднее — тем худший прогноз, чем сильнее меняется $y$.
• Высокий $R^2$ не означает правильную модель. Прямая, проведённая через точки на параболе, может иметь $R^2$ около 1 и остатки, которые все положительны на краях и отрицательны в середине. Смотрите на остатки.

Квадратичная ошибка и выбросы#

Квадратичная ошибка оценивает условное среднее, а среднее следует за самыми большими значениями. Другие функции потерь оценивают другие величины (другие функции потерь — другие цели): лучшая константа при абсолютной ошибке — медиана. Для пяти значений $1, 2, 4, 8, 16$ лучшая константа при квадратичной ошибке — их среднее $6.2$, при абсолютной ошибке — их медиана $4$. Поэтому одно неверное значение может сдвинуть решение метода наименьших квадратов сколь угодно далеко, а решение на основе медианы оно сдвигает, только пока не пересечёт другие точки.

Рисунок

Один выброс и три функции потерь

У выброса на рисунке остаток по методу наименьших квадратов равен $6.07$; его квадрат, $36.9$, больше суммы квадратов остатков остальных девяти точек. Метод наименьших квадратов наклоняет прямую к нему: наклон растёт с $0.48$ без этой точки до $0.82$ с ней. Влияние больших остатков ограничивают две функции потерь:

• Абсолютная ошибка, $\rho(r) = \lvert r\rvert$, — метод наименьших модулей (least absolute deviations). Каждая точка тянет с одной и той же силой независимо от расстояния. Подобранная прямая здесь — $1.27 + 0.47x$; точка минимума не всегда единственна, а функция потерь недифференцируема в 0.
• Функция потерь Хьюбера (Huber loss; Huber, 1964): $\rho_\delta(r) = r^2/2$ при $\lvert r\rvert \le \delta$ и $\delta(\lvert r\rvert - \delta/2)$ дальше. Квадратичная для малых остатков, поэтому сохраняет эффективность метода наименьших квадратов на чистых гауссовских данных; линейная для больших, поэтому сила, с которой тянет выброс, ограничена величиной $\delta$. При $\delta = 1$ прямая — $0.99 + 0.53x$.

Ограничена производная по остатку. Градиент по весам умножает эту производную на $x_i$, поэтому точка с крайним значением входа, точка с большим рычагом (high leverage), по-прежнему тянет сильно. Десять точек на прямой $y = x$ и одна грубая ошибка в $(100, 0)$ одинаково выравнивают все три подгонки: наклоны наименьших квадратов, Хьюбера и наименьших модулей равны $-0.04$. Эти функции потерь защищают от выбросов по $y$, а не по $x$.

Ни у одной из них нет замкнутой формулы. Обе можно подобрать, многократно решая задачу взвешенных наименьших квадратов с весами, которые уменьшаются для больших остатков, — так их вычисляет рисунок, — или общим методом из следующей заметки.

Что изменилось#

	Ближайшие соседи	Метод наименьших квадратов
Модель	среднее целевых значений соседей	глобальная линейная функция с $d$ коэффициентами
Обучение	запомнить данные	решить $X^\top X w = X^\top y$, $O(nd^2 + d^3)$
Прогноз	поиск по обучающей выборке	одно скалярное произведение, $O(d)$
Диабет, $R^2$ на кросс-валидации	$0.454$ в лучшем случае ($k = 20$)	$0.489$
Неопределённость	не встроена	$\operatorname{Var}\hat w = \sigma^2 (X^\top X)^{-1}$, стандартные ошибки
Интерпретация	нет	по весу на признак, если признаки не коллинеарны
Слабое место	большая размерность, масштаб	нелинейные зависимости, коллинеарность, выбросы, $d$ близко к $n$

Ограничения#

• Линейность по выбранным признакам. Модель видит взаимодействия и кривизну, только если их кто-то заложил в $X$; раздел о деревьях и раздел о нейросетях выучивают их сами.
• Неустойчивость при коррелированных признаках или $d$, близком к $n$. Дисперсия коэффициентов взрывается вдоль почти постоянных направлений; при $d \ge n$ и $\operatorname{rank} X = n$ ошибка на обучении равна нулю, так что подгонка на обучающих данных ничего не говорит об обобщении. Какое из многих интерполирующих решений выбрано (например, с наименьшей нормой) и геометрия данных всё равно решают, полезны ли его прогнозы, как показывает оговорка о двойном спуске.
• Чувствительность к выбросам. Одна грубая ошибка может повернуть прямую, потому что квадратичная функция потерь взвешивает остатки по их величине; устойчивые функции потерь помогают против выбросов по $y$, но не против крайних значений входов.
• Формула дисперсии предполагает некоррелированный шум с одинаковой дисперсией. При гетероскедастичных или коррелированных ошибках (временные ряды, сгруппированные данные) коэффициенты остаются несмещёнными, но стандартные ошибки неверны.
• Замкнутая формула стоит $O(nd^2 + d^3)$ и существует только для квадратичной ошибки.

Что нужно запомнить#

• Метод наименьших квадратов минимизирует $\lVert y - Xw\rVert^2$; нормальные уравнения $X^\top X\hat w = X^\top y$ имеют единственное решение $(X^\top X)^{-1}X^\top y$, когда столбцы $X$ независимы.
• Геометрически $\hat y$ — ортогональная проекция $y$ на пространство столбцов; остаток ортогонален каждому признаку, а при наличии свободного члена остатки в сумме дают ноль.
• При гауссовском шуме метод наименьших квадратов — это метод максимального правдоподобия, и $\hat w \sim N\big(w, \sigma^2 (X^\top X)^{-1}\big)$; $\sigma^2$ оценивают как $\mathrm{RSS}/(n - d)$.
• Почти коллинеарные признаки делают отдельные коэффициенты бессмысленными, хотя прогнозы остаются устойчивыми; здесь вступает регуляризация.
• Никогда не решайте нормальные уравнения численно для плохо обусловленных признаков: построение $X^\top X$ возводит число обусловленности в квадрат. Используйте QR или SVD (lstsq).
• $R^2$ сравнивает модель с прогнозом средним; на обучающих данных он не убывает при добавлении признаков, поэтому сообщайте его на валидационных данных.
• Квадратичная ошибка оценивает условное среднее и следует за выбросами; абсолютная ошибка оценивает медиану, а функция потерь Хьюбера — промежуточный вариант между ними.

Код#

Разобранные примеры, выборочное распределение коэффициентов, сравнение с ближайшими соседями на данных о диабете, ловушки $R^2$, коллинеарность, обусловленность и устойчивые подгонки:

import numpy as np
from scipy import optimize
from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold, cross_val_score
from sklearn.neighbors import KNeighborsRegressor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

rng = np.random.default_rng(0)
r2 = lambda y, p: 1 - ((y - p) ** 2).sum() / ((y - y.mean()) ** 2).sum()

# разобранный пример: четыре точки, нормальные уравнения, остатки ортогональны столбцам
x = np.array([0., 1, 2, 3]); y = np.array([1., 3, 2, 5])
X = np.column_stack([np.ones(4), x])
w = np.linalg.solve(X.T @ X, X.T @ y)
r = y - X @ w
H = X @ np.linalg.inv(X.T @ X) @ X.T
print(w, r.round(2), (X.T @ r).round(12) + 0, round(r @ r, 2), round(r2(y, X @ w), 3), np.diag(H).round(2))
# [1.1 1.1] [-0.1  0.8 -1.3  0.6] [0. 0.] 2.7 0.691 [0.7 0.3 0.3 0.7]
s2 = r @ r / (4 - 2)                                  # несмещённая оценка дисперсии шума, n - d = 2
print(round(s2, 2), np.sqrt(s2 * np.diag(np.linalg.inv(X.T @ X))).round(3))   # 1.35 [0.972 0.52 ]

# выборочное распределение w-hat: 30 точек на [0, 1], w = (1, 2), sigma = 0.5, 20 000 реализаций шума
Xs = np.column_stack([np.ones(30), np.linspace(0, 1, 30)])
W = np.array([np.linalg.lstsq(Xs, Xs @ [1, 2] + 0.5 * rng.normal(size=30), rcond=None)[0] for _ in range(20_000)])
print(W.mean(0).round(3), W.std(0).round(3), (0.5 * np.sqrt(np.diag(np.linalg.inv(Xs.T @ Xs)))).round(3))
# [1.001 2.   ] [0.177 0.303] [0.178 0.306]: средние по моделированию, стандартные отклонения по моделированию и по теории

# диабет (442 пациента, 10 признаков): scikit-learn вычисляет именно замкнутую формулу; линейная модель против kNN
Xd, yd = load_diabetes(return_X_y=True)
A = np.column_stack([np.ones(len(yd)), Xd])
wd = np.linalg.lstsq(A, yd, rcond=None)[0]
lr = LinearRegression().fit(Xd, yd)
print(f'{np.abs(wd - np.r_[lr.intercept_, lr.coef_]).max():.0e}', round(r2(yd, A @ wd), 3))   # 2e-12 0.518
cv = KFold(5, shuffle=True, random_state=0)
print(round(cross_val_score(LinearRegression(), Xd, yd, cv=cv, scoring='r2').mean(), 3),
      [round(float(cross_val_score(make_pipeline(StandardScaler(), KNeighborsRegressor(k)), Xd, yd, cv=cv, scoring="r2").mean()), 3) for k in (5, 10, 20, 40)])
# 0.489 [0.382, 0.444, 0.454, 0.428]

# R² на обучении не убывает с бесполезными признаками; на новых данных он рушится (50 + 50 точек)
z = rng.normal(size=(100, 1)); t = z[:, 0] + rng.normal(size=100)
for k in (0, 10, 20, 40):
    Z = np.column_stack([np.ones(100), z, rng.normal(size=(100, k))])
    v = np.linalg.lstsq(Z[:50], t[:50], rcond=None)[0]
    print(k, round(r2(t[:50], Z[:50] @ v), 3), round(r2(t[50:], Z[50:] @ v), 3))
# 0 0.572 0.325 | 10 0.649 0.154 | 20 0.747 -0.158 | 40 0.921 -0.71

# два почти одинаковых признака: коэффициенты скачут, а их сумма и прогнозы — нет
x1 = rng.normal(size=100); x2 = x1 + 0.01 * rng.normal(size=100); X2 = np.column_stack([x1, x2])
print(round(np.corrcoef(x1, x2)[0, 1], 5), f'{np.linalg.cond(X2.T @ X2):.1e}')   # 0.99994 3.5e+04
for _ in range(3):
    v = np.linalg.lstsq(X2, x1 + x2 + rng.normal(size=100), rcond=None)[0]
    print(v.round(2), round(v.sum(), 2))
# [-5.    6.88] 1.88 | [11.88 -9.9 ] 1.98 | [-8.83 10.94] 2.11
print(np.sqrt(np.diag(np.linalg.inv(X2.T @ X2))).round(2), round(1 / np.sqrt(x1 @ x1), 2))   # [9.74 9.73] 0.1: стандартные ошибки при sigma = 1

# обусловленность: полиномиальные признаки степени 10 на 50 точках; нормальные уравнения против QR/SVD (lstsq)
V = np.vander(np.linspace(0, 1, 50), 11, increasing=True); yv = V @ np.ones(11)
print(f'{np.linalg.cond(V):.1e} {np.linalg.cond(V.T @ V):.1e}',
      f'{np.abs(np.linalg.solve(V.T @ V, V.T @ yv) - 1).max():.0e} {np.abs(np.linalg.lstsq(V, yv, rcond=None)[0] - 1).max():.0e}')
# 2.0e+07 4.1e+14 4e-03 4e-10

# один выброс: наименьшие квадраты, Хьюбер (delta = 1) и наименьшие модули
xo = np.arange(10.); yo = np.array([1.3, 1.1, 2.2, 3.0, 2.7, 3.6, 3.5, 4.9, 13.0, 5.3]); Xo = np.column_stack([np.ones(10), xo])
huber = lambda r: np.where(np.abs(r) <= 1, 0.5 * r ** 2, np.abs(r) - 0.5)
fits = {'ls': np.linalg.lstsq(Xo, yo, rcond=None)[0],
        'ls without x=8': np.linalg.lstsq(np.delete(Xo, 8, 0), np.delete(yo, 8), rcond=None)[0],
        'huber': optimize.minimize(lambda v: huber(yo - Xo @ v).sum(), [0, 0], method='BFGS', options={'gtol': 1e-10}).x,
        'lad': optimize.linprog(np.r_[0, 0, np.ones(10)], A_ub=np.block([[Xo, -np.eye(10)], [-Xo, -np.eye(10)]]),
                                b_ub=np.r_[yo, -yo], bounds=[(None, None)] * 2 + [(0, None)] * 10).x[:2]}
for k, v in fits.items(): print(k, v.round(3))
# ls [0.373 0.819] | ls without x=8 [1.107 0.477] | huber [0.986 0.533] | lad [1.267 0.467]

# грубая ошибка при крайнем значении входа: десять точек на y = x и одна в (100, 0)
xl = np.r_[np.arange(10.), 100]; yl = np.r_[np.arange(10.), 0]; Xl = np.column_stack([np.ones(11), xl])
print(round(float(np.linalg.lstsq(Xl, yl, rcond=None)[0][1]), 2),
      round(float(optimize.minimize(lambda v: huber(yl - Xl @ v).sum(), [0, 1], method='BFGS', options={'gtol': 1e-10}).x[1]), 2),
      round(float(optimize.linprog(np.r_[0, 0, np.ones(11)], A_ub=np.block([[Xl, -np.eye(11)], [-Xl, -np.eye(11)]]),
                                   b_ub=np.r_[yl, -yl], bounds=[(None, None)] * 2 + [(0, None)] * 11).x[1]), 2))
# -0.04 -0.04 -0.04: наклоны наименьших квадратов, Хьюбера и наименьших модулей

Дальше#

Нормальным уравнениям нужна вся матрица $X^\top X$ и решение системы размера $d \times d$, что недостижимо для миллиона признаков, и они существуют только для квадратичной ошибки: у функций потерь Хьюбера и абсолютной ошибки из последнего раздела, как и у логистической функции потерь из раздела о классификации, замкнутой формулы нет вовсе. Следующая заметка заменяет формулу итерацией, которой нужен только градиент, здесь $2X^\top(Xw - y)$, и показывает, что её скорость задаёт то же число обусловленности, из-за которого нормальные уравнения были неточны.

Метод наименьших квадратов встречается во всех курсах: векторы слов GloVe подбираются взвешенным методом наименьших квадратов по логарифмам числа совместных вхождений, GAN с функцией потерь наименьших квадратов заменяет кросс-энтропию дискриминатора квадратичной ошибкой, а метод чередующихся наименьших квадратов подбирает матричную факторизацию в рекомендательных системах.

Литература#

• E. Sokolov et al. Машинное обучение (HSE), конспекты лекций, лекция 2: линейные модели, функции потерь и метрики качества регрессии и лекция 3: переобучение, обучение линейной регрессии.
• T. Hastie, R. Tibshirani, J. Friedman. The Elements of Statistical Learning, 2nd ed. Springer, 2009. Раздел 3.2 (метод наименьших квадратов, его геометрия и выборочное распределение коэффициентов), раздел 3.2.2 (теорема Гаусса–Маркова).
• C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006. Раздел 3.1.1 (максимальное правдоподобие и наименьшие квадраты), раздел 3.1.2 (геометрия метода наименьших квадратов).
• M. P. Deisenroth, A. A. Faisal, C. S. Ong. Mathematics for Machine Learning. Cambridge University Press, 2020. Глава 9 (линейная регрессия), раздел 3.8 (ортогональные проекции).
• L. N. Trefethen, D. Bau III. Numerical Linear Algebra. SIAM, 1997. Лекция 11 (задачи наименьших квадратов), лекция 19 (устойчивость алгоритмов наименьших квадратов: нормальные уравнения против QR и SVD).
• P. J. Huber. Robust estimation of a location parameter. Annals of Mathematical Statistics, 35(1), 1964.