МНК и нормальные уравнения: линейная алгебра подгонки прямой

Задача, в которой нет ни одного лишнего числа: по точкам $(1,1)$, $(2,3)$, $(3,4)$ провести прямую $y=c+mx$. Попробуйте — и быстро упрётесь в невозможное: прямая по двум точкам определяется однозначно, третья почти всегда мимо. Пусть измерения содержат погрешность — тогда честная постановка другая: найти прямую, для которой сумма квадратов вертикальных отклонений минимальна. Это и есть метод наименьших квадратов (МНК), и решается он не перебором, а двумя матричными строками — нормальными уравнениями. Заодно увидите, что весь метод — это знакомая ортогональная проекция в новой одежде.

Переопределённая система: решения нет, задача есть#

Запишем требование «точка лежит на прямой» для каждой точки: $c+m\cdot 1=1$, $c+m\cdot 2=3$, $c+m\cdot 3=4$. Три уравнения, две неизвестных — система переопределённая. В матричном виде:

Точного решения у такой системы, как правило, нет: столбец $\vec{b}$ не лежит в плоскости, натянутой на столбцы $A$. Зато можно искать наилучшее приближение — вектор $A{\vec{x}}^{\ast }$, ближайший к $\vec{b}$, то есть минимизирующий длину невязки $\vec{r}=\vec{b}-A\vec{x}$:

Почему квадраты, а не модули или просто сумма? Три причины: квадрат сглаживает функцию (можно дифференцировать), штрафует крупные отклонения сильнее мелких и — что приятнее всего — задача распадается в чистую линейную алгебру, без численных итераций. Ответ будет точным, за один проход.

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

Геометрический вывод. Вектор $A\vec{x}$ при переборе всех $\vec{x}$ пробегает пространство столбцов $A$ — плоскость. Минимизируемая величина — квадрат расстояния от $\vec{b}$ до этой плоскости, а достигается минимум в проекции $\vec{b}$ на неё. Проекция узнаётся по невязке: она обязана быть перпендикулярна всей плоскости, то есть обоим столбцам $A$. Два условия перпендикулярности $(A\vec{x}-\vec{b},A{\vec{a}}_1)=0$ и $(A\vec{x}-\vec{b},A{\vec{a}}_2)=0$ собираются в одно матричное равенство:

Аналитический вывод даёт то же самое честным дифференцированием: берём частные производные $Q$ по $c$ и по $m$, приравниваем к нулю — и после группировки слагаемых получаем в точности $A^{\top }A\vec{x}=A^{\top }\vec{b}$. Совпадение двух выводов — красивое подтверждение: минимум суммы квадратов — это проекция, просто увиденная со стороны вычислений. Свойства ортогональности невязки разобраны в уроке про евклидовы пространства.

Техническое замечание, которое экономит баллы: система $A\vec{x}=\vec{b}$ может не иметь решений, но нормальная система $A^{\top }A\vec{x}=A^{\top }\vec{b}$ имеет решение всегда. А если столбцы $A$ линейно независимы, то $A^{\top }A$ обратима и решение единственно — вот почему для МНК важно, чтобы признаки не дублировали друг друга. Проверяется это через ранг или определитель матрицы $A^{\top }A$.

Пример по точкам: считаем до конца#

Возвращаемся к точкам $(1,1)$, $(2,3)$, $(3,4)$. Матрицы уже выписаны выше. Считаем обе части нормальных уравнений — тут потребуется только умножение матриц:

Нормальная система: $3c+6m=8$, $6c+14m=19$. Определитель $3\cdot 14-6\cdot 6=42-36=6\ne 0$ — решение единственно. Пройдите расчёт по шагам, каждая дробь здесь короткая:

Шаг второй — самое место для практики на аналогичной системе: механизм исключения тот же, что и в методе Гаусса, просто матрица двух на два. Нажимайте «Шаг» и следите за ведущим элементом:

Мост в статистику#

Решение нормальных уравнений можно выписать формулами, знакомыми по школьной «статистике». Для подгонки с свободным членом: наклон $m=\frac{\sum (x_i-\bar{x})(y_i-\bar{y})}{\sum (x_i-\bar{x}{)}^2}$ и сдвиг $c=\bar{y}-m\bar{x}$, где черта — среднее. Подставьте наши точки: $\bar{x}=2$, $\bar{y}=\frac{8}{3}$, числитель $(-1)(-\frac{5}{3})+0+(1)(\frac{4}{3})=3$, знаменатель $2$, итого $m=\frac{3}{2}$ и $c=\frac{8}{3}-3=-\frac{1}{3}$ — те же значения, третьим независимым способом. Никакой магии: это та же линейная система, свёрнутая в формулы.

Отсюда же — язык корреляции и регрессии: наклон МНК-прямой выражается через коэффициент корреляции и выборочные дисперсии, а уравнение регрессии — та самая $y=c+mx$. Разница в интерпретации: линал говорит «проекция», статистика — «оценка связи»; математика одна.

Типовые ошибки#

• Минимизировать сумму отклонений без квадратов: положительные и отрицательные невязки взаимно гасятся, и «лучшей» окажется любая прямая с нулевой суммой. Квадраты обязательны.
• Перепутать порядок: нормальные уравнения — это $A^{\top }A\vec{x}=A^{\top }\vec{b}$, а не $A{A}^{\top }$. Размеры подскажут: $A^{\top }A$ квадратная $2\times 2$ при двух неизвестных, $A{A}^{\top }$ — $3\times 3$, не то.
• Забыть столбец единиц в $A$: без него модель получится без свободного члена, $y=mx$, и ответ будет систематически сдвинут. Столбец единиц — часть конструкции.
• После решения не проверить невязки. Одна строка — сумма невязок ноль, сумма их квадратов минимальна — ловит арифметические промахи мгновенно.

Числовой сквозной пример: подгоним прямую $y=kx+b$ по трём точкам $(0,1,1)$, $(1,2,0)$, $(2,3,4)$. Матрица плана $A$ состоит из строк $(1,x_i)$, вектор $b$ — из $y_i$. Нормальные уравнения $A^{T}A\vec{z}=A^{T}b$ дают: $A^{T}A=\left(\begin{array}{cc}3& 3\\ 3& 5\end{array}\right)$, $A^{T}b=(6,5;10,9)$, откуда $k\approx 1,15$, $b\approx 1,02$. Проверка качества: прогнозы $1,02$, $2,17$, $3,32$ против фактических $1,1$, $2,0$, $3,4$ — каждая точка объяснена с ошибкой меньше $0,17$. Именно эта механика, записанная один раз, переезжает в эконометрику и статистику практически без изменений.

Куда дальше: симметричный сосед по курсу — графический метод линейного программирования, где тоже минимизируется функция на плоскости, но с ограничениями-неравенствами вместо квадратов. Систематическое применение МНК к данным — в уроке корреляция и регрессия; решать нормальные системы механически тренирует решатель Гаусса.