PCA数学原理推导及python实现

1. PCA介绍

PCA（Principal components analysis），又称主成分分析，是一种统计分析、简化数据集的方法。 它利用正交变换来对一系列可能相关的变量的观测值进行线性变换，从而投影为一系列线性不相关变量的值，这些不相关变量称为主成分（Principal Components）。

具体地，主成分可以看做一个线性方程，其包含一系列线性系数来指示投影方向。PCA对原始数据的正则化或预处理敏感（相对缩放）。

基本思想

（1）将坐标轴中心移到数据的中心，然后旋转坐标轴，使得数据在C1轴上的方差最大，即全部n个数据个体在该方向上的投影最为分散。意味着更多的信息被保留下来。C1成为第一主成分。

（2）C2第二主成分：找一个C2，使得C2与C1的协方差（相关系数）为0，以免与C1信息重叠，并且使数据在该方向的方差尽量最大。

（3）以此类推，找到第三主成分，第四主成分……第p个主成分。p个随机变量可以有p个主成分。

主成分分析经常用于减少数据集的维数，同时保留数据集当中对方差贡献最大的特征。这是通过保留低维主成分，忽略高维主成分做到的。这样低维成分往往能够保留住数据的最重要部分。

2. PCA的数学原理

2.1 PCA的降维思路

定义一个$n\times m$ 的矩阵 $\mathbf{X}$, $n$ 代表样本个数，$m$ 代表变量个数。假设矩阵已经零均值化，比如每列的值都减去该列的均值。

于是，$m\times m$ 维的协方差矩阵 $\mathbf{C}$ 可以得到：

$$\mathbf{C}={\mathbf{X}}^T\mathbf{X}/(n-1)$$

协方差矩阵是一个对称矩阵，可以进行对角化：

$$\mathbf{C}=\mathbf{V}\mathbf{L}{\mathbf{V}}^T$$

其中，矩阵$\mathbf{V}$ 是特征向量矩阵，每一列都是一个特征向量；矩阵$\mathbf{L}$ 是一个对角矩阵，在对角处特征值${\lambda }_i$ 以递减的顺序排列。 特征向量也称为主元的轴或者是主元方向。数据在主元方向的投影称为主成分（principal components, PC scores）。

第$j$个主成分是矩阵$\mathbf{X}\mathbf{V}$的第 $j$ 列。第 $i$ 个数据点在新的主成分空间的坐标是矩阵$\mathbf{X}\mathbf{V}$的第 $i$ 行。

如果将原始数据样本矩阵 $\mathbf{X}$ 进行奇异值分解，可以得到：

$$\mathbf{X}=\mathbf{U}\mathbf{S}{\mathbf{V}}^T$$

其中，$\mathbf{U}$是酉矩阵（酉矩阵和它的共轭转置 = 单位矩阵）；$\mathbf{S}$是对角矩阵，对角元素是奇异值$s_i$, 因此，可以得到：

$$\mathbf{C}=\mathbf{V}\mathbf{S}{\mathbf{U}}^T\mathbf{U}\mathbf{S}{\mathbf{V}}^T/(n-1)=\mathbf{V}\frac{{\mathbf{S}}^2}{n-1}{\mathbf{V}}^T$$

矩阵$\mathbf{V}$的每列（右奇异向量）代表了主元方向；每个奇异值和协方差矩阵的特征值相关，并且有：

$${\lambda }_i=\frac{s_i^2}{n-1}$$

主成分可以通过以下获得：

$$\mathbf{Y}=\mathbf{X}\mathbf{V}=\mathbf{U}\mathbf{S}{\mathbf{V}}^T\mathbf{V}=\mathbf{U}\mathbf{S}$$

其中，$\mathbf{Y}$ 是降维后的样本。

当$n-1>m$时，降维矩阵$\mathbf{Y}$是唯一确定的。

为了实现降维到 $k$维，其中，$k<m$, 从 $\mathbf{U}$矩阵中选择前$k$列，以及矩阵$\mathbf{S}$的左上$k\times k$部分.
降维的数据表示为：

$$\mathbf{Y}={\mathbf{U}}_k{\mathbf{S}}_k$$

其中， $\mathbf{Y}\in R^{n\times k}$.

注意的是，上述对协方差矩阵除了SVD分解，还可以直接求特征值Eigenvalues和特征向量Eigenvectors，

定义载荷矩阵loading ${\mathbf{V}}^T$为：

$$\text{Loadings}=\text{Eigenvectors}\cdot \sqrt{\text{Eigenvalues}}.$$

根据矩阵分解，还可以得到得分矩阵$\mathbf{U}$。

$\mathbf{V}$是loading 矩阵（载荷矩阵），$\mathbf{U}$是scoring矩阵（得分矩阵）。

每个样本乘于载荷矩阵，就可以得到了降维后的新样本，载荷矩阵的含义是每个主成分由各个维度变量构成的权重（PCA是线性变换模型）。

2.2 PCA的推导过程

PCA的目标是找到一组新的正交基 { u 1 , u 2 , ⋯   , u k } \left\{u_{1}, u_{2}, \cdots, u_{k}\right\} {u1​,u2​,⋯,uk​}（从$m$维下降到$k$维），使得数据点在该正交基构成的平面上投影后，数据间的方差最大。

所有数据在该基底上投影的方差为：

$$J_j=\frac{1}{n}\sum _{i=1}^n{\left({\mathbf{X}}_i{u}_j-{\mathbf{X}}_{\text{center}}{u}_j\right)}^2$$

其中，$n$ 为样本数量，如果事先对数据 $\mathbf{X}$ 进行0均值初始化，则${\mathbf{X}}_{\text{center}}$是零向量。

那么，

$$J_j=\frac{1}{n}\sum _{i=1}^n{\left({\mathbf{X}}_i{u}_j\right)}^2=\frac{1}{n}\sum _{i=1}^n\left(u_j^T{\mathbf{X}}_i^T\cdot {\mathbf{X}}_i{u}_j\right)=u_j^T\cdot \frac{1}{n}\sum _{i=1}^n\left({\mathbf{X}}_i^T{\mathbf{X}}_i\right)\cdot u_j$$

所以,

$$J_j=\frac{1}{n}{u}_j^T{\mathbf{X}}^T\mathbf{X}{u}_j$$

记：

$$\mathbf{C}={\mathbf{X}}^T\mathbf{X}/(n-1)$$

根据拉格朗日算子（求极值）求解：

$$\begin{gathered} J_j\approx u_j^T\mathbf{C}{u}_j \\ \text{ s. t .}{u}_j^T{u}_j=1 \end{gathered}$$

求解 $\frac{\partial F}{\partial u_j}=0$, 得：

$$2\mathbf{C}\cdot u_j-2{\lambda }_j\cdot u_j=0\Rightarrow \mathbf{C}{u}_j={\lambda }_j{u}_j$$

当 $u_j,{\lambda }_j$ 分别为 $\mathbf{C}$ 矩阵的特征向量、特征值时，$J_j$ 有极值:

$$J_j=u_j^T{\lambda }_j{u}_j={\lambda }_j$$

所以对于任意满足条件的正交基，对应的数据在上面投影后的方差值为S矩阵的特征向量，从而：

$$J_{\max }=\sum _{j=1}^k{\lambda }_j,$$

将 ${\lambda }_j$ 进行从大到小排序，获得对应的特征向量 $u_j$, 就可以得到期待的正交基 { u 1 , u 2 , ⋯   , u k } \left\{u_{1}, u_{2}, \cdots, u_{k}\right\} {u1​,u2​,⋯,uk​}。

2.3 PCA 求解步骤

步骤：

1.用X表示原有数据；
2.零均值化；
3.求协方差矩阵；
4.求特征值和特征向量；
5.根据相应的特征值把特征向量从大到小排序，从组成的矩阵选取K行代表降维的基（K维）；
6.降维的基和原有数据X相乘，即为降维后的数据Y

PCA主成分选择

假如协方差矩阵$\mathbf{C}$ 有$p$个，即有$p$个主成分，为了合理选择$k$个实现最终数据降维，引入方差贡献率：

$${\mathbf{C}}_i=\frac{{\lambda }_i}{{\sum }_{k=1}^p{\lambda }_k}$$

方差贡献率反映了信息量的大小.

一般通过累计贡献率达到足够大的阈值来确定$k$的大小（一般情况选取0.85）

$$G(k)=\frac{{\sum }_{i=1}^k{\lambda }_i}{{\sum }_{i=1}^p{\lambda }_i}$$

PCA 特点：

优点：

(1) 以方差衡量信息的无监督学习，不受样本标签限制；
(2) 各成分之间正交，消除原始数据成分间可能存在的相关性；
(3) 降维；

缺点：

(1) 主成分解释具有一定的模糊性；
(2) 贡献率小的主成分可能含有一定信息量；
(3) PCA对变量的缩放很敏感。如果我们只有两个变量，而且它们具有相同的样本方差，并且成正相关，那么PCA将涉及两个变量的主成分的旋转。

但是，如果把第一个变量的所有值都乘以100，那么第一主成分就几乎和这个变量一样，另一个变量只提供了很小的贡献，第二主成分也将和第二个原始变量几乎一致。这就意味着当不同的变量代表不同的单位（如温度和质量）时，PCA是一种比较武断的分析方法。

3. python实现PCA

#pca
import numpy as np


dataMat= np.loadtxt(open(r'G:\训练小样本2.csv',"rb"),delimiter=",",skiprows=0)  

#零均值化
def zeroMean(dataMat):        
    meanVal=np.mean(dataMat,axis=0)     #按列求均值，即求各个特征的均值  
    newData=dataMat-meanVal  
    return newData,meanVal  
    

def pca(dataMat,n):  
    #求协方差矩阵
    newData,meanVal=zeroMean(dataMat)  
    covMat=np.cov(newData,rowvar=0)    #求协方差矩阵,return ndarray；若rowvar非0，一列代表一个样本，为0，一行代表一个样本  
    #求特征值、特征矩阵 
    eigVals,eigVects=np.linalg.eig(np.mat(covMat))#求特征值和特征向量,特征向量是按列放的，即一列代表一个特征向量  
    # print(eigVals)
    # print(eigVects)
    #保留主要的成分[即保留值比较大的前n个特征]
    eigValIndice=np.argsort(eigVals)            #对特征值从小到大排序  
    # print(eigValIndice)
    n_eigValIndice=eigValIndice[-1:-(n+1):-1]   #最大的n个特征值的下标  
    # print(n_eigValIndice)
    n_eigVect=eigVects[:,n_eigValIndice]        #最大的n个特征值对应的特征向量  
    print(n_eigVect)
    lowDDataMat=newData*n_eigVect               #低维特征空间的数据  
    reconMat=(lowDDataMat*n_eigVect.T)+meanVal  #重构数据  
    # print(lowDDataMat)
    # print(reconMat)
    return lowDDataMat,reconMat  

# [newData,meanVal]=zeroMean(dataMat)
[a,b]=pca(dataMat,2)
# print(a)
# print(b)
# print(newData)
# print(meanVal)

代码解释：

1.newData=dataMat-meanVal ，返回零均值化结果

2.eigVals,eigVects=np.linalg.eig(np.mat(covMat)),返回矩阵的特征值、特征向量

#特征值
[  4.27623823e+01   4.39882621e+00   1.00919422e+00   1.79586824e-03
   1.29584233e-03]
   
#特征向量
[[ -2.49078536e-01   9.59959063e-01   1.28022644e-01  -6.97723720e-03
   -2.95702270e-05]
 [  1.33074844e-02   3.12351351e-02  -1.56277491e-01   9.53854600e-01
    2.54137809e-01]
 [ -8.09197847e-03  -4.41528939e-03   3.14396428e-02   2.62359608e-01
   -9.64413817e-01]
 [  4.91864942e-02   1.41753127e-01  -9.75180514e-01  -1.45780429e-01
   -7.25104921e-02]
 [ -9.67108061e-01  -2.39561098e-01  -8.49827190e-02   5.31259008e-03
    7.88617062e-03]]

3.eigValIndice=np.argsort(eigVals)特征值从小到大排序

[4 3 2 1 0]

4、n_eigValIndice=eigValIndice[-1:-(n+1):-1]最大n个特征的索引值（下标）

[0 1]#2维

5.n_eigVect=eigVects[:,n_eigValIndice]，最大n个特征值的特征向量（取n列，避免了转置步骤）

[[-0.24907854  0.95995906]
 [ 0.01330748  0.03123514]
 [-0.00809198 -0.00441529]
 [ 0.04918649  0.14175313]
 [-0.96710806 -0.2395611 ]]

6.reconMat=(lowDDataMat*n_eigVect.T)+meanVal #重构数据

3.2 sklearn API实现PCA

#pca降维2
#%matplotlib inline
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.datasets import load_iris

data = load_iris()
y = data.target
X = data.data
pca = PCA(n_components=2)
reduced_X = pca.fit_transform(X)

red_x, red_y = [], []
blue_x, blue_y = [], []
green_x, green_y = [], []

for i in range(len(reduced_X)):
    if y[i] == 0:
        red_x.append(reduced_X[i][0])
        red_y.append(reduced_X[i][1])
    elif y[i] == 1:
        blue_x.append(reduced_X[i][0])
        blue_y.append(reduced_X[i][1])
    else:
        green_x.append(reduced_X[i][0])
        green_y.append(reduced_X[i][1])

plt.scatter(red_x, red_y, c='r', marker='x')
plt.scatter(blue_x, blue_y, c='b', marker='D')
plt.scatter(green_x, green_y, c='g', marker='.')
plt.show()

降维效果：

---

参考：

1. PCA的数学原理；
2. python实现PCA；
3. PCA 理论推导；
4. 详细推导PCA算法;
5. Relationship between SVD and PCA. How to use SVD to perform PCA?;
6. PCA - Loadings and Scores;
7. Loadings vs eigenvectors in PCA: when to use one or another?