Всем привет! Сразу признаюсь, что я немного обманываю вас в заголовке: взглядов будет не пятьдесят, а всего пять. Но, надеюсь, каждый из них окажется достаточно полезным, чтобы по-новому посмотреть на один и тот же метод.

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

PCA (Principal Component Analysis, метод главных компонент) - как раз такой случай. Это один из старейших и наиболее известных методов анализа данных. Его раннюю формулировку предложил Карл Пирсон в 1901 году, а позднее метод был развит Гарольдом Хотеллингом.

Статья вдохновлена курсом Geometrical Methods of Machine Learning, который читает в Сколтехе Александр Бернштейн.

Немного теории

Здесь распишу достаточно большой блок по всему математическому аппарату, который понадобится для понимания статьи. Те, кто еще ~~не успел забыть~~ помнит курсы вышмата, можно смело пропускать этот раздел, а при необходимости возвращаться к нему по ходу чтения.

Собственные значения и собственные векторы

Начнем с самого базового аппарата линейной алгебры, на который в итоге сведется практически все остальное в этой статье.

Пусть $A\in\mathbb{R}^{p\times p}$ - квадратная матрица. Ненулевой вектор называется собственным вектором матрицы , если существует число $\lambda$ такое, что

$Ae=\lambda e.$

Число $\lambda$ называется соответствующим собственным значением.

Геометрически это означает, что действие матрицы на собственный вектор не меняет его направление: вектор только растягивается, сжимается или, в общем случае, меняет ориентацию при отрицательном собственном значении.

Нас будут интересовать симметричные матрицы, $A^{\top}=A$ . По спектральной теореме такая матрица всегда допускает ортонормированный базис из собственных векторов, то есть допускает разложение

$A=E\Lambda E^{\top},$

где

$E= \begin{pmatrix} e_1 & e_2 & \cdots & e_p \end{pmatrix}, \qquad E^{\top}E=I_p,$

а

$\Lambda= \operatorname{diag}(\lambda_1,\lambda_2,\ldots,\lambda_p).$

Договоримся всегда упорядочивать собственные значения по убыванию:

$\lambda_1\geq \lambda_2\geq \cdots\geq \lambda_p.$

Ниже в роли почти везде будет выступать ковариационная матрица $\Sigma$ - мы определим ее чуть позже. Забегая вперед: она не только симметрична, но и положительно полуопределена, поэтому все её собственные значения дополнительно неотрицательны:

$\lambda_1\geq \lambda_2\geq \cdots\geq \lambda_p\geq 0.$

Ведущее собственное подпространство

Небольшая напоминалка на будущее, поскольку словосочетание «ведущее собственное подпространство» размерности будет несколько раз всплывать по ходу статьи как решение самых разных оптимизационных задач.

Ведущим собственным подпространством размерности называется

$\operatorname{span}(e_1,\ldots,e_q),$

то есть подпространство, натянутое на собственные векторы, отвечающие наибольшим собственным значениям $\lambda_1\geq\cdots\geq\lambda_q$ .

Слово «ведущее» здесь важно: у матрицы есть много разных -мерных подпространств, натянутых на всевозможные наборы из собственных векторов (например, $\operatorname{span}(e_2,e_5,e_7)$ - такое же законное инвариантное подпространство). Но именно набор из первых векторов, отвечающих самым большим собственным значениям, максимизирует сумму $\sum_{k}\lambda_k$ по выбранным направлениям - и именно поэтому оно и оказывается решением всех вариационных задач ниже.

Также несколько нюансов, которые стоит иметь в виду:

если $\lambda_q>\lambda_{q+1}$ (есть спектральный зазор), подпространство определено однозначно как множество - сами базисные векторы внутри него определены лишь с точностью до знака;
если $\lambda_q=\lambda_{q+1}$ , однозначного разделения между «взять» и «не взять» уже нет: граница проходит внутри общего собственного подпространства, отвечающего повторяющемуся собственному значению, и конкретный ортонормированный базис в нем можно выбирать по-разному.

След матрицы и норма Фробениуса

След квадратной матрицы - сумма ее диагональных элементов:

$\operatorname{Tr}(A)=\sum_i A_{ii}.$

Полезные свойства следа:

$\operatorname{Tr}(AB)=\operatorname{Tr}(BA),$ $\operatorname{Tr}(ABC)=\operatorname{Tr}(BCA)=\operatorname{Tr}(CAB).$

Квадрат нормы Фробениуса равен

$\lVert A\rVert_F^2 = \sum_{i,j}A_{ij}^2 = \operatorname{Tr}(A^{\top}A).$

Ортогональные проекторы

Пусть столбцы матрицы $E_q\in\mathbb{R}^{p\times q}$ ортонормированы, то есть $E_q^{\top}E_q=I_q$ . Тогда

$P_q=E_qE_q^{\top}$

является матрицей ортогональной проекции на подпространство, натянутое на столбцы .

Для ортогонального проектора выполняются свойства

$P_q^{\top}=P_q,$

Вектор раскладывается на ортогональные компоненты:

поэтому по теореме Пифагора

$\lVert x\rVert^2 = \lVert P_qx\rVert^2 + \lVert(I-P_q)x\rVert^2.$

Эта формула будет часто появляться в геометрических интерпретациях PCA.

Аффинное преобразование и аффинное подпространство

~~Так давайте же наконец ответим на вопрос, что такое афинное преобразование~~

Отображение $f:\mathbb{R}^p\to\mathbb{R}^q$ называется аффинным, если его можно записать в виде

где $A\in\mathbb{R}^{q\times p}$ - линейное отображение (матрица), а $b\in\mathbb{R}^q$ - фиксированный вектор сдвига. Другими словами, аффинное отображение - это линейное отображение плюс постоянный сдвиг ; при оно превращается в обычное линейное отображение.

Аналогично, аффинным подпространством размерности в $\mathbb{R}^p$ называется множество вида

$L=x_0+V=\{x_0+v:v\in V\},$

где $V\subset\mathbb{R}^p$ - линейное подпространство размерности , а $x_0\in\mathbb{R}^p$ - произвольная точка сдвига. В отличие от линейного подпространства, аффинное подпространство не обязано проходить через начало координат - это просто линейное подпространство , сдвинутое на вектор .

Если $V=\operatorname{span}(e_1,\ldots,e_q)$ задано ортонормированным базисом $E=(e_1\;\cdots\;e_q)$ , $E^{\top}E=I_q$ , то точки удобно параметризовать как

$L=\{x_0+Et:t\in\mathbb{R}^q\}.$

Данные и центрирование

Пусть у нас есть выборка из объектов:

$X_1, X_2, \ldots, X_n \in \mathbb{R}^p.$

Каждый объект описывается числовыми признаками. Например, если объект - изображение размером $32 \times 32$ пикселя, то после вытягивания изображения в вектор получаем точку в пространстве размерности .

Выборочное среднее равно

$\overline{X}=\frac{1}{n}\sum_{i=1}^{n}X_i.$

Центрированными наблюдениями назовем векторы

$\widetilde{X}_i=X_i-\overline{X}.$

После центрирования среднее выборки становится нулевым:

$\frac{1}{n}\sum_{i=1}^{n}\widetilde{X}_i=0.$

Центрирование принципиально важно. Обычный PCA ищет линейные направления для центрированных данных, а в исходном пространстве этим направлениям соответствует аффинное подпространство, проходящее через среднее выборки.

Соберем центрированные наблюдения в матрицу по столбцам:

$\widetilde{X}= \begin{pmatrix} \widetilde{X}_1 & \widetilde{X}_2 & \cdots & \widetilde{X}_n \end{pmatrix} \in \mathbb{R}^{p\times n}.$

Иногда в машинном обучении объекты записывают по строкам. Тогда матрица данных имеет размер $n\times p$ . Обе конвенции эквивалентны, но формулы транспонируются. В этой статье объекты будут записаны по столбцам.

Ковариационная матрица

Выборочная ковариационная матрица в принятой здесь нормировке имеет вид

$\Sigma = \frac{1}{n}\sum_{i=1}^{n} \widetilde{X}_i\widetilde{X}_i^{\top} = \frac{1}{n}\widetilde{X}\widetilde{X}^{\top}.$

В несмещенном случае вместо имеем :

$\widehat{\Sigma}_{\mathrm{unbiased}} = \frac{1}{n-1}\widetilde{X}\widetilde{X}^{\top}.$

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

Матрица $\Sigma$ симметрична:

$\Sigma^{\top}=\Sigma,$

и положительно полуопределена, поскольку для любого вектора $a\in\mathbb{R}^p$

$a^{\top}\Sigma a = \frac{1}{n}\sum_{i=1}^{n} \left(a^{\top}\widetilde{X}_i\right)^2 \geq 0.$

Отсюда следует, что все собственные значения $\Sigma$ вещественны и неотрицательны - это в точности тот случай матрицы , о котором шла речь в самом первом разделе. Значит, для $\Sigma$ верно спектральное разложение

$\Sigma=E\Lambda E^{\top}, \qquad \lambda_1\geq\cdots\geq\lambda_p\geq 0,$

а первые собственных векторов образуют матрицу

$E_q= \begin{pmatrix} e_1 & e_2 & \cdots & e_q \end{pmatrix} \in\mathbb{R}^{p\times q}.$

Векторы $e_1,\ldots,e_p$ задают главные оси данных. Координаты объекта в пространстве главных компонент равны

$y_i=E_q^{\top}(X_i-\overline{X}) \in\mathbb{R}^q,$

а обратное приближенное восстановление имеет вид

$\widehat{X}_i = \overline{X}+E_qy_i = \overline{X}+E_qE_q^{\top}(X_i-\overline{X}).$

Наконец, полная дисперсия центрированных данных равна следу ковариационной матрицы:

$\frac{1}{n}\sum_{i=1}^{n} \lVert X_i-\overline{X}\rVert^2 = \operatorname{Tr}(\Sigma) = \sum_{k=1}^{p}\lambda_k.$

Квадратичная форма и отношение Рэлея

Для единичного вектора величина

$e^{\top}\Sigma e$

называется квадратичной формой ковариационной матрицы. В контексте PCA она равна дисперсии данных после проекции на направление .

Отношение Рэлея определяется как

$R_{\Sigma}(e) = \frac{e^{\top}\Sigma e}{e^{\top}e}.$

Для симметричной матрицы оно удовлетворяет неравенству

$\lambda_p \leq R_{\Sigma}(e) \leq \lambda_1.$

Максимум достигается на собственном векторе , а минимум - на . Именно это наблюдение лежит в основе формулировки PCA через максимизацию дисперсии.

Whitening-преобразование

Whitening (отбеливание) - это линейное преобразование, которое переводит центрированные данные в новые координаты с единичной ковариационной матрицей, то есть выравнивает дисперсию по всем направлениям.

Пусть $\Sigma=E\Lambda E^{\top}$ - введенное выше спектральное разложение, и все

$\lambda_k>0$ (в общем случае можно ограничиться первыми направлениями с положительными собственными значениями и работать с $E_q,\Lambda_q$ ). Определим

$\Lambda^{-1/2}=\operatorname{diag}\!\left(\lambda_1^{-1/2},\ldots,\lambda_p^{-1/2}\right)$

и преобразование

$Z=\Lambda^{-1/2}E^{\top}(X-\overline{X}).$

Несложно проверить, что $\operatorname{Cov}(Z)=I_p$ . Действительно,

$\operatorname{Cov}(Z) = \Lambda^{-1/2}E^{\top}\Sigma E\Lambda^{-1/2} = \Lambda^{-1/2}E^{\top}E\Lambda E^{\top}E\Lambda^{-1/2} = \Lambda^{-1/2}\Lambda\Lambda^{-1/2} = I_p,$

где мы воспользовались тем, что $E^{\top}E=I_p$ , а диагональная матрица $\Lambda$ коммутирует с $\Lambda^{-1/2}$ .

Whitening удобно представлять как композицию двух шагов:

поворот $E^{\top}(X-\overline{X})$ - это в точности переход в координаты главных компонент, где ковариация уже диагональна и равна $\Lambda$ ;
покомпонентное масштабирование $\Lambda^{-1/2}(\cdot)$ , которое растягивает каждую координату так, чтобы ее дисперсия стала равна единице.

Whitening легко комбинируется с понижением размерности - если оставить только первые компонент,

$Z=\Lambda_q^{-1/2}E_q^{\top}(X-\overline{X})\in\mathbb{R}^q,$

то $\operatorname{Cov}(Z)=I_q$ по тем же самым вычислениям, только с урезанными и $\Lambda_q$ .

PCA как максимизация дисперсии и декорреляция

Начнем с наиболее распространенного взгляда: PCA как метод поиска направлений, вдоль которых данные имеют максимальную дисперсию - многие видели известную иллюстрацию с эллипсами.

Первая главная компонента

Рассмотрим центрированный случай. Спроецируем случайный вектор на единичное направление :

$Y_e=e^{\top}X, \qquad e^{\top}e=1.$

Дисперсия проекции равна

$\operatorname{Var}(Y_e) = \operatorname{Var}(e^{\top}X) = e^{\top}\Sigma e.$

Первая главная компонента определяется как решение задачи

$\max_{e\in\mathbb{R}^p} e^{\top}\Sigma e \quad \text{при условии} \quad e^{\top}e=1.$

Используем метод множителей Лагранжа:

$\mathcal{L}(e,\mu) = e^{\top}\Sigma e - \mu(e^{\top}e-1).$

Условие стационарности по дает

$\nabla_e\mathcal{L} = 2\Sigma e-2\mu e=0.$

Следовательно,

$\Sigma e=\mu e.$

То есть оптимальное направление должно быть собственным вектором ковариационной матрицы. Для собственного вектора

$e^{\top}\Sigma e = e^{\top}(\lambda e) = \lambda e^{\top}e = \lambda.$

Чтобы получить максимальную дисперсию, нужно выбрать наибольшее собственное значение. Значит,

$e_{\text{PC1}}=e_1,$ $\operatorname{Var}(e_1^{\top}X)=\lambda_1.$

Последующие компоненты

Вторая главная компонента должна иметь максимально возможную дисперсию среди направлений, ортогональных первой:

$\max_e e^{\top}\Sigma e$

при ограничениях

$e^{\top}e=1, \qquad e^{\top}e_1=0.$

Решением будет , соответствующий второму по величине собственному значению $\lambda_2$ . Аналогично, для компоненты с номером добавляются ограничения ортогональности ко всем предыдущим направлениям:

$e_k^{\top}e_j=0, \qquad j<k.$

В результате получаем последовательность

$e_1,e_2,\ldots,e_q,$

где дисперсия компоненты с номером равна

$\operatorname{Var}(e_k^{\top}X)=\lambda_k.$

Одновременная оптимизация нескольких направлений

Вместо последовательного поиска можно сразу искать матрицу с ортонормированными столбцами:

$E_q^{\top}E_q=I_q.$

Суммарная дисперсия проекции равна

$\operatorname{Tr}(E_q^{\top}\Sigma E_q).$

Задача принимает вид

$\max_{E_q^{\top}E_q=I_q} \operatorname{Tr}(E_q^{\top}\Sigma E_q).$

По принципу Куранта - Фишера или теореме Кая Фана максимум равен сумме первых собственных значений:

$\max_{E_q^{\top}E_q=I_q} \operatorname{Tr}(E_q^{\top}\Sigma E_q) = \sum_{k=1}^{q}\lambda_k.$

Он достигается, когда столбцы образуют базис ведущего собственного подпространства ковариационной матрицы.

Про сам принцип Куранта - Фишера, теорему Кая Фана и связанное с ними неравенство Вейля можно почитать здесь. А если хочется углубиться в тему собственных значений эрмитовых матриц еще подробнее - есть отличный блог самого Теренса Тао: Eigenvalues and sums of Hermitian matrices.

Если собственные значения различны, отдельные направления определены с точностью до знака. Если некоторые собственные значения совпадают, то единственным является соответствующее собственное подпространство, а ортонормированный базис внутри него можно выбирать по-разному.

Почему компоненты декоррелированы

Запишем все главные компоненты в виде

$Y=E^{\top}X.$

Их ковариационная матрица равна

$\operatorname{Cov}(Y) = E^{\top}\Sigma E.$

Используя спектральное разложение $\Sigma=E\Lambda E^{\top}$ , получаем

$\operatorname{Cov}(Y) = E^{\top}E\Lambda E^{\top}E = \Lambda.$

Матрица $\Lambda$ диагональна, поэтому разные главные компоненты некоррелированы:

$\operatorname{Cov}(Y_j,Y_k)=0, \qquad j\neq k.$

При этом

$\operatorname{Var}(Y_k)=\lambda_k.$

Важно: некоррелированность не означает независимость в общем случае. Для совместно гауссовского распределения некоррелированные компоненты действительно независимы, но для произвольных данных зависимость более высокого порядка может сохраняться.

Доля объясненной дисперсии

Доля дисперсии, объясняемая компонентой , равна

$\mathrm{EVR}_k = \frac{\lambda_k}{\sum_{j=1}^{p}\lambda_j}.$

Суммарная доля дисперсии, сохраняемая первыми компонентами:

$\mathrm{CEVR}(q) = \frac{\sum_{k=1}^{q}\lambda_k} {\sum_{k=1}^{p}\lambda_k}.$

Например, можно выбрать минимальное , для которого

$\mathrm{CEVR}(q)\geq 0.95.$

PCA и whitening

PCA сам по себе только поворачивает координаты и, при необходимости, отбрасывает часть направлений. Дисперсии сохраненных компонент остаются равными $\lambda_k$ .

Чтобы дополнительно привести дисперсию каждой компоненты к единице, применяют whitening:

$Z=\Lambda_q^{-1/2}E_q^{\top}(X-\overline{X}).$

Тогда

$\operatorname{Cov}(Z)=I_q.$

Whitening и PCA - связанные, но не тождественные преобразования. Whitening включает дополнительное масштабирование и может сильно усиливать шум в направлениях с малыми собственными значениями.

PCA как частный случай SVD

Второй взгляд на PCA - с точки зрения вычислительной линейной алгебры. Вместо явного построения ковариационной матрицы можно применить сингулярное разложение непосредственно к центрированной матрице данных.

Что такое SVD

Для любой матрицы

$A\in\mathbb{R}^{m\times n}$

существует сингулярное разложение

$A=U\Gamma V^{\top},$

где столбцы и ортонормированы, а $\Gamma$ - диагональная или прямоугольная диагональная матрица с неотрицательными сингулярными числами:

$\sigma_1\geq\sigma_2\geq\cdots\geq 0.$

Для компактного SVD матрицы ранга можно записать

$A=U_r\Gamma_rV_r^{\top},$

где

$U_r\in\mathbb{R}^{m\times r}, \qquad \Gamma_r\in\mathbb{R}^{r\times r}, \qquad V_r\in\mathbb{R}^{n\times r}.$

SVD связано со спектральными разложениями матриц $AA^{\top}$ и $A^{\top}A$ :

$AA^{\top}=U\Gamma^2U^{\top},$ $A^{\top}A=V\Gamma^2V^{\top}.$

Левые сингулярные векторы являются собственными векторами $AA^{\top}$ , правые сингулярные векторы - собственными векторами $A^{\top}A$ , а квадраты сингулярных чисел являются их ненулевыми собственными значениями.

Связь SVD и ковариационной матрицы

Применим SVD к центрированной матрице данных:

$\widetilde{X}=U\Gamma V^{\top}.$

Ковариационная матрица равна

$\Sigma = \frac{1}{n}\widetilde{X}\widetilde{X}^{\top}.$

Подставляя SVD, получаем

$\Sigma = \frac{1}{n} U\Gamma V^{\top}V\Gamma U^{\top} = U\frac{\Gamma^2}{n}U^{\top}.$

Следовательно, столбцы - собственные векторы ковариационной матрицы, а собственные значения равны

$\lambda_k=\frac{\sigma_k^2}{n}.$

При нормировке ковариации на формула меняется на

$\lambda_k=\frac{\sigma_k^2}{n-1}.$

Таким образом, PCA можно вычислить без отдельного спектрального разложения $\Sigma$ : достаточно выполнить SVD центрированной матрицы данных.

Координаты объектов в PCA-пространстве

Матрица координат всех объектов в пространстве главных компонент равна

$Y=U_q^{\top}\widetilde{X}.$

Подставим усеченное SVD:

$\widetilde{X}=U\Gamma V^{\top}.$

Тогда

$Y = U_q^{\top}U\Gamma V^{\top} = \Gamma_qV_q^{\top}.$

То есть координаты объектов можно получить двумя эквивалентными способами:

$Y=U_q^{\top}\widetilde{X},$

или

$Y=\Gamma_qV_q^{\top}.$

При записи объектов по строкам соответствующая привычная формула имеет вид

Усеченное SVD и приближение матрицы данных

Оставим только первые сингулярных компонент:

$\widetilde{X}_q=U_q\Gamma_qV_q^{\top}.$

Теорема Эккарта - Янга - Мирского утверждает, что это наилучшее приближение матрицы $\widetilde{X}$ матрицей ранга не выше в норме Фробениуса:

$\widetilde{X}_q = \underset{\operatorname{rank}(A)\leq q}{\operatorname{argmin}} \lVert\widetilde{X}-A\rVert_F.$

Доказательство этого факта можно почитать, например, здесь.

Минимальная ошибка равна

$\lVert\widetilde{X}-\widetilde{X}_q\rVert_F^2 = \sum_{k=q+1}^{r}\sigma_k^2.$

После деления на получаем ошибку PCA-реконструкции:

$\frac{1}{n} \lVert\widetilde{X}-\widetilde{X}_q\rVert_F^2 = \sum_{k=q+1}^{r}\lambda_k.$

Эта формула уже связывает SVD с геометрическим приближением данных подпространством.

Почему PCA обычно считают через SVD

Наивный алгоритм мог бы сначала построить

$\Sigma=\frac{1}{n}\widetilde{X}\widetilde{X}^{\top},$

а затем найти ее собственные значения и собственные векторы. Но явное вычисление $\widetilde{X}\widetilde{X}^{\top}$ может ухудшать численную устойчивость, поскольку число обусловленности фактически возводится в квадрат.

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

Выбор вычислительной схемы зависит от размеров задачи:

если число признаков умеренно, удобно работать с $p\times p$ ковариационной матрицей;
если $p\gg n$ , ненулевых компонент не больше после центрирования, и иногда выгоднее использовать матрицу Грама $\widetilde{X}^{\top}\widetilde{X}$ размера $n\times n$ ;
для разреженных и очень больших данных применяют truncated SVD или итерационные методы, не формируя плотные матрицы.

Нужно учитывать, что TruncatedSVD в некоторых библиотеках по умолчанию не центрирует данные. Поэтому он совпадает с PCA только после корректного центрирования или в специально оговоренных задачах.

PCA как классическое метрическое многомерное шкалирование

Теперь посмотрим на PCA через попарные расстояния между объектами.

Задача многомерного шкалирования

Пусть исходные точки находятся в пространстве высокой размерности:

$X_1,\ldots,X_n\in\mathbb{R}^p.$

Мы хотим найти точки

$y_1,\ldots,y_n\in\mathbb{R}^q, \qquad q<p,$

так, чтобы расстояния между ними были максимально похожи на исходные.

Обозначим через матрицу квадратов исходных расстояний:

$(D_X)_{ij}=\lVert X_i-X_j\rVert^2.$

Аналогично,

$(D_Y)_{ij}=\lVert y_i-y_j\rVert^2.$

В общей постановке метрического MDS можно минимизировать некоторую функцию расхождения расстояний $\operatorname{Disc}(Y)$ . Самый известный вариант такой функции - критерий stress, предложенный Крускалом:

$\operatorname{Disc}(Y) = \sum_{i<j} \left( \lVert X_i-X_j\rVert - \lVert y_i-y_j\rVert \right)^2.$

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

Разница здесь важна: минимизация stress по расстояниям и минимизация strain по матрице Грама - в общем случае разные оптимизационные задачи с разными решениями. Эквивалентность с PCA относится именно к классическому MDS, то есть к минимизации strain для евклидовых расстояний.

От расстояний к скалярным произведениям

Для центрированных точек определим матрицу Грама

$S_X=\widetilde{X}^{\top}\widetilde{X} \in\mathbb{R}^{n\times n}.$

Ее элементы равны попарным скалярным произведениям:

$(S_X)_{ij} = \langle X_i-\overline{X},X_j-\overline{X}\rangle.$

Квадрат расстояния выражается через матрицу Грама:

$\lVert X_i-X_j\rVert^2 = (S_X)_{ii}+(S_X)_{jj}-2(S_X)_{ij}.$

Введем матрицу центрирования

$H=I_n-\frac{1}{n}\mathbf{1}\mathbf{1}^{\top}.$

Несложно заметить, что она симметрична и идемпотентна:

$H^{\top}=H,$

Симметричность следует прямо из симметричности $\mathbf{1}\mathbf{1}^{\top}$ . Идемпотентность проверяется коротким вычислением: поскольку $\mathbf{1}^{\top}\mathbf{1}=n$ ,

$H^2 = \left(I_n-\frac{1}{n}\mathbf{1}\mathbf{1}^{\top}\right)^2 = I_n-\frac{2}{n}\mathbf{1}\mathbf{1}^{\top} + \frac{1}{n^2}\mathbf{1}\left(\mathbf{1}^{\top}\mathbf{1}\right)\mathbf{1}^{\top} = I_n-\frac{2}{n}\mathbf{1}\mathbf{1}^{\top}+\frac{1}{n}\mathbf{1}\mathbf{1}^{\top} = H.$

Двойное центрирование матрицы квадратов расстояний дает матрицу Грама:

$S_X=-\frac{1}{2}HD_XH.$

Таким образом, если известны только попарные евклидовы расстояния, мы можем восстановить центрированные скалярные произведения, не зная исходных координат.

Спектральное решение классического MDS

Классическое MDS ищет координаты $Y\in\mathbb{R}^{q\times n}$ , для которых

$Y^{\top}Y$

как можно лучше приближает исходную матрицу Грама:

$\min_{Y\in\mathbb{R}^{q\times n}} \lVert S_X-Y^{\top}Y\rVert_F^2.$

Матрица $Y^{\top}Y$ имеет ранг не выше и положительно полуопределена. Поэтому решение получается усеченным спектральным разложением .

Пусть

$S_X=V\Gamma^2V^{\top}.$

Тогда оптимальные координаты можно выбрать как

$Y=\Gamma_qV_q^{\top}.$

Но из SVD центрированной матрицы данных

$\widetilde{X}=U\Gamma V^{\top}$

следует

$\widetilde{X}^{\top}\widetilde{X} = V\Gamma^2V^{\top}.$

Следовательно, координаты классического MDS равны

$Y_{\mathrm{MDS}}=\Gamma_qV_q^{\top}.$

Ранее для PCA мы получили

$Y_{\mathrm{PCA}}=U_q^{\top}\widetilde{X}=\Gamma_qV_q^{\top}.$

Значит,

$Y_{\mathrm{MDS}}=Y_{\mathrm{PCA}}$

с точностью до ортогонального преобразования внутри пространства одинаковых собственных значений и до соглашений о знаках компонент.

Геометрический смысл

PCA обычно формулируется через признаки и ковариацию, а классическое MDS - через расстояния между объектами. Но для евклидовых данных эти описания содержат одну и ту же информацию после центрирования:

$\text{квадраты расстояний} \quad\longleftrightarrow\quad \text{центрированная матрица Грама}.$

Поэтому PCA и классическое MDS приходят к одному спектральному разложению.

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

PCA как наилучшее линейное приближение

Теперь рассмотрим чисто геометрическую задачу: какое -мерное аффинное подпространство проходит ближе всего к облаку данных?

Аффинное подпространство

Линейное подпространство обязательно проходит через начало координат. Аффинное подпространство может быть сдвинуто и задается как

$L^{(q)} = \left\{ X_0+Et \;:\; t\in\mathbb{R}^q \right\},$

где

$E= \begin{pmatrix} e_1 & \cdots & e_q \end{pmatrix}, \qquad E^{\top}E=I_q.$

Точка задает сдвиг, а столбцы - ортонормированный базис направляющего линейного подпространства.

Ортогональная проекция точки на $L^{(q)}$ равна

$\operatorname{Pr}_{L^{(q)}}(X) = X_0+EE^{\top}(X-X_0).$

Или покомпонентно:

$\operatorname{Pr}_{L^{(q)}}(X) = X_0+ \sum_{k=1}^{q} \langle X-X_0,e_k\rangle e_k.$

Наша цель - минимизировать средний квадрат ортогонального расстояния от данных до подпространства:

$J(X_0,E) = \frac{1}{n} \sum_{i=1}^{n} \left\lVert X_i- \operatorname{Pr}_{L^{(q)}}(X_i) \right\rVert^2.$

Почему оптимальная плоскость проходит через среднее

Обозначим ортогональный проектор

$P=EE^{\top}.$

Остаток после проекции равен

$X_i- \operatorname{Pr}_{L^{(q)}}(X_i) = (I-P)(X_i-X_0).$

Разложим

$X_i-X_0 = (X_i-\overline{X})+(\overline{X}-X_0).$

Тогда средняя ошибка равна

$J(X_0,E) = \frac{1}{n} \sum_{i=1}^{n} \left\lVert (I-P) \left[ (X_i-\overline{X})+(\overline{X}-X_0) \right] \right\rVert^2.$

После раскрытия квадрата смешанный член исчезает, потому что

$\sum_{i=1}^{n}(X_i-\overline{X})=0.$

Получаем

$J(X_0,E) = \frac{1}{n} \sum_{i=1}^{n} \lVert(I-P)(X_i-\overline{X})\rVert^2 + \lVert(I-P)(\overline{X}-X_0)\rVert^2.$

Второе слагаемое неотрицательно и обращается в ноль, если

$(I-P)(\overline{X}-X_0)=0.$

Это означает, что среднее должно лежать в аффинном подпространстве. В частности, можно выбрать

$X_0=\overline{X}.$

После этого задача сводится к поиску линейного подпространства для центрированных данных.

От минимизации расстояния к максимизации проекции

При $X_0=\overline{X}$ имеем

$\operatorname{Pr}_{L^{(q)}}(X_i) = \overline{X}+P(X_i-\overline{X}).$

По теореме Пифагора

$\lVert X_i-\overline{X}\rVert^2 = \lVert P(X_i-\overline{X})\rVert^2 + \lVert(I-P)(X_i-\overline{X})\rVert^2.$

Следовательно,

$J(E) = \frac{1}{n} \sum_{i=1}^{n} \lVert X_i-\overline{X}\rVert^2 - \frac{1}{n} \sum_{i=1}^{n} \lVert P(X_i-\overline{X})\rVert^2.$

Первое слагаемое не зависит от . Поэтому минимизация ортогональной ошибки эквивалентна максимизации энергии проекции:

$\max_{E^{\top}E=I_q} \frac{1}{n} \sum_{i=1}^{n} \lVert EE^{\top}(X_i-\overline{X})\rVert^2.$

Так как

$\lVert E^{\top}(X_i-\overline{X})\rVert^2 = \sum_{k=1}^{q} \langle X_i-\overline{X},e_k\rangle^2,$

получаем

$\frac{1}{n} \sum_{i=1}^{n} \lVert E^{\top}(X_i-\overline{X})\rVert^2 = \operatorname{Tr}(E^{\top}\Sigma E).$

Значит, задача имеет вид

$\max_{E^{\top}E=I_q} \operatorname{Tr}(E^{\top}\Sigma E),$

и ее решением снова являются первые собственных векторов $\Sigma$ .

Минимальная средняя ошибка равна сумме отброшенных собственных значений:

$J_{\min} = \sum_{k=q+1}^{p}\lambda_k.$

PCA-линия и линейная регрессия - не одно и то же

Для двумерных данных первая главная компонента задает прямую, минимизирующую сумму квадратов ортогональных расстояний от точек до прямой.

Обычная линейная регрессия

минимизирует вертикальные остатки:

$\sum_{i=1}^{n} \left(y_i-(ax_i+b)\right)^2.$

PCA симметрично относится ко всем координатам, а стандартная регрессия явно разделяет зависимую и независимую переменные. Поэтому наклон PCA-прямой и регрессионной прямой в общем случае различается.

Почему здесь возникает именно аффинное преобразование

Отображение из исходного пространства в пространство главных компонент имеет вид

$f(X)=E_q^{\top}(X-\overline{X}).$

Это не чисто линейное отображение по переменной из-за вычитания среднего. Его можно переписать как

$f(X)=E_q^{\top}X-E_q^{\top}\overline{X},$

то есть это аффинное отображение.

Восстановление также аффинно:

$g(y)=\overline{X}+E_qy.$

Композиция дает

$g(f(X)) = \overline{X}+E_qE_q^{\top}(X-\overline{X}),$

то есть ортогональную проекцию на аффинное подпространство

$\overline{X}+\operatorname{span}(e_1,\ldots,e_q).$

Именно поэтому корректнее говорить, что PCA ищет наилучшее аффинное подпространство в исходных координатах и наилучшее линейное подпространство после центрирования.

PCA как наилучший метод снижения размерности

Пятый взгляд формулирует PCA как задачу кодирования и восстановления данных.

Пусть мы хотим сжать центрированный вектор из $\mathbb{R}^p$ в вектор из $\mathbb{R}^q$ , а затем восстановить исходный объект.

Линейный энкодер и декодер

Линейный энкодер задается матрицей

$W\in\mathbb{R}^{q\times p}:$ $y_i=W(X_i-\overline{X}).$

Линейный декодер задается матрицей

$V\in\mathbb{R}^{p\times q}:$ $\widehat{X}_i = \overline{X}+Vy_i.$

Композиция имеет вид

$\widehat{X}_i = \overline{X}+VW(X_i-\overline{X}).$

Средняя ошибка восстановления:

$\mathcal{E}(V,W) = \frac{1}{n} \sum_{i=1}^{n} \left\lVert X_i- \overline{X}-VW(X_i-\overline{X}) \right\rVert^2.$

Ортогональный вариант

Сначала рассмотрим естественное ограничение

$WW^{\top}=I_q$

и связанный декодер

$V=W^{\top}.$

Тогда

$P=W^{\top}W$

является ортогональным проектором на пространство строк , а восстановление равно

$\widehat{X}_i = \overline{X}+P(X_i-\overline{X}).$

Ошибка принимает вид

$\mathcal{E}(W) = \frac{1}{n} \sum_{i=1}^{n} \lVert(I-P)(X_i-\overline{X})\rVert^2.$

По теореме Пифагора

$\mathcal{E}(W) = \frac{1}{n} \sum_{i=1}^{n} \lVert X_i-\overline{X}\rVert^2 - \frac{1}{n} \sum_{i=1}^{n} \lVert W(X_i-\overline{X})\rVert^2.$

Второй член равен

$\operatorname{Tr}(W\Sigma W^{\top}).$

Поэтому

$\min_{WW^{\top}=I_q}\mathcal{E}(W)$

эквивалентно

$\max_{WW^{\top}=I_q} \operatorname{Tr}(W\Sigma W^{\top}).$

Оптимальный энкодер:

$W_{\mathrm{PCA}}=E_q^{\top}.$

Оптимальный декодер:

$V_{\mathrm{PCA}}=E_q.$

Минимальная ошибка восстановления равна

$\mathcal{E}_{\min} = \sum_{k=q+1}^{p}\lambda_k.$

Общий линейный автоэнкодер и неединственность параметров

Теперь снимем ограничения $V=W^{\top}$ и ортонормированности. Будем оптимизировать произвольные матрицы и при размерности скрытого представления .

Композиция имеет ранг не выше . Поэтому задача сводится к поиску наилучшего линейного оператора ранга не выше для восстановления центрированных данных:

$\min_{\operatorname{rank}(A)\leq q} \lVert\widetilde{X}-A\widetilde{X}\rVert_F^2.$

Оптимальное восстановленное подпространство совпадает с главным собственным подпространством PCA. Однако сами матрицы энкодера и декодера не обязаны быть равны $E_q^{\top}$ и .

Для любой обратимой матрицы

$R\in\mathbb{R}^{q\times q}$

можно взять

$W=RE_q^{\top},$ $V=E_qR^{-1}.$

Тогда

$VW = E_qR^{-1}RE_q^{\top} = E_qE_q^{\top}.$

Реконструкция останется той же, хотя скрытые координаты будут преобразованы матрицей . Следовательно, общий линейный автоэнкодер восстанавливает то же оптимальное подпространство, что и PCA, но его латентные координаты не обязаны быть ортогональными, декоррелированными или упорядоченными по дисперсии.

Чтобы получить именно канонические PCA-координаты, нужны дополнительные ограничения: ортонормированность, связанные веса или последующая ортогонализация.

Линейные и нелинейные методы снижения размерности

PCA использует аффинное кодирование

$y=E_q^{\top}(X-\overline{X})$

и аффинное восстановление

$\widehat{X}=\overline{X}+E_qy.$

Следовательно, все восстанавливаемые точки лежат в одном -мерном аффинном подпространстве. Это сильное ограничение.

Если данные расположены около нелинейного многообразия, линейная плоскость может описывать их плохо. Классический пример - «швейцарский рулет»: двумерная поверхность свернута в трехмерном пространстве. Глобальная линейная проекция не может развернуть такую структуру без сильных искажений.

Нелинейные методы допускают отображение вида

где нелинейно. В зависимости от цели используют разные подходы:

kernel PCA выполняет линейный PCA в неявном пространстве признаков, заданном ядром;
нелинейные автоэнкодеры обучают нейронные сети для кодирования и декодирования;
Isomap пытается сохранять геодезические расстояния на многообразии;
LLE сохраняет локальные линейные отношения между соседями;
t-SNE и UMAP в первую очередь предназначены для визуализации локальной структуры, а не для точного глобального восстановления данных.

Нелинейность дает большую выразительность, но имеет цену:

задача оптимизации обычно становится невыпуклой;
результат сильнее зависит от гиперпараметров;
может отсутствовать простое обратное преобразование;
новые точки не всегда легко отображать без дополнительной модели;
глобальные расстояния и дисперсии могут искажаться;
интерпретация компонент обычно сложнее, чем в PCA.

PCA остается полезной базовой моделью именно потому, что имеет закрытое решение, устойчивую геометрическую интерпретацию и прозрачную связь ошибки со спектром собственных значений.

Связь с линейной регрессией и моделями низкого ранга

Линейный автоэнкодер можно рассматривать как модель низкого ранга для всей матрицы данных:

$\widetilde{X} \approx U_q\Gamma_qV_q^{\top}.$

При этом

$U_q\Gamma_q$

описывает базис и масштабы вариации в пространстве признаков, а

$V_q^{\top}$

содержит коэффициенты, специфичные для объектов.

Как выбрать число компонент

Все рассмотренные формулировки упорядочивают компоненты по собственным значениям

$\lambda_1\geq\lambda_2\geq\cdots\geq\lambda_p.$

Число компонент можно выбирать несколькими способами.

Порог объясненной дисперсии

Выбираем минимальное , для которого сохраняется заданная доля дисперсии :

$q(P) = \min\left\{ q: \frac{\sum_{k=1}^{q}\lambda_k} {\sum_{k=1}^{p}\lambda_k} \geq P \right\}.$

Часто используют

$P=0.90, \qquad P=0.95, \qquad P=0.99.$

Но универсального правильного порога не существует: он зависит от задачи и стоимости потери информации.

Метод локтя

На scree plot изображают собственные значения в порядке убывания. Выбирают точку, после которой спектр заметно выравнивается. Формально можно смотреть на спектральные разрывы

$\lambda_q-\lambda_{q+1}.$

Большой разрыв указывает на относительно устойчивое отделение ведущего подпространства от остальных направлений.

Кросс-валидация по целевой задаче

Если PCA используется как часть supervised-пайплайна, разумнее выбирать по качеству конечной модели на валидации, а не только по объясненной дисперсии.

Компонента с большой дисперсией не обязательно наиболее полезна для предсказания целевой переменной. И наоборот, направление с небольшой дисперсией может содержать важный для классификации сигнал.

Ошибка восстановления

Средняя квадратичная ошибка при сохранении компонент равна

$\mathcal{E}_q = \sum_{k=q+1}^{p}\lambda_k.$

Поэтому выбор можно формулировать через допустимый бюджет ошибки:

$\sum_{k=q+1}^{p}\lambda_k \leq \varepsilon.$

Некоторая суммаризация

Дефакто все эти задачи сводятся к одной и той же задаче оптимизации:

$\max_{E^{\top}E=I_q} \operatorname{Tr}(E^{\top}\Sigma E).$

Ее оптимальное значение равно

$\sum_{k=1}^{q}\lambda_k,$

а решение задается ведущим собственным подпространством ковариационной матрицы.

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

$\min_{\operatorname{rank}(A)\leq q} \lVert M-A\rVert_F^2.$

Разница между взглядами - это только разница в том, что мы подставляем в роли и какие дополнительные ограничения накладываем на :

если $M=\widetilde{X}$ - сама центрированная матрица данных, а ищется без ограничений на структуру, получаем теорему Эккарта - Янга - Мирского и PCA как частный случай SVD;
если - матрица Грама попарных скалярных произведений, получаем классическое метрическое MDS;
если $A=P\widetilde{X}$ , где - ортогональный проектор, а $M=\widetilde{X}$ , получаем задачу о наилучшем линейном (аффинном, если вернуться к нецентрированным данным) подпространстве;
если $A=VW\widetilde{X}$ для произвольных матриц и подходящих размеров, а $M=\widetilde{X}$ , получаем общий линейный автоэнкодер;
наконец, максимизация $\operatorname{Tr}(E^{\top}\Sigma E)$ - это тот же самый вопрос, только заданный не через ошибку приближения, а через энергию, которая в этом приближении сохраняется; по теореме Пифагора одно спокойно пересчитывается в другое.

То есть внешне разные постановки - через дисперсию, через SVD, через расстояния, через геометрию подпространства, через кодирование и декодирование - на самом деле являются разными подходами к одному и тому же математическому факту:

Усечение спектрального разложения симметричной положительно полуопределенной матрицы - или SVD произвольной матрицы - до ведущих компонент дает наилучшее приближение ранга в смысле теоремы Эккарта-Янга-Мирского и одновременно выделяет ведущее -мерное подпространство, характеризуемое вариационными принципами Куранта-Фишера и Кая Фана..

Если хотите потыкать интерактивные визуализации, то можно посмотреть и потыкать мой англоязычный бложик. Также если знаете какие-то другие интересные взгляды на PCA, буду рад комментариям и идеям для расширения материала!

50 оттенков PCA