Метод главных компонент

У каждого из 469 городов мира 24 числа: средняя температура и сумма осадков каждого месяца за 2001–2025 годы. В 24-мерном пространстве города не разглядеть, а из 24 чисел не видно, кто на кого похож. Метод главных компонент поворачивает это пространство так, чтобы первая ось шла вдоль самого большого разброса городов, вторая — вдоль следующего, и так далее. На плоскости первых двух осей помещается 68,8 % всего разброса, и соседи на ней — города с похожим климатом. Цвет города — его место на этой плоскости. На карте те же цвета складываются в климатические пояса, хотя ни широт, ни долгот алгоритм не видел.

Что означают оси, видно по их весам под картой. У первой все 24 веса положительны: она отличает тепло и влажно от холодно и сухо. На одном её краю — Норильск, Якутск и Мирный, на другом — Арава, Бинтулу и Тимика. Вторая противопоставляет дожди жаркому лету: на её отрицательном краю пустыни — Донгола, Джидда, Эс-Сулайиль, на положительном — дождливые города с прохладным летом, Берген и Джуно, и снова Тимика. Переключите оси на 1–3: третья компонента — время дождей, летом северного полушария или летом южного; на её краях Дакка и Мумбаи против Дарвина и Антананариву. Коснитесь любого города: пунктир — его год, собранный обратно всего из k чисел.

Идея

Сдвинем облако точек так, чтобы его центр масс оказался в нуле: из каждого столбца вычтем его среднее. Возьмём единичное направление u и спроецируем на него все точки. Разброс проекций — их средний квадрат — зависит от u, и это квадратичная форма u^TCu, где C — матрица ковариаций: на её диагонали дисперсии столбцов, вне диагонали — ковариации пар. Направление с наибольшим разбросом — собственный вектор C с наибольшим собственным значением, а само значение и есть этот разброс. Почему так, подробно доказано в «Царице наук»; там же показано, что это направление — прямая, к которой точки ближе всего по перпендикуляру.

Главные компоненты — собственные векторы матрицы ковариаций, упорядоченные по убыванию собственных значений. Первая — направление наибольшего разброса данных, каждая следующая — наибольшего разброса среди направлений, перпендикулярных предыдущим. Собственное значение, делённое на сумму всех, — доля разброса, которую объясняет компонента. Оставив первые k, мы описываем каждую точку k числами вместо d и теряем как можно меньше.

Главная компонента — это направление, а координата города вдоль неё — счёт. Сумма собственных значений равна сумме дисперсий всех столбцов: поворот пространства не меняет общего разброса, а только перераспределяет его между осями. Поэтому доля первых двух компонент — 68,8 % — говорит, сколько всего разброса осталось на плоскости сцены.

Как найти главную компоненту

Формулы для собственных векторов матрицы 24 × 24 нет: собственные значения — корни многочлена 24-й степени. Их ищут последовательными приближениями, и самый простой способ — степенной метод. Возьмём любой вектор, умножим на C и нормируем, снова умножим и снова нормируем. Вектор поворачивается к главной компоненте.

Для стандартизованного климата λ₂/λ₁ = 0,487: каждый шаг уменьшает угол примерно вдвое, и за 35 шагов вектор перестаёт поворачиваться — изменение меньше 10⁻¹². Пунктир на графике — обещание теоремы, и сплошная линия идёт вдоль него. Скорость зависит только от отношения двух собственных значений: чем ближе λ₂ к λ₁, тем дольше спор между двумя направлениями.

Чтобы найти вторую компоненту, из матрицы вычитают уже найденную: C ← C − λ₁u₁u₁^T. Это исчерпывание: у новой матрицы те же собственные векторы, но на месте λ₁ стоит ноль, и наибольшим становится λ₂. Для второй компоненты климата решающее отношение — λ₃/λ₂ = 15,8/22,5 ≈ 0,70, и степенному методу нужно 72 шага; для третьей — λ₄/λ₃ ≈ 0,45 и снова 34. Ошибка каждой найденной компоненты переходит в следующие, поэтому исчерпыванием берут немного первых компонент, а весь спектр считают другими методами — например, вращениями Якоби, как во второй задаче статьи.

Код

Листинг делает всё сам, без библиотек: центрирует столбцы (и, если нужно, делит на их стандартные отклонения), собирает матрицу ковариаций, находит четыре компоненты степенным методом с исчерпыванием и печатает, что они значат. Знак собственного вектора произволен — u и −u задают одно направление, — поэтому договоримся, что самая большая по модулю координата положительна; иначе Python, Go и сцена могли бы показать зеркальные картинки.

standardized columns
  PC1:  46.2% of the variance,  46.2% together,   35 steps; heaviest column t10
  PC2:  22.5% of the variance,  68.8% together,   72 steps; heaviest column p11
  PC3:  15.8% of the variance,  84.6% together,   34 steps; heaviest column p8
  PC4:   7.1% of the variance,  91.6% together,   34 steps; heaviest column t7

Четыре числа вместо двадцати четырёх держат 91,6 % разброса. Четвёртая компонента — тоже знакомая вещь: её веса температуры положительны летом северного полушария и отрицательны зимой. На одном её краю Якутск и Багдад с их огромной разницей между июлем и январём, на другом — Неукен, Комодоро-Ривадавия и Сантьяго, где тёплый месяц — январь.

Градусы против миллиметров

Переключите в сцене «как есть: °C и мм». Картина меняется целиком:

raw columns (°C and mm)
  PC1:  52.1% of the variance,  52.1% together,   44 steps; heaviest column p12
  PC2:  32.1% of the variance,  84.2% together,   19 steps; heaviest column p7

Обе первые компоненты теперь про осадки — декабрьские и июльские. Дело не в климате, а в единицах. Стандартное отклонение осадков месяца по городам — от 75 до 100 миллиметров, температуры — от 6 до 16 градусов, а метод главных компонент сравнивает дисперсии просто как числа: у осадков они в десятки раз больше. Запиши мы осадки в сантиметрах, главными стали бы температуры.

Поэтому столбцы в разных единицах сначала приводят к одному масштабу: делят каждый на его стандартное отклонение. Это стандартизация, и после неё матрица ковариаций превращается в матрицу корреляций. Когда столбцы измерены в одном и том же и их разброс сам по себе важен — пиксели картинки, координаты точек, — стандартизовать не нужно: это исказило бы данные. Решение принимает человек, а не алгоритм, и популярные библиотеки сами его не принимают: PCA в scikit-learn центрирует столбцы, но не масштабирует.

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

Оставим от каждого города k счетов и соберём его обратно: среднее плюс счета, умноженные на компоненты. Насколько он изменится?

k = 1: typical miss 5.20 °C and 77.2 mm
k = 2: typical miss 4.22 °C and 60.5 mm
k = 3: typical miss 3.77 °C and 35.5 mm
k = 4: typical miss 1.12 °C and 33.3 mm

Третья компонента почти вдвое уменьшает промах в осадках — она отвечает за время дождей. Четвёртая втрое уменьшает промах в температуре: без неё модель не знает, в каком полушарии город и какой месяц у него летний. Это общая черта: компонента с маленькой долей разброса может нести то, что для человека важно.

Это не случайность, а теорема, которую в 1936 году доказали Карл Эккарт и Гейл Янг. Среди всех k-мерных плоскостей наименьшая сумма квадратов расстояний от точек — у плоскости, проходящей через центр масс вдоль первых k главных компонент, и эта сумма, делённая на число точек, равна сумме отброшенных собственных значений λ_k+1 + … + λ_d (если матрицу ковариаций делить на n). Доказали они её в форме для матриц: лучшее приближение матрицы матрицей меньшего ранга даёт усечённое сингулярное разложение. Доказательство той же теоремы для наибольшего растяжения разобрано в «Царице наук»; для суммы квадратов оно длиннее, и здесь мы его не повторяем. Отсюда правило выбора k: смотрят на накопленную долю разброса (оранжевая линия под сценой) и останавливаются, когда она достаточна для задачи, — 90 % или 99 % — или там, где столбики резко мельчают. Для наглядной картинки берут две-три компоненты, сколько бы их ни требовалось для точности.

Сколько это стоит

Матрица ковариаций для n точек в d измерениях стоит O(n·d²): каждая точка добавляет все попарные произведения своих координат. Все собственные векторы матрицы d × d — O(d³), одна компонента степенным методом — O(d²) на шаг. Для климата это мгновения: 469 × 24² ≈ 270 тысяч умножений и 35 шагов по 576.

Наклоны обеих линий на логарифмических осях близки к двум — и матрица, и шаг степенного метода растут как d²; при малых d их немного прижимают накладные расходы. Когда d — тысячи и миллионы (пиксели фотографий, гены, слова), матрица d × d перестаёт помещаться в память, и её не строят. Шаг степенного метода можно сделать прямо по данным: Cv = X^T(Xv)/(n − 1) — два умножения на таблицу, O(n·d). А если точек мало, а измерений много, собственные векторы ищут у маленькой матрицы XX^T размера n × n и пересчитывают обратно; так в 1991 году Мэтью Тёрк и Алекс Пентленд нашли главные компоненты изображений лиц — «собственные лица». Для больших таблиц есть рандомизированные методы, которые описали Натан Халко, Пер-Гуннар Мартинссон и Джоэл Тропп: сначала умножить данные на несколько случайных векторов, чтобы поймать подпространство главных направлений, и только потом решать маленькую задачу. PCA в scikit-learn выбирает способ сам: при числе столбцов меньше тысячи и строк больше чем вдесятеро — собственные векторы матрицы ковариаций, как у нас, для больших таблиц — рандомизированный способ, иначе — полное сингулярное разложение.

Где метод ошибается

• Только прямые направления. Главные компоненты — линейные комбинации столбцов. Если данные лежат на изогнутой поверхности, например на свёрнутом рулоном листе, никакая плоскость её не развернёт. Для наглядных картинок таких данных есть нелинейные методы вроде t-SNE и UMAP, но их оси уже не имеют простого смысла.
• Квадраты и выбросы. Разброс считается через квадраты, и одна далёкая точка может развернуть первую компоненту на себя — та же слабость, что у метода наименьших квадратов.
• Разброс — не важность. Направление, вдоль которого точки различаются сильнее всего, не обязательно отделяет то, что нужно отличать. Две группы могут различаться как раз вдоль малой компоненты, и, выбросив её, мы выбросим ответ.
• Компоненты — не причины. «Тепло и влажно» — удобное имя для первой оси климата, но ось — смесь всех 24 столбцов, и имена ей придумывает человек. На других данных компоненты бывает трудно назвать вовсе.
• Это не кластеризация. Метод поворачивает облако, но не делит его на группы. Делить удобно уже в нескольких главных компонентах — k-средних или DBSCAN.

История

Карл Пирсон в 1901 году искал прямые и плоскости, ближайшие к системе точек, когда ошибаются все координаты сразу, — и пришёл к главным осям облака. Гарольд Хотеллинг в 1933 году развил метод в статье о разложении набора связанных статистических переменных и дал ему имя «главные компоненты». Степенной метод для собственных векторов обычно связывают с Рихардом фон Мизесом и Хильдой Поллачек-Гейрингер, описавшими его в 1929 году, хотя его идеи встречаются и раньше. В 1936 году Эккарт и Янг доказали теорему о наилучшем приближении матрицы меньшего ранга. С появлением компьютеров метод стал первым, что делают с многомерными данными: в 1987 году Лоуренс Сирович и Майкл Кирби описали изображения лиц немногими главными компонентами, а в 1991 году Тёрк и Пентленд построили на этом распознавание лиц.

Задачи

Первая задача — главная ось облака на плоскости, у которой есть формула; на её тестах — пояса землетрясений. Вторая — сколько потеряет климат, если сжать его до k чисел, и тут нужен весь спектр матрицы.