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

Формальное описание

Пусть имеется матрица переменных $\mathbf{X}$ размерностью ($I×J$), где $I$ – число образцов (строк), а $J$ – это число независимых переменных (столбцов), которых, как правило, много ($J>>1$). В методе главных компонент используются новые, формальные переменные $\mathbf{t}_a$ ($a=1,\dots,A$), являющиеся линейной комбинацией исходных переменных $\mathbf{x}_j$ ($j=1,\dots,J$):

$$ \mathbf{t}_a = \mathbf{p}_{a1} \mathbf{x}_1 + \dots + \mathbf{p}_{aJ} \mathbf{x}_J $$

С помощью этих новых переменных матрица X разлагается в произведение двух матриц $\mathbf{T}$ и $\mathbf{P}$:

$$ \mathbf{X} = \mathbf{TP}^\mathrm{t} + \mathbf{E} = \sum_{a=1}^A \mathbf{t}_a\mathbf{p}_a^\mathrm{t} + \mathbf{E} $$

Матрица $\mathbf{T}$ называется матрицей счетов (scores). Ее размерность ($I×A$).

Матрица $\mathbf{P}$ называется матрицей нагрузок (loadings). Ее размерность ($J×A$).

$\mathbf{E}$ – это матрица остатков, размерностью ($I×J$).

Новые переменные $\mathbf{t}_a$ называются главными компонентами (Principal Components), поэтому и сам метод называется методом главных компонент (PCA). Число столбцов $\mathbf{t}_a$ в матрице $\mathbf{T}$, и $\mathbf{p}_a$ в матрице $\mathbf{P}$, равно $A$, которое называется числом главных компонент (PC). Эта величина заведомо меньше числа переменных $J$ и числа образцов $I$.

Важным свойством PCA является ортогональность (независимость) главных компонент. Поэтому матрица счетов $\mathbf{T}$ не перестраивается при увеличении числа компонент, а к ней просто прибавляется еще один столбец – соответствующий новому направлению. Тоже происходит и с матрицей нагрузок $\mathbf{P}$.

Алгоритм NIPALS

Чаще всего для построения PCA счетов и нагрузок, используется рекуррентный алгоритм NIPALS, который на каждом шагу вычисляет одну компоненту. Сначала исходная матрица $\mathbf{X}$ преобразуется (как минимум – центрируется) и превращается в матрицу $\mathbf{E}_0$, $a=0$. Далее применяют следующий алгоритм.

1. Выбрать начальный вектор $\mathbf{t}$
2. Расчитать вектор нагрузок: $\mathbf{p}^\mathrm{t} = \mathbf{t}^\mathrm{t} \mathbf{E}_a / \mathbf{t}^\mathrm{t}\mathbf{t}$
3. Сделать нормировку: $\mathbf{p} = \mathbf{p} / (\mathbf{p}^\mathrm{t}\mathbf{p})^½$
4. Расчитать вектор счетов: $\mathbf{t} = \mathbf{E}_a \mathbf{p} / \mathbf{p}^\mathrm{t}\mathbf{p}$
5. Проверить сходимость, если нет, то идти на 2

После вычисления очередной ($a$-ой) компоненты, полагаем $\mathbf{t}_a=\mathbf{t}$ и $\mathbf{p}_a=\mathbf{p}$. Для получения следующей компоненты надо вычислить остатки $\mathbf{E}_{a+1} = \mathbf{E}_a – \mathbf{t}\mathbf{p}^\mathrm{t}$ и применить к ним тот же алгоритм, заменив индекс $a$ на $a+1$.

В этом пособии для построения PCA используется специальная надстройка для программы Excel (Add–In) Chemometrics.xla. Она дополняет список стандартных функций Excel и позволяет проводить PCA разложение на листах рабочей книги.

После того, как построено пространство из главных компонент, новые образцы $\mathbf{X}_{new}$ могут быть на него спроецированы, иными словами – определены матрицы их счетов $\mathbf{T}_{new}$. В методе PCA это делается очень просто

$$\mathbf{T}_{new} = \mathbf{X}_{new}\mathbf{P}$$

PCA и SVD

Метод главных компонент тесно связан с другим разложением – по сингулярным значениям, SVD. В последнем случае исходная матрица $\mathbf{X}$ разлагается в произведение трех матриц

$$\mathbf{X}=\mathbf{U}\mathbf{S}\mathbf{V}^\mathrm{t}$$

Здесь $\mathbf{U}$ – матрица, образованная ортонормированными собственными векторами $\mathbf{u}_r$ матрицы $\mathbf{XX}^\mathrm{t}$, соответствующим значениям $λ_r$:

$$\mathbf{XX}^\mathrm{t}\mathbf{u}_r = λ_r\mathbf{u}_r$$

$\mathbf{V}$ – матрица, образованная ортонормированными собственными векторами $\mathbf{v}_r$ матрицы $\mathbf{X}^\mathrm{t}\mathbf{X}$:

$$\mathbf{X}^\mathrm{t}\mathbf{X}\mathbf{v}_r = λ_r\mathbf{v}_r$$

$\mathbf{S}$ – положительно определенная диагональная матрица, элементами которой являются сингулярные значения $\sigma_1 ≥ \dots ≥ \sigma_R ≥ 0$ равные квадратным корням из собственных значений $λ_r$: $\sigma_r = \sqrt{\lambda_r}$.

Связь между PCA и SVD определяется следующими простыми соотношениями

$$\mathbf{T} = \mathbf{US}$$ $$\mathbf{P} = \mathbf{V}$$

Счета

Матрица счетов $\mathbf{T}$ дает нам проекции исходных образцов ($J$–мерных векторов $\mathbf{x}_1,\dots,\mathbf{x}_I$) на подпространство главных компонент ($A$-мерное). Строки $\mathbf{t}_1,\dots,\mathbf{t}_I$ матрицы $\mathbf{T}$ – это координаты образцов в новой системе координат. Столбцы $\mathbf{t}_1,\dots,\mathbf{t}_A$ матрицы $\mathbf{T}$ – ортогональны и представляют проекции всех образцов на одну новую координатную ось.

При исследовании данных методом PCA, особое внимание уделяется графикам счетов. Они несут в себе информацию, полезную для понимания того, как устроены данные. На графике счетов каждый образец изображается в координатах ($\mathbf{t}_i, \mathbf{t}_j$), чаще всего – ($\mathbf{t}_1, \mathbf{t}_2$), обозначаемых PC1 и PC2. Близость двух точек означает их схожесть, т.е. положительную корреляцию. Точки, расположенные под прямым углом, являются некоррелироваными, а расположенные диаметрально противоположно – имеют отрицательную корреляцию.

Подробнее о том, как из графиков счетов извлекается полезная информация, будет рассказано в примере.

Для матрицы счетов имеют место следующие соотношения:

$$\mathbf{T}^\mathrm{t}\mathbf{T} = \mathbf{\Lambda} = \mathrm{diag}\{λ_1,\dots,λ_A\}$$

где величины $λ_1≥\dots≥λ_A≥0$ – это собственные значения. Они характеризуют важность каждой компоненты

Нулевое собственное значение $λ_0$ определяется как сумма всех собственных значений, т.е.

Для вычисления PCA-счетов в надстройке Chemometrics Add-In используется функция ScoresPCA, которая имеет следующий синтаксис:

ScoresPCA (X [, PC] [,CentWeightX] [, Xnew])

• X – множество значений $\mathbf{X}$ (калибровочный набор)
• PC – необязательный аргумент. Целая величина, определяющая число главных компонент ($A$), используемых в методе главных компонент.
• CentWeightX – необязательный аргумент. Целая величина, определяющая проводится ли центрирование и/или шкалирование X переменных.
• Xnew – необязательный аргумент. Множество новых значений $\mathbf{X}_{new}$ (проверочный набор) для которых вычисляются и выводятся значения счетов $\mathbf{T}$.

Матрица новых значений Xnew должна иметь столько же столбцов (переменных), как и матрица X. Если аргумент Xnew опущен, то предполагается, что он совпадает с массивом X. Результат является массивом (матрицей), в которой число строк равно числу строк (образцов, $I$) в массиве Xnew, а число столбцов равно числу ГК ($A$).

Аналогичная стандартная функция листа: ТЕНДЕНЦИЯ.

Функция ScoresPCA является формулой массива и ее ввод должен завершаться нажатием комбинации CTRL+SHIFT+ENTER.

Нагрузки

Матрица нагрузок $\mathbf{}P$ – это матрица перехода из исходного пространства переменных $\mathbf{x}_1,\dots,\mathbf{x}_J$ ($J$-мерного) в пространство главных компонент ($A$-мерное). Каждая строка матрицы $\mathbf{P}$ состоит из коэффициентов, связывающих переменные $\mathbf{t}$ и $\mathbf{x}$. Например, $a$-я строка – это проекция всех переменных $\mathbf{x}_1,\dots,\mathbf{x}_J$ на $a$-ю ось главных компонент. Каждый столбец $\mathbf{P}$ – это проекция соответствующей переменной $\mathbf{x}_j$ на новую систему координат.

График нагрузок применяется для исследования роли переменных. На этом графике каждая переменная $\mathbf{x}_j$ отображается точкой в координатах ($\mathbf{p}_i, \mathbf{p}_j$), например ($\mathbf{p}_1, \mathbf{p}_2$). Анализируя его аналогично графику счетов, можно понять, какие переменные связаны, а какие независимы. Совместное исследование парных графиков счетов и нагрузок, также может дать много полезной информации о данных.

В методе главных компонент нагрузки – это ортогональные нормированные вектора, т.е.

$$\mathbf{P}^\mathrm{t}\mathbf{P} = \mathbf{I}$$

Для вычисления PCA-нагрузок в надстройке Chemometrics Add-In используется функция LoadingsPCA:

LoadingsPCA (X [, PC] [, CentWeightX])

• X – множество значений X (калибровочный набор)
• PC – необязательный аргумент. Целая величина, определяющая число главных компонент ($A$), используемых в методе главных компонент.
• CentWeightX – необязательный аргумент. Целая величина, определяющая проводится ли центрирование и/или шкалирование X переменных.

Результат является массивом (матрицей), в которой число столбцов равно числу ГК ($A$), а число строк равно числу столбцов (переменных, $J$) в массиве X. Аналогичная стандартная функция листа: МОБР

Функция LoadingsPCA является формулой массива и ее ввод должен завершаться нажатием комбинации CTRL+SHIFT+ENTER.

Данные специального вида

Результат моделирования методом главных компонент не зависит от порядка, в котором следуют образцы и/или переменные. Иными словами строки и столбцы в исходной матрице $\mathbf{X}$ можно переставить, но ничего принципиально не изменится. Однако, в некоторых случаях, сохранять и отслеживать этот порядок очень полезно – это позволяет лучше понять устройство моделируемых данных.

Рассмотрим простой пример – моделирование данных, полученных методом высокоэффективной жидкостной хроматографией с детектированием на диодной матрице (ВЭЖХ–ДДМ). Данные представляются матрицей, размерностью 30 образцов ($I$) на 28 переменных ($J$). Образцы соответствуют временам удерживания от 0 до 30 с, а переменные – длинам волн от 220 до 350 нм, на которых происходит детектирование. Данные ВЭЖХ–ДДМ представлены на Рис 10.

Эти данные хорошо моделируются методом PCA с двумя главными компонентами. Ясно, что в этом примере нам важен порядок, в котором идут образцы и переменные – он задается естественным ходом времени и спектральным диапазоном. Полученные счета и нагрузки полезно изобразить на графиках в зависимости от соответствующего параметра – счета от времени, а нагрузки от длины волны. (см. Рис 11)

Подробнее этот пример разобран в пособии Разрешение многомерных кривых.

Погрешности

PCA декомпозиция матрицы $\mathbf{X}$ является последовательным, итеративным процессом, который можно оборвать на любом шаге $a=A$. Получившаяся матрица

$$\hat{\mathbf{X}} = \mathbf{TP}^\mathrm{t}$$

вообще говоря, отличается от матрицы $\mathbf{X}$. Разница между ними

$$\mathbf{E} = \mathbf{X} - \hat{\mathbf{X}}$$

называется матрицей остатков.

Рассмотрим геометрическую интерпретацию остатков. Каждый исходный образец $\mathbf{x}_i$ (строка в матрице $\mathbf{X}$) можно представить как вектор в $J$–мерном пространстве с координатами

$$\mathbf{x}_i = (x_{i1}, x_{i1}, \dots, x_{iJ})$$

PCA проецирует его в вектор, лежащий в пространстве главных компонент, $\mathbf{t}_i=(t_{i1}, t_{i2},\dots,t_{iA})$ размерностью $A$. В исходном пространстве этот же вектор $\mathbf{t}_i$ имеет координаты

$$\hat{\mathbf{x}_i} = (\hat{x}_{i1}, \hat{x}_{i1}, \dots, \hat{x}_{iJ})$$

Разница между исходным вектором и его проекцией является вектором остатков

$$\mathbf{e} = \mathbf{x} - \hat{\mathbf{x}_i} = (e_{i1}, e_{i1}, \dots, e_{iJ})$$

Он образует $i$–ю строку в матрице остатков $\mathbf{E}$.

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

Величина –

$$v_i = \sum_{j = 1}^J e^2_{ij}$$

определяет квадрат отклонения исходного вектора $\mathbf{x}_i$ от его проекции на пространство PC. Чем оно меньше, тем лучше приближается $i$–ый образец. Для вычисления отклонений можно использовать стандартные функции листа или специальную пользовательскую функцию.

Эта же величина, деленная на число переменных

$$d_i = \frac{1}{J} \sum_{j = 1}^J e^2_{ij}$$

дает оценку дисперсии (вариации) $i$–го образца.

Среднее (для всех образцов) расстояние $v_0$ вычисляется как

$$v_0 = \frac{1}{I} \sum_{i = 1}^I v_i$$

Оценка общая (для всех образцов) дисперсии вычисляется так –

$$D = \frac{1}{I} \sum_{i = 1}^I d_i$$

Проверка

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

В некоторых случаях объем данных слишком мал для такой проверки. Тогда применяют другой метод – перекрестной проверки (кросс–валидация). Используется также проверка методом коррекции размахом, суть которой предлагается изучить самостоятельно.

"Качество" декомпозиции

Результатом PCA моделирования являются величины $\hat{\mathbf{X}}_c$ – оценки, найденные по модели, построенной на обучающем наборе $\mathbf{X}_c$. Результатом проверки служат величины $\hat{\mathbf{X}}_t$ – оценки проверочных значений $\mathbf{X}_t$, вычисленные по той же модели, но как новые образцы. Отклонение оценки от проверочного значения вычисляют как матрицу остатков:

в обучении: $\mathbf{E}_c = \mathbf{X}_c - \hat{\mathbf{X}}_c$

и в проверке: $\mathbf{E}_t = \mathbf{X}_t - \hat{\mathbf{X}}_t$

Следующие величины характеризуют "качество" моделирования в среднем.

Полная дисперсия остатков в обучении ($TRVC$) и в проверке ($TRVP$) –

$$TRVC = \frac{1}{I_c J} \sum_{i=1}^{I_c} \sum_{j=1}^J e_{ij}^2 = D_c$$ $$TRVP = \frac{1}{I_t J} \sum_{i=1}^{I_t} \sum_{j=1}^J e_{ij}^2 = D_t$$

Полная дисперсия выражается в тех же единицах (точнее их квадратах), что и исходные величины X.

Объясненная дисперсия остатков в обучении ($ERVC$) и в проверке ($ERVP$)

$$ERVC = 1 - \sum_{i=1}^{I_c} \sum_{j=1}^J e_{ij}^2 \bigg/ \sum_{i=1}^{I_c} \sum_{j=1}^J x_{ij}^2$$ $$ERVP = 1 - \sum_{i=1}^{I_t} \sum_{j=1}^J e_{ij}^2 \bigg/ \sum_{i=1}^{I_t} \sum_{j=1}^J x_{ij}^2$$

Объясненная дисперсия – это относительная величина. При ее вычислении используется естественная нормировка – сумма квадратов всех исходных величин $x_{ij}$. Обычно она выражается в процентах или в долях единицы. Во всех этих формулах величины $e_{ij}$ – это элементы матриц $\mathbf{E}_c$ или $\mathbf{E}_t$. Для характеристик, наименование которых оканчивается на $C$ (например, $TRVC$), используется матрица $\mathbf{E}_c$ (обучение), а для тех, которые оканчиваются на $P$ (например, $TRVP$), берется матрица $\mathbf{E}_t$ (проверка).

Выбор числа главных компонент

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

Для выбора значения числа главных компонент обычно используется график, на котором объясненная дисперсия ($ERV$) изображается в зависимости от числа PC. Пример такого графика приведен на Рис. 15.

Из этого графика видно, что правильное число PC – это 3 или 4. Три компоненты объясняют 95%, а четыре 98% исходной вариации. Окончательное решение о величине $A$ можно принять только после содержательного анализа данных.

Другим полезным инструментом является график, на котором изображаются собственные значения в зависимости от числа PC. Пример показан на Рис.16.

Из этого рисунка опять видно, что для $a=3$ происходит резкое изменение формы графика – излом. Поэтому верное число PC – это три или четыре.

Неединственность PCA

Разложение по методу главных компонент

$$\mathbf{X} = \mathbf{TP}^\mathrm{t}+\mathbf{E}$$

не является единственным. Вместо матриц $\mathbf{T}$ и $\mathbf{P}$ можно использовать другие матрицы $\tilde{\mathbf{T}}$ и $\tilde{\mathbf{P}}$, которые дадут аналогичную декомпозицию:

$$\mathbf{X} = \tilde{\mathbf{T}}\tilde{\mathbf{P}}^\mathrm{t}+\mathbf{E}$$

с той же матрицей ошибок $\mathbf{E}$. Простейший пример – это одновременное изменение знаков у соответствующих компонент векторов $\mathbf{t}_a$ и $\mathbf{p}_a$, при котором произведение

$$\tilde{\mathbf{t}}_a\tilde{\mathbf{p}}_a^\mathrm{t} = \mathbf{t}_a\mathbf{p}_a^\mathrm{t}$$

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

Более сложный случай – это одновременное вращение матриц $\mathbf{T}$ и $\mathbf{P}$. Пусть $\mathbf{R}$ – это ортогональная матрица вращения размерностью $A×A$ , т.е такая матрица, что $\mathbf{R}^\mathrm{t}=\mathbf{R}^{–1}$. Тогда

$$ \mathbf{TP}^\mathrm{t} = \mathbf{TRR}^{-1}\mathbf{P}^\mathrm{t} = (\mathbf{TR})(\mathbf{PR})^\mathrm{t} = \tilde{\mathbf{T}}\tilde{\mathbf{P}}^\mathrm{t} $$

Заметим, что новые матрицы счетов и нагрузок сохраняют все свойства старых,

$$\tilde{\mathbf{T}}^\mathrm{t}\tilde{\mathbf{T}} = \mathbf{T}^\mathrm{t}\mathbf{T} = \mathbf{\Lambda}$$ $$\tilde{\mathbf{P}}^\mathrm{t}\tilde{\mathbf{P}} = \mathbf{P}^\mathrm{t}\mathbf{P} = \mathbf{I}$$

Это свойство PCA называется вращательной неопределенностью. Оно интенсивно используется при решении задач разделения кривых, в частности методом прокрустова вращения. Если отказаться от условий ортогональности главных компонент, то декомпозиция матрицы станет еще более общей. Пусть теперь $\mathbf{R}$ – это произвольная невырожденная матрица размерностью $A×A$. Тогда

$$\tilde{\mathbf{T}} = \mathbf{TR}$$ $$\tilde{\mathbf{P}} = \mathbf{P}(\mathbf{R}^\mathrm{t})^{-1}$$

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

Подготовка данных

Во многих случаях, перед применением PCA, исходные данные нужно предварительно подготовить: отцентрировать и/или отнормировать. Эти преобразования проводятся по столбцам – переменным.

Центрирование – это вычитание из каждого столбца $\mathbf{x}_j$ среднего (по столбцу) значения:

$$m_j = (x_{1j} + \dots + x_{Ij}) / I$$

Центрирование необходимо потому, что оригинальная PCA модель, $\mathbf{X} = \mathbf{TP}^\mathrm{t}+\mathbf{E}$, не содержит свободного члена.

Второе простейшее преобразование данных – это нормирование. Это преобразование выравнивает вклад разных переменных в PCA модель. При этом преобразовании каждый столбец $\mathbf{x}_j$ делится на свое стандартное отклонение.

$$s_j = \sqrt{\frac{1}{I} \sum_{i=1}^I (x_{ij} - m_j)^2}$$

Комбинация центрирования и нормирования по столбцам называется автошкалированием.

$$\tilde{x}_{ij} = (x_{ij} - m_j) / s_j$$

Любое преобразование данных – центрирование, нормирование, и т.п. – всегда делается сначала на обучающем наборе. По этому набору вычисляются значения $m_j$ и $s_j$, которые затем применяются и к обучающему, и к проверочному набору.

В надстройке Chemometrics Add In подготовка данных проводится автоматически. Если подготовку нужно провести вручную, то для нее можно использовать стандартные функции листа или специальную пользовательскую функцию.

В задачах, где структура исходных данных $\mathbf{X}$ априори предполагает однородность и гомоскедастичность, подготовка данных не только не нужна, но и вредна. Именно такой случай представляют ВЭЖХ–ДДМ данные

Размах и отклонение

При заданном числе главных компонент $A$, величина

$$ h_i = \mathbf{t}_i^\mathrm{t}(\mathbf{T}^\mathrm{t}_A\mathbf{T}_A)^{-1}\mathbf{t}_i = \sum_{a=1}^A \frac{t_{ia}^2}{\lambda_a} $$

называется размахом (leverage). Эта величина равна квадрату расстояния Махаланобиса от центра модели до $i$–го образца в пространстве счетов, поэтому размах характеризует как далеко находится каждый образец в гиперплоскости главных компонент.

Для размаха имеет место соотношение

$$h_0 = \frac{1}{I} \sum_{i=1}^I h_i \equiv \frac{A}{I}$$

которое выполняется тождественно – по построению PCA.

Для вычисления размахов можно использовать стандартные функции листа или специальную пользовательскую функцию.

Другой важной характеристикой PCA модели является отклонение $v_i$, которое вычисляется как сумма квадратов остатков – квадрат эвклидова расстояния от плоскости главных компонент до объекта $i$.

Для вычисления отклонений можно использовать стандартные функции листа или специальную пользовательскую функцию.

Две эти величины: $h_i$ и $v_i$ определяют положение объекта (образца) относительно имеющейся PCA модели. Слишком большие значения размаха и/или отклонения свидетельствуют об особенности такого объекта, который может быть экстремальным или выпадающим образцом.

Анализ величин $h_i$ и $v_i$ составляет основу SIMCA – метода классификации с обучением.