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

Ме́тод головни́х компоне́нтів (МГК, англ. principal component analysis, PCA) — метод факторного аналізу в статистиці, який використовує ортогональне перетворення множини спостережень з можливо пов'язаними змінними (сутностями, кожна з яких набуває різних числових значень) у множину змінних без лінійної кореляції, які називаються головними компонентами.

Метод головних компонент — один з основних способів зменшити розмірність даних, втративши найменшу кількість інформації. Винайдений Карлом Пірсоном у та доповнений і розширений Гарольдом Готелінґом в 1933 р. Застосовується в багатьох галузях, зокрема, в економетриці, біоінформатиці, обробці зображень, для стиснення даних, у суспільних науках.

Обчислення головних компонент може бути зведене до обчислення сингулярного розкладу матриці даних або до обчислення власних векторів і власних чисел коваріаційної матриці початкових даних. Іноді метод головних компонент називають перетворенням Кархунена — Лоева або перетворенням Хотеллінга (англ. Hotelling transform).

Загальна характеристика

Схема математичних перетворень. На малюнку позначення: X — матриця вихідних даних розмірністю n×m (n — число об'єктів спостереження, m — число елементарних аналітичних ознак); Z — матриця центрованих і нормованих значень ознак; R — матриця парних кореляцій. Λ — діагональна матриця власних (характеристичних) чисел. A — матриця факторного відображення, її елементи a_rj — вагові коефіцієнти. Спочатку A має розмірність m*m — за кількістю елементарних ознак Xj, потім в аналізі залишається r найвагоміших компонент, r ≤ m. Обчислюють матрицю A за відомими даними матриці власних чисел Λ і нормованими власними векторами V за формулою A = VΛ^1/2. F — матриця значень головних компонент розмірністю r×n.

Метод головних компонент — один із найпоширеніших методів факторного аналізу.

Серед інших подібних методів, що дозволяють узагальнювати значення елементарних ознак, МГК виділяється простою логічною конструкцією, і, разом з тим, на його прикладі стають зрозумілими загальна ідея й цілі чисельних методів факторного аналізу.

Метод головних компонент дає можливість за m — числом початкових ознак виділити r головних компонент, або узагальнених ознак. Простір головних компонент ортогональний. Математична модель методу головних компонент базується на логічному припущенні, що значення множини взаємозалежних ознак породжують деякий загальний результат.

Розв'язування задачі методом головних компонент зводиться до поетапного перетворення матриці початкових даних X.

Формальна постановка задачі

Задача аналізу головних компонент має щонайменше чотири базових версії:

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

Перші три версії оперують скінченними множинами даних. Вони еквівалентні і не використовують жодної гіпотези про статистичне породження даних. Четверта версія оперує випадковими величинами. Скінченні множини з'являються тут як вибірки з даного розподілу, а розв'язання трьох перших завдань — як наближення до розкладу за (теоремою Кархунена — Лоева) («істинного перетворення Кархунена — Лоева»). При цьому виникає додаткове і не цілком тривіальне питання про точність цього наближення.

Апроксимація даних лінійними многовидами

Ілюстрація до знаменитої роботи Пірсона (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}+...+\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].}$

У лівій частині цієї тотожності стоїть середньоквадратична відстань між точками, а в квадратних дужках праворуч — вибіркова дисперсія. Таким чином, у методі головних компонент шукаються підпростори, в проєкції на які середньоквадратична відстань між точками максимальна (або, що те ж саме, її спотворення внаслідок проєкції мінімальне). Таке переформулювання дозволяє будувати узагальнення зі зважуванням різних парних відстаней (а не тільки точок).

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

Для заданої ${\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 ...\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}$.

Хоча формально завдання сингулярного розкладу матриці даних і спектрального розкладу коваріаційної матриці збігаються, алгоритми обчислення сингулярного розкладу безпосередньо, без обчислення коваріаційної матриці і її спектра, більш ефективні і стійкі.

Теорія сингулярного розкладання була створена Джеймсом Джозефом Сильвестром у і викладена в усіх докладних посібниках з теорії матриць.

Простий ітераційний алгоритм сингулярного розкладання

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

${\displaystyle F(b,a)={\frac {1}{2}}\sum _{i=1}^{m}\sum _{j=1}^{n}(x_{ij}-b_{i}a_{j})^{2}\to \min }$

Розв'язок цієї задачі дається послідовними ітераціями за явними формулами. При фіксованому векторі ${\displaystyle a=(a_{j})}$ значення ${\displaystyle b=(b_{i})}$, що надають мінімум формі ${\displaystyle F(b,a)}$однозначно і явно визначаються з рівностей ${\displaystyle \partial F/\partial b_{i}=0}$ :

${\displaystyle {\frac {\partial F}{\partial b_{i}}}=-\sum _{j=1}^{n}(x_{ij}-b_{i}a_{j})a_{j}=0;\;\;b_{i}={\frac {\sum _{j=1}^{n}x_{ij}a_{j}}{\sum _{j=1}^{n}a_{j}^{2}}}\,.}$

Аналогічно, при фіксованому векторі ${\displaystyle b=(b_{i})}$ визначаються значення ${\displaystyle a=(a_{j})}$:

${\displaystyle a_{j}={\frac {\sum _{i=1}^{m}b_{i}x_{ij}}{\sum _{i=1}^{m}b_{i}^{2}}}\,.}$

Як початкове наближення вектора ${\displaystyle a}$ береться випадковий вектор одиничної довжини, обчислюємо вектор ${\displaystyle b}$, далі для цього вектора ${\displaystyle b}$ обчислюємо вектор ${\displaystyle a}$ і т. д. Кожен крок зменшує значення ${\displaystyle F(b,a)}$. Як критерій зупинки використовується малість відносного зменшення значення функціоналу ${\displaystyle F(b,a)}$, що мінімізується, за крок ітерації (${\displaystyle \Delta F/F}$) або малість самого значення ${\displaystyle F}$.

Внаслідок цього для матриці ${\displaystyle X=(x_{ij})}$ отримується найкраще наближення матрицею ${\displaystyle P_{1}}$ виду ${\displaystyle b^{1}\otimes a^{1}=(b_{i}^{1}a_{j}^{1})}$ (тут верхнім індексом позначено номер наближення). Далі, з матриці ${\displaystyle X}$ віднімається отримана матриця ${\displaystyle P_{1}}$, і для отриманої матриці ухилень ${\displaystyle X_{1}=X-P_{1}}$ знову шукається найкраще наближення ${\displaystyle P_{2}}$ цього ж виду і т. д., поки, наприклад, норма ${\displaystyle X_{k}}$ не стане достатньо малою. В результаті отримали ітераційну процедуру розкладання матриці ${\displaystyle X}$ у вигляді суми матриць рангу 1, тобто ${\displaystyle X=P_{1}+P_{2}+...+P_{q}\;(P_{l}=b^{l}\otimes a^{l})}$ . Приймаємо ${\displaystyle \sigma _{l}=\|a^{l}\|\|b^{l}\|}$ і нормуємо вектори ${\displaystyle a^{l}\,,\,b^{l}}$: ${\displaystyle a^{l}:=a^{l}/\|a^{l}\|;\,\,b^{l}:=b^{l}/\|b^{l}\|.}$ В результаті отримано апроксимацію сингулярних чисел ${\displaystyle \sigma _{l}}$ і сингулярних векторів (правих — ${\displaystyle a^{l}}$ і лівих — ${\displaystyle b^{l}}$).

До переваг цього алгоритму відноситься його виняткова простота і можливість майже без змін перенести його на дані з прогалинами, а також зважені дані.

Існують різні модифікації базового алгоритму, що поліпшують точність і стійкість. Наприклад, вектори головних компонент ${\displaystyle a^{l}}$ при різних ${\displaystyle l}$ повинні бути ортогональні «за побудовою», однак за великого числа ітерацій (велика розмірність, багато компонент) малі відхилення від ортогональності накопичуються і може знадобитися спеціальна корекція ${\displaystyle a^{l}}$ на кожному кроці, що забезпечує його ортогональність раніше знайденим головним компонентам.

Для квадратних симетричних додатно визначених матриць описаний алгоритм перетворюється на метод прямих ітерацій для пошуку власних векторів (див. статтю Власні вектори та власні числа).

Сингулярне розкладання тензорів і тензорний метод головних компонент

Часто вектор даних має додаткову структуру прямокутної таблиці (наприклад, плоске зображення) або навіть багатовимірної таблиці — тобто тензора: ${\displaystyle x_{i_{1}i_{2}...i_{q}}}$, ${\displaystyle 1\leq i_{j}\leq n_{j}}$. У цьому випадку також ефективно застосовувати сингулярне розкладання. Визначення, основні формули та алгоритми переносяться практично без змін: замість матриці даних маємо ${\displaystyle q+1}$-індексну величину ${\displaystyle \operatorname {X} =(x_{i_{0}i_{1}i_{2}...i_{q}})}$, де перший індекс ${\displaystyle i_{0}}$-номер точки (тензора) даних.

Основна процедура — пошук найкращого наближення тензора ${\displaystyle x_{i_{0}i_{1}i_{2}...i_{q}}}$ тензором виду ${\displaystyle a_{i_{0}}^{0}a_{i_{1}}^{1}a_{i_{2}}^{2}...a_{i_{q}}^{q}}$ (де ${\displaystyle a^{0}=(a_{i_{0}}^{0})}$ — ${\displaystyle m}$-вимірний вектор (${\displaystyle m}$ — кількість точок даних), ${\displaystyle a^{l}=(a_{i_{l}}^{l})}$ — вектор розмірності ${\displaystyle n_{l}}$ при ${\displaystyle l>0}$) методом найменших квадратів:

${\displaystyle F={\frac {1}{2}}\sum _{i_{0}=1}^{m}\sum _{i_{1}=1}^{n_{1}}...\sum _{i_{q}=1}^{n_{q}}(x_{i_{0}i_{1}...i_{q}}-a_{i_{0}}^{0}a_{i_{1}}^{1}...a_{i_{q}}^{q})^{2}\to \min }$

Розв'язок цієї задачі отримується послідовними ітераціями за явними формулами. Якщо задано всі вектори-співмножники крім одного ${\displaystyle a_{i_{k}}^{k}}$, то він визначається явно з достатніх умов мінімуму.

${\displaystyle a_{i_{k}}^{k}={\frac {\sum _{i_{0}=1}^{m}\sum _{i_{1}=1}^{n_{1}}...\sum _{i_{k-1}=1}^{n_{k-1}}\sum _{i_{k+1}=1}^{n_{k+1}}...\sum _{i_{q}=1}^{n_{q}}x_{i_{0}i_{1}...i_{k-1}i_{k}i_{k+1}...i_{q}}a_{i_{0}}^{0}a_{i_{k-1}}^{k-1}a_{i_{k+1}}^{k+1}...a_{i_{q}}^{q}}{\prod _{j\neq k}\|a^{j}\|^{2}}}\,.}$

Як початкове наближення векторів ${\displaystyle a^{l}=(a_{i_{l}}^{l})}$ (${\displaystyle l>0}$) беруться випадкові вектори одиничної довжини, обчислимо вектор ${\displaystyle a^{0}}$далі для цього вектора ${\displaystyle a^{0}}$ і даних векторів ${\displaystyle a^{2},a^{3},...}$ обчислюється вектор ${\displaystyle a^{1}}$ і так далі (циклічно перебираючи індекси). Кожен крок зменшує значення ${\displaystyle F(b,a)}$. Алгоритм, очевидно, збігається. Як критерій зупинки використовується малість відносного зменшення значення функціоналу ${\displaystyle F}$, що мінімізується, за цикл або малість самого значення ${\displaystyle F}$. Далі, від тензора ${\displaystyle \operatorname {X} }$ віднімається отримане наближення ${\displaystyle a_{i_{0}}^{0}a_{i_{1}}^{1}a_{i_{2}}^{2}...a_{i_{q}}^{q}}$ і для залишку знову шукається найкраще наближення цього ж вигляду і т. д., поки, наприклад, норма чергового залишку не стане достатньо малою.

Це багатокомпонентне сингулярне розкладання (тензорний метод головних компонент) успішно застосовується при обробці зображень, відеосигналів, і, ширше, будь-яких даних, що мають табличну або тензорну структуру.

Матриця перетворення до головних компонентів

Матриця ${\displaystyle A}$ перетворення даних до головних компонент складається з векторів головних компонент, розташованих у порядку спадання власних значень:

${\displaystyle A=\left\{a_{1},...,a_{n}\right\}^{T}}$ (${\displaystyle {\,}^{T}}$ означає транспонування), причому

${\displaystyle AA^{T}=1.}$

Тобто, матриця ${\displaystyle A}$ є ортогональною.

Велика частина варіації даних буде зосереджена в перших координатах, що дозволяє перейти до простору меншої розмірності.

Залишкова дисперсія

Нехай дані центровані, ${\displaystyle {\overline {X}}=0}$. При заміні векторів даних ${\displaystyle x_{i}}$ їх проєкцією на перші ${\displaystyle k}$ головних компонент ${\displaystyle x_{i}\mapsto \sum _{j=1}^{k}a_{j}(a_{j},x_{i})}$ вноситься середній квадрат помилки в розрахунку на один вектор даних:

${\displaystyle {\frac {1}{m}}\sum _{i=1}^{m}\left\Vert x_{i}-\sum _{j=1}^{k}a_{j}(a_{j},x_{i})\right\Vert ^{2}=\sum _{l=k+1}^{n}\lambda _{l},}$

де ${\displaystyle \lambda _{1}\geq \lambda _{2}\geq ...\geq \lambda _{n}\geq 0}$ — власні значення емпіричної коваріаційної матриці ${\displaystyle C}$, розташовані в порядку спадання, з урахуванням кратності.

Ця величина називається залишковою дисперсією. Величина

${\displaystyle {\frac {1}{m}}\sum _{i=1}^{m}\left\Vert \sum _{j=1}^{k}a_{j}(a_{j},x_{i})\right\Vert ^{2}={\frac {1}{m}}\sum _{i=1}^{m}\sum _{j=1}^{k}(a_{j},x_{i})^{2}=\sum _{l=1}^{k}\lambda _{l}}$

називається поясненою дисперсією. Їх сума дорівнює вибірковій дисперсії. Відповідний квадрат відносної помилки — це відношення залишкової дисперсії до вибіркової дисперсії (тобто частка непоясненої дисперсії):

${\displaystyle \delta _{k}^{2}={\frac {\lambda _{k+1}+\lambda _{k+2}+...+\lambda _{n}}{\lambda _{1}+\lambda _{2}+...+\lambda _{n}}}.}$

За відносною помилкою ${\displaystyle \delta _{k}}$ оцінюється придатність методу головних компонент з проєціюванням на перші ${\displaystyle k}$ компонент.

Зауваження: в більшості обчислювальних алгоритмів власні числа ${\displaystyle \lambda _{i}}$ з відповідними власними векторами головними компонентами ${\displaystyle a_{i}}$ обчислюються в порядку «від великих ${\displaystyle \lambda _{i}}$ — до менших». Для обчислення ${\displaystyle \delta _{k}}$ достатньо обчислити перші ${\displaystyle k}$ власних чисел і слід емпіричної коваріаційної матриці ${\displaystyle C}$, ${\displaystyle \operatorname {tr} C}$ (суму діагональних елементів ${\displaystyle C}$, тобто дисперсій за осями). Тоді

${\displaystyle \delta _{k}^{2}={\frac {1}{\operatorname {tr} C}}\left(\operatorname {tr} C-\sum _{i=1}^{k}\lambda _{i}\right).}$

Відбір головних компонент за правилом Кайзера

Цільовий підхід до оцінки числа головних компонент з необхідною часткою поясненої дисперсії формально застосовується завжди, однак неявно він припускає, що немає поділу на «сигнал» і «шум», і будь-яка наперед задана точність має сенс. Тому часто більш продуктивна інша евристика, яка ґрунтується на гіпотезі про наявність «сигналу» (порівняно мала розмірність, відносно велика амплітуда) і «шуму» (велика розмірність, відносно мала амплітуда). З цієї точки зору метод головних компонент працює як фільтр: сигнал міститься, в основному, в проєкції на перші головні компоненти, а в інших компонентах пропорція шуму значно вища.

Питання: як оцінити число необхідних головних компонент, якщо відношення «сигнал/шум» заздалегідь не відоме?

Найпростіший і найстаріший метод відбору головних компонент дає правило Кайзера (англ. Kaiser's rule): значущі ті головні компоненти, для яких

${\displaystyle \lambda _{i}>{\frac {1}{n}}\operatorname {tr} C,}$

тобто ${\displaystyle \lambda _{i}}$ перевищує середнє значення ${\displaystyle \lambda }$ (середню вибіркову дисперсію координат вектора даних). Правило Кайзера добре працює в простих випадках, коли є декілька головних компонент з ${\displaystyle \lambda _{i}}$, які значно перевищують середнє значення, а інші власні числа, менші за нього. У більш складних випадках воно може давати занадто багато значущих головних компонент. Якщо дані нормовані на одиничну вибіркову дисперсію за осями, то правило Кайзера набуває особливо простого вигляду: значущі тільки ті головні компоненти, для яких ${\displaystyle \lambda _{i}>1}$.

Оцінка числа головних компонент за правилом зламаної тростини

Приклад: оцінка числа головних компонент за правилом зламаної тростини в розмірності 5.

Одним з найбільш популярних евристичних підходів до оцінки числа необхідних головних компонент є правило зламаної тростини (англ. Broken stick model). Набір нормованих на одиничну суму власних чисел (${\displaystyle \lambda _{i}/\operatorname {tr} C}$, ${\displaystyle i=1,...n}$) порівнюється з розподілом довжин уламків тростини одиничної довжини, зламаної в ${\displaystyle n-1}$-й випадково вибраній точці (точки зламу вибираються незалежно і рівномірно розподілені вздовж тростини). Нехай ${\displaystyle L_{i}}$ (${\displaystyle i=1,...n}$) — довжини отриманих шматків тростини, занумеровані в порядку зменшення довжини${\displaystyle L_{1}\geq L_{2}\geq ...L_{n}}$. Нескладно знайти математичне сподівання ${\displaystyle L_{i}}$:

${\displaystyle l_{i}=\operatorname {E} (L_{i})={\frac {1}{n}}\sum _{j=i}^{n}{\frac {1}{j}}.}$

За правилом зламаної тростини ${\displaystyle k}$-й власний вектор (у порядку зменшення власних чисел ${\displaystyle \lambda _{i}}$) зберігається в списку головних компонент, якщо

${\displaystyle {\frac {\lambda _{1}}{\operatorname {tr} C}}>l_{1}\;and\;{\frac {\lambda _{2}}{\operatorname {tr} C}}>l_{2}\;and\;...{\frac {\lambda _{k}}{\operatorname {tr} C}}>l_{k}.}$

На рисунку наведено приклад для 5-вимірного випадку:

${\displaystyle l_{1}}$=(1+1/2+1/3+1/4+1/5)/5; ${\displaystyle l_{2}}$=(1/2+1/3+1/4+1/5)/5; ${\displaystyle l_{3}}$=(1/3+1/4+1/5)/5; ${\displaystyle l_{4}}$=(1/4+1/5)/5; ${\displaystyle l_{5}}$=(1/5)/5.

Для прикладу вибрано

${\displaystyle {\frac {\lambda _{1}}{\operatorname {tr} C}}}$=0.5; ${\displaystyle {\frac {\lambda _{2}}{\operatorname {tr} C}}}$=0.3; ${\displaystyle {\frac {\lambda _{3}}{\operatorname {tr} C}}}$=0.1; ${\displaystyle {\frac {\lambda _{4}}{\operatorname {tr} C}}}$=0.06; ${\displaystyle {\frac {\lambda _{5}}{\operatorname {tr} C}}}$=0.04.

За правилом зламаної тростини в цьому прикладі слід залишати 2 головні компоненти:

${\displaystyle {\frac {\lambda _{1}}{\operatorname {tr} C}}>l_{1}\;;\;{\frac {\lambda _{2}}{\operatorname {tr} C}}>l_{2}\;;\;{\frac {\lambda _{3}}{\operatorname {tr} C}}<l_{3}\;.}$

За оцінками користувачів, правило зламаної тростини має тенденцію занижувати кількість значущих головних компонент.

Нормування

Нормування після зведення до головних компонент

Після проєціювання на перші ${\displaystyle k}$ головних компонент з ${\displaystyle \lambda _{1}\geq \lambda _{2}\geq ...\geq \lambda _{k}>0}$ зручно провести нормування на одиничну (вибіркову) дисперсію за осями. Дисперсія вздовж ${\displaystyle i}$-ї головної компоненти дорівнює ${\displaystyle \lambda _{i}>0\;(1\leq i\leq k}$), тому для нормування треба поділити відповідну координату на ${\displaystyle {\sqrt {\lambda _{i}}}}$. Це перетворення не є ортогональним і не зберігає скалярного добутку. Коваріаційна матриця проєкції даних після нормування стає одиничною, проєкції на будь-які два ортогональні напрямки стають незалежними величинами, а будь-який ортонормований базис стає базисом головних компонент (нагадаємо, що нормування змінює відношення ортогональності векторів). Відображення простору початкових даних на перші ${\displaystyle k}$ головних компонент разом з нормуванням задається матрицею

${\displaystyle K=\left\{{\frac {a_{1}}{\sqrt {\lambda _{1}}}},{\frac {a_{2}}{\sqrt {\lambda _{2}}}},...,{\frac {a_{k}}{\sqrt {\lambda _{k}}}}\right\}^{T}}$.

Саме це перетворення найчастіше називається перетворенням Кархунена — Лоева. Тут ${\displaystyle a_{i}}$ — вектори-стовпці, а верхній індекс ${\displaystyle T}$ означає транспонування.

Нормування до обчислення головних компонент

Попередження: не слід плутати нормування, що проводиться після перетворення до головних компонент, з нормуванням і «знерозмірюванням» при передобробці даних, що проводиться до обчислення головних компонент. Попереднє нормування потрібне для обґрунтованого вибору метрики, в якій буде обчислюватися найкраща апроксимація даних, або будуть шукатися напрямки найбільшого розкиду (що еквівалентно). Наприклад, якщо дані являють собою просторові вектори з «метрів, літрів і кілограмів», то при використанні стандартної евклідової відстані різниця на 1 метр за першою координатою робитиме такий самий внесок, що й різниця на 1 літр за другою, або на 1 кг за третьою. Зазвичай системи одиниць, в яких подані початкові дані, недостатньо точно відображають наші уявлення про природні масштаби за осями, і проводиться «знерозмірювання»: кожна координата ділиться на певний масштаб, який визначається даними, цілями їх обробки і процесами вимірювання і збору даних.

Є три істотно різних стандартних підходи до такого нормування: на одиничну дисперсію за осями (масштаби за осями дорівнюють середнім квадратичним ухиленням — після цього перетворення коваріаційна матриця збігається з матрицею коефіцієнтів кореляції), на рівну точність вимірювання (масштаб за віссю пропорційний точності вимірювання даної величини) і на рівні вимоги в задачі (масштаб за віссю визначається необхідною точністю прогнозу даної величини або допустимим її спотворенням — рівнем толерантності). На вибір передобробки впливають змістовна постановка задачі, а також умови збору даних (наприклад, якщо колекція даних принципово не завершена і дані будуть ще надходити, то нераціонально вибирати нормування строго на одиничну дисперсію, навіть якщо це відповідає змісту завдання, оскільки це передбачає ренормалізацію всіх даних після отримання нової порції; розумніше вибрати певний масштаб, що грубо оцінює стандартне відхилення, і далі його не змінювати).

Попереднє нормування на одиничну дисперсію за осями руйнується поворотом системи координат, якщо осі не є головними компонентами, і нормування при передобробці даних не замінює нормування після зведення до головних компонент.

Механічна аналогія і метод головних компонент для зважених даних

Якщо зіставити кожному вектору даних одиничну масу, то емпірична коваріаційна матриця ${\displaystyle C}$ збіжиться з тензором інерції цієї системи точкових мас (поділеним на повну масу ${\displaystyle m}$), а задача про головні компоненти — із задачею зведення тензора інерції до головних осей. Можна використовувати додаткову свободу у виборі значень мас для врахування важливості точок даних або надійності їхніх значень (важливим даними або даними з більш надійних джерел приписуються великі маси). Якщо вектору даних ${\displaystyle x_{l}}$ надається маса ${\displaystyle w_{l}}$, то замість емпіричної коваріаційної матриці ${\displaystyle C}$ отримаємо

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

Всі подальші операції зі зведення до головних компонент виконуються так само, як і в основній версії методу: шукається ортонормований власний базис ${\displaystyle C^{w}}$, впорядковується за спаданням власних значень, оцінюється середньозважена помилка апроксимації даних першими ${\displaystyle k}$ компонентами (за сумами власних чисел ${\displaystyle C^{w}}$), проводиться нормування і так далі.

Більш загальний спосіб зважування дає максимізація зваженої суми попарних відстаней між проєкціями. Для кожних двох точок даних, ${\displaystyle x_{l},\ x_{q}}$ вводиться вага ${\displaystyle d_{lq}}$; ${\displaystyle d_{lq}=d_{ql}}$ і ${\displaystyle d_{l}=\sum _{q=1}^{m}d_{lq}}$. Замість емпіричної коваріаційної матриці ${\displaystyle C}$ використовується

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

При ${\displaystyle d_{lq}>0}$ симетрична матриця ${\displaystyle C^{d}}$ додатно визначена, оскільки додатна квадратична форма:

${\displaystyle \sum _{ij}c_{ij}^{d}a_{i}a_{j}={\frac {1}{2}}\sum _{lq}d_{lq}\left(\sum _{i}a_{i}(x_{li}-x_{qi})\right)^{2}.}$

Далі шукаємо ортонормований власний базис ${\displaystyle C^{d}}$, впорядковуємо його за спаданням власних значень, оцінюємо середню помилку апроксимації даних першими ${\displaystyle k}$ компонентами і т. д. — так само, як і в основному алгоритмі.

Цей спосіб застосовується за наявності класів: для ${\displaystyle x_{l},\ x_{q}}$ з різних класів вага ${\displaystyle d_{lq}}$ вага вибирається більшою, ніж для точок одного класу. Як наслідок, у проєкції на зважені головні компоненти різні класи «розсуваються» на більшу відстань.

Інше застосування — зниження впливу великих ухилень, так званих викидів (en.:outlier), які можуть спотворювати картину через використання середньоквадратичної відстані: якщо вибрати ${\displaystyle d_{lq}=1/\|x_{l}-x_{q}\|}$, то вплив великих ухилень буде зменшено. Таким чином, описана модифікація методу головних компонент є більш робастною, ніж класична.

Спеціальна термінологія

В статистиці при використанні методу головних компонент використовують кілька спеціальних термінів.

• Матриця даних — ${\displaystyle \mathbf {X} =\{x_{1},...x_{m}\}^{T}}$; кожен рядок — вектор передоброблених даних (центрованих і правильно нормованих), число рядків — ${\displaystyle m}$ (кількість векторів даних), число стовпців — ${\displaystyle n}$ (розмірність простору даних);
• Матриця навантажень (англ. loadings) — ${\displaystyle \mathbf {P} =\{a_{1},...a_{k}\}}$; кожен стовпець — вектор головних компонент, число рядків — ${\displaystyle n}$ (розмірність простору даних), число стовпців — ${\displaystyle k}$ (кількість векторів головних компонент, вибраних для проєціювання);
• Матриця рахунків (англ. scores) — ${\displaystyle \mathbf {T} =[t_{ij}];\;t_{ij}=(x_{i},a_{j})}$; кожен рядок — проєкція вектора даних на ${\displaystyle k}$ головних компонент; число рядків — ${\displaystyle m}$ (кількість векторів даних), число стовпців — ${\displaystyle k}$ (кількість векторів головних компонент, вибраних для проєціювання);
• Матриця ${\displaystyle Z}$-рахунків (англ. ${\displaystyle Z}$-scores) — ${\displaystyle \mathbf {Z} =[z_{ij}];\;z_{ij}={\frac {(x_{i},a_{j})}{\sqrt {\lambda _{j}}}}}$; кожен рядок — проєкція вектора даних на ${\displaystyle k}$ головних компонент, нормована на одиничну вибіркову дисперсію; число рядків — ${\displaystyle m}$ (кількість векторів даних), число стовпців — ${\displaystyle k}$ (кількість векторів головних компонент, вибраних для проєктування);
• Матриця помилок (або залишків) (англ. errors або residuals) — ${\displaystyle \mathbf {E} =\mathbf {X} -\mathbf {T} \mathbf {P} ^{T}}$.
• Основна формула: ${\displaystyle \mathbf {X} =\mathbf {T} \mathbf {P} ^{T}+\mathbf {E} }$.

Межі застосування та обмеження ефективності методу

Метод головних компонент застосовний завжди. Розповсюджене твердження про те, що він застосовний тільки до нормально розподілених даних (або для розподілів близьких до нормальних) хибне: у початковому формулюванні Пірсона ставиться задача про наближення скінченної множини даних та відсутня навіть гіпотеза про їх статистичне породження, не кажучи вже про розподіл.

Однак метод не завжди ефективно знижує розмірність за заданих обмежень на точність ${\displaystyle \delta _{k}}$. Прямі і площини не завжди забезпечують гарну апроксимацію. Наприклад, дані можуть з хорошою точністю дотримуватися якоїсь кривої, а ця крива може бути складно розташована у просторі даних.