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

---

Введение

Часто n-мерное пространство признаков оказывается избыточным. Объекты (точки) концентрируются в нём на вложенных m-мерных поверхностях, где m существенно меньше n. Выявление таких поверхностей приводит к снижению размерности задачи и упрощению дальнейшего анализа данных (в задачах классификации и т.п.). Метод главных компонент (principal component analysis) является простейшим методом, в котором, в качестве m-мерной поверхности выступает m-плоскость.

Напомним, что $m$ векторов $\mathbf{e}_k=\{\mathbf{e}_1,...,\mathbf{e}_m\}$ задают уравнение $m$-плоскости в $n$-мерном пространстве: $$\begin{equation}\label{plane_eq} \mathbf{x}= \mathbf{x}_0 + \sum^m_{k=1} t_k\mathbf{e}_k. \end{equation}$$ Параметры $\{t_1,...,t_m\}$ выступают декартовыми координатами в базисе $\{\mathbf{e}_1,...,\mathbf{e}_m\}$ (если он ортогонален). Вектор $\mathbf{x}_0$ задаёт положение начала координат (точка в которой все координаты на плоскости $t_k$ равны нулю). Если $m=1$, то ($\ref{plane_eq}$) будет уравнением прямой (одномерный объект). При $m=2$ мы имеем двумерную плоскость и т.д.

Квадрат кратчайшего евклидового расстояния $d^2=(\mathbf{X}-\mathbf{x})^2$ от произвольной точки $\mathbf{X}$ до $m$-плоскости ($\ref{plane_eq}$) равно $$\begin{equation}\label{d2_to_plane} d^2 = \bigr(\mathbf{X}-\mathbf{x}_0\bigr)^2 - \sum^m_{k=1} \bigr[(\mathbf{X}-\mathbf{x}_0)\,\mathbf{e}_k \bigr]^2, \end{equation}$$ а проекция точки $\mathbf{X}$ на $m$-плоскости имеет следующие координаты (в $n$- и $m$- пространствах): $$\begin{equation}\label{projec_coord} \mathbf{X} ~\mapsto~\mathbf{x}=\mathbf{x}_0 + \sum^m_{k=1} \bigr[(\mathbf{X}-\mathbf{x}_0)\,\mathbf{e}_k\bigr]\mathbf{e}_k, ~~~~~~~~~~~~~ t_k = (\mathbf{X}-\mathbf{x}_0)\,\mathbf{e}_k. \end{equation}$$

---

Центр тяжести

В методе главных компонент ищут $m$-плоскость, сумма квадратов расстояний от которой до множества точек (обучающий набор объектов) $\mathbf{X}_1,...,\mathbf{X}_N$ минимальна (нижний индекс - номер объекта). Будем в дальнейшем обозначать угловыми скобками усреднение по всем точкам. Например, среднее значение (центр тяжести $N$ объектов) равен: $$\langle \mathbf{X} \rangle = \frac{1}{N}\sum^N_{i=1} \mathbf{X}_i.$$ Таким образом, мы требуем, чтобы $$\langle d^2 \rangle = \langle\bigr(\mathbf{X}-\mathbf{x}_0\bigr)^2\bigr\rangle - \sum^m_{k=1} \bigr\langle\bigr[(\mathbf{X}-\mathbf{x}_0)\,\mathbf{e}_k \bigr]^2\bigr\rangle = min.$$ Определим сначала $\mathbf{x}_0$. Беря производную по $\mathbf{x}_0$ и приравнивая её нулю, имеем: $$\langle\mathbf{X}-\mathbf{x}_0\rangle - \sum^m_{k=1} \bigr[(\langle\mathbf{X}\rangle-\mathbf{x}_0)\,\mathbf{e}_k \bigr]\,\mathbf{e}_k = 0.$$ Это соотношение тождественно выполняется, если $\mathbf{x}_0=\langle\mathbf{X}\rangle$, т.е. плоскость проходит через центр тяжести. Поэтому вычтем его из каждой точки: $\mathbf{X}_i \mapsto \mathbf{X}_i - \langle \mathbf{X} \rangle$ и далее будем считать, что $\langle \mathbf{X} \rangle=0$ и $\mathbf{x}_0=0$.

---

Первый компонент

Для определения направляющих векторов $\mathbf{e}_1,...,\mathbf{e}_m$ (базис в плоскости) необходимо найти максимум (знак минус в ($\ref{d2_to_plane}$)!) выражения: $$\begin{equation}\label{maxD} D = \sum^m_{k=1} \bigr\langle\bigr(\mathbf{X}\,\mathbf{e}_k \bigr)^2\bigr\rangle = \bigr\langle\bigr(\mathbf{X}\,\mathbf{e}_1 \bigr)^2\bigr\rangle+...+\bigr\langle\bigr(\mathbf{X}\,\mathbf{e}_m \bigr)^2\bigr\rangle = max,~~~~~~~~~~~\mathbf{e}_k\mathbf{e}_p = \delta_{kp}, \end{equation}$$ при сохранении ортонормированности векторов. Максимальным должно быть каждое слагаемое в сумме. Чем больше значение слагаемого, тем более значимым является $k$-й вектор. Пусть эти слагаемые отсортированы по убыванию значимости. Определим сначала $\mathbf{e}_1$, который называют первым компонентом. В плоскости он может иметь произвольное направление (произвол ориентации базиса в плоскости). Поэтому стоит задача поиска максимума $\bigr\langle\bigr(\mathbf{X}\,\mathbf{e}_1 \bigr)^2\bigr\rangle$ при условии, что $\mathbf{e}^2_1=1.$ Эту задачу можно решать методом множителей Лагранжа: $$\begin{equation}\label{max_eqX} \bigr\langle\bigr(\mathbf{X}\,\mathbf{e}_1 \bigr)^2\bigr\rangle + \lambda\,(1-\mathbf{e}^2_1) = max. \end{equation}$$ Беря производные по $\mathbf{e}_1$, $\lambda$ и приравнивая их нулю, имеем: $$\begin{equation}\label{e1} \bigr\langle\bigr(\mathbf{X}\,\mathbf{e}_1 \bigr)\mathbf{X}\bigr\rangle = \lambda\,\mathbf{e}_1,~~~~~~~~~~~~~\mathbf{e}^2_1 = 1. \end{equation}$$ Обозначим в круглых скобках номер компоненты вектора в $n$-мерном пространстве: $\mathbf{e}_1=\{e^{(1)}_1,...,e^{(n)}_1\}$. Первое уравнение переписывается в форме уравнения на собственные значения: $$\begin{equation}\label{eq_Ce} \sum^n_{j=1}C_{ij}\,e^{(j)}_1 = \lambda\,e^{(i)}_1 \end{equation}$$ ковариационной матрицы: $$\begin{equation}\label{eq_cov} C_{ij} = \langle X^{(i)}X^{(j)}\bigr\rangle = \frac{1}{N-1}\sum^N_{k=1} X^{(i)}_kX^{(j)}_k. \end{equation}$$

Из ($\ref{e1}$) следует (умножаем первое соотношение на $\mathbf{e}_1$), что $\lambda=\bigr\langle\bigr(\mathbf{X}\,\mathbf{e}_1 \bigr)^2\bigr\rangle$, поэтому из всех собственных значений необходимо выбрать максимальное. Впрочем, это не обязательно. Если будут найдены все собственные значения матрицы $C_{ij}$ (см. ниже), то их можно отсортировать в порядке убывания, оставив только $m$ самых больших. Соответствующие им собственные векторы ковариационной матрицы и являются искомыми направляющими векторами $m$-плоскости.

---

Второй компонент

Пусть первый компонент $\mathbf{e}_1$ и первое собственное значение $\lambda_1=\lambda$ матрицы $C_{ij}$ уже найдены. Изменим координаты точек следующущим образом: $$\begin{equation}\label{tranform} \mathbf{X}' = \mathbf{X} - (\mathbf{X}\mathbf{e}_1)\,\mathbf{e}_1,~~~~~~~~~\mathbf{X}'\mathbf{e}_1=0. \end{equation}$$ Среднее значение останется нулевым $\langle\mathbf{X}'\rangle=0$. Остаток максимизируемой суммы ($\ref{maxD}$) также не поменяется: $$D' =\bigr\langle\bigr(\mathbf{X}\,\mathbf{e}_2 \bigr)^2\bigr\rangle+...+\bigr\langle\bigr(\mathbf{X}\,\mathbf{e}_m \bigr)^2\bigr\rangle =\bigr\langle\bigr(\mathbf{X}'\,\mathbf{e}_2 \bigr)^2\bigr\rangle+...+\bigr\langle\bigr(\mathbf{X}'\,\mathbf{e}_m \bigr)^2\bigr\rangle =max.$$ Его необходимо продолжать максимизировать. Заметим, что пересчитывать ковариационную матрицу ($\ref{eq_cov}$) по всем точкам $\mathbf{X}'$ не обязательно, так после преобразования ($\ref{tranform}$) она меняется следующим образом: $$C'_{ij} = C_{ij} - \sum^{n}_{k=1} C_{ki}\, e^{(k)}_1e^{(j)}_1 - \sum^{n}_{k=1} C_{kj}\, e^{(k)}_1e^{(i)}_1 + \sum^{n}_{k=1} C_{kp}\, e^{(k)}_1e^{(p)}_1 \,e^{(i)}_1e^{(j)}_1.$$ Учитывая уравнение на собственные значения ($\ref{eq_Ce}$) это соотношение можно упростить: $$\begin{equation} C'_{ij} = C_{ij} - \lambda_1 \, e^{(i)}_1e^{(j)}_1. \end{equation}$$ Теперь процесс поиска максимума можно повторять до тех пор, пока все векторы базиса плоскости не будут найдены. При этом, в силу построения по формуле ($\ref{tranform}$), все эти векторы окажутся ортогональными.

Несложно видеть, собственные значения матрицы $C'_{ij}$ одновременно являются собственными значениями исходной матрицы $C_{ij}$. Поэтому $m$-плоскость строиться по $m$ собственным векторами ковариационной матрицы, соответствующим максимальным собственным значениям.

Ещё одна интерпретация собственных векторов ковариационной матрицы - это оси эллипсоида характерного "разброса" точек в пространстве. Если эллипсоид вытянут вдоль одной из осей - это главная (наиболее значимая) компонента. Следующая по "вытянутости" ось - вторая компонента и т.д.

---

Пример 1: ирисы Фишера

Рассмотрим в качестве примера набор из 150 измерений длины и ширины лепестков цветков ириса (ирисы Фишера). В данном случае пространство признаков 4-мерно и в задаче классификации присутствуют 3 класса.

Статистики по всем признакам имеют вид:

xi	 aver	 sigma	  min	  max
x1	 5.843	 0.828	 4.30	 7.90
x2	 3.054	 0.434	 2.00	 4.40
x3	 3.759	 1.764	 1.00  6.90
x4	 1.199	 0.763	 0.10  2.50

После нормировки $x_i \mapsto (x_i-\langle x_i\rangle)/\sigma_i$ ковариационная матрица $C_{ij}$ равна:

 1.000	-0.109	 0.872	 0.818
-0.109	 1.000	-0.421	-0.357
 0.872	-0.421	 1.000	 0.963
 0.818	-0.357	 0.963	 1.000

Её собственные значения и собственные векторы:

λ1=2.911;   e1=(-0.522,  0.263,	-0.581,	-0.566)
λ2=0.921;   e2=(-0.372, -0.926,	-0.021,	-0.065)
λ3=0.147;   e3=( 0.721, -0.242,	-0.141,	-0.634)
λ4=0.021;   e4=(-0.262,  0.124,  0.801,	-0.524)

Главная компонента (наибольшее собственное значение λ₁=2.911) cущественно отличается от остальных. Это свидетельствует о том, что точки в 4-мерном пространстве в основном вытянуты вдоль прямой. Следующее по значимости собственное значение λ₂=0.921 и соответствующий ему собственный вектор e₂, вместе с вектором e₁ образует двумерную плоскость. Справа на рисунке приведена проекция данных на эту плоскость в диапазоне [-4,..,4] по каждой координате. Вектор $\mathbf{e}_1$ направлен вправо, а $\mathbf{e}_2$ - вниз. Статистики координат на этой плоскости имеют вид:

ti	 aver	sigma	   min	  max
t1	0.000	1.706	-3.298	2.765
t2	0.000	0.960	-2.713	2.649

Синие круги - первый класс, красные - второй и зелёные - третий.

---

Пример 2: рукописные цифры MNIST

База рукописных цифр MNIST содержит содержит 60000 обучающих 10000 тестовых изображений 28x28 пикселей в градациях серого. Все картинки отцентрованы по центру тяжести и более или менее приведены к единому масштабу.

Уменьшим цифры в 2 раза (картинки 14×14) и посчитаем ковариационную матрицу (в её графическом представлении красный цвет +1, синий -1). Собственные значения ковариационной матрицы приведены справа от неё. Видно, что размерность пространства признаков существенно меньше чем 14·14=196

PCA 2-мерная ($m=2$) плоскость для MNIST имеет вид (цветом помечены различные классы=цифры, расшифровка которых приведена в верхней части рисунка):

---

Проблемы PCA в задачах классификации

Метод главных компонент позволяет выяснить характерную размерность данных, т.е. размерность поверхности, вложенной в пространство признаков. Однако, проектирование данных на эту $m$-мерную плоскость может существенно ухудшить классификацию классов. Справа на рисунке приведены два класса (круги и крестики). Они являются хорошо (линейно) разделимыми. Однако при проектировании на линию (одна главная компонента), объекты различных класов перемешаются и перестанут быть разделимыми.

Например, в задаче ирисов Фишера метод главных компонент свидетельствует о "двумерности" исходных данных. Однако, первая главная компонента направлена не в направлении наиболее значимого (четвёртого) признака. Ниже, слева от вертикальной линии приведены исходные данные ирисов Фишера во всех парах координат. Справа те-же данные после устранения корреляции (т.е. PCA c $m=4$):

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

---