Понижение размерности

Смеси распределений из предыдущей заметки объясняли каждую точку дискретной латентной переменной — кластером, из которого она пришла. Многие данные лучше объясняются несколькими непрерывными. Изображение рукописной цифры $8 \times 8$ состоит из 64 значений пикселей, но изображения меняются вдоль гораздо меньшего числа направлений: какая это цифра, её наклон, ширина, толщина штриха. Точки заполняют тонкую область в $\mathbb R^{64}$, близкую к подпространству или поверхности малой размерности. Найти эту область — значит получить координаты, которые сжимают данные, убирают шум, делают расстояния осмысленнее и позволяют посмотреть на данные в двух измерениях. Эта заметка выводит линейный ответ, метод главных компонент, с двух сторон, которые сходятся в одной задаче на собственные значения, даёт ему вероятностное прочтение как модели с латентными переменными, а затем переходит к нелинейным методам, которые сохраняют окрестности, а не дисперсию.

Обозначения. $x_1, \dots, x_n \in \mathbb R^d$ — точки данных, центрированные так, что $\sum_i x_i = 0$, а $X \in \mathbb R^{n \times d}$ — матрица со строками $x_i^\top$. $S = \frac1n X^\top X$ — выборочная ковариационная матрица с собственными значениями $\lambda_1 \ge \dots \ge \lambda_d \ge 0$ и ортонормированными собственными векторами $v_1, \dots, v_d$. $k < d$ — целевая размерность, $W \in \mathbb R^{d \times k}$ — матрица с ортонормированными столбцами ($W^\top W = I_k$), а $z_i = W^\top x_i \in \mathbb R^k$ — новые координаты $x_i$. $X = U \Sigma V^\top$ — сингулярное разложение с сингулярными числами $\sigma_1 \ge \sigma_2 \ge \dots$.

Зачем понижать размерность#

Четыре причины, часто все сразу:

• Посмотреть. У графиков две оси. Отображение в две размерности — единственный способ увидеть все 64 признака набора данных на одной диаграмме рассеяния.
• Сжатие. Хранение $k$ координат вместо $d$ экономит память и ускоряет всё, что идёт дальше.
• Шумоподавление. Если сигнал живёт в нескольких направлениях, а шум распределён по всем, отбрасывание слабых направлений убирает в основном шум.
• Лучшие расстояния. В высокой размерности расстояния между точками концентрируются вокруг общего значения (сферический слой гауссианы), что вредит ближайшим соседям и методам кластеризации этого раздела. Расстояния в хорошо выбранном пространстве малой размерности информативнее.

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

Спроецируем данные на одно направление — единичный вектор $w$. Проекции $w^\top x_i$ имеют нулевое среднее, потому что данные центрированы, и дисперсию

$$\frac1n \sum_{i=1}^n (w^\top x_i)^2 = w^\top \Big(\frac1n \sum_{i=1}^n x_i x_i^\top\Big) w = w^\top S w .$$

Направление, вдоль которого точки остаются разнесёнными, сохраняет различия между ними, поэтому ищем направление наибольшей дисперсии: максимизируем $w^\top S w$ при условии $w^\top w = 1$. У лагранжиана $w^\top S w - \lambda (w^\top w - 1)$ градиент $2Sw - 2\lambda w$, и он обращается в ноль, когда

$$S w = \lambda w .$$

Значит, $w$ должен быть собственным вектором $S$, и тогда дисперсия равна $w^\top S w = \lambda w^\top w = \lambda$: лучшее направление — собственный вектор с наибольшим собственным значением, $v_1$, первая главная компонента (first principal component). Вторая компонента — направление наибольшей дисперсии, ортогональное первой. Записав любой единичный $w$ в собственном базисе, $w = \sum_j c_j v_j$ с $\sum_j c_j^2 = 1$, получаем $w^\top S w = \sum_j \lambda_j c_j^2$; при $c_1 = 0$, которого требует ортогональность, максимум равен $\lambda_2$ и достигается при $w = v_2$. Продолжая, получаем, что $k$ старших собственных векторов — это $k$ главных компонент, а дисперсия данных вдоль $v_j$ в точности равна $\lambda_j$.

В заметке о многомерном нормальном распределении те же собственные векторы встречались как оси эллипса ковариации. Метод главных компонент (principal component analysis, PCA) — это та же геометрия, оценённая по данным: главные компоненты — оси эллипса, лучше всего описывающего облако, упорядоченные по длине.

PCA как минимизация ошибки восстановления#

Теперь зададим другой вопрос. Представим каждую точку $k$ числами $z_i = W^\top x_i$ и восстановим её как $\hat x_i = W z_i = W W^\top x_i$. При каком $W$ восстановленные точки ближе всего к данным? Поскольку $W W^\top$ — ортогональный проектор, остаток $x_i - \hat x_i$ ортогонален $\hat x_i$, и по теореме Пифагора

$$\lVert x_i - W W^\top x_i \rVert^2 = \lVert x_i \rVert^2 - \lVert W^\top x_i \rVert^2 .$$

Усредняя по точкам, получаем

$$\frac1n \sum_{i=1}^n \lVert x_i - \hat x_i \rVert^2 = \operatorname{tr} S - \operatorname{tr}(W^\top S W) .$$

Первое слагаемое — полная дисперсия, заданная данными; второе — дисперсия, сохранённая проекцией. Поэтому минимизация ошибки восстановления — та же задача, что максимизация сохранённой дисперсии, и решение у неё то же: $W = (v_1, \dots, v_k)$, с минимальной ошибкой

$$\frac1n \sum_{i=1}^n \lVert x_i - \hat x_i \rVert^2 = \sum_{j > k} \lambda_j ,$$

суммой отброшенных собственных значений. (То, что $k$ старших собственных векторов максимизируют $\operatorname{tr}(W^\top S W)$ по всем ортонормированным $W$, следует из того же рассуждения в собственном базисе, применённого к каждому столбцу, или из теоремы Ки Фаня.) Два вывода — это два исторических истока PCA: Pearson (1901) искал прямую и плоскость, ближе всего подходящие к точкам, Hotelling (1933) — компоненты наибольшей дисперсии.

Рисунок

Две стороны одной задачи: дисперсия и ошибка восстановления

Разобранный пример. Для ковариационной матрицы $\Sigma = \begin{pmatrix} 2 & 1.2 \\ 1.2 & 1 \end{pmatrix}$ из курса теории вероятностей собственные значения равны $2.8$ и $0.2$, а первый собственный вектор — $v_1 = (0.832, 0.555)$. Первая компонента сохраняет $2.8 / 3 = 93.3\%$ дисперсии, а восстановление только по ней оставляет среднеквадратичную ошибку $0.2$. Точка $x = (1, 1)$ получает координату $z = v_1^\top x = 1.387$ и восстановление $\hat x = z v_1 = (1.154, 0.769)$; квадрат длины остатка $(-0.154, 0.231)$ равен $0.077$, и действительно $\lVert x \rVert^2 = 2 = 1.387^2 + 0.077$.

Через SVD#

У любой матрицы есть сингулярное разложение (singular value decomposition, SVD) $X = U \Sigma V^\top$, где у $U \in \mathbb R^{n \times r}$ и $V \in \mathbb R^{d \times r}$ ортонормированные столбцы, а $\Sigma$ диагональна с $\sigma_1 \ge \dots \ge \sigma_r > 0$ ($r$ — ранг). Тогда

$$S = \frac1n X^\top X = \frac1n V \Sigma U^\top U \Sigma V^\top = V \frac{\Sigma^2}{n} V^\top ,$$

— спектральное разложение: главные компоненты — правые сингулярные векторы $V$, собственные значения равны $\lambda_j = \sigma_j^2 / n$, а координаты всех точек — $XV = U\Sigma$. Есть две причины вычислять PCA именно так. Построение $X^\top X$ возводит в квадрат число обусловленности, поэтому малые собственные значения теряют точность, а SVD работает с $X$ напрямую. И для больших данных усечённое или рандомизированное SVD вычисляет только $k$ старших компонент, без полной матрицы $d \times d$.

SVD показывает и то, что PCA оптимален в более сильном смысле. По теореме Эккарта–Янга (Eckart–Young theorem; Eckart and Young, 1936) лучшее приближение $X$ матрицей ранга $k$ по норме Фробениуса — усечённое SVD $X_k = U_k \Sigma_k V_k^\top$ с ошибкой $\lVert X - X_k \rVert_F^2 = \sum_{j > k} \sigma_j^2$. Восстановление PCA $X W W^\top$ при $W = V_k$ — это в точности $X_k$: никакое линейное описание данных ранга $k$ не лучше. Усечённое SVD матрицы счётчиков «слово–документ» — это латентно-семантический анализ, базовый метод из заметки о GloVe.

Объяснённая дисперсия и выбор k#

Доля дисперсии, сохраняемая первыми $k$ компонентами, — это доля объяснённой дисперсии (explained variance ratio)

$$\frac{\lambda_1 + \dots + \lambda_k}{\lambda_1 + \dots + \lambda_d} ,$$

а единица минус она — относительная ошибка восстановления. На 1797 изображениях цифр из scikit-learn, по 64 пикселя в каждом, первые две компоненты объясняют 14.9% и 13.6% дисперсии; 5 компонент объясняют половину, 21 — 90%, 29 — 95%. При 21 компоненте среднеквадратичная ошибка восстановления на изображение равна 116.3 — сумме 43 отброшенных собственных значений, как и предсказывает вывод выше.

Единого правила для $k$ нет. Распространённые варианты:

• Доля дисперсии, например 90% или 95%. Просто, но порог произволен.
• Излом на графике каменистой осыпи (scree plot), то есть собственных значений в порядке убывания: оставить компоненты до того места, где кривая выходит на «пол» из малых похожих собственных значений, — так выглядит изотропный шум.
• Последующая задача: выбрать $k$ по ошибке на валидации модели, которая использует компоненты.
• Правдоподобие вероятностного PCA (ниже), на отложенных данных или со штрафом вроде BIC.

Рисунок

Сколько компонент оставить

На рисунке первые два правила расходятся. У данных три направления сигнала и шум во всех десяти; график осыпи показывает три, но они объясняют только 89.9% дисперсии, поэтому правило «90%» оставляет четвёртую компоненту, которая является чистым шумом.

Стандартизация и отбеливание#

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

Отбеливание (whitening) идёт на шаг дальше и масштабирует каждую компоненту к единичной дисперсии:

$$z_i = \Lambda_k^{-1/2} W^\top x_i, \qquad \Lambda_k = \operatorname{diag}(\lambda_1, \dots, \lambda_k) .$$

Новые координаты некоррелированы и имеют единичную дисперсию, $\frac1n \sum_i z_i z_i^\top = I_k$, — выборочный вариант отбеливающего преобразования из заметки о многомерном нормальном распределении. Отбеливание помогает методам, которые предполагают сравнимые некоррелированные входы, например анализу независимых компонент или методам, основанным на расстояниях. Его риск — деление на $\sqrt{\lambda_j}$: направления с крошечными собственными значениями, часто чистый шум, усиливаются до масштаба сигнала, поэтому отбеливание сочетают с усечением до $k$ компонент или с добавлением небольшой константы к каждому $\lambda_j$.

Вероятностный PCA#

PCA в описанном виде — проекция, а не модель: он не приписывает точке вероятности. Вероятностный PCA (probabilistic PCA; Tipping and Bishop, 1999) — непрерывный аналог гауссовской смеси. Латентный вектор $z \sim \mathcal N(0, I_k)$ линейно отображается в пространство данных и размывается изотропным шумом, $x = Wz + \mu + \varepsilon$ с $\varepsilon \sim \mathcal N(0, \sigma^2 I_d)$; по правилу для аффинных отображений $x \sim \mathcal N(\mu, W W^\top + \sigma^2 I_d)$ — гауссиана, ковариация которой — «ранг $k$ плюс шум». Её решение методом максимального правдоподобия — замаскированный PCA: столбцы $W$ порождают подпространство $k$ старших собственных векторов, $W = V_k (\Lambda_k - \sigma^2 I)^{1/2} R$ для любого поворота $R$, а дисперсия шума равна среднему отброшенных собственных значений, $\sigma^2 = \frac{1}{d - k}\sum_{j > k} \lambda_j$ (Bishop, 2006, раздел 12.2.1). На рисунке с графиком осыпи это среднее равно $0.250$, истинной дисперсии шума; для цифр при $k = 10$ оно равно $5.83$. Апостериорное среднее латентного вектора равно $\mathbb E[z \mid x] = (W^\top W + \sigma^2 I_k)^{-1} W^\top (x - \mu)$. Это не вектор главных координат $V_k^\top (x - \mu)$: в $W$ лежат масштабированные и повёрнутые нагрузки, а не ортонормированные собственные векторы, и при $\sigma^2 \to 0$ среднее стремится к $R^\top \Lambda_k^{-1/2} V_k^\top (x - \mu)$ — главным координатам, делённым на свои стандартные отклонения и повёрнутым на $R$. Проекцией PCA становится восстановление в пространстве данных, $W \, \mathbb E[z \mid x] + \mu \to V_k V_k^\top (x - \mu) + \mu$. Модель можно обучать и EM-алгоритмом, который обходится без ковариационной матрицы $d \times d$ и справляется с пропущенными значениями, а поскольку это плотность, она даёт правдоподобие для выбора $k$ и оценку необычности точек — отправную точку следующей заметки.

Ядерный PCA#

PCA находит только линейную структуру. У двух концентрических окружностей нет направления большой дисперсии, которое бы их разделяло: любая проекция на прямую смешивает кольца. Ядерный PCA (kernel PCA; Schölkopf, Smola and Müller, 1998) выполняет PCA в пространстве признаков $\varphi(x)$, доступном через ядро $k(x, x') = \varphi(x)^\top \varphi(x')$, как в заметке о ядрах. Главные направления в пространстве признаков — комбинации отображённых точек, $\sum_i \alpha_i \varphi(x_i)$, и задача на собственные значения превращается в задачу для матрицы ядра $n \times n$, $K_{ij} = k(x_i, x_j)$, центрированной в пространстве признаков:

$$\tilde K = K - \mathbf 1_n K - K \mathbf 1_n + \mathbf 1_n K \mathbf 1_n , \qquad \tilde K \alpha = n \lambda \alpha ,$$

где $\mathbf 1_n$ — матрица $n \times n$, все элементы которой равны $1/n$, а $\alpha$ нормирован так, чтобы направление имело единичную длину, $n\lambda\, \alpha^\top \alpha = 1$. Координата точки на компоненте равна $\sum_i \alpha_i \tilde k(x_i, x)$ и вычисляется только через ядра. При $k(x, x') = (x^\top x')^2$ пространство признаков содержит $x_1^2$, $x_2^2$ и $x_1 x_2$, поэтому квадрат радиуса $x_1^2 + x_2^2$, разделяющий окружности, там линейная функция. В разделе «Код» логистическая регрессия на двух обычных главных компонентах окружностей даёт точность на обучающих данных 0.502 — уровень случайного угадывания, а на трёх ядерных компонентах — 1.0. Обе модели обучаются и оцениваются на одних и тех же 400 точках, поэтому числа показывают, что каждое представление может выразить линейно, а не качество на тесте. Только третья ядерная компонента разделяет кольца сама по себе: ядерный PCA упорядочивает направления по дисперсии в пространстве признаков, а не по полезности. Цена — задача на собственные значения размера $n \times n$, а у компонент нет простой интерпретации во входном пространстве.

t-SNE и UMAP: окрестности вместо дисперсии#

Для визуализации сохранение больших расстояний — часто неверная цель. Двумерная линейная проекция цифр сохраняет направления наибольшего разброса, и десять классов цифр сильно перекрываются. t-SNE (van der Maaten and Hinton, 2008) требует другого: точки, близкие в данных, должны остаться близкими на плоскости.

1. Превратить расстояния в данных в вероятности соседства. Для каждой точки $x_i$ положим $p_{j \mid i} \propto \exp\big(-\lVert x_i - x_j \rVert^2 / 2\sigma_i^2\big)$ с нормировкой по $j \ne i$, где каждая ширина $\sigma_i$ подбирается так, чтобы у распределения была заданная перплексия (perplexity), $2^{H(p_{\cdot \mid i})}$, — примерно число соседей, которые учитываются; типичные значения от 5 до 50. Симметризовать: $p_{ij} = (p_{j \mid i} + p_{i \mid j}) / 2n$.
2. Определить такие же вероятности для точек $y_i$ на плоскости, с ядром $t$-распределения Стьюдента (Student-t) с тяжёлыми хвостами: $q_{ij} \propto (1 + \lVert y_i - y_j \rVert^2)^{-1}$ с нормировкой по всем парам.
3. Двигать $y_i$ градиентным спуском, минимизируя $D_{\mathrm{KL}}(P \,\|\, Q) = \sum_{i \ne j} p_{ij} \log (p_{ij} / q_{ij})$.

Что сохраняется, решает асимметрия KL-дивергенции. Пара с большим $p_{ij}$, настоящие соседи, поставленная далеко (малое $q_{ij}$), стоит дорого; пара с малым $p_{ij}$, поставленная близко, стоит дёшево. t-SNE сохраняет окрестности и свободно обращается со всем остальным. Тяжёлый хвост отвечает на проблему скученности (crowding problem): у точки в 64 измерениях может быть много равноудалённых соседей, у точки на плоскости — нет, и ядро $t$ позволяет отодвинуть умеренно далёкие точки далеко, а не сваливать их в середину. В разделе «Код» 58.4% из 10 ближайших соседей каждого изображения цифры в 64 измерениях входят в число его 10 ближайших соседей на карте t-SNE, против 11.8% в двумерной проекции PCA.

UMAP (McInnes, Healy and Melville, 2018) следует тому же плану с другими решениями: взвешенный граф $k$ ближайших соседей в данных, похожий граф на плоскости и кросс-энтропия между весами их рёбер, оптимизируемая стохастическим градиентным спуском. Он гораздо быстрее, масштабируется на миллионы точек и умеет вкладывать новые точки, и часто говорят, что он лучше сохраняет глобальную структуру; правила чтения ниже относятся и к нему.

Как читать такие карты (Wattenberg, Viégas and Johnson, 2016):

• Гиперпараметры действительно важны. Перплексия (у UMAP — число соседей) меняет картину; стоит посмотреть на несколько значений.
• Размеры кластеров ничего не значат. t-SNE расширяет плотные кластеры и сжимает разреженные, поэтому площадь кластера ничего не говорит о его разбросе в данных.
• Расстояния между кластерами могут ничего не значить. Два кластера, далёкие на карте, могут быть близки в данных, и наоборот.
• Случайный шум может выглядеть кластеризованным, особенно при малой перплексии.
• Для форм и топологии одного графика мало. Кольцо или линия на карте может быть артефактом одного запуска.

Надёжно карта показывает, какие точки — соседи. Она для того, чтобы смотреть, а не измерять: у её осей нет смысла, в обычном t-SNE нет отображения для новых точек, а её координаты — плохие признаки для модели. Заметка об автоэнкодерах читает графики t-SNE выученных кодов именно с этими оговорками, а заметка о векторных представлениях слов сравнивает проекции PCA и t-SNE векторов слов.

Автоэнкодеры: нелинейный PCA#

PCA — лучшая линейная пара «энкодер–декодер»: $z = W^\top x$, $\hat x = Wz$, обученная минимизировать квадратичную ошибку восстановления. Обучение линейной сети с узким местом (bottleneck) из $k$ нейронов на той же функции потерь градиентным спуском приходит к тому же подпространству: её минимум порождает главное подпространство, и других локальных минимумов у неё нет, только седловые точки (Baldi and Hornik, 1989). Замена двух линейных отображений нейросетями, энкодером $z = f_\varphi(x)$ и декодером $\hat x = g_\theta(z)$, даёт автоэнкодер (autoencoder), который может следовать изогнутой поверхности, а не плоскому подпространству. Заметка об автоэнкодерах курса генеративных моделей строит на этой идее свёрточные автоэнкодеры для изображений, а вероятностная версия, в которой энкодер — приближённое апостериорное распределение, — это вариационный автоэнкодер, вероятностный PCA, ставший нелинейным.

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

	PCA	Вероятностный PCA	Ядерный PCA	t-SNE / UMAP	Автоэнкодер
Отображение	линейное, $z = W^\top x$	линейное, с плотностью	линейное в пространстве признаков	нет: оптимизируются положения	нейросетевые энкодер и декодер
Сохраняет	дисперсию, большие расстояния	то же, плюс правдоподобие	дисперсию в пространстве признаков	окрестности	то, что нужно для восстановления
Новые точки	проекция	проекция	ядра с обучающей выборкой	нет в обычном t-SNE	кодирование
Решение	SVD, точное	замкнутый вид или EM	задача на собственные значения $n \times n$	невыпуклое, случайное	невыпуклое, градиентный спуск
Применение	сжатие, шумоподавление, признаки	то же, плотность, пропущенные данные	нелинейные признаки	только визуализация	представления, генерация

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

• Линейность и ориентация на дисперсию. PCA сохраняет направления наибольшей дисперсии, а они не обязаны быть важными для задачи; направление с малой дисперсией может нести всю информацию о классе.
• Чувствительность к масштабу и выбросам. Единицы измерения определяют компоненты, и одна далёкая точка может забрать себе компоненту; для этого существуют робастные варианты PCA.
• Компоненты трудно интерпретировать. Каждая — смесь всех признаков, а знаки определены лишь с точностью до $\pm 1$.
• Нелинейные карты — чтобы смотреть. t-SNE и UMAP искажают размеры и расстояния, зависят от гиперпараметров и начальных значений генератора а обычный t-SNE не даёт функции для новых точек. У UMAP она есть (transform в библиотеке umap-learn), но новые положения наследуют те же искажения; конвейер, который передаёт вложение модели, должен обучать его только на обучающей части.
• Внутренняя размерность — выбор при моделировании. Каждое правило для $k$ закладывает предположение о том, что сигнал, а что шум.

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

• PCA проецирует центрированные данные на старшие собственные векторы ковариационной матрицы $S$; дисперсия вдоль $v_j$ равна $\lambda_j$.
• Максимизация сохранённой дисперсии и минимизация квадратичной ошибки восстановления — одна задача, поскольку $\frac1n\sum \lVert x_i - \hat x_i\rVert^2 = \operatorname{tr} S - \operatorname{tr}(W^\top S W)$; минимальная ошибка равна $\sum_{j > k} \lambda_j$.
• PCA вычисляют через SVD $X = U\Sigma V^\top$: компоненты $V$, собственные значения $\sigma_j^2/n$, координаты $U\Sigma$; по теореме Эккарта–Янга никакое приближение ранга $k$ не лучше.
• $k$ выбирают по объяснённой дисперсии, излому графика осыпи, последующей задаче или правдоподобию; правила могут расходиться.
• Признаки с разными единицами измерения нужно стандартизовать; отбеливание масштабирует компоненты к единичной дисперсии и усиливает шумовые направления.
• У вероятностного PCA, $x \sim \mathcal N(\mu, WW^\top + \sigma^2 I)$, решение методом максимального правдоподобия — PCA, а $\sigma^2$ — среднее отброшенных собственных значений.
• Ядерный PCA — это PCA центрированной матрицы ядра: нелинейные компоненты ценой задачи на собственные значения размера $n \times n$.
• t-SNE и UMAP сохраняют окрестности; размеры кластеров и промежутки между ними на их картах несут мало информации.
• Линейный автоэнкодер выучивает подпространство PCA; нелинейный его обобщает.

Код#

PCA через SVD, объяснённая дисперсия и ошибка восстановления на цифрах, вероятностный PCA, отбеливание, ядерный PCA и окрестности на картах PCA и t-SNE:

import numpy as np
from sklearn.datasets import load_digits, make_circles
from sklearn.decomposition import PCA, KernelPCA
from sklearn.manifold import TSNE
from sklearn.neighbors import NearestNeighbors
from sklearn.linear_model import LogisticRegression

# распределение: собственные значения ковариационной матрицы и доля первой компоненты
S = np.array([[2.0, 1.2], [1.2, 1.0]])
lam = np.linalg.eigvalsh(S)[::-1]
print(lam.round(3), round(lam[0] / lam.sum(), 3))                  # [2.8 0.2] 0.933

# PCA выборки через SVD центрированной матрицы данных: то же, что собственные векторы ковариационной матрицы
rng = np.random.default_rng(0)
X = rng.multivariate_normal([0, 0], S, size=2000)
Xc = X - X.mean(axis=0)
U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
ev = np.linalg.eigvalsh(Xc.T @ Xc / len(X))[::-1]
print((s ** 2 / len(X)).round(3), ev.round(3), np.abs(Vt[0]).round(3))   # [2.812 0.198] [2.812 0.198] [0.826 0.563]

# цифры: 1797 изображений 8 x 8 пикселей; доля первых двух компонент, число компонент для 50%, 90%, 95%
D = load_digits().data
p = PCA().fit(D)
cum = np.cumsum(p.explained_variance_ratio_)
print(p.explained_variance_ratio_[:2].round(3), [int(np.searchsorted(cum, q) + 1) for q in (0.5, 0.9, 0.95)])
# [0.149 0.136] [5, 21, 29]

# ошибка восстановления по k компонентам = сумма отброшенных собственных значений (на изображение, ковариация с 1/n)
for k in (2, 10, 21):
    pk = PCA(k).fit(D)
    err = ((D - pk.inverse_transform(pk.transform(D))) ** 2).sum(axis=1).mean()
    lam_n = p.explained_variance_ * (len(D) - 1) / len(D)
    print(k, round(err, 1), round(lam_n[k:].sum(), 1), round(cum[k - 1], 3))
# 2 858.9 858.9 0.285 | 10 314.5 314.5 0.738 | 21 116.3 116.3 0.903

# вероятностный PCA: дисперсия шума — среднее отброшенных собственных значений
p10 = PCA(10).fit(D)
print(round(p10.noise_variance_, 3), round(p.explained_variance_[10:].mean(), 3))   # 5.828 5.828
# отбеливание: у масштабированных координат единичная ковариационная матрица
Z = PCA(10, whiten=True).fit_transform(D)
print(np.allclose(np.cov(Z.T), np.eye(10)))                         # True

# ядерный PCA на двух концентрических окружностях: точность линейного классификатора на компонентах на обучающих данных
Xr, yr = make_circles(n_samples=400, factor=0.4, noise=0.05, random_state=0)
acc = lambda z: round(LogisticRegression().fit(z, yr).score(z, yr), 3)
z_lin = PCA(2).fit_transform(Xr)
z_ker = KernelPCA(3, kernel='poly', degree=2, gamma=1, coef0=0).fit_transform(Xr)   # k(x, y) = (x . y)^2
print(acc(z_lin), acc(z_ker), [acc(z_ker[:, [j]]) for j in range(3)])   # 0.502 1.0 [0.53, 0.535, 1.0]

# окрестности: доля 10 ближайших соседей каждого изображения в 64-D, оставшихся среди 10 ближайших в 2-D
def kept(A, B, k=10):
    ia = NearestNeighbors(n_neighbors=k + 1).fit(A).kneighbors(A, return_distance=False)[:, 1:]
    ib = NearestNeighbors(n_neighbors=k + 1).fit(B).kneighbors(B, return_distance=False)[:, 1:]
    return np.mean([len(set(a) & set(b)) / k for a, b in zip(ia, ib)])
for name, E in (('pca', PCA(2).fit_transform(D)), ('t-sne', TSNE(2, perplexity=30, random_state=0).fit_transform(D))):
    print(name, round(kept(D, E), 3))                               # pca 0.118 | t-sne 0.584

Дальше#

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

Литература#

• K. Pearson. On lines and planes of closest fit to systems of points in space. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(11), 559–572, 1901.
• H. Hotelling. Analysis of a complex of statistical variables into principal components. Journal of Educational Psychology, 24(6), 417–441, 1933.
• C. Eckart, G. Young. The approximation of one matrix by another of lower rank. Psychometrika, 1(3), 211–218, 1936.
• M. E. Tipping, C. M. Bishop. Probabilistic Principal Component Analysis. Journal of the Royal Statistical Society B, 61(3), 611–622, 1999.
• B. Schölkopf, A. Smola, K.-R. Müller. Nonlinear Component Analysis as a Kernel Eigenvalue Problem. Neural Computation, 10(5), 1299–1319, 1998.
• L. van der Maaten, G. Hinton. Visualizing Data using t-SNE. JMLR, 9, 2579–2605, 2008.
• L. McInnes, J. Healy, J. Melville. UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. arXiv:1802.03426, 2018.
• M. Wattenberg, F. Viégas, I. Johnson. How to Use t-SNE Effectively. Distill, 2016.
• P. Baldi, K. Hornik. Neural networks and principal component analysis: Learning from examples without local minima. Neural Networks, 2(1), 53–58, 1989.
• C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006. Раздел 12.1 (PCA, обе постановки), раздел 12.2 (вероятностный PCA; 12.2.1 — решение методом максимального правдоподобия), раздел 12.3 (ядерный PCA).
• M. P. Deisenroth, A. A. Faisal, C. S. Ong. Mathematics for Machine Learning. Cambridge University Press, 2020. Глава 10 (понижение размерности методом PCA), раздел 4.5 (SVD), раздел 4.6 (приближение матриц).
• E. Sokolov et al. Машинное обучение (HSE): конспекты лекций, лекция 12: обучение без учителя (2021); задачи по PCA (2025).