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

Метод главных компонент (англ. principal component analysis, PCA) — один из основных способов уменьшить размерность данных, потеряв наименьшее количество информации. Изобретён Карлом Пирсоном в 1901 году. Применяется во многих областях, в том числе в эконометрике, биоинформатике, обработке изображений, для сжатия данных, в общественных науках^[⇨].

Вычисление главных компонент может быть сведено к вычислению сингулярного разложения матрицы данных^[⇨] или к вычислению собственных векторов и собственных значений ковариационной матрицы исходных данных^[⇨]. Иногда метод главных компонент называют преобразованием Кархунена — Лоэва^[1] или преобразованием Хотеллинга (англ. Hotelling transform).

Задача анализа главных компонент имеет, как минимум, четыре базовые версии:

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

Первые три версии оперируют конечными множествами данных. Они эквивалентны и не используют никакой гипотезы о статистическом порождении данных. Четвёртая версия оперирует случайными величинами. Конечные множества появляются здесь как выборки из данного распределения, а решение трёх первых задач — как приближение к разложению по теореме Кархунена — Лоэва («истинному преобразованию Кархунена — Лоэва»). При этом возникает дополнительный и не вполне тривиальный вопрос о точности этого приближения.

Аппроксимация данных линейными многообразиями

Иллюстрация к знаменитой работе Пирсона (1901): даны точки ${\displaystyle P_{i}}$ на плоскости, ${\displaystyle p_{i}}$ — расстояние от ${\displaystyle P_{i}}$ до прямой ${\displaystyle AB}$. Ищется прямая ${\displaystyle AB}$, минимизирующая сумму ${\displaystyle \sum _{i}p_{i}^{2}}$

Метод главных компонент начинался с задачи наилучшей аппроксимации конечного множества точек прямыми и плоскостями (Пирсон, 1901). Дано конечное множество векторов ${\displaystyle x_{1},x_{2},\dots ,x_{m}\in \mathbb {R} ^{n}}$, для каждого ${\displaystyle k=0,1,\dots ,n-1}$ среди всех ${\displaystyle k}$-мерных линейных многообразий в ${\displaystyle \mathbb {R} ^{n}}$ найти такое ${\displaystyle L_{k}\subset \mathbb {R} ^{n}}$, что сумма квадратов уклонений ${\displaystyle x_{i}}$ от ${\displaystyle L_{k}}$ минимальна:

${\displaystyle \sum _{i=1}^{m}\operatorname {dist} ^{2}(x_{i},L_{k})\to \min }$,

где ${\displaystyle \operatorname {dist} (x_{i},L_{k})}$ — евклидово расстояние от точки до линейного многообразия. Всякое ${\displaystyle k}$-мерное линейное многообразие в ${\displaystyle \mathbb {R} ^{n}}$ может быть задано как множество линейных комбинаций ${\displaystyle L_{k}=\{a_{0}+\beta _{1}a_{1}+\dots +\beta _{k}a_{k}|\beta _{i}\in \mathbb {R} \}}$, где параметры ${\displaystyle \beta _{i}}$ пробегают вещественную прямую ${\displaystyle \mathbb {R} }$, ${\displaystyle a_{0}\in \mathbb {R} ^{n}}$ а ${\displaystyle \left\{a_{1},\dots ,a_{k}\right\}\subset \mathbb {R} ^{n}}$ — ортонормированный набор векторов

${\displaystyle \operatorname {dist} ^{2}(x_{i},L_{k})=\Vert x_{i}-a_{0}-\sum _{j=1}^{k}a_{j}(a_{j},x_{i}-a_{0})\Vert ^{2}}$,

где ${\displaystyle \Vert \cdot \Vert }$ евклидова норма, ${\displaystyle \left(a_{j},x_{i}\right)}$ — евклидово скалярное произведение, или в координатной форме:

${\displaystyle \operatorname {dist} ^{2}(x_{i},L_{k})=\sum _{l=1}^{n}\left(x_{il}-a_{0l}-\sum _{j=1}^{k}a_{jl}\sum _{q=1}^{n}a_{jq}(x_{iq}-a_{0q})\right)^{2}}$.

Решение задачи аппроксимации для ${\displaystyle k=0,1,\dots ,n-1}$ даётся набором вложенных линейных многообразий ${\displaystyle L_{0}\subset L_{1}\subset \dots L_{n-1}}$, ${\displaystyle L_{k}=\{a_{0}+\beta _{1}a_{1}+\ldots +\beta _{k}a_{k}|\beta _{i}\in \mathbb {R} \}}$. Эти линейные многообразия определяются ортонормированным набором векторов ${\displaystyle \left\{a_{1},...,a_{n-1}\right\}}$ (векторами главных компонент) и вектором ${\displaystyle a_{0}}$. Вектор ${\displaystyle a_{0}}$ ищется как решение задачи минимизации для ${\displaystyle L_{0}}$:

${\displaystyle a_{0}={\underset {a_{0}\in \mathbb {R} ^{n}}{\operatorname {argmin} }}\left(\sum _{i=1}^{m}\operatorname {dist} ^{2}(x_{i},L_{0})\right)}$,

то есть

${\displaystyle a_{0}={\underset {a_{0}\in \mathbb {R} ^{n}}{\operatorname {argmin} }}\left(\sum _{i=1}^{m}\Vert x_{i}-a_{0}\Vert ^{2}\right)}$.

Это — выборочное среднее: ${\displaystyle a_{0}={\frac {1}{m}}\sum _{i=1}^{m}x_{i}={\overline {X}}}$.

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

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

1. Централизуются данные (вычитанием среднего): ${\displaystyle x_{i}:=x_{i}-{\overline {X}}}$. Теперь ${\displaystyle \sum _{i=1}^{m}x_{i}=0}$;
2. Отыскивается первая главная компонента как решение задачи:

   ${\displaystyle a_{1}={\underset {\Vert a_{1}\Vert =1}{\operatorname {argmin} }}\left(\sum _{i=1}^{m}\Vert x_{i}-a_{1}(a_{1},x_{i})\Vert ^{2}\right)}$.

   если решение не единственно, то осуществляется выбор одного из них.

3. Из данных вычитается проекция на первую главную компоненту:

   ${\displaystyle x_{i}:=x_{i}-a_{1}\left(a_{1},x_{i}\right)}$;

4. Отыскивается вторая главная компонента как решение задачи:

   ${\displaystyle a_{2}={\underset {\Vert a_{2}\Vert =1}{\operatorname {argmin} }}\left(\sum _{i=1}^{m}\Vert x_{i}-a_{2}(a_{2},x_{i})\Vert ^{2}\right)}$.

   Если решение не единственно, то выбирается одно из них.

Далее процесс продолжается, то есть на шаге ${\displaystyle 2k-1}$ вычитается проекция на ${\displaystyle (k-1)}$-ю главную компоненту (к этому моменту проекции на предшествующие ${\displaystyle (k-2)}$ главные компоненты уже вычтены):

${\displaystyle x_{i}:=x_{i}-a_{k-1}\left(a_{k-1},x_{i}\right)}$;

и на шаге ${\displaystyle 2k}$ определяется ${\displaystyle k}$-я главная компонента как решение задачи:

${\displaystyle a_{k}={\underset {\Vert a_{k}\Vert =1}{\operatorname {argmin} }}\left(\sum _{i=1}^{m}\Vert x_{i}-a_{k}(a_{k},x_{i})\Vert ^{2}\right)}$ (если решение не единственно, то выбирается одно из них).

На каждом подготовительном шаге ${\displaystyle (2k-1)}$ вычитается проекция на предшествующую главную компоненту. Найденные векторы ${\displaystyle \left\{a_{1},...,a_{n-1}\right\}}$ ортонормированы просто в результате решения описанной задачи оптимизации, однако чтобы не дать ошибкам вычисления нарушить взаимную ортогональность векторов главных компонент, можно включать ${\displaystyle a_{k}\bot \{a_{1},...,a_{k-1}\}}$ в условия задачи оптимизации.

Неединственность в определении ${\displaystyle a_{k}}$ помимо тривиального произвола в выборе знака (${\displaystyle a_{k}}$ и ${\displaystyle -a_{k}}$ решают ту же задачу) может быть более существенной и происходить, например, из условий симметрии данных. Последняя главная компонента ${\displaystyle a_{n}}$ — единичный вектор, ортогональный всем предыдущим ${\displaystyle a_{k}}$.

Поиск ортогональных проекций с наибольшим рассеянием

Первая главная компонента максимизирует выборочную дисперсию проекции данных

Пусть нам дан центрированный набор векторов данных ${\displaystyle x_{i}\in \mathbb {R} ^{n}\;(i=1,...,m)}$ (среднее арифметическое значение ${\displaystyle x_{i}}$ равно нулю). Задача — найти такое ортогональное преобразование в новую систему координат, для которого были бы верны следующие условия:

• Выборочная дисперсия данных вдоль первой координаты максимальна (эту координату называют первой главной компонентой);
• Выборочная дисперсия данных вдоль второй координаты максимальна при условии ортогональности первой координате (вторая главная компонента);
• …
• Выборочная дисперсия данных вдоль значений ${\displaystyle k}$-ой координаты максимальна при условии ортогональности первым ${\displaystyle k-1}$ координатам;
• …

Выборочная дисперсия данных вдоль направления, заданного нормированным вектором ${\displaystyle a_{k}}$, это

${\displaystyle S_{m}^{2}\left[(X,a_{k})\right]={\frac {1}{m}}\sum \limits _{i=1}^{m}(a_{k},x_{i})^{2}={\frac {1}{m}}\sum \limits _{i=1}^{m}\left(\sum \limits _{j=1}^{n}x_{ij}a_{kj}\right)^{2}}$

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

Решение задачи о наилучшей аппроксимации даёт то же множество главных компонент ${\displaystyle \left\{a_{i}\right\}}$, что и поиск ортогональных проекций с наибольшим рассеянием, по очень простой причине: ${\displaystyle \Vert x_{i}-a_{k}(a_{k},x_{i})\Vert ^{2}=\Vert x_{i}\Vert ^{2}-(a_{k},x_{i})^{2},}$ и первое слагаемое не зависит от ${\displaystyle a_{k}}$.

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

Ещё одна эквивалентная формулировка следует из очевидного тождества, верного для любых ${\displaystyle m}$ векторов ${\displaystyle x_{i}}$:

${\displaystyle {\frac {1}{m(m-1)}}\sum _{i,j=1}^{m}(x_{i}-x_{j})^{2}={\frac {2m^{2}}{m(m-1)}}\left[{\frac {1}{m}}\sum _{i=1}^{m}x_{i}^{2}-\left({\frac {1}{m}}\sum _{i}^{m}x_{i}\right)^{2}\right].}$

В левой части этого тождества стоит среднеквадратичное расстояние между точками, а в квадратных скобках справа — выборочная дисперсия. Таким образом, в методе главных компонент ищутся подпространства, в проекции на которые среднеквадратичное расстояние между точками максимально (или, что то же самое, его искажение в результате проекции минимально)^[2]. Такая переформулировка позволяет строить обобщения с взвешиванием различных парных расстояний (а не только точек).

Аннулирование корреляций между координатами

Для заданной ${\displaystyle n}$-мерной случайной величины ${\displaystyle X}$ найти такой ортонормированный базис, ${\displaystyle \left\{a_{1},...,a_{n}\right\}}$, в котором коэффициент ковариации между различными координатами равен нулю. После преобразования к этому базису

${\displaystyle \operatorname {cov} (X_{i},X_{j})=0}$ для ${\displaystyle i\neq j}$.

Здесь ${\displaystyle \operatorname {cov} (X_{i},X_{j})=\operatorname {E} [(X_{i}-\operatorname {E} [X_{i}])(X_{j}-\operatorname {E} [X_{j}])]}$ — коэффициент ковариации, где ${\displaystyle \operatorname {E} }$ — математическое ожидание.

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

${\displaystyle C=[c_{ij}],\ c_{ij}={\frac {1}{m-1}}\sum _{l=1}^{m}(x_{li}-{\overline {X_{i}}})(x_{lj}-{\overline {X_{j}}}).}$

Ковариационная матрица многомерной случайной величины ${\displaystyle X}$, это

${\displaystyle \Sigma =[\sigma _{ij}],\ \sigma _{ij}=\operatorname {cov} (X_{i},X_{j})=\operatorname {E} [(X_{i}-\operatorname {E} [X_{i}])(X_{j}-\operatorname {E} [X_{j}])].}$

Векторы главных компонент для задач о наилучшей аппроксимации и о поиске ортогональных проекций с наибольшим рассеянием — это ортонормированный набор ${\displaystyle \left\{a_{1},...,a_{n}\right\}}$ собственных векторов эмпирической ковариационной матрицы ${\displaystyle C}$, расположенных в порядке убывания собственных значений ${\displaystyle \lambda$ :\lambda _{1}\geq \lambda _{2}\geq \ldots \geq \lambda _{n}\geq 0.} Эти векторы служат оценкой для собственных векторов ковариационной матрицы ${\displaystyle \operatorname {cov} (X_{i},X_{j})}$. В базисе из собственных векторов ковариационной матрицы она, естественно, диагональна, и в этом базисе коэффициент ковариации между различными координатами равен нулю.

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

Основная статья: Сингулярное разложение

Идея сингулярного разложения

Математическое содержание метода главных компонент — это спектральное разложение ковариационной матрицы ${\displaystyle C}$, то есть представление пространства данных в виде суммы взаимно ортогональных собственных подпространств ${\displaystyle C}$, а самой матрицы ${\displaystyle C}$ — в виде линейной комбинации ортогональных проекторов на эти подпространства с коэффициентами ${\displaystyle \lambda _{i}}$. Если ${\displaystyle \operatorname {X} =\left\{x_{1},...,x_{m}\right\}^{T}}$ — матрица, составленная из векторов-строк (размерности ${\displaystyle n}$) центрированных данных, то ${\displaystyle C={\frac {1}{m-1}}\operatorname {X} ^{T}\operatorname {X} }$ и задача о спектральном разложении ковариационной матрицы ${\displaystyle C}$ превращается в задачу о сингулярном разложении матрицы данных ${\displaystyle \operatorname {X} }$.

Число ${\displaystyle \sigma \geq 0}$ называется сингулярным числом матрицы ${\displaystyle \operatorname {X} }$ тогда и только тогда, когда существуют правый и левый сингулярные векторы: такие ${\displaystyle m}$-мерный вектор-строка ${\displaystyle b_{\sigma }}$ и ${\displaystyle n}$-мерный вектор-столбец ${\displaystyle a_{\sigma }}$ (оба единичной длины), что выполнено два равенства:

${\displaystyle \operatorname {X} a_{\sigma }=\sigma b_{\sigma }^{T};\,\,b_{\sigma }\operatorname {X} =\sigma a_{\sigma }^{T}.}$

Пусть ${\displaystyle p=\operatorname {rang} \operatorname {X} \leq \min\{n,m\}}$ — ранг матрицы данных. Сингулярное разложение матрицы данных ${\displaystyle \operatorname {X} }$ — это её представление в виде

${\displaystyle \operatorname {X} =\sum _{l=1}^{p}\sigma _{l}b_{l}^{T}a_{l}^{T};\;\operatorname {X} ^{T}=\sum _{l=1}^{p}\sigma _{l}a_{l}b_{l}\;\left(x_{ij}=\sum _{l=1}^{p}\sigma _{l}b_{li}a_{lj}\right),}$

где ${\displaystyle \sigma _{l}>0}$ — сингулярное число, ${\displaystyle a_{l}=(a_{lj}),\,j=1,...n}$ — соответствующий правый сингулярный вектор-столбец, а ${\displaystyle b_{l}=(b_{li}),\,i=1,...m}$ — соответствующий левый сингулярный вектор-строка (${\displaystyle l=1,...p}$). Правые сингулярные векторы-столбцы ${\displaystyle a_{l}}$, участвующие в этом разложении, являются векторами главных компонент и собственными векторами эмпирической ковариационной матрицы ${\displaystyle C={\frac {1}{m-1}}\operatorname {X} ^{T}\operatorname {X} }$, отвечающими положительным собственным числам ${\displaystyle \lambda _{l}={\frac {1}{m-1}}\sigma _{l}^{2}>0}$.

Хотя формально задачи сингулярного разложения матрицы данных и спектрального разложения ковариационной матрицы совпадают, алгоритмы вычисления сингулярного разложения напрямую, без вычисления ковариационной матрицы и её спектра, более эффективны и устойчивы^[3].

Теория сингулярного разложения была создана Джеймсом Джозефом Сильвестром в 1889 году и изложена во всех подробных руководствах по теории матриц^[4].

Простой итерационный алгоритм сингулярного разложения

Основная процедура — поиск наилучшего приближения произвольной ${\displaystyle m\times n}$ матрицы ${\displaystyle X=(x_{ij})}$ матрицей вида ${\displaystyle b\otimes a=(b_{i}a_{j})}$ (где ${\displaystyle b}$