主成分分析是在做特征筛选时的重要手段，这个方法在大部分的书中都只是介绍了步骤方法，并没有从头到尾把这个事情给说清楚。本文的目的是把PCA和KPCA给说清楚。主要参考了YouTube上李政轩的Principal Component Analysis and Kernel Principal Component Analysis这个视频（强烈推荐看一下）。
主成分分析所做的工作就是将数据集从高维投影到低维，从而用极少的几个特征来涵盖大部分的数据集信息。
所谓的投影，就是下图所示的这样。
x j x_j xj投影到 v v v上的向量为
x j ′ = ( ∣ ∣ x j ∣ ∣ c o s θ ) v ∣ ∣ v ∣ ∣ x_j'=(||x_j||cos\theta)\dfrac{v}{||v||} xj′=(∣∣xj∣∣cosθ)∣∣v∣∣v
其中，
θ
\theta
θ为
x
j
x_j
xj与
v
v
v的夹角。
由于向量之间的内积为
< x j , v > = ∣ ∣ x j ∣ ∣ ⋅ ∣ ∣ v ∣ ∣ ⋅ c o s θ <x_j, v>=||x_j|| \cdot ||v|| \cdot cos\theta <xj,v>=∣∣xj∣∣⋅∣∣v∣∣⋅cosθ
故有
x j ′ = < x j , v > ∣ ∣ v ∣ ∣ 2 v x_j' = \dfrac{<x_j,v>}{||v||^2}v xj′=∣∣v∣∣2<xj,v>v
如果我们把 v v v设置成单位向量的话（即 ∣ ∣ v ∣ ∣ = 1 ||v||=1 ∣∣v∣∣=1）就有
x j ′ = < x j , v > v x_j' = <x_j,v>v xj′=<xj,v>v
也就是说我们只要求出
<
x
j
,
v
>
<x_j,v>
<xj,v>就可以知道
x
j
x_j
xj投影到
v
v
v上的大小了。
又由于在坐标当中，内积可以表示为
< x j , v > = x j T ⋅ v = v T ⋅ x j <x_j, v>=x_j^T \cdot v=v^T \cdot x_j <xj,v>=xjT⋅v=vT⋅xj
故可以用 v T ⋅ x j v^T \cdot x_j vT⋅xj来表示投影后的数值大小。
主成分分析认为，沿某特征分布的数据的方差越大，则该特征所包含的信息越多，也就是所谓的主成分。
我们已经知道了可以用
v
T
⋅
x
j
v^T \cdot x_j
vT⋅xj来表示投影后的数值大小，那么我们现在就可以算出投影后的方差大小了。注意我们么已经把数据标准化过了，所以
v
T
⋅
x
v^T \cdot x
vT⋅x的均值为
v
T
⋅
0
=
0
v^T \cdot 0=0
vT⋅0=0。
σ 2 = 1 N − 1 ∑ i = 1 N ( v T x i − 0 ) 2 = 1 N − 1 ∑ i = 1 N ( v T x i ) ( v T x i ) \sigma^2 = \dfrac{1}{N - 1}\sum_{i=1}^{N}(v^Tx_i-0)^2=\dfrac{1}{N - 1}\sum_{i=1}^{N }(v^Tx_i)(v^Tx_i) σ2=N−11i=1∑N(vTxi−0)2=N−11i=1∑N(vTxi)(vTxi)
注意到 v T x i v^Tx_i vTxi是一个数值，不是向量，故有 v T x i = ( v T x i ) T v^Tx_i=(v^Tx_i)^T vTxi=(vTxi)T于是
σ 2 = 1 N − 1 ∑ i = 1 N v T x i x i T v = v T ( 1 N − 1 ∑ i = 1 N x i x i T ) v = v T C v \sigma^2=\dfrac{1}{N - 1}\sum_{i=1}^{N}v^Tx_ix_i^Tv=v^T(\dfrac{1}{N- 1}\sum_{i=1}^{N}x_ix_i^T)v=v^TCv σ2=N−11i=1∑NvTxixiTv=vT(N−11i=1∑NxixiT)v=vTCv
其中，
C
=
1
N
−
1
∑
i
=
1
N
x
i
x
i
T
C=\dfrac{1}{N - 1}\sum_{i=1}^{N}x_ix_i^T
C=N−11∑i=1NxixiT是一个
m
×
m
m \times m
m×m的矩阵，
m
m
m为特征的个数。
好了，如果我们要找到最大的方差，也就是要找到一个向量
v
v
v使得方差最大。
我们可以将求最大方差的问题写成
m a x v T C v s . t . ∣ ∣ v ∣ ∣ = 1 max \quad v^TCv \\ s.t. \quad ||v||=1 maxvTCvs.t.∣∣v∣∣=1
又由于 ∣ ∣ v ∣ ∣ = v T v ||v||=v^Tv ∣∣v∣∣=vTv ，故上式即
m a x v T C v s . t . v T v = 1 max \quad v^TCv \\ s.t. \quad v^Tv=1 maxvTCvs.t.vTv=1
利用拉格朗日乘子法可以将上述问题转化为
f ( v , λ ) = v T C v − λ ( v T v − 1 ) f(v,\lambda)=v^TCv-\lambda (v^Tv-1) f(v,λ)=vTCv−λ(vTv−1)
其中， f ( v , λ ) f(v, \lambda) f(v,λ)的平稳点，和我们所要求的最大方差问题是等价的，即求下述方程式的解
{ ∂ f ∂ v = 2 C v − 2 λ v = 0 ∂ f ∂ λ = v T v − 1 = 0 \begin{cases}\dfrac{\partial f}{\partial v}=2Cv-2\lambda v=0 \\ \dfrac{\partial f}{\partial \lambda}=v^Tv-1=0 \end{cases} ⎩⎪⎨⎪⎧∂v∂f=2Cv−2λv=0∂λ∂f=vTv−1=0
上述方程组等价于
{ C v = λ v ∣ ∣ v ∣ ∣ = 1 \begin{cases}Cv=\lambda v \\ ||v|| =1\end{cases} {Cv=λv∣∣v∣∣=1
看到了没， C v = λ v Cv=\lambda v Cv=λv不就是求特征值和特征向量的方程吗！更神奇的地方在下面，我们再回到最初求最大方差的问题
v T C v = v T λ v = λ v T v = λ v^TCv=v^T\lambda v=\lambda v^Tv=\lambda vTCv=vTλv=λvTv=λ
是不是很神奇！要求的方差就是我们这里的特征值！所以我们只需要把 C v = λ v Cv=\lambda v Cv=λv的特征值求出来，然后按大小排个序就，选出最大的几个特征值，并求出对应的特征向量，最后用这几个特征向量来完成数据集在其上的投影 v T x v^Tx vTx，这样就完成了特征的筛选！
值得注意的是， C C C是一个 m × m m \times m m×m的矩阵
C = 1 N − 1 ∑ i = 1 N x i x i T = 1 N − 1 [ x 1 , x 2 , . . . , x N ] [ x 1 T x 2 T . . . x N T ] C=\dfrac{1}{N - 1}\sum_{i=1}^{N}x_ix_i^T=\dfrac{1}{N - 1}[x_1,x_2,...,x_N] \begin{bmatrix}x_1^T \\ x_2^T \\... \\ x_N^T \end{bmatrix} C=N−11i=1∑NxixiT=N−11[x1,x2,...,xN]⎣⎢⎢⎡x1Tx2T...xNT⎦⎥⎥⎤
其中，每个 x i x_i xi为一个列向量
x i = [ x i ( 1 ) x i ( 2 ) . . . x i ( m ) ] x_i=\begin{bmatrix}x_i^{(1)} \\ x_i^{(2)} \\... \\x_i^{(m)} \end{bmatrix} xi=⎣⎢⎢⎢⎡xi(1)xi(2)...xi(m)⎦⎥⎥⎥⎤
其中，
m
m
m为特征的个数。
为了方便表示，我们作出如下定义
X T = [ x 1 , x 2 , . . . , x N ] X^T=[x_1,x_2,...,x_N] XT=[x1,x2,...,xN]
于是， C C C可以表示为
C = 1 N − 1 X T X C=\dfrac{1}{N - 1}X^TX C=N−11XTX
基于核函数的主成分分析和主成分分析的步骤是一样的，只不过用核函数替代了原来的数据。这里对什么是核函数不作说明，请参考其它文章。
对于线性不可分的数据集，我们可以将其映射到高维上，再进行划分。
C = 1 N − 1 ∑ i = 1 N ϕ ( x i ) ϕ ( x i ) T = 1 N [ ϕ ( x 1 ) , . . . , ϕ ( x N ) ] [ ϕ ( x 1 ) T . . . ϕ ( x N ) T ] C=\dfrac{1}{N - 1}\sum_{i=1}^{N}\phi (x_i)\phi(x_i)^T=\dfrac{1}{N}[\phi(x_1),...,\phi(x_N)]\begin{bmatrix}\phi(x_1)^T \\ ... \\ \phi(x_N)^T \end{bmatrix} C=N−11i=1∑Nϕ(xi)ϕ(xi)T=N1[ϕ(x1),...,ϕ(xN)]⎣⎡ϕ(x1)T...ϕ(xN)T⎦⎤
我们令
X T = [ ϕ ( x 1 ) , . . . , ϕ ( x N ) ] X^T=[\phi(x_1),...,\phi(x_N)] XT=[ϕ(x1),...,ϕ(xN)]
那么
C = 1 N − 1 X T X C=\dfrac{1}{N - 1}X^TX C=N−11XTX
在这里， ϕ ( x ) \phi(x) ϕ(x)我们是不知道的，所以上式是没法算的。就算知道了，计算成本也太大了。故引入核函数，我们知道核函数有
K = X X T = [ ϕ ( x 1 ) T . . . ϕ ( x N ) T ] [ ϕ ( x 1 ) , ⋯ , ϕ ( x N ) ] = [ κ ( x 1 , x 1 ) . . . κ ( x 1 , x N ) ⋮ ⋱ ⋮ κ ( x N , x 1 ) ⋯ κ ( x N , x N ) ] K=XX^T=\begin{bmatrix} \phi(x_1)^T \\...\\ \phi(x_N)^T \end{bmatrix} [ \phi(x_1) , \cdots ,\phi(x_N)]=\begin{bmatrix} \kappa(x_1,x_1) & ... & \kappa(x_1,x_N) \\ \vdots & \ddots & \vdots \\ \kappa(x_N, x_1) & \cdots & \kappa(x_N,x_N) \end{bmatrix} K=XXT=⎣⎡ϕ(x1)T...ϕ(xN)T⎦⎤[ϕ(x1),⋯,ϕ(xN)]=⎣⎢⎡κ(x1,x1)⋮κ(xN,x1)...⋱⋯κ(x1,xN)⋮κ(xN,xN)⎦⎥⎤
上述的
K
K
K我们根据核函数的性质是可以算出来的，现在来看看
K
K
K和
C
C
C之间有没有关系。
如果要求
K
K
K的特征值和特征向量的话，我们有下式
( X X T ) u = λ u (XX^T)u=\lambda u (XXT)u=λu
其中，
u
u
u为矩阵
K
K
K的特征向量，
λ
\lambda
λ为矩阵
K
K
K的特征值。
我们对左右两边同时左乘一个
X
T
X^T
XT有
X T ( X X T ) u = λ X T u X^T(XX^T)u=\lambda X^Tu XT(XXT)u=λXTu
即
( X T X ) ( X T u ) = λ ( X T u ) (X^TX)(X^Tu)=\lambda (X^Tu) (XTX)(XTu)=λ(XTu)
又由于
(
N
−
1
)
⋅
C
=
X
T
X
(N - 1) \cdot C=X^TX
(N−1)⋅C=XTX，所以我们发现矩阵
K
K
K和
C
C
C的特征值是相同的，都为
λ
\lambda
λ，
C
C
C的特征向量为
X
T
u
X^Tu
XTu。
由于我们希望特征向量是单位向量，所以我们对其做一下单位化
v = 1 ∣ ∣ X T u ∣ ∣ X T u = 1 u T X X T u X T u = 1 u T K u X T u = 1 u T λ u X T u = 1 λ X T u v=\dfrac{1}{||X^Tu||}X^Tu=\dfrac{1}{\sqrt{u^TXX^Tu}}X^Tu=\dfrac{1}{\sqrt{u^TKu}}X^Tu=\dfrac{1}{\sqrt{u^T\lambda u}}X^Tu=\dfrac{1}{\sqrt{\lambda}}X^Tu v=∣∣XTu∣∣1XTu=uTXXTu 1XTu=uTKu 1XTu=uTλu 1XTu=λ 1XTu
在上式中， λ \lambda λ和 u u u可以通过矩阵 K K K求得，但是 X T X^T XT仍旧是不可知的。那么 C C C的特征向量还是算不出来，难道费了这么大的劲，我们白算了？不急，我们接着往下看。虽然求不出 v v v，但是 v v v并不是我们的最终目标，我们只要知道 x x x在 v v v上的投影就可以了
v T ϕ ( x j ) = ( 1 λ X T u ) T ϕ ( x j ) = 1 λ u T X ϕ ( x j ) = 1 λ u T [ ϕ ( x 1 ) T ⋮ ϕ ( x N ) T ] ϕ ( x j ) = 1 λ u T [ κ ( x 1 , x j ) ⋮ κ ( x N , x j ) ] v^T\phi(x_j)=(\dfrac{1}{\sqrt{\lambda}}X^Tu)^T\phi(x_j)=\dfrac{1}{\sqrt{\lambda}}u^TX\phi(x_j)=\dfrac{1}{\sqrt{\lambda}}u^T\begin{bmatrix} \phi(x_1)^T \\ \vdots \\ \phi(x_N)^T \end{bmatrix} \phi(x_j)=\dfrac{1}{\sqrt{\lambda}}u^T\begin{bmatrix} \kappa(x_1, x_j) \\ \vdots \\ \kappa(x_N, x_j) \end{bmatrix} vTϕ(xj)=(λ 1XTu)Tϕ(xj)=λ 1uTXϕ(xj)=λ 1uT⎣⎢⎡ϕ(x1)T⋮ϕ(xN)T⎦⎥⎤ϕ(xj)=λ 1uT⎣⎢⎡κ(x1,xj)⋮κ(xN,xj)⎦⎥⎤
上式中所有的量都是可以求得的，也就说我们在没有求出特征向量的情况下，直接算出了样本在特征向量上的投影！
这样一来问题就解决了！是不是很神奇！
在python的sklearn包中，已经对PCA和KPCA进行了实现，我们只需要调用函数即可，非常方便。
我们用的数据集是UCI上关于葡萄酒的数据集，得到数据集后对其进行预处理，使得其均值为0。
import pandas as pd
from sklearn.preprocessing import StandardScaler
df = pd.read_csv('http://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data', header=None)
x, y = df.iloc[:, 1:].values, df.iloc[:, 0].values
sc = StandardScaler()
x = sc.fit_transform(x)
这个时候得到的
x
x
x是一个
178
×
13
178 \times 13
178×13规模的数据集，也就是说有
13
13
13个特征，每个特征下有
178
178
178个数据。
我们用主成分分析法将
13
13
13个特征通过线性组合得到一个
2
2
2个特征的数据集。
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
x_pca = pca.fit_transform(x)
然后我们来看下效果
import matplotlib.pyplot as plt
plt.scatter(x_pca[y==1, 0], x_pca[y==1, 1], color='red', marker='^', alpha=0.5)
plt.scatter(x_pca[y==2, 0], x_pca[y==2, 1], color='blue', marker='o', alpha=0.5)
plt.scatter(x_pca[y==3, 0], x_pca[y==3, 1], color='lightgreen', marker='s', alpha=0.5)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.show()
可以得到
很显然，此时已经可以看成是线性可分的数据集了，效果不错。
PCA的使用是有局限性的，如果遇到了，一个像下面这样的线性不可分的数据集，就比较麻烦了。
from sklearn.datasets import make_moons
x2, y2 = make_moons(n_samples=100, random_state=123)
plt.scatter(x2_std[y2==0, 0], x2_std[y2==0, 1], color='red', marker='^', alpha=0.5)
plt.scatter(x2_std[y2==1, 0], x2_std[y2==1, 1], color='blue', marker='o', alpha=0.5)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.show()
不相信的话我们可以用PCA先试下看
x2_std = sc.fit_transform(x2)
x_spca = pca.fit_transform(x2_std)
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(14,6))
ax[0].scatter(x_spca[y2==0, 0], x_spca[y2==0, 1], color='red', marker='^', alpha=0.5)
ax[0].scatter(x_spca[y2==1, 0], x_spca[y2==1, 1], color='blue', marker='o', alpha=0.5)
ax[1].scatter(x_spca[y2==0, 0], np.zeros((50,1))+0.02, color='red', marker='^', alpha=0.5)
ax[1].scatter(x_spca[y2==1, 0], np.zeros((50,1))+0.02, color='blue', marker='o', alpha=0.5)
ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')
plt.show()
从图中可以看出，经过主成分分析之后，数据仍旧是线性不可分的。接下来，我们用基于核函数的主成分分析来试下看。
from sklearn.decomposition import KernelPCA
kpca = KernelPCA(n_components=2, kernel='rbf', gamma=15)
x_kpca = kpca.fit_transform(x2_std)
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(14,6))
ax[0].scatter(x_kpca[y2==0, 0], x_kpca[y2==0, 1], color='red', marker='^', alpha=0.5)
ax[0].scatter(x_kpca[y2==1, 0], x_kpca[y2==1, 1], color='blue', marker='o', alpha=0.5)
ax[1].scatter(x_kpca[y2==0, 0], np.zeros((50,1))+0.02, color='red', marker='^', alpha=0.5)
ax[1].scatter(x_kpca[y2==1, 0], np.zeros((50,1))+0.02, color='blue', marker='o', alpha=0.5)
ax[0].set_xlabel('PC1')
ax[0].set_ylabel('PC2')
ax[1].set_ylim([-1, 1])
ax[1].set_yticks([])
ax[1].set_xlabel('PC1')
plt.show()
由图可知，只需要把数据集投影到经变换后的特征 P C 1 PC1 PC1上就可以实现线性划分了，这个时候只需要一个特征 P C 1 PC1 PC1就够了。
[1] https://www.youtube.com/watch?v=G2NRnh7W4NQ&t=1s
[2] Raschka S. Python Machine Learning[M]. Packt Publishing, 2015.