背景與原理：

PCA（主成分分析）是將一個數據的特徵數量減少的同時儘可能保留最多資訊的方法。所謂降維，就是在說對於一個$n$維資料集，其可以看做一個$n$維空間中的點集（或者向量集），而我們要把這個向量集投影到一個$k<n$維空間中，這樣當然會導致資訊損失，但是如果這個$k$維空間的基底選取的足夠好，那麼我們可以在投影過程中儘可能多地保留原資料集的資訊。

資料降維的目的在於使得資料更直觀、更易讀、降低演算法的計算開銷、去除噪聲。

接下來我們討論下如何選取$k$維子空間：

假設原資料集有$m$條資料，每條資料有$n$維，那麼可以將其拼成一個$n*m$的矩陣M，而我們想投影到的$k$維空間的一個單位正交基底為$(p_{1},...,p_{k})$，那麼我們想把這$m$維向量投影到這個空間中實際上就是進行一次矩陣乘法$\begin{pmatrix} p_{1}\\p_{2}\\...\\p_{k} \end{pmatrix} M$

這個道理是簡單易懂的。對於一個向量$\alpha$，其在另一個向量$\beta$方向上的投影是$\dfrac{\alpha \cdot \beta }{|\beta|}$（高中數學）

如果$|\beta|$是一個單位向量，那麼這個投影即為$\alpha \cdot \beta=\beta^{T} \alpha$，於是投影到$k$個單位向量為基底的空間中的情況即如上述所示。

因此我們要找到這$k$個單位向量作為基底，然後拼出$P=\begin{pmatrix} p_{1}\\p_{2}\\...\\p_{k} \end{pmatrix}$即可。

那麼怎麼找呢？我們考慮降維之後我們需要什麼，由於我們在降維之後要儘可能多地保留原始資訊，因此降維之後的資料要提供最大的資訊量，那麼這個資訊量在這裡可以用資料的方差來反映，方差越大，資料的離散程度越高，那麼資料的自身特徵保留的就越好（個人理解：PCA降維的目的在於突出資料的個體特徵，減少資訊損失，而如果降維之後資料離散程度低，意味著這些資料全都堆在一起，那資料的特徵體現的就不明顯了——所有資料全都差不多，這樣個體資訊保留的就不好了）。

對於一個特徵，在$m$組資料中的方差為$\sigma^{2}=\dfrac{1}{m}\sum_{i=1}^{m}(x_{i}-\overline{x})^{2}$，為了便於討論，我們對所有特徵零均值化（即把每個$x_{i}$預先減去$\overline{x}$），這樣一個特徵的方差即為$\sigma^{2}=\dfrac{1}{m}\sum_{i=1}^{m}x_{i}^{2}$

但是降維過程中只考慮方差是不夠的——如果我們發現兩個特徵之間有很強的線性相關性，那麼這兩個特徵其實差別就不大了，我們當然不需要同時保留這兩個特徵，因此我們還希望降維之後任意兩個特徵的協方差（$cov(a,b)=\dfrac{1}{m}\sum_{i=1}^{m}a_{i}b_{i}$，因為已經進行過零均值化了）為零，也就是在說我們選取的子空間的基底一定是正交的。

那麼現在的問題就轉化成了：對於一個$n$維$m$組資料的$n*m$資料矩陣$X$，我們希望將其投影到$n$維空間的一個$k$維子空間中，因此我們要找到$k$個單位正交基$(p_{1},...,p_{k})$，而如果這$k$個單位正交基構成的矩陣$P=\begin{pmatrix} p_{1}\\p_{2}\\...\\p_{k} \end{pmatrix}$，那麼投影過程即為$Y=PX$，$Y$即為降維後所得的$k$維資料集

而結合上述討論，我們希望$Y$各個特徵的方差最大，同時$Y$的兩特徵的協方差為零，這怎麼操作呢？

對於一個$n$維有$m$組資料的$n*m$矩陣$X$，我們考察$C=\dfrac{1}{m}XX^{T}$，那麼我們看到如果$X=\begin{pmatrix} x_{11} & x_{12} & ... & x_{1m}\\...\\x_{n1} & x_{n2} &... & x_{nm}\end{pmatrix}$，我們有：

$XX^{T}=\begin{pmatrix} \sum_{i=1}^{m}x_{1i}^{2} & \sum_{i=1}^{m} x_{1i}x_{2i} &...& \sum_{i=1}^{m}x_{1i}x_{ni}\\...\\ \sum_{i=1}^{m}x_{ni}x_{1i} & \sum_{i=1}^{m} x_{ni}x_{2i} &...& \sum_{i=1}^{m}x_{ni}^{2}\end{pmatrix}$

我們稱$\dfrac{1}{m}XX^{T}$為協方差矩陣，因為我們看到按照我們上面的解釋，這個矩陣是一個實對稱矩陣，其主對角線上的元素是一個特徵維度的方差，而其餘位置上的元素是兩個對應特徵的協方差！

那麼我們的目的是要最大化主對角線上的元素，同時讓其餘位置上的元素為$0$，那麼我們進行的不就是實對稱矩陣的正交相似對角化嘛！

形式化地解釋一下：我們設$Y$的協方差矩陣為$D$，那麼我們希望$D$是一個對角矩陣，同時$D$的主對角線上的元素要儘可能大，那麼我們有：

$D=\dfrac{1}{m}YY^{T}=\dfrac{1}{m}(PX)(X^{T}P^{T})=P(\dfrac{1}{m}XX^{T})P^{T}$

那麼我們實際進行的不就是把$C=\dfrac{1}{m}XX^{T}$這個協方差矩陣正交相似對角化嘛！

至於我們希望主對角線元素儘可能大，那我們就選取$C$的前$k$大的特徵值組成$D$就好了嘛，而此時的$P$就對應於前$k$大的特徵值對應的$k$個正交的特徵向量構成的矩陣。

那麼我們的演算法步驟如下：

對於$n$行$m$列的矩陣$X$，我們解釋成其有$m$組資料，每組資料有$n$個特徵，現在我們欲將其變成$k*m$的矩陣$Y$，表示降維後每組資料只有$k$個特徵。

（1）零均值化：對$X$的每個元素，減去自己所在行的均值（即我們是逐特徵操作，一行對應於同一個特徵，不要搞錯這一點）

（2）計算協方差矩陣$C=\dfrac{1}{\textbf{m}}XX^{T}$

（3）對協方差矩陣對角化$C=P\Sigma P^{T}$，找到其單位正交的特徵向量$e_{1},...,e_{n}$

（4）選取最大的$k$個特徵值對應的特徵向量$e_{i_{1}},...,e_{i_{k}}$，拼成一個變換矩陣$P_{k}=\begin{pmatrix} e_{i_{1}}\\e_{i_{2}}\\...\\e_{i_{k}} \end{pmatrix}$

（5）降維後的資料即為$Y=P_{k}X$

（6）如果希望根據降維後的資料集$Y$近似還原原資料集$\hat{X}$，我們有$\hat{X}=P_{k}^{T}Y$（這裡的邏輯是如果我們不降維，那麼$P_{k}=P$就是一個正交矩陣，那麼$P^{T}P=I$，相當於此時資料集沒有損失，那麼類比這個過程就能匯出近似還原方法$\hat{X}=P_{k}^{T}Y$）

程式碼實現：

import numpy as np
from sympy.matrices import Matrix,GramSchmidt
np.random.seed(1)

x = 7*np.random.rand(100)
y = 0.5*x + 1 + 3*np.random.rand(100)

X = np.hstack([x.reshape(100, 1), x.reshape(100, 1), y.reshape(100, 1), x.reshape(100, 1)])

def centerData(X):
    X = X.copy()
    X -= np.mean(X, axis=0)
    return X

X = centerData(X)print(X[7][2])
C= (np.transpose(X)@X)/100
val,fea=np.linalg.eig(C)
dic=dict()
for i in range(0,4):
    dic[val[i]]=fea[:,i]
val=abs(np.sort(-val))
P=np.vstack([dic[val[0]],dic[val[1]]])
Y=[email protected]
reconstruct_X=Y@P
print(reconstruct_X[7][2])

值得注意的問題是這段程式碼中生成的資料每組資料是一個行向量，每列對應於一個特徵，因此所有的計算和上面的推導都構成一個轉置。

此外這裡使用了numpy裡面的linalg.eig方法用來求一個實對稱矩陣的特徵值和特徵向量，返回的val是特徵值，fea是特徵向量，要特別說明的是不出意外的情況下這裡的val都是按從大到小排序的，而fea實際上是一個矩陣，這個矩陣每個列向量對應於一個特徵值，因此一定要注意選取的方法。上面程式碼使用前兩個特徵向量作為主成分恢復原矩陣，可以看到恢復效果還是不錯的。

import numpy as np
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris

data=load_iris()
X=data.data
Y=data.target

pca=PCA(n_components=2)

X_2d=pca.fit_transform(X)

plt.scatter(X_2d[:,0],X_2d[:,1],c=Y)
plt.show()

當然，PCA也可以直接使用sklearn裡面的包，上述程式碼載入了經典的鳶尾花資料集，然後進行PCA降維（降成二維，這個n_components引數給出了要降到的維度），然後能清楚看到三個鳶尾花的類別，效果很好。