主成分分析 Principal component analysis(PCA)
===
###### tags: `Dimension reduction`

## 前言：
主成分分析 Principal component analysis(PCA)的核心概念在於以更有價值的基底描述原始數據。

## 核心概念
1. 數據正規化
在PCA中，我們會先將原始數據做正規化處理以利後續的分析，正規化就是將數據減去平均再除以標準差：$$ z = (x - \mu)/\sigma \tag{1} $$
正規化過後的數據會有平均為0、標準差為1的特性。
2. 共變異數
描述兩個變數之間聯合變化的程度，其計算公式為：
$$ Cov(f_1, f_2)=\frac {1} {N-1}\sum_{i=1}^{N} {(f_{i1}-\tilde{f_1})(f_{i2}-\tilde{f_2})} \tag{2}$$
其中，$f_1$、$f_2$ 為兩個變數，$\tilde{f_1}$、$\tilde{f_2}$ 為平均數。

3. 線性轉換
若某一線性轉換會將基底$\begin{bmatrix} 1 \\ 0 \\ \end{bmatrix}$變成$\begin{bmatrix} 2 \\ 1 \\ \end{bmatrix}$，將基底$\begin{bmatrix} 0 \\ 1 \\ \end{bmatrix}$變成$\begin{bmatrix} 3 \\ 2 \\ \end{bmatrix}$，那我們就可以定義一個線性轉換矩陣$\begin{bmatrix} 2 & 3 \\ 1 & 2 \\ \end{bmatrix}$，如此一來透過矩陣乘法我們就可以找出任一向量被轉換後的座標，例如，$\begin{bmatrix} 5 \\ 4 \\ \end{bmatrix}$被轉換後的結果為$\begin{bmatrix} 2 & 3 \\ 1 & 2 \\ \end{bmatrix}$$\begin{bmatrix} 5 \\ 4 \\ \end{bmatrix}=$ $\begin{bmatrix} 22 \\ 13 \\ \end{bmatrix}$

4. 特徵值與特徵向量
特徵向量的意義在於對於給定的線性轉換(任一矩陣都是一線性轉換)，有某個向量轉換後方向不變，只是長度改變（以特徵值的量做倍數伸縮）。
舉例，假設我們有一個二維線性轉換為$\begin{bmatrix} 3 & 1 \\ 0 & 2 \\ \end{bmatrix}$，特徵向量為$\begin{bmatrix} 1 \\ 0 \\ \end{bmatrix}$及$\begin{bmatrix} -1 \\ 1 \\ \end{bmatrix}$，對應特徵值3和2。
你可以發現，當你以$\begin{bmatrix} 3 & 1 \\ 0 & 2 \\ \end{bmatrix}$$\begin{bmatrix} 1 \\ 0 \end{bmatrix}$計算$\begin{bmatrix} 3 & 1 \\ 0 & 2 \\ \end{bmatrix}$這個線性轉換作用在$\begin{bmatrix} 1 \\ 0 \end{bmatrix}$的結果為$\begin{bmatrix} 3 \\ 0 \end{bmatrix}$，
以$\begin{bmatrix} 3 & 1 \\ 0 & 2 \\ \end{bmatrix}$$\begin{bmatrix} -1 \\ 1 \end{bmatrix}$計算$\begin{bmatrix} 3 & 1 \\ 0 & 2 \\ \end{bmatrix}$這個線性轉換作用在$\begin{bmatrix} -1 \\ 1 \end{bmatrix}$的結果為$\begin{bmatrix} -2 \\ 2 \end{bmatrix}$，特徵向量們只是變成原本的"特徵值倍"
## PCA運作流程

在以下文章，本文將帶入實際數據操作展示PCA的流程，以利讀者直觀了解PCA流程，為此，首先假設我們有一個二維（2個特徵）且樣本數為20的數據，你可以在Google Colab以下列程式碼產生這些數據並正規化。

import numpy as np
from numpy.testing import assert_almost_equal

# 設定np只顯示小數點以下兩位，保持畫面整潔
np.set_printoptions(precision=2)
# 給定RandomState一個random seed，同樣的seed可以得到同樣的random結果，因此方便你核對程式碼結果
# 初始化數據 ###
rng = np.random.RandomState(1)
w = rng.rand(2, 2) 
x_normal = rng.normal(scale=5, size=(2, 20)) #由scale=5的高斯分布創(2, 20)的隨機值陣列
x_orig = w @ x_normal
x_mean = x_orig.mean(axis=1)[:, np.newaxis]
x_std = x_orig.std(axis=1)[:, np.newaxis]
x = (x_orig - x_mean) / x_std
mean = x.mean(axis=1)

# # 檢查樣本平均為0
assert_almost_equal(0, mean)

print('x.shape', x.shape, '\n')
print('x',x, '\n')
print('x mean', x.sum(axis=1)/20)
print('x std', x.std(axis=1))

import matplotlib.pyplot as plt
plt.scatter(x[0], x[1])
plt.xlabel('f1')
plt.ylabel('f2')
plt.grid()
plt.show()

如此一來你可以將我們的原始數據畫出來如下圖：<center></center>

1. 計算兩變數之間的共變異數矩陣($K_{f_1 f_2}$)及其特徵向量與特徵值
接著我們計算這筆數據的共變異數矩陣，為什麼要計算共變異矩陣？因為它代表著原始數據的變異以及特徵之間的線性關係。$$K_{f_1 f_2} = \begin{bmatrix} C_{11} & C_{12} \\ C_{21} & C_{22} \\ \end{bmatrix}\tag{3}$$其中 $C_{11}$ 為 $f_1$ 和 $f_1$ 的共變異數，$C_{12}$ 為 $f_1$ 和 $f_2$ 的共變異數，以此類推。
而其實已經正規化數據會簡化共變異數的公式，因為平均數為0，所以你會得到公式：$$Cov(f_1, f_2)=\frac {1} {N-1}\sum_{i=1}^{N} {(f_{i1})(f_{i2})} \tag{4}$$以程式實作計算共變異數矩陣：

# 產生共變異數矩陣

f1_bar = x[0].mean()
f2_bar = x[1].mean()
#證實兩個feature經過數據前處理後平均數為0
assert_almost_equal(0, f1_bar)
assert_almost_equal(0, f2_bar)

# 共變異矩陣為2 * 2矩陣，C11代表f1和f1自己的變異，C12代表f1和f2之間的變異，以此類推
# C11
cov_c11 = 0
for x_ in x.T:
  f1 = x_[0]
  cov_c11 += (f1 - f1_bar) * (f1 - f1_bar)
n = x.shape[1]
cov_c11 /= (n - 1)
print('cov_c11: {:.2f}'.format(cov_c11))

# C12 = C21
cov_c12 = 0
for x_ in x.T:
  f1 = x_[0]
  f2 = x_[1]
  cov_c12 += (f1 - f1_bar) * (f2 - f2_bar)
n = x.shape[1]
cov_c12 /= (n - 1)
print('cov_c12: {:.2f}'.format(cov_c12))
cov_c21 = cov_c12

# C22
cov_c22 = 0
for x_ in x.T:
  f2 = x_[1]
  cov_c22 += (f2 - f2_bar) * (f2 - f2_bar)
n = x.shape[1]
cov_c22 /= (n - 1)
print('cov_c22: {:.2f}'.format(cov_c22))

print('-' * 40)
Cov = np.array([[cov_c11, cov_c12], [cov_c21, cov_c22]])
print("Cov Matrix: \n", Cov)

# 也可以以numpy計算共變異矩陣
Cov_numpy = np.cov(x)
print('Cov_numpy: \n', Cov_numpy)

你可以發現 $K_{f_1 f_2}$ 是一個對稱矩陣，且如上述所說，這個矩陣也是一個線性轉換，若將此線性轉換作用在原始數據上，將可以增強原始數據的變異傾向，使得每個樣本依照兩特徵 $f_1$ 與 $f_2$ 的線性關係做延伸。<br/>
那 $K_{f_1 f_2}$ 的特徵向量又代表什麼意思呢？<br/><br/><center><strong>==在這個增強變異傾向的線性轉換下，方向不變並且延伸的向量==</strong></center><br/>也就是說，這個向量其實就是原始數據的變異方向才會發生沿線性關係增強變異後方向不變還得到延伸的現象！這就是我們所要找的"主成分"！<br/>
於是我們可以著手計算 $K_{f_1 f_2}$ 的特徵向量，你可以依照下列程式得到特徵向量以及特徵值

# 解共變異數矩陣的特徵向量
eig_vals, eig_vecs = np.linalg.eig(Cov)
print(np.linalg.eig(Cov))

# eig_vecs 的每一行（column）為一個 Eigenvectors
print(f"eig_vecs.shape:", eig_vecs.shape)
print(eig_vecs)
print('-' * 40)
# eig_vals 裡則是對應的 Eigenvalues
print("eig_vals.shape:", eig_vals.shape)
print(eig_vals)

請注意！向量的表示法是行向量，接著我們將特徵向量標示在原始數據上

import matplotlib.pyplot as plt
plt.figure()
plt.scatter(x[0], x[1])
plt.scatter(-0.71, 0.71)
plt.scatter(-0.71, -0.71)
plt.arrow(0, 0, -.71, 0.71,length_includes_head=True, head_width=0.1)
plt.arrow(0, 0, -0.71, -0.71,length_includes_head=True, head_width=0.1)
plt.xlabel('f1')
plt.ylabel('f2')
plt.grid()
plt.show()

<center><img src="https://i.imgur.com/vikHJLG.png"></center>

2. 去關聯、降維
有了以上代表 $f_1$、 $f_2$ 變化方向的向量後，我們就可以對原始數據做降維或者去關聯。
- 去關聯
首先，我們要先了解何謂"關聯"，在原始數據中 $f_1$ 和 $f_2$ "有關聯"，因為隨 $f_1$ 變動 $f_2$ 也會跟著變動，那是因為我們正使用最普遍的基底$\begin{bmatrix} 1 \\ 0 \end{bmatrix}$、$\begin{bmatrix} 0 \\ 1 \end{bmatrix}$描述這些數據，而原始數據的變化方向並不平行這兩個基底，所以會有"隨 $f_1$ 變動的 $f_2$ "、"隨 $f_2$ 變動的 $f_1$ "，所以若我們以兩變化方向(共變異矩陣的特徵向量，它們正交)做為基底便可以做到兩個方向上的變動互不相干，即"去關聯"。<br/>為此，我們需要做基底變換，也就是說，以新的基底描述數據使得數據更"乾淨"。<br/>
於是我們拿出剛剛算出來的特徵向量矩陣$\begin{bmatrix} -0.71 & -0.71 \\ 0.71 & -0.71 \\ \end{bmatrix}$，注意！行向量才我們表示向量的方式，這個以$\begin{bmatrix} 1 \\ 0 \end{bmatrix}$、$\begin{bmatrix} 0 \\ 1 \end{bmatrix}$描述的兩個向量形成的矩陣可以做為基底變換矩陣將"原始數據的原始基底座標"轉換成"在$\begin{bmatrix} -0.71 \\ 0.71 \end{bmatrix}$、$\begin{bmatrix} -0.71 \\ -0.71 \end{bmatrix}$為基底的空間座標"，只需要簡單的一次矩陣相乘。<br/>於是我們可以實作並畫出比較圖：

plt.figure()
plt.scatter(x[0], x[1])
eigen_vec_matrix = np.array([[-0.71, -0.71], 
                [0.71, -0.71]])
Change_base = eigen_vec_matrix @ x
plt.scatter(Change_base[0], Change_base[1])
plt.grid()
plt.show()

<center><img src="https://i.imgur.com/3VCEwSz.png"></center>

圖中的箭頭就是$\begin{bmatrix} -0.71 \\ 0.71 \end{bmatrix}$、$\begin{bmatrix} -0.71 \\ -0.71 \end{bmatrix}$在以$\begin{bmatrix} -0.71 \\ 0.71 \end{bmatrix}$、$\begin{bmatrix} -0.71 \\ -0.71 \end{bmatrix}$為基底的世界的樣子，你可以看到以變異方向為基底(橘)和舊的基底(藍)描述數據的差別，原始數據的變異方向在新基底的描述下只是水平或鉛直的方向，並不互相影響。
- 降維
首先我們介紹降維的理想的降維結果應該要有哪些特性：<div><br/>
- 最大變異：降維後所得到的新 K 維特徵 L 具有最大的變異（Variance）
- 最小錯誤：用 K 維的新特徵 L 重新建構回 N 維數據能得到最小的重建錯誤</div>
我們需要讓資料投影後有最大變異，到這裡，應該就知道要選擇哪一個向量了吧！沒錯，就是對應特徵值較大的那個特徵向量，回憶前面，我們提到特徵值是特徵向量伸縮的倍數，而事實上，這兩個條件會同時達成。
現在我們可以計算投影過後的數據變成什麼樣子：

v = np.array([-0.71, -0.71])
print("v   :",v)

# 使用v建立投影矩陣p1，順便轉置v
p1 = v[np.newaxis, :]
print("p1   :",p1)

# 利用p1將數據x投影到v所在之子空間，做點積等於將v.T矩陣相乘x
L = p1 @ x

print("L[:, :]", L[:, :])
plt.figure()
plt.scatter(L[0], [0] * 20) # 沒有y座標，補上0以呈現數線
plt.grid()
plt.show()

<center><img src="https://i.imgur.com/FHoBJTn.png"></center>

那這些一維數據到底解釋了多少原始數據的變異量？於是，我們可以計算變異解釋率（Explained Variance Ratio）$$\frac{\lambda_i}{\Sigma{\lambda_i}}\mid K\vec{v}=\lambda_i\vec{v}$$白話一點就是特徵值所佔的權重，因此由前段程式結果我們可以計算投影在$\begin{bmatrix} -0.71 \\ -0.71 \end{bmatrix}$的變異解釋率，$\frac{1.87}{0.24+1.87} =0.886$，這樣的一維數據其實就已經解釋了原始變異的88%，一般來說85%以上就合格。

## PCA的應用實例
接下來我們以sklearn的PCA套件實作PCA在數據分析上應用的例子，這邊以從Kaggle取得的"Valorant-stats"資料集為例，這是FPS遊戲"Valorant"的槍枝數據，而程式碼撰寫環境為Google Colab。<br/><br/>首先，我們先讀入資料集：

from google.colab import files
import pandas as pd
from IPython.display import display

# 因為資料集為本地端的csv檔案，所以上傳至colab的workspace
valorant_data = files.upload()
# 以csv檔案建立一個dataframe，並以Name欄位做為樣本標籤
df = pd.read_csv("/content/valorant-stats.csv", sep=",", index_col="Name")
# 代換Wall Penetrarion欄位的值，因為要做運算所以把'Low', 'Medium', 'High'代換成1, 2, 3
df['Wall Penetration'] = df['Wall Penetration'].replace(['Low', 'Medium', 'High'], [1, 2, 3])

# 顯示我們建立的dataframe前5列
print("df.shape:", df.shape)
display(df.head(5))

# 顯示各特徵的平均與標準差
df.describe().loc[['mean', 'std']].style.format("{:.2f}")

<center><img src="https://i.imgur.com/U7th3KW.png" width=700></center>

<br/><br/>今天假設我們沒有關於槍枝的領域知識(Domain Knowledge)（狙擊槍射速慢、衝鋒槍射速快...等），我們使用PCA做主成分分析來告訴我們槍枝的特性，但在實際使用PCA之前，別忘了我們要把數據正規化：

from sklearn.preprocessing import StandardScaler
from numpy.testing import assert_almost_equal

# 將Name, Weapon以外的特徵全取出
X = df.iloc[:, 2:]  # (n_samples, n_features)
# 使用 scikit-learn 內建的 API 正規化
scaler = StandardScaler()
Z_sk = scaler.fit_transform(X)  # 注意維度

df.iloc[:, 2:] = Z_sk

# 展示前5列
display(df.head(5).style\
        .format("{:.2f}", subset=df.columns[1:]))

# 顯示各特徵的平均與標準差
df_stats = df.describe().loc[['mean', 'std']]
df_stats.style.format("{:.2f}")

<center><img src="https://i.imgur.com/lPR4xrN.png" width=700></center>

<br/><br/>接著我們實際對處理過後的數據做主成分分析，並降成2維：

from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

# 我們只要最大的兩個主成分。scikit-learn 會自動幫我們
# 依照 eigenvalue 的大小排序共變異數矩陣的 eigenvectors
n_components = 2
random_state = 9527

pca = PCA(n_components=n_components, 
          random_state=random_state)

# 注意我們是對正規化後的特徵 Z_sk 做 PCA
L = pca.fit_transform(Z_sk)  # (n_samples, n_components)

#　以'Weapon Type'的5個類作圖，將投影到第一主成分的數值顯示在 x 軸，第二主成分在 y 軸
label = list(set(df['Weapon Type']))
for i in label:
  plt.scatter(L[df['Weapon Type']==i, 0], L[df['Weapon Type']==i, 1], label=i)

plt.legend()
plt.axis('equal')
plt.grid()
plt.show()

import numpy as np
pcs = np.array(pca.components_) # (n_comp, n_features)

df_pc = pd.DataFrame(pcs, columns=df.columns[2:])
df_pc.index = [f"第{c}主成分" for c in['一', '二']]
# 以梯度變化的程度用'Greys_r'色條為背景上色，色條的兩端為黑與白
# 也就是說黑、白代表梯度變化最高
df_pc.style\
    .background_gradient(cmap='Greys_r', axis=None)\
    .format("{:.2}")

<center><img src="https://i.imgur.com/f5ix9OU.png" width=500></center>
<center><img src="https://i.imgur.com/dtVQY5r.png" width=800></center>

<br/><br/>接下來我們實際對得到的結果做觀察及分析，首先我們可以看到第一主成分的梯度變化最大的是"Fire Rate"欄位且為負值，所以我們可以想見在散佈圖的x軸上，越右邊的點"Fire Rate"越小，第二主成分的梯度變化則在"Magazine Capacity"欄位最大且為正值，代表在y軸上越上面的點有越大的"Magazine Capacity"。

所以，我們可以以主成分分析的結果做出以下歸納：

- Heavy(重型武器)：射速中等偏快、彈匣容量最高
- Sidearm(隨身武器)：射速中等偏快、彈匣容量少
- SMG(衝鋒槍)：射速快、彈匣容量中等
- Sniper(狙擊槍)：射速最慢、彈匣容量少
- Shotgun(霰彈槍)：射速最快、彈匣容量最少
- Rifle(步槍)：射速中等、彈匣容量中等

最後，我們可以試計算降到2維的資料實際有多少變異量：

pca_10d = PCA(10, random_state=random_state)
# fit只算出主成分，不做投影降維
pca_10d.fit(Z_sk)
# round為取近似值的函式，這裡指定到小數點後第二位
np.round(pca_10d.explained_variance_ratio_, 2)

\begin{array}{c|lcr}
主成分 & \text{1} & \text{2} & \text{3} & \text{4} & \text{5} & \text{6} & \text{7} & \text{8} & \text{9} & \text{10} \\
\hline
變異解釋率 & 0.75 & 0.14 & 0.06 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
\end{array}
其實降到2維的資料就已經擁有原始數據89%的變異量了！
<br/>以上就是PCA在資料分析上的實例運用，學會之後不彷自己參考並實作一份感興趣的資料，相信動手操作會讓你收穫更多！
<br/>參考資料及資源：
<a href="https://leemeng.tw/essence-of-principal-component-analysis.html">世上最生動的PCA：直觀理解並應用主成分分析</a>
<a href="https://www.kaggle.com/aadhavvignesh/valorant-weapon-stats#">Kaggle：Valorant 資料集</a>
<a href="https://www.youtube.com/playlist?list=PLZHQObOWTQDPD3MizzM2xVFitgF8hE_ab">3Blue1Brown：線性代數系列影片</a>