Як працює метод головних компонент (PCA) на простому прикладі



У цій статті я б хотів розповісти про те, як саме працює метод аналізу головних компонент (PCA – principal component analysis) з точки зору інтуїції, що стоїть за її математичним апаратом. Максимально просто, але докладно.

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

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

Наприклад, витрата палива у нас вимірюється в літрах на 100 км, а в США миль на галон. На перший погляд, різні величини, але насправді вони строго залежать один від одного. Милі 1600м, а в галоні 3.8 л. Одна ознака строго залежить від іншого, знаючи один, знаємо і іншого.

Але набагато частіше буває так, що ознаки залежать один від одного не так суворо і (що важливо!) не так явно. Об'єм двигуна в цілому позитивно впливає на розгін до 100 км/год, але це вірно не завжди. А ще може виявитися, що з урахуванням не помітних на перший погляд факторів (типу поліпшення якості палива, використання більш легких матеріалів та інших сучасних досягнень), рік автомобіля не сильно, але теж впливає на його розгін.

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

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

Крок 1. Підготовка даних
Тут для простоти прикладу я не буду брати реальні навчальні датасеты на десятки ознак і сотні спостережень, а зроблю свій, максимально простий іграшковий приклад. 2 ознаки і 10 спостережень буде цілком достатньо для опису того, що, а головне – навіщо, відбувається в надрах алгоритму.

Згенеруємо вибірку:

x = np.arange(1,11)
y = 2 * x + np.random.randn(10)*2
X = np.vstack((x,y))
print X

OUT:
[[ 1. 2. 3. 4. 5. 
6. 7. 8. 9. 10. ]
[ 2.73446908 4.35122722 7.21132988 11.24872601 9.58103444 
12.09865079 13.78706794 13.85301221 15.29003911 18.0998018 ]]

У даній вибірці у нас є дві ознаки, сильно корелюють один з одним. За допомогою алгоритму PCA ми зможемо легко знайти ознака-комбінацію і, ціною частини інформації, висловити обидва ці ознаки одним новим. Отже, давайте розбиратися!

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

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

Xcentered = (X[0] - x.mean(), X[1] - y.mean())
m = (x.mean(), y.mean())
print Xcentered
print "Mean vector: ", m

OUT:
(array([-4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5]),
array([-8.44644233, -8.32845585, -4.93314426, -2.56723136, 1.01013247,
0.58413394, 1.86599939, 7.00558491, 4.21440647, 9.59501658]))
Mean vector: (5.5, 10.314393916)

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

Крок 2. Ковариационная матриця
У випадку з багатовимірної випадкової величиною (випадковим вектором) положення центру все так само буде мат. очікуваннями її проекцій на осі. А ось для опису її форми вже недостатньо толькое її дисперсій по осях. Подивіться на ці графіки, у всіх трьох випадкових величин однакові мат.очікування і дисперсії, а їх проекції на осі в цілому виявляться однакові!



Для опису форми випадкового вектора необхідна ковариационная матриця.

Це матриця, у якої (i,j)-елемент є кореляцією ознак (Xі, Xj). Згадаємо формулу коваріації:



У нашому випадку вона спрощується, так як E(Xі) = E(Xj) = 0:



Зауважимо, що коли Xі = Xj:



і це справедливо для будь-яких випадкових величин.

Таким чином, у нашій матриці по діагоналі будуть дисперсії ознак (т. к. i = j), а в інших комірках – коваріації відповідних пар ознак. А в силу симетричності матриця коваріації теж буде симетрична.

Зауваження: Ковариационная матриця є узагальненням дисперсії на випадок багатовимірних випадкових величин – вона так само описує форму (розкид) випадкової величини, як і дисперсія.
І дійсно, дисперсія одновимірної випадкової величини – це ковариационная матриця розміру 1x1, в якій її єдиний член задано формулою Cov(X,X) = Var(X).

Отже, сформуємо ковариационную матрицю Σ для нашої вибірки. Для цього обчислимо дисперсії Xі Xj, а також їх ковариацию. Можна скористатися вышенаписанной формулою, але раз вже ми озброїлися Python'ом, то гріх не скористатися функцією numpy.cov(X). Вона приймає на вхід список всіх ознак випадкової величини і повертає її ковариационную матрицю і де X – n-вимірний випадковий вектор (n-кількість рядків). Функція відмінно підходить і для розрахунку незміщеної дисперсії і коваріації для двох величин, і для складання ковариационной матриці.
(Нагадаю, що в Python матриця представляється масивом-стовпцем масивів-рядків.

covmat = np.cov(Xcentered)
print covmat, "\n"
print "Variance of X: ", np.cov(Xcentered)[0,0]
print "Variance of Y: ", np.cov(Xcentered)[1,1]
print "Covariance X and Y: ", np.cov(Xcentered)[0,1]

OUT:
[[ 9.16666667 17.93002811]
[ 17.93002811 37.26438587]] 

Variance of X: 9.16666666667
Variance of Y: 37.2643858743
Covariance X and Y: 17.9300281124

Крок 3. Власні вектори й значення (айгенпары)
О'кей, ми отримали матрицю, яка описує форму нашої випадкової величини, з якої ми можемо отримати її розміри по x і y (тобто X1 X2), а також приблизну форму на площині. Тепер треба знайти такий вектор (в нашому випадку тільки один), при якому максимизировался б розмір (дисперсія) проекції нашої вибірки на нього.

Зауваження: Узагальнення дисперсії на вищі розмірності — ковариационная матриця, і ці два поняття еквівалентні. При проекції на вектор максимізується дисперсія проекції, при проекції на простір великих порядків – вся її ковариационная матриця.
Отже, візьмемо одиничний вектор на який будемо проектувати наш випадковий вектор X. Тоді проекція на нього буде дорівнює vTX. Дисперсія проекції на вектор буде відповідно дорівнює Var(vTX). У загальному вигляді у векторній формі (для центрированных величин) дисперсія виражається так:



Відповідно, дисперсія проекції:



Легко помітити, що дисперсія максимізується при максимальному значенні vT Σv. Тут нам допоможе відношення Релея. Не вдаючись занадто глибоко в математику, просто скажу, що у відношення Релея є спеціальний випадок для коваріаційних матриць:


і


Остання формула повинна бути знайома по темі розкладання матриці на власні вектори та значення. x є власним вектором, а λ – власним значенням. Кількість власних векторів і значень дорівнюють розміру матриці (та значення можуть повторюватися).

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

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

Розмірність нашої вибірки дорівнює двом і кількість айгенвекторов у неї, відповідно, 2. Знайдемо їх.

В бібліотеці numpy реалізована функція numpy.linalg.eig(X), де X – квадратна матриця. Вона повертає 2 масиви – масив айгензначений і масив айгенвекторов (вектори-стовпці). І вектори нормовані — їх довжина дорівнює 1. Якраз те, що треба. Ці 2 вектора задають новий базис для вибірки, такий що його осі збігаються з півосями апроксимуючого еліпса нашої вибірки.


На цьому графіку ми апроксимировали нашу вибірку еліпсом з радіусами в 2 сігми (тобто він повинен містити в собі 95% всіх спостережень – що в принципі ми тут і спостерігаємо). Я инвертировал більший вектор (функція eig(X) направляла його у зворотний бік) – нам важливо напрямок, а не орієнтація вектора.

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

Зауваження: діагональні елементи ковариационной матриці показують дисперсії за початкового базису, а її власні значення – за новим (по головних компонентів).
Часто потрібно оцінити обсяг втраченої (та збереженої інформації. Найзручніше представити у відсотках. Ми беремо дисперсії по кожній з осей і ділимо на загальну суму дисперсій по осях (тобто суму всіх власних чисел ковариационной матриці).
Таким чином, наш більший вектор описує 45.994 / 46.431 * 100% = 99.06%, а менший, відповідно, приблизно 0.94%. Відкинувши менший вектор і спроектувавши дані на більший, ми втратимо менше 1% інформації! Відмінний результат!

Зауваження: На практиці, в більшості випадків, якщо сумарна втрата інформації складає не більше 10-20%, то можна спокійно знижувати розмірність.
Для проведення проекції, як вже згадувалося раніше на кроці 3, треба провести операцію vTX (вектор повинен бути довжини 1). Або, якщо у нас не один вектор, а гіперплощина, то замість вектора vT беремо матрицю базисних векторів VT. Отриманий вектор (або матриця) буде масивом проекцій наших спостережень.

v = (-vecs[0][1], -vecs[1][1])
Xnew = dot(v,Xcentered)
print Xnew

OUT:
[ -9.56404107 -9.02021624 -5.52974822 -2.96481262 0.68933859
0.74406645 2.33433492 7.39307974 5.3212742 10.59672425]

dot(X,Y) — почленное твір (так ми перемножуємо вектори та матриці в Python)

Неважко помітити, що значення проекцій відповідають картині на попередньому графіку.

Крок 5. Відновлення даних
З проекцією зручно працювати, будувати на її основі гіпотези і розробляти моделі. Але не завжди отримані головні компоненти будуть мати очевидний, зрозумілий сторонній людині, сенс. Іноді корисно розкодувати, приміром, виявлені викиди, щоб подивитися, що за спостереження за ними стоять.

Це дуже просто. У нас є вся необхідна інформація, а саме координати базисних векторів в заданому базисі (вектори, на які ми проектували) і вектор середніх (для скасування центрування). Візьмемо, приміром, найбільше значення: 10.596… і раскодируем його. Для цього помножимо його праворуч на транспонований вектор і додамо вектор середніх, або в загальному вигляді для всієї выбоки: XTvT+m

Xrestored = dot(Xnew[9],v) + m
print 'Down: ', Xrestored
print 'Original: ', X[:,9]

OUT:
Restored: [ 10.13864361 19.84190935]
Original: [ 10. 19.9094105]

Різниця невелика, але вона є. Адже втрачена інформація не відновлюється. Тим не менш, якщо простота важливіше точності, відновлене значення відмінно апроксимує вихідну.

Замість висновку – перевірка алгоритму
Отже, ми розібрали алгоритм, показали як він працює на іграшковому прикладі, тепер залишилося тільки порівняти його з PCA, реалізованим в sklearn – адже будемо користуватися саме ним.

from sklearn.decomposition import PCA
pca = PCA(n_components = 1)
XPCAreduced = pca.fit_transform(transpose(X))

Параметр n_components вказує на кількість вимірювань, на які буде проводитися проекція, тобто до скількох вимірювань ми хочемо знизити наш датасет. Іншими словами – це n айгенвекторов з самими великими власними числами. Перевіримо результат зниження розмірності:

print 'Our reduced X: \n', Xnew
print 'Sklearn reduced X: \n', XPCAreduced

OUT:
Our reduced X: 
[ -9.56404106 -9.02021625 -5.52974822 -2.96481262 0.68933859
0.74406645 2.33433492 7.39307974 5.3212742 10.59672425]
Sklearn reduced X: 
[[ -9.56404106]
[ -9.02021625]
[ -5.52974822]
[ -2.96481262]
[ 0.68933859]
[ 0.74406645]
[ 2.33433492]
[ 7.39307974]
[ 5.3212742 ]
[ 10.59672425]]

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

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

Перевіримо та інші параметри моделі – функція має ряд атрибутів, що дозволяють отримати доступ до проміжних змінних:

— Вектор середніх: mean_
— Вектор(матриця) проекції: components_
— Дисперсії осей проекції (вибіркова): explained_variance_
— Частка інформації (частка від загальної дисперсії): explained_variance_ratio_

Зауваження: explained_variance_ показує вибіркову дисперсію, тоді як функція cov() для побудови ковариационной матриці розраховує незміщені дисперсії!
Порівняємо отримані нами значення зі значеннями бібліотечної функції.

print 'Mean vector: ', pca.mean_, m
print 'Projection: ', pca.components_, v
print 'Explained variance ratio: ', pca.explained_variance_ratio_, l[1]/sum(l)

OUT: 
Mean vector: [ 5.5 10.31439392] (5.5, 10.314393916)
Projection: [[ 0.43774316 0.89910006]] (0.43774316434772387, 0.89910006232167594)
Explained variance: [ 41.39455058] 45.9939450918
Explained variance ratio: [ 0.99058588] 0.990585881238

Єдина відмінність – в дисперсіях, але як уже згадувалося, ми використовували функцію cov(), яка використовує незміщену дисперсію, тоді як атрибут explained_variance_ повертає вибіркову. Вони відрізняються тільки тим, що перша для отримання мат.очікування ділить на (n-1), а друга – на n. Легко перевірити, що 45.99 ∙ (10 — 1) / 10 = 41.39.

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



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


Отже, ми розглянули принципи роботи алгоритму PCA та його реалізації у sklearn. Я сподіваюся, ця стаття була досить зрозуміла тим, хто тільки починає знайомство з аналізом даних, а також хоч трохи інформації для тих, хто добре знає даний алгоритм. Інтуїтивне уявлення вкрай корисно для розуміння того, як працює метод, а розуміння дуже важливо для правильного налаштування обраної моделі. Спасибі за увагу!

P. S.: Прохання не лаяти автора за можливі неточності. Автор сам у процесі знайомства з дата-аналізом і хоче допомогти таким же, як він у процесі освоєння цієї дивовижної галузі знань! Але конструктивна критика і різноманітний досвід всіляко вітаються!
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.