主成分分析法（PCA）详解与python实现

PCA原理

主成分分析（principal component analysis,PCA）是一种数据分析方法，出发点是从一组特征中计算出一组按重要性从大到小排列的新特征，它们是原有特征的线性组合，并且相互之间是不相关的。主要用来数据降维、可视化、去噪等。
以样本数据有两个特征为例，如图：

我们将样本点向x轴投影，得到图中红色点列；向y轴投影，得到图中绿色点列，可以看到红色点列比绿色的更好地保留了原数据样本的距离特点，因此若在x、y轴这两个中选一个，我们应选择x轴方向作为降维后的特征。
实际上，我们可以选择的最优方向是图中的黑色虚线方向，当样本数据点向黑色虚线投影后得到的点列比红色点列更优。这里为了进行定量计算，我们引用方差作为判别函数，也就是投影后新点列的方差越大，则该投影方向越优。
将二维推广至n维：设$x_1,x_2,\cdots ,x_n$为数据的原始n个特征，${\xi }_1,{\xi }_2,\cdots ,{\xi }_n$为变换后的n个新特征，而变换方法是对原始特征作线性组合，即：$${\xi }_i=\sum _{j=1}^n{w}_{ij}{x}_j\cdots \cdots （\ast ）$$
而对变换的要求是新特征之间互不相关，（另外为了统一${\xi }_i$的尺度，要求线性组合系数的模为1），选择函数则是新特征的方差，我们会选择使得新特征的方差达到最大值的一组系数作为主成分的变换系数，进而根据方差依次递减求出n个新特征。

利用梯度上升法求解主成分

现在来看主成分的一种求解方法。
将样本集如下表示成mxn矩阵：
$$X=\left(\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1n}\\ x_{21}& x_{22}& \cdots & x_{2n}\\ ⋮& ⋮& \cdots & ⋮\\ x_{m1}& x_{m2}& \cdots & x_{mn}\end{array}\right)$$
每行表示一个样本，共m个样本；每列表示一个特征，共n个特征。
变换向量：$w=\left(\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_n\end{array}\right)$
m个样本求得的主成分：${\xi }_i=x_{i}w$，其中$x_i=(x_{i1}{x}_{i2}\cdots x_{in}),i=1,2,\cdots ,m$
于是有$$Var(\xi )=\frac{{\sum }_{i=1}^m({\xi }_i-\bar{\xi }{)}^2}{m-1}$$

公式推导

为了方便计算，我们将样本每个维度的特征值减去每个维度的均值。即：$X=X-\bar{X}$，容易计算出现在的$\bar{\xi }=0$，从而$$Var(\xi )=\frac{{\sum }_{i=1}^m({\xi }_i-\bar{\xi }{)}^2}{m-1}=\frac{{\sum }_{i=1}^m{\xi }_i^2}{m-1}=\frac{{\sum }_{i=1}^m(x_{i}w{)}^2}{m-1}$$
可求得梯度：$$\begin{gathered} \frac{\partial Var}{\partial w_k}=\frac{2}{m-1}\sum _{i=1}^m(x_{i}w)\ast x_{ik}= \\ \frac{2}{m-1}\ast (x_{1}w{x}_{2}w\cdots x_{m}w)\left(\begin{array}{c}{x}_{1k}\\ x_{2k}\\ ⋮\\ x_{mk}\end{array}\right)= \\ \frac{2}{m-1}\ast (Xw{)}^T\left(\begin{array}{c}{x}_{1k}\\ x_{2k}\\ ⋮\\ x_{mk}\end{array}\right) \end{gathered}$$
因此$$\nabla Var=\left(\begin{array}{c}\frac{\partial Var}{\partial w_1}\\ \frac{\partial Var}{\partial w_2}\\ ⋮\\ \frac{\partial Var}{\partial w_n}\end{array}\right)=\frac{2}{m-1}{X}^T(Xw)$$
求得方差函数的梯度后，我们便可按照梯度上升法求其最大值

python实现

仍以两个特征的样本为例
导入模块：

import numpy as np
import matplotlib.pyplot as plt

创建样本集：

X = np.empty((100,2))
X[:,0] = np.random.uniform(0,10,size=100)
X[:,1] = 3*X[:,0]+2+np.random.normal(0,1,size=100)
m = 100#样本数

初始样本集如图：

定义归零函数并将样本集归零：

def Zero(X):
    return X-X.mean(axis=0)
X = Zero(X)

归零后的样本集如图：

定义归一函数，初始化w：

def One(w):
    return w/np.linalg.norm(w)
w = One(np.ones((2,1)))

定义方差函数：

def Var(X,w):
    return (1/m-1)*np.sum(np.dot(X,w)**2)

下面进入梯度上升法循环：

alpha = 0.1#学习步长
#进入循环：
t = 0#循环次数
while t<1000:
    nabla = (2/m-1)*np.dot(X.T,np.dot(X,w))#梯度
    w0 = w#保存上次循环w值
    w += alpha*nabla#更新w
    w = One(w)#归一
    if np.abs(Var(X,w)-Var(X,w0))<0.001:#如果更新后的w使得方差变化不明显，退出循环
        break
    t += 1#更新循环次数

得出结果如图：

w = [[0.32036449]
[0.94729435]]
可知$\xi =0.32036449x_1+0.94729435x_2$即为降维后的新特征，从图中来看，将其作为主成分是合理的。

计算前k个主成分

之前我们计算了主成分，现在我们来看如何确定得出n个新特征并且选出k个最重要的

公式推导

不同的新特征对应的$w$不同，由此会产生一个特征变换矩阵$W$，有以下关系：
$$\xi =W^{T}x$$其中字母含义同公式$(\ast )$
现在我们要求的是最优的正交变换$W$，使新特征${\xi }_i$的方差达到极值（正交的原因是使得新特征之间两两不相关）
考虑第一个新特征${\xi }_1$，$${\xi }_1=\sum _{j=1}^n{w}_{1j}{x}_j=w_1^{T}x$$
它的方差是$$Var({\xi }_1)=E({\xi }_1^2)-(E({\xi }_1){)}^2=E(w_1^{T}x{x}^T{w}_1)-E(w_1^{T}x)E(x^T{w}_1)=w_1^{T}S{w}_1$$其中$S$为样本特征的协方差矩阵，可以用样本来估计，E是数学期望。又考虑到约束条件$w_1^T{w}_1=1$，因此可以用拉格朗日乘子法确定最大值
$$f(w_1)=w_1^{T}S{w}_1-\lambda (w_1^T{w}_1-1)$$将上述函数对$w_1$求导并令其等于0，解得：
$$S{w}_1=\lambda w_1$$
这说明$w_1$是矩阵$S$的特征向量，而$\lambda$是对应的特征值。将此结果代入之前的方差表达式中
$$Var({\xi }_1)=w_1^{T}S{w}_1=\lambda w_1^T{w}_1=\lambda$$这说明最大的方差就等于最大的特征值，从而第一主成分对应的$w$即为矩阵$S$最大特征值对应的特征向量。
进一步按照不相关以及模为1的要求，可以得出第二主成分对应的向量为矩阵$S$的第二大特征值对应的特征向量，以此类推

python实现

仍以二维特征为例

import numpy as np
import matplotlib.pyplot as plt

#创建样本集：
X = np.random.uniform(0,10,size=100)
Y = 3*X+2+np.random.normal(0,1,size=100)
m = 100#样本数

#定义归零函数并将样本集归零：
def Zero(X):
    return X-X.mean(axis=0)
X = Zero(X)
Y = Zero(Y)

#协方差矩阵：
S = np.cov(X,Y)

#特征值与特征向量：
l,W = np.linalg.eig(S)#l为特征值集合，W为特征向量集合（W的列向量表示每个特征向量）
n = np.argmax(l)#取将l排序后最大值的下标
w = W.T[n]#第一主成分的向量
print(w)
print(l[n]/np.sum(l))
'''
#可视化：
plt.rcParams['font.sans-serif'] = [u'SimHei']
plt.rcParams['axes.unicode_minus'] = False

x = np.arange(-5,5,0.01)
y = w[1]*x/w[0]
plt.plot(x,y,c='r')
plt.text(1,0,'主成分（新特征）',c='r',fontsize=16)
plt.scatter(X,Y,c='lightskyblue')
plt.xlabel('原特征1')
plt.ylabel('原特征2')
plt.title('显示主成分')
plt.grid()
plt.show()
'''

结果：
[-0.3106425 -0.95052682]
0.9986514522294836
这说明该主成分所代表的的数据占全部方差的比例是0.9986514522294836，非常接近1，从而可以用其作为降维后的特征。
如果是多维的特征，可以用下式计算前k个主成分所占比例：
$$\frac{{\sum }_{i=1}^k{\lambda }_i}{{\sum }_{i=1}^n{\lambda }_i}$$进而根据需要选择使用几个主成分来表示数据，如使得比例达到80%或90%等

小结

使用PCA可以实现对特征的变换和降维，这种变化呢是非监督的，也就是没有考虑样本的类别信息。在监督模式识别情况下，这种以方差为目标函数的方法并不一定总有利于后续分类。
另外，除了降维，我们可以把特征值很小的成分视为噪声，将其去掉可以理解为对样本进行了去噪处理。