Ортогональность, МНК и SVD

Прошлая глава оставила нас перед облаком точек, через которое не провести прямую. Возьмём облако поменьше — из трёх точек. К пружине подвешивают груз в $1$, $2$ и $3$ кг, и стрелка на шкале встаёт на деления $1$, $2$ и $2$. По закону Гука показание растёт линейно: $y = C + Dt$, где $t$ — масса груза, $C$ — где стояла бы стрелка без груза, а $D$ — на сколько делений её сдвигает каждый килограмм. Три замера дают три уравнения с двумя неизвестными:

$$C + D = 1, \qquad C + 2D = 2, \qquad C + 3D = 2.$$

Из первых двух $D = 1$ и $C = 0$, а третье тогда требует $3 = 2$. Прямой через три точки нет. Но жёсткость у пружины есть, просто дрогнула рука или шкала нанесена неровно. Задача меняется: найти прямую, которая ошибается меньше всех. Что значит «меньше всех», придётся договориться, и договорённость двухсотлетней давности оказалась удачнее, чем можно было ожидать. Из неё вырастет метод наименьших квадратов, а из него — разложение любой матрицы на поворот, растяжение и ещё один поворот.

Глава устроена как лабораторный журнал. В нём шесть опытов с данными, а между опытами — инструменты, без которых опыт не поставить.

Опыт 1. Три замера и ни одной прямой

Запишем систему матрицей, как в главе о методе Гаусса:

$$\underbrace{\begin{pmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \end{pmatrix}}_{A}\underbrace{\begin{pmatrix} C \\ D \end{pmatrix}}_{\mathbf x} = \underbrace{\begin{pmatrix} 1 \\ 2 \\ 2 \end{pmatrix}}_{\mathbf b}.$$

Теперь посмотрим на неё по столбцам. Произведение $A\mathbf x$ — это $C\,\mathbf a_1 + D\,\mathbf a_2$, где $\mathbf a_1 = (1;\,1;\,1)$ и $\mathbf a_2 = (1;\,2;\,3)$ — столбцы $A$. Все такие комбинации заполняют плоскость в трёхмерном пространстве, проходящую через начало координат: это линейная оболочка столбцов, или образ матрицы (глава 37). Решить систему — значит добраться до точки $\mathbf b = (1;\,2;\,2)$, не сходя с этой плоскости. А $\mathbf b$ лежит вне её, поэтому решения нет.

Раз до $\mathbf b$ не дойти, дойдём до ближайшей к ней точки плоскости. Расстояние в пространстве — корень из суммы квадратов разностей координат, как у диагонали коробки из главы о теореме Пифагора. А координаты разности $\mathbf b - A\mathbf x$ — это в точности промахи прямой в трёх замерах: $y_i - (C + Dt_i)$.

Получается, что «ближайшая точка плоскости» и «прямая с наименьшей суммой квадратов промахов» — одна и та же задача, сказанная на двух языках. Геометрию мы знаем лучше, поэтому решать будем на её языке. Правда, замеров обычно не три, а сотни, и вектор $\mathbf b$ живёт в пространстве $\mathbb R^{100}$ или $\mathbb R^{1000}$. Нарисовать его не получится, так что нужна линейка, которой всё равно, сколько у пространства измерений.

Линейка для любого числа измерений

В главе о векторах скалярное произведение на плоскости считалось по координатам: $a_1b_1 + a_2b_2$. Для векторов из $\mathbb R^n$ формула та же, только слагаемых $n$:

$$\langle \mathbf u, \mathbf v\rangle = u_1v_1 + u_2v_2 + \dots + u_nv_n = \mathbf u^{\mathsf T}\mathbf v.$$

Через него определяются длина $|\mathbf u| = \sqrt{\langle \mathbf u, \mathbf u\rangle}$ (теорема Пифагора, применённая $n - 1$ раз) и расстояние $|\mathbf u - \mathbf v|$. На плоскости скалярное произведение было ещё и равно $|\mathbf u|\,|\mathbf v|\cos\varphi$. В $\mathbb R^n$ транспортира нет, и этим равенством угол как раз определяют.

Угол между ненулевыми векторами, от $0^\circ$ до $180^\circ$. Скалярное произведение: сумма попарных произведений координат. Произведение длин. Деление на него убирает масштаб: удлините любой вектор вдвое, и угол не изменится. Пример: для столбцов $\mathbf a_1 = (1;\,1;\,1)$ и $\mathbf a_2 = (1;\,2;\,3)$ получаем $\langle \mathbf a_1, \mathbf a_2\rangle = 6$, $|\mathbf a_1| = \sqrt3$, $|\mathbf a_2| = \sqrt{14}$ и $\cos\theta = \frac{6}{\sqrt{42}} \approx 0{,}926$, то есть $\theta \approx 22{,}2^\circ$. Столбцы нашей матрицы смотрят почти в одну сторону.

Определение честное, только если дробь справа никогда не выходит за пределы отрезка $[-1;\,1]$, иначе такого угла просто нет. Это и утверждает главное неравенство линейной алгебры.

Хитрость в том, что два вектора всегда лежат в одной плоскости, даже в $\mathbb R^{1000}$. Опустим перпендикуляр из конца $\mathbf u$ на прямую вектора $\mathbf v$: тень не может оказаться длиннее самого вектора. Все шаги ниже — вычисления со скалярным произведением, так что чертёж ничего не подсказывает сверх выкладок.

Если $\mathbf v = \mathbf 0$, обе части равны нулю. Иначе рассмотрим вектор $\mathbf p = \frac{\langle \mathbf u, \mathbf v\rangle}{\langle \mathbf v, \mathbf v\rangle}\,\mathbf v$ на прямой вектора $\mathbf v$. Его длина равна $\frac{|\langle \mathbf u, \mathbf v\rangle|}{|\mathbf v|^2}\,|\mathbf v| = \frac{|\langle \mathbf u, \mathbf v\rangle|}{|\mathbf v|}$. Остаток $\mathbf e = \mathbf u - \mathbf p$ перпендикулярен $\mathbf v$. Скалярное произведение линейно по каждому аргументу, поэтому $\langle \mathbf e, \mathbf v\rangle = \langle \mathbf u, \mathbf v\rangle - \frac{\langle \mathbf u, \mathbf v\rangle}{\langle \mathbf v, \mathbf v\rangle}\langle \mathbf v, \mathbf v\rangle = 0$. Вектор $\mathbf p$ кратен $\mathbf v$, так что и $\langle \mathbf p, \mathbf e\rangle = 0$. Теорема Пифагора в $\mathbb R^n$: раскроем скобки в $|\mathbf u|^2 = \langle \mathbf p + \mathbf e, \mathbf p + \mathbf e\rangle = |\mathbf p|^2 + 2\langle \mathbf p, \mathbf e\rangle + |\mathbf e|^2$. Среднее слагаемое равно нулю по предыдущему шагу, и остаётся $|\mathbf u|^2 = |\mathbf p|^2 + |\mathbf e|^2$. Квадрат $|\mathbf e|^2$ не отрицателен, поэтому $|\mathbf p| \le |\mathbf u|$, то есть $\frac{|\langle \mathbf u, \mathbf v\rangle|}{|\mathbf v|} \le |\mathbf u|$. Умножив на $|\mathbf v|$, получаем неравенство. Равенство означает $\mathbf e = \mathbf 0$: вектор $\mathbf u$ совпал со своей тенью $\mathbf p$ и кратен $\mathbf v$.

Для сумм это неравенство опубликовал Огюстен Коши в 1821 году, для интегралов — Виктор Буняковский в 1859-м и независимо Герман Шварц в 1888-м; в западных книгах его обычно называют неравенством Коши — Шварца. Интегралы здесь не случайны: в главе о рядах Фурье скалярным произведением функций служил $\int f g\,dx$, и всё, что мы докажем о векторах, верно и для сигналов.

Самый полезный угол — прямой.

При каком $a$ векторы $(1;\,2;\,a;\,1)$ и $(3;\,-1;\,2;\,1)$ из $\mathbb R^4$ ортогональны?

Скалярное произведение равно $1 \cdot 3 + 2 \cdot (-1) + a \cdot 2 + 1 \cdot 1 = 2 + 2a$. Оно обращается в ноль при $a = -1$.

Тень на плоскость

Доказательство неравенства попутно решило задачу о ближайшей точке прямой. Тень $\mathbf p = \frac{\langle \mathbf b, \mathbf a\rangle}{\langle \mathbf a, \mathbf a\rangle}\,\mathbf a$ вектора $\mathbf b$ на прямую вектора $\mathbf a$ — вектор, длина которого со знаком и есть проекция из главы о векторах: та полезная часть усилия бурлаков, что тянула баржу вдоль реки. Теперь нам нужна тень не на прямую, а на плоскость и вообще на любое подпространство. Главное свойство тени от этого не меняется.

Идея: соедините $\mathbf b$ с любой другой точкой $\mathbf y$ подпространства — получится прямоугольный треугольник, и перпендикуляр в нём короче гипотенузы. Чертёж изображает плоскость, в которой лежат $\mathbf b$, $\mathbf p$ и $\mathbf y$; в $\mathbb R^n$ такая плоскость найдётся всегда.

Пусть $\mathbf p \in W$ и $\mathbf b - \mathbf p$ ортогонален $W$. Возьмём любую другую точку $\mathbf y \in W$. Разность $\mathbf p - \mathbf y$ тоже лежит в $W$: подпространство замкнуто относительно вычитания. Значит, $\langle \mathbf b - \mathbf p, \mathbf p - \mathbf y\rangle = 0$, и угол треугольника при вершине $\mathbf p$ прямой. Так как $\mathbf b - \mathbf y = (\mathbf b - \mathbf p) + (\mathbf p - \mathbf y)$, раскрытие скобок, как в доказательстве неравенства Коши — Буняковского, даёт теорему Пифагора: $|\mathbf b - \mathbf y|^2 = |\mathbf b - \mathbf p|^2 + |\mathbf p - \mathbf y|^2$. Второе слагаемое положительно при $\mathbf y \ne \mathbf p$, поэтому $\mathbf y$ дальше от $\mathbf b$, чем $\mathbf p$. Заодно доказана единственность. Обратно: пусть $\mathbf b - \mathbf q$ не ортогонален какому-то вектору $\mathbf w \in W$. Двигаясь из $\mathbf q$ вдоль $\mathbf w$, мы идём по хорде окружности с центром $\mathbf b$ и радиусом $|\mathbf b - \mathbf q|$, которая проходит через $\mathbf q$. Хорда — не касательная, потому что радиус в $\mathbf q$ к ней не перпендикулярен, и внутри круга есть точки $W$, которые ближе к $\mathbf b$. В числах: $|\mathbf b - \mathbf q - s\mathbf w|^2 = |\mathbf b - \mathbf q|^2 - 2s\langle \mathbf b - \mathbf q, \mathbf w\rangle + s^2|\mathbf w|^2$, и при малом $s$ нужного знака это меньше $|\mathbf b - \mathbf q|^2$. Ближайшая точка подпространства — основание перпендикуляра, опущенного из $\mathbf b$, и только оно.

Вернёмся к пружине. Подпространство $W$ — образ матрицы $A$, его векторы имеют вид $A\mathbf x$. Разность $\mathbf b - A\hat{\mathbf x}$ ортогональна $W$ ровно тогда, когда она ортогональна каждому столбцу: $\mathbf a_j^{\mathsf T}(\mathbf b - A\hat{\mathbf x}) = 0$. Строки $\mathbf a_j^{\mathsf T}$ — это строки транспонированной матрицы, так что все условия разом записываются как $A^{\mathsf T}(\mathbf b - A\hat{\mathbf x}) = \mathbf 0$.

Квадратная симметричная матрица, составленная из скалярных произведений столбцов $A$ друг на друга: в клетке $(i, j)$ стоит $\langle \mathbf a_i, \mathbf a_j\rangle$. Её размер равен числу неизвестных, а не числу замеров. Лучшее приближённое решение: при нём сумма квадратов невязок $|\mathbf b - A\mathbf x|^2$ наименьшая. Скалярные произведения столбцов на правую часть. Пример: у пружины $A^{\mathsf T}A = \left(\begin{smallmatrix} 3 & 6 \\ 6 & 14 \end{smallmatrix}\right)$ и $A^{\mathsf T}\mathbf b = (5;\,11)$. Система $3C + 6D = 5$, $6C + 14D = 11$ даёт $D = \frac12$, $C = \frac23$. Прямая $y = \frac23 + \frac t2$ предсказывает $\frac76$, $\frac53$, $\frac{13}{6}$, невязки равны $-\frac16$, $\frac13$, $-\frac16$, а сумма их квадратов — $\frac16$. Проверка перпендикулярности: $-\frac16 + \frac13 - \frac16 = 0$ и $-\frac16 + \frac23 - \frac12 = 0$.

Остался вопрос, всегда ли нормальные уравнения решаются однозначно.

Одна и та же задача живёт в двух картинках сразу. На плоскости замеров мы двигаем прямую и смотрим на промахи. В пространстве, где каждая ось — это один замер, мы опускаем перпендикуляр из $\mathbf b$ на плоскость столбцов. Сравните их.

Опыт 2. Лучшая прямая

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

Для $n$ точек $(t_i;\,y_i)$ у матрицы $A$ две колонки: единицы и значения $t_i$. Нормальные уравнения превращаются в систему двух уравнений, которую можно решить раз и навсегда. Обозначим средние чертой: $\bar t = \frac1n\sum t_i$, $\bar y = \frac1n\sum y_i$.

Прямая $y = C + Dt$ проведена через $50$ точек по методу наименьших квадратов. Чему равна сумма невязок $e_1 + e_2 + \dots + e_{50}$?

Ровно нулю, при любых точках Наименьшему возможному числу Сумме квадратов невязок Зависит от расположения точек

Невязка ортогональна каждому столбцу $A$, в том числе столбцу из одних единиц: $\langle (1;\,\dots;\,1), \mathbf e\rangle = e_1 + \dots + e_n = 0$. Промахи вверх и вниз уравновешивают друг друга. Если прямую заставить проходить через начало координат (без $C$), столбца единиц нет, и это свойство теряется.

У квадратов есть оборотная сторона. Промах в $10$ делений весит столько же, сколько сто промахов в одно деление, поэтому единственная грубая ошибка — сбитый датчик, опечатка в таблице — может развернуть прямую, построенную по тысяче хороших точек. Попробуйте это в виджете. Можно минимизировать сумму модулей невязок: такую прямую выброс сдвигает гораздо меньше, но у неё нет формулы, и она бывает не единственной.

И ещё одно ограничение. Мы считали промахи по вертикали, потому что масса груза известна точно, а ошибается шкала. Если неточны обе координаты, честнее мерить расстояние до прямой по перпендикуляру, и тогда ответ другой — к нему приведёт опыт 5.

Зато слово «линейный» в названии метода относится к неизвестным, а не к графику. Параболу $y = C + Dt + Et^2$ подгоняют теми же нормальными уравнениями, добавив в $A$ третий столбец из чисел $t_i^2$: неизвестные $C$, $D$, $E$ по-прежнему входят в каждое уравнение в первой степени. Экспоненту $y = ke^{rt}$ сводят к прямой логарифмом: $\ln y = \ln k + rt$.

Четыре замера: $(0;\,1)$, $(1;\,3)$, $(2;\,4)$, $(3;\,4)$. Какое значение $y$ при $t = 4$ предсказывает прямая, проведённая по методу наименьших квадратов? Ответ дайте числом.

Средние: $\bar t = 1{,}5$, $\bar y = 3$. Отклонения $t$: $-1{,}5$; $-0{,}5$; $0{,}5$; $1{,}5$, отклонения $y$: $-2$; $0$; $1$; $1$. Числитель $3 + 0 + 0{,}5 + 1{,}5 = 5$, знаменатель $2{,}25 + 0{,}25 + 0{,}25 + 2{,}25 = 5$, так что $D = 1$ и $C = 3 - 1{,}5 = 1{,}5$. При $t = 4$ прямая $y = 1{,}5 + t$ даёт $5{,}5$.

Проекции, прямые и нормальные уравнения можно потренировать.

Инструмент: прямые углы вместо уравнений

Нормальные уравнения приходится решать, потому что столбцы $\mathbf a_1$ и $\mathbf a_2$ не перпендикулярны: между ними $22^\circ$, и тень на один столбец «заходит» на другой. Если бы они были перпендикулярны и единичной длины, матрица $A^{\mathsf T}A$ стала бы единичной и решать было бы нечего: $\hat{\mathbf x} = A^{\mathsf T}\mathbf b$.

Координаты в таком базисе считаются без всяких систем. Если $\mathbf v = c_1\mathbf q_1 + \dots + c_k\mathbf q_k$, умножим обе части скалярно на $\mathbf q_j$: все слагаемые, кроме $j$-го, обнулятся, и останется $c_j = \langle \mathbf v, \mathbf q_j\rangle$. По той же причине проекция любого $\mathbf b$ на оболочку равна $\langle \mathbf b, \mathbf q_1\rangle\mathbf q_1 + \dots + \langle \mathbf b, \mathbf q_k\rangle\mathbf q_k$: разность $\mathbf b$ минус эта сумма ортогональна каждому $\mathbf q_j$. Коэффициенты Фурье из главы 33 вычислялись точно так же. Осталось научиться выпрямлять базис.

Строим векторы по одному: из каждого нового вычитаем его тени на уже построенные, а остаток нормируем. На чертеже — первые два шага; дальше всё повторяется.

Вектор $\mathbf a_1$ ненулевой: иначе набор был бы зависимым. Положим $\mathbf q_1 = \mathbf a_1 / |\mathbf a_1|$ — единичный вектор с той же оболочкой. Вычтем из $\mathbf a_2$ его тень на $\mathbf q_1$: $\mathbf a_2' = \mathbf a_2 - \langle \mathbf a_2, \mathbf q_1\rangle\mathbf q_1$. Остаток ортогонален $\mathbf q_1$: $\langle \mathbf a_2', \mathbf q_1\rangle = \langle \mathbf a_2, \mathbf q_1\rangle - \langle \mathbf a_2, \mathbf q_1\rangle\langle \mathbf q_1, \mathbf q_1\rangle = 0$, ведь $\langle \mathbf q_1, \mathbf q_1\rangle = 1$. Остаток не нулевой: иначе $\mathbf a_2$ был бы кратен $\mathbf a_1$. Положим $\mathbf q_2 = \mathbf a_2' / |\mathbf a_2'|$. Оболочка не изменилась: $\mathbf q_1$, $\mathbf q_2$ выражаются через $\mathbf a_1$, $\mathbf a_2$ и наоборот, так как $\mathbf a_2 = |\mathbf a_2'|\,\mathbf q_2 + \langle \mathbf a_2, \mathbf q_1\rangle\mathbf q_1$. На шаге $j$ вычитаем тени на все построенные векторы: $\mathbf a_j' = \mathbf a_j - \sum_{i < j}\langle \mathbf a_j, \mathbf q_i\rangle\mathbf q_i$. Умножив на $\mathbf q_l$ при $l < j$, получим $\langle \mathbf a_j, \mathbf q_l\rangle - \langle \mathbf a_j, \mathbf q_l\rangle = 0$: из суммы выживает одно слагаемое, остальные гасит ортогональность $\mathbf q_i$. Остаток не нулевой, потому что $\mathbf a_j$ не лежит в оболочке предыдущих $\mathbf a$, которая совпадает с оболочкой предыдущих $\mathbf q$. Нормируем: $\mathbf q_j = \mathbf a_j' / |\mathbf a_j'|$. Индукция по $j$ даёт весь набор, и на каждом шаге оболочка сохраняется по тем же причинам, что и на втором.

Выпрямим столбцы пружины. Первый даёт $\mathbf q_1 = \frac{1}{\sqrt3}(1;\,1;\,1)$. Тень второго на него — $\frac{6}{3}(1;\,1;\,1) = (2;\,2;\,2)$, и остаток $(1;\,2;\,3) - (2;\,2;\,2) = (-1;\,0;\,1)$. Вычесть тень на столбец единиц — значит вычесть из $t_i$ их среднее $\bar t = 2$; вот откуда отклонения в формуле наклона. Теперь $\mathbf q_2 = \frac{1}{\sqrt2}(-1;\,0;\,1)$, и проекция $\mathbf b = (1;\,2;\,2)$ получается без единого уравнения: $\frac53(1;\,1;\,1) + \frac12(-1;\,0;\,1) = \left(\frac76;\,\frac53;\,\frac{13}{6}\right)$ — те же предсказания, что и раньше.

Процесс носит имена Йоргена Грама (1883) и Эрхарда Шмидта (1907). Если записать его матрицами, получится разложение $A = QR$: столбцы $Q$ ортонормированы, а $R$ — верхнетреугольная, потому что $\mathbf a_j$ выражается только через $\mathbf q_1, \dots, \mathbf q_j$.

На плоскости ортогональных матриц совсем немного. Первый столбец — единичный вектор $(\cos\varphi;\,\sin\varphi)$, второй — единичный и перпендикулярный ему, а таких два: $(-\sin\varphi;\,\cos\varphi)$ и $(\sin\varphi;\,-\cos\varphi)$. В первом случае получается поворот на угол $\varphi$ с определителем $1$, во втором — отражение с определителем $-1$. Это знакомые движения из главы о симметрии, те, что оставляют начало координат на месте. Например, $\frac15\left(\begin{smallmatrix} 3 & -4 \\ 4 & 3 \end{smallmatrix}\right)$ поворачивает плоскость на $\operatorname{arctg}\frac43 \approx 53{,}13^\circ$. Ортогональна и любая матрица перестановки, у которой в каждой строке и каждом столбце одна единица: она лишь переставляет координаты.

Опыт 3. Круг превращается в эллипс

У прямоугольной матрицы нет собственных векторов: она переводит векторы из одного пространства в другое, и сравнить $A\mathbf v$ с $\mathbf v$ нельзя. Даже у квадратной матрицы, например у сдвига, собственных направлений может быть мало. Но спросим иначе: во что матрица превращает единичную окружность?

Какую бы матрицу вы ни задали, окружность становится эллипсом (или отрезком), и у эллипса есть две перпендикулярные полуоси. Самое неожиданное — их прообразы: на окружности всегда найдутся два перпендикулярных вектора $\mathbf v_1$, $\mathbf v_2$, которые матрица переводит в эти полуоси, то есть снова в перпендикулярные векторы. Даже сдвиг, у которого всего одно собственное направление, одну пару прямых углов сохраняет. Это верно для любой матрицы.

Собственные векторы есть не у $A$, а у симметричной матрицы $A^{\mathsf T}A$, и для неё охота из прошлой главы всегда удачна. Её собственные векторы и будут нужными перпендикулярами. Слева на чертеже — пространство, где живут $\mathbf v$, справа — где живут $A\mathbf v$.

Матрица $A^{\mathsf T}A$ размера $n \times n$ симметрична: $(A^{\mathsf T}A)^{\mathsf T} = A^{\mathsf T}(A^{\mathsf T})^{\mathsf T} = A^{\mathsf T}A$. По спектральной теореме у неё есть ортонормированный базис из собственных векторов $\mathbf v_1, \dots, \mathbf v_n$ с собственными значениями $\lambda_1, \dots, \lambda_n$. Пронумеруем их по убыванию. Все $\lambda_i$ неотрицательны: $\lambda_i = \lambda_i\langle \mathbf v_i, \mathbf v_i\rangle = \langle \mathbf v_i, A^{\mathsf T}A\mathbf v_i\rangle = \langle A\mathbf v_i, A\mathbf v_i\rangle = |A\mathbf v_i|^2$. Положим $\sigma_i = \sqrt{\lambda_i} = |A\mathbf v_i|$ — это длины образов. Образы перпендикулярны: при $i \ne j$ имеем $\langle A\mathbf v_i, A\mathbf v_j\rangle = \langle \mathbf v_i, A^{\mathsf T}A\mathbf v_j\rangle = \lambda_j\langle \mathbf v_i, \mathbf v_j\rangle = 0$. Пусть $\sigma_1, \dots, \sigma_r$ — ненулевые из этих чисел. Векторы $\mathbf u_i = A\mathbf v_i / \sigma_i$ при $i \le r$ единичны и попарно ортогональны; дополним их до ортонормированного базиса $\mathbb R^m$ процессом Грама — Шмидта. При $i > r$ имеем $|A\mathbf v_i| = 0$, то есть $A\mathbf v_i = \mathbf 0$. Векторы $A\mathbf v_1, \dots, A\mathbf v_n$ порождают весь образ $A$, ненулевые среди них — первые $r$, а ненулевые попарно ортогональные векторы линейно независимы. Поэтому $r$ совпадает с рангом. Равенства $A\mathbf v_i = \sigma_i\mathbf u_i$ разом записываются как $AV = U\Sigma$, где столбцы $V$ и $U$ — векторы $\mathbf v_i$ и $\mathbf u_i$, а $\Sigma$ — прямоугольная матрица с числами $\sigma_i$ на диагонали и нулями вне её. Матрица $V$ ортогональна, $V^{-1} = V^{\mathsf T}$, и $A = U\Sigma V^{\mathsf T}$. Точка окружности $\cos s\,\mathbf v_1 + \sin s\,\mathbf v_2$ переходит в $\sigma_1\cos s\,\mathbf u_1 + \sigma_2\sin s\,\mathbf u_2$ — это и есть эллипс с полуосями $\sigma_1$ и $\sigma_2$. Ортогональная матрица $n \times n$ — поворот, возможно с отражением. Она переводит $\mathbf v_i$ в $\mathbf e_i$: сначала пространство разворачивают так, чтобы нужные направления легли на оси. Прямоугольная диагональная матрица $m \times n$: растягивает $i$-ю ось в $\sigma_i$ раз, а лишние координаты обнуляет или дописывает нулями, если $m \ne n$. Ортогональная матрица $m \times m$: поворачивает оси $\mathbf e_i$ в направления $\mathbf u_i$. Читаем, как всегда, справа налево. Пример: для $A = \left(\begin{smallmatrix} 3 & 0 \\ 4 & 5 \end{smallmatrix}\right)$ матрица $A^{\mathsf T}A = \left(\begin{smallmatrix} 25 & 20 \\ 20 & 25 \end{smallmatrix}\right)$ имеет собственные значения $45$ и $5$ с векторами $\mathbf v_1 = \frac{1}{\sqrt2}(1;\,1)$, $\mathbf v_2 = \frac{1}{\sqrt2}(1;\,-1)$. Значит, $\sigma_1 = 3\sqrt5 \approx 6{,}708$ и $\sigma_2 = \sqrt5 \approx 2{,}236$. Образы: $A\mathbf v_1 = \frac{1}{\sqrt2}(3;\,9)$ и $A\mathbf v_2 = \frac{1}{\sqrt2}(3;\,-1)$, они перпендикулярны, а $\mathbf u_1 = \frac{1}{\sqrt{10}}(1;\,3)$, $\mathbf u_2 = \frac{1}{\sqrt{10}}(3;\,-1)$. Проверка: $\sigma_1\sigma_2 = 15 = |\det A|$.

Из доказательства видно, как считать: сингулярные числа — корни из собственных значений $A^{\mathsf T}A$. Для квадратной матрицы их произведение равно $|\det A|$: повороты площадь не меняют, а $\Sigma$ растягивает её в $\sigma_1\sigma_2\cdots\sigma_n$ раз. У симметричной матрицы с неотрицательными собственными значениями сингулярное разложение совпадает со спектральным. А у сдвига $\left(\begin{smallmatrix} 1 & 1 \\ 0 & 1 \end{smallmatrix}\right)$, от которого охотник за собственными векторами ушёл почти ни с чем, сингулярные числа — золотое сечение $\varphi \approx 1{,}618$ и $1/\varphi \approx 0{,}618$: $A^{\mathsf T}A = \left(\begin{smallmatrix} 1 & 1 \\ 1 & 2 \end{smallmatrix}\right)$, и его характеристическое уравнение $\lambda^2 - 3\lambda + 1 = 0$ имеет корни $\varphi^2$ и $\varphi^{-2}$.

Найдите сингулярные числа матрицы $\left(\begin{smallmatrix} 3 & -8 \\ 4 & 6 \end{smallmatrix}\right)$.

Столбцы $(3;\,4)$ и $(-8;\,6)$ ортогональны: $-24 + 24 = 0$. Поэтому $A^{\mathsf T}A = \left(\begin{smallmatrix} 25 & 0 \\ 0 & 100 \end{smallmatrix}\right)$, собственные значения $100$ и $25$, а сингулярные числа — $10$ и $5$, длины столбцов. Правые сингулярные векторы здесь — просто $\mathbf e_2$ и $\mathbf e_1$ (первым идёт тот, что растягивается сильнее): окружность растягивается вдоль осей координат, а $U$ поворачивает результат. Проверка: $10 \cdot 5 = 50 = |3 \cdot 6 - (-8) \cdot 4|$.

Опыт 4. Фотография по слоям

Чёрно-белая фотография — это таблица яркостей, то есть матрица. Её сингулярное разложение можно переписать как сумму: столбец $\mathbf u_i$, умноженный на строку $\mathbf v_i^{\mathsf T}$, даёт матрицу того же размера, что и $A$, и

$$A = \sigma_1\mathbf u_1\mathbf v_1^{\mathsf T} + \sigma_2\mathbf u_2\mathbf v_2^{\mathsf T} + \dots + \sigma_r\mathbf u_r\mathbf v_r^{\mathsf T}.$$

Каждое слагаемое — картинка ранга $1$. В ней все строки пропорциональны одной строке $\mathbf v_i^{\mathsf T}$, а яркость пикселя равна произведению «веса строки» на «вес столбца». Так выглядит шотландка: вертикальные полосы, умноженные на горизонтальные. Фотография — сумма шотландок, и слагаемые идут по убыванию громкости $\sigma_i$.

Самый громкий слой. Для фотографии это обычно размытое «среднее» изображение: общая освещённость по строкам и столбцам. Последний оставленный слой. Всё, что после него, выбрасываем. Пример: чтобы хранить $A_k$ для фотографии $m \times n$, нужны $k$ столбцов $\mathbf u_i$ по $m$ чисел, $k$ строк $\mathbf v_i^{\mathsf T}$ по $n$ чисел и $k$ чисел $\sigma_i$, всего $k(m + n + 1)$. Для снимка $1000 \times 1000$ и $k = 50$ это $100\,050$ чисел вместо миллиона — около $10$ %.

Почему выбрасывать надо именно тихие слои, а не какие-то другие? Чтобы ответить, нужно договориться, как мерить ошибку приближения.

Хитрость в подсчёте размерностей: у матрицы малого ранга большое ядро, и оно обязательно пересекается с направлениями, которые $A$ растягивает сильнее всего. На чертеже случай $2 \times 2$ и $k = 1$; в общем случае рассуждение то же.

Сначала сам $A_k$. Разность $A - A_k = \sigma_{k+1}\mathbf u_{k+1}\mathbf v_{k+1}^{\mathsf T} + \dots$ уже записана в виде сингулярного разложения, и её наибольшее растяжение равно $\sigma_{k+1}$. На чертеже $A - A_1$ сплющивает окружность в отрезок длиной $2\sigma_2$ вдоль $\mathbf u_2$. Теперь любая $B$ ранга не больше $k$. По теореме о ранге и дефекте (глава 37) её ядро — векторы, которые $B$ отправляет в ноль, — имеет размерность не меньше $n - k$. На чертеже ядро — прямая; поворачивайте её, это и есть выбор $B$. Оболочка векторов $\mathbf v_1, \dots, \mathbf v_{k+1}$ имеет размерность $k + 1$, а $(n - k) + (k + 1) > n$. Сложим базисы ядра и оболочки: это не меньше $n + 1$ векторов в $\mathbb R^n$, и они линейно зависимы. Перенесём в зависимости слагаемые из базиса ядра налево, остальные направо: слева вектор ядра, справа такой же вектор оболочки. Он не нулевой, иначе обе части были бы нулевыми комбинациями независимых векторов и все коэффициенты оказались бы нулями. Возьмём этот общий вектор единичной длины и назовём $\mathbf x$. На плоскости $\mathbf v_1$ и $\mathbf v_2$ порождают всё, и $\mathbf x$ — просто единичный вектор ядра. Так как $B\mathbf x = \mathbf 0$, получаем $(A - B)\mathbf x = A\mathbf x$. Разложим $\mathbf x = c_1\mathbf v_1 + \dots + c_{k+1}\mathbf v_{k+1}$ с $c_1^2 + \dots + c_{k+1}^2 = 1$. Тогда $A\mathbf x = \sum c_i\sigma_i\mathbf u_i$ и $|A\mathbf x|^2 = \sum c_i^2\sigma_i^2 \ge \sigma_{k+1}^2$. На чертеже $A\mathbf x$ лежит на эллипсе, а эллипс не заходит внутрь окружности радиуса $\sigma_2$. Значит, $\|A - B\| \ge |(A - B)\mathbf x| \ge \sigma_{k+1} = \|A - A_k\|$: ни одна матрица ранга $k$ не приближает $A$ лучше, чем сумма первых $k$ слоёв.

А если мерить ошибку суммой квадратов всех пикселей

Для картинки естественнее складывать квадраты ошибок во всех пикселях: $\|M\|_F^2 = \sum_{i,j} m_{ij}^2$ (норма Фробениуса). Умножение на ортогональную матрицу слева не меняет длин столбцов, справа — длин строк, поэтому $\|M\|_F^2 = \|\Sigma\|_F^2 = \sigma_1^2(M) + \sigma_2^2(M) + \dots$. Для $A_k$ ошибка равна $\sigma_{k+1}^2 + \dots + \sigma_r^2$, и это тоже наилучший результат.

Нам понадобится оценка: у подпространств $U$ и $W$ пространства $\mathbb R^n$ пересечение имеет размерность не меньше $\dim U + \dim W - n$. Она следует из теоремы о ранге и дефекте, применённой к отображению пар $(\mathbf y, \mathbf z) \mapsto \mathbf y - \mathbf z$ из $U \times W$ (размерность $\dim U + \dim W$) в $\mathbb R^n$: его ядро состоит из пар $(\mathbf y, \mathbf y)$ с $\mathbf y \in U \cap W$, а образ не больше $\mathbb R^n$.

Докажем, что для любой $B$ ранга не больше $k$ и любого $i \ge 1$ выполнено $\sigma_{k+i}(A) \le \sigma_i(A - B)$. Возьмём три подпространства: оболочку $\mathbf v_1, \dots, \mathbf v_{k+i}$ (размерность $k + i$), ядро $B$ (не меньше $n - k$) и оболочку правых сингулярных векторов матрицы $A - B$ с номерами от $i$ до $n$ (размерность $n - i + 1$). По оценке пересечение первых двух имеет размерность не меньше $(k + i) + (n - k) - n = i$, а его пересечение с третьим — не меньше $i + (n - i + 1) - n = 1$. Возьмём в общем пересечении единичный вектор $\mathbf x$. Как в основном доказательстве, $|A\mathbf x| \ge \sigma_{k+i}(A)$. С другой стороны, $A\mathbf x = (A - B)\mathbf x$, а $\mathbf x$ раскладывается по сингулярным векторам $A - B$ с номерами от $i$, которые растягиваются не больше чем в $\sigma_i(A - B)$ раз, поэтому $|(A - B)\mathbf x| \le \sigma_i(A - B)$. Остаётся сложить квадраты: $\|A - B\|_F^2 = \sum_i \sigma_i^2(A - B) \ge \sum_i \sigma_{k+i}^2(A) = \|A - A_k\|_F^2$.

Теорему опубликовали в 1936 году физик Карл Эккарт и математик Гейл Янг, причём в журнале «Психометрика»: психологам нужно было описывать большие таблицы результатов тестов немногими общими «способностями». В компьютерном сжатии фотографий SVD почти не используют: хранить приходится и сами векторы $\mathbf u_i$, $\mathbf v_i$, свои для каждой картинки. Формат JPEG раскладывает блоки $8 \times 8$ пикселей по косинусам из главы о рядах Фурье, одинаковым для всех изображений. Но идея та же: оставить громкие слагаемые и выбросить тихие, а теорема Эккарта — Янга говорит, что для одной конкретной матрицы лучше SVD не сделать.

Картинка размером $400 \times 300$ пикселей. Сколько чисел нужно сохранить для её приближения ранга $20$?

По формуле $k(m + n + 1) = 20 \cdot (400 + 300 + 1) = 14\,020$. Это около $11{,}7$ % от $120\,000$ пикселей исходной картинки.

Опыт 5. Облако и его главная ось

Вернёмся к облаку из конца прошлой главы: рост и вес тысячи человек. Здесь неточны обе координаты, и ни одну из них не назначишь «правильной». Промахи естественно мерить по перпендикуляру к прямой. Такую задачу поставил Карл Пирсон в 1901 году в статье «О прямых и плоскостях, наиболее близких к системам точек в пространстве».

Сдвинем начало координат в центр масс облака, вычтя из каждой координаты её среднее, и запишем точки строками матрицы $X$ размера $n \times 2$ (или $n \times d$, если признаков $d$).

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

Пусть $\mathbf w$ — единичный вектор направления прямой, а $\mathbf x_1, \dots, \mathbf x_n$ — точки после сдвига в центр масс. Тень точки на прямую имеет длину $|\langle \mathbf x_i, \mathbf w\rangle|$, промах $d_i$ — расстояние от точки до прямой. По теореме Пифагора $|\mathbf x_i|^2 = \langle \mathbf x_i, \mathbf w\rangle^2 + d_i^2$. Складываем по всем точкам: $\sum |\mathbf x_i|^2 = \sum \langle \mathbf x_i, \mathbf w\rangle^2 + \sum d_i^2$. Левая часть от прямой не зависит — столбики на чертеже всегда одной высоты. Меньше промахов ровно тогда, когда больше сумма квадратов теней. Сумма квадратов теней — квадратичная форма: $\sum \langle \mathbf x_i, \mathbf w\rangle^2 = |X\mathbf w|^2 = \mathbf w^{\mathsf T}(X^{\mathsf T}X)\mathbf w$. Матрица $S = X^{\mathsf T}X$ симметрична, и по спектральной теореме у неё есть ортонормированные собственные векторы $\mathbf v_1, \mathbf v_2$ с $\lambda_1 \ge \lambda_2 \ge 0$. Разложим $\mathbf w = c_1\mathbf v_1 + c_2\mathbf v_2$, где $c_1^2 + c_2^2 = 1$. Тогда $\mathbf w^{\mathsf T}S\mathbf w = \lambda_1c_1^2 + \lambda_2c_2^2 \le \lambda_1(c_1^2 + c_2^2) = \lambda_1$, и равенство достигается при $\mathbf w = \pm\mathbf v_1$. В $d$ измерениях выкладка та же, только слагаемых $d$. Лучшая прямая идёт вдоль $\mathbf v_1$; сумма квадратов теней на ней равна $\lambda_1 = \sigma_1^2$, а сумма квадратов промахов — $\lambda_2 = \sigma_2^2$ (в $d$ измерениях — сумме остальных $\lambda_i$).

По облаку «рост — вес» провели прямую МНК, предсказывающую вес по росту, и вторую — предсказывающую рост по весу. Как они расположены?

Это одна и та же прямая Обе проходят через центр масс, а главная ось лежит между ними Они параллельны Они перпендикулярны

Каждая прямая МНК проходит через центр масс. Прямая «вес по росту» наклонена положе главной оси, а «рост по весу», если нарисовать её в тех же осях, — круче. Вертикальные промахи «прижимают» прямую к горизонтали, горизонтальные — к вертикали, а перпендикулярные никакой оси не предпочитают. Чем ближе облако к прямой, тем ближе все три.

Термин «главные компоненты» ввёл Гарольд Хотеллинг в 1933 году, и с тех пор метод стал первым, что делают с многомерными данными. Если признаков сто, а главных осей, которые держат почти весь разброс, три-четыре, данные можно рассматривать в трёх-четырёх измерениях. В 2008 году группа генетиков во главе с Джоном Новембре изучила около полумиллиона генетических маркеров в выборке из трёх тысяч европейцев и нарисовала людей точками по первым двум главным компонентам. Получилась узнаваемая карта Европы: соседи по карте оказались соседями и по генам.

И здесь есть где ошибиться. Главные оси зависят от единиц измерения: переведите рост из метров в сантиметры, и разброс по росту вырастет в $10\,000$ раз, а главная ось повернётся к нему. Поэтому признаки обычно сначала приводят к одному масштабу. А если облако почти круглое, $\lambda_1 \approx \lambda_2$, главная ось неустойчива: сдвиньте одну точку — и она развернётся. В виджете это видно сразу.

С сингулярным разложением метод связан напрямую. Если $X = U\Sigma V^{\mathsf T}$, то $X^{\mathsf T}X = V\Sigma^{\mathsf T}\Sigma V^{\mathsf T}$: главные оси — столбцы $V$, а разброс вдоль $i$-й оси равен $\sigma_i^2 / n$. Приближение $X_k$ по Эккарту — Янгу — это облако, спроектированное на первые $k$ главных осей.

Опыт 6. Кто что посмотрит

Последний опыт — с таблицей, в которой почти ничего не записано. В октябре 2006 года Netflix выложил $100\,480\,507$ оценок, которые $480\,189$ зрителей поставили $17\,770$ фильмам, и пообещал миллион долларов тому, кто предскажет оценки на $10$ % точнее, чем их собственная программа. Таблица «зрители × фильмы» заполнена чуть больше чем на $1$ %, и задача — угадать остальные клетки.

Помогает предположение о малом ранге. Пусть вкус зрителя описывают несколько чисел — сколько он любит комедии, ужасы, авторское кино, — и фильм описывают столько же чисел. Оценка тогда примерно равна скалярному произведению вектора вкуса на вектор фильма, а вся таблица — произведению двух узких матриц, то есть матрица малого ранга. Сингулярное разложение по пустой таблице не посчитать, но можно подобрать узкие матрицы по методу наименьших квадратов только по известным клеткам. Зафиксируем векторы фильмов — вектор каждого зрителя находится обычной задачей МНК из опыта 2. Зафиксируем зрителей — так же находятся фильмы. Чередуя, быстро приходим к хорошему приближению.

В декабре 2006 года программист Брэндин Уэбб под псевдонимом Саймон Фанк описал в блоге такую модель: узкие множители подбираются только по известным оценкам (он подбирал их градиентным спуском, а не чередованием). Он назвал её SVD, и метод быстро разошёлся среди участников. Приз вручили в сентябре 2009 года команде BellKor's Pragmatic Chaos, которая улучшила результат на $10{,}06$ %. Команда The Ensemble добилась того же результата, но отправила решение на $20$ минут позже. Победная программа смешивала больше сотни моделей, и разложение матрицы на узкие множители было среди главных.

Куда дальше

Ортогональные матрицы $2 \times 2$ оказались поворотами и отражениями. Возьмите из них те восемь, что переводят в себя квадрат с центром в начале координат: четыре поворота (на $0^\circ$, $90^\circ$, $180^\circ$, $270^\circ$) и четыре отражения. Перемножьте любые две — снова получится одна из восьми. У этих матриц своя таблица умножения $8 \times 8$, в которой не встречается ничего постороннего. Похожие таблицы есть у поворотов граней кубика Рубика, у перетасовок колоды и у сложения часов на циферблате, хотя ни чисел, ни матриц там нет. Что у всех них общего и что из этого общего следует, расскажет глава 40.