PCA——主成分分析

PCA

　　主成分分析（Principal Components Analysis, PCA）是一种降维方法。假设数据集$X\in R^{n\times m}$包含$n$条$m$维数据，PCA即实现线性映射$Y=XD\in R^{n\times k}$。其中矩阵$D\in R^{m\times k},k<m$，其各列向量之间单位正交。PCA目标就是找到一个最优矩阵$D$，使降维映射产生的数据信息丢失越小越好。

　　下面介绍两种PCA的解释及相应的推导，分别从方差与数据压缩的角度理解PCA。

最大化映射方差

　　从方差的角度解释应该是PCA的原始定义。将$D$中的$k$个列向量看做映射向量，则$XD$就是把每个数据映射到各个向量上。而PCA的目标就是使数据在这些映射向量上的映射的方差之和最大。从信息论的角度看，方差越大表明数据在该方向上携带信息越多，因此这个目标是符合信息丢失越小越好的。

　　要计算数据在各方向上映射后的方差，先计算均值如下：

$\displaystyle \mu = \frac{1}{n}ee^TXD$

　　得到与$XD$同尺寸的列方向上的均值$\mu$，从而能进行后续计算。其中$e$表示$n$维全为1的列向量。则我们可以计算数据在$k$个映射向量上的协方差矩阵为

\begin{align}\Sigma &= \frac{1}{n-1}(XD-\mu)^T(XD-\mu) \\&=   \frac{1}{n-1}\left(XD-\frac{1}{n}ee^TXD\right)^T\left(XD-\frac{1}{n}ee^TXD\right)  \\  \end{align}

　　则目标最大化$k$个方向上的方差，就是最大化协方差矩阵的迹（矩阵对角线之和），即

\begin{align}\max\limits_D  \text{Tr}\left(\Sigma \right)& =      \text{Tr}\left[\frac{1}{n-1}\left(XD-\frac{1}{n}ee^TXD\right)^T\left(XD-\frac{1}{n}ee^TXD\right)\right]  \\    \end{align}

　　将两个括号的$D$分别提取到括号左右得

$\displaystyle\max\limits_D   \text{Tr}\left[D^T\frac{1}{n-1}\left( X^T-\frac{1}{n}ee^TX\right)^T\left( X^T-\frac{1}{n}ee^TX\right)D\right]\\$

　　可以看出$D^T$与$D$中间的式子就是计算数据在原始$m$个维度上的协方差矩阵。根据线性代数，下列优化式

$\displaystyle \max\limits_D D^TAD$

在$D$的列向量取对称正定矩阵$A$的前$k$个最大特征值的特征向量时最优。说明，PCA就是计算数据$m$个维度的协方差矩阵的前$k$大特征值对应的$k$个单位正交特征向量，而这些特征值就是数据映射到相应特征向量（主成分）上的方差值。根据特征值的大小排序，我们把这些特征向量称为第一主成分、第二主成分等。

最小化重建损失

　　该方法是花书《深度学习》中的推导方式，表示为降维后重建损失的最小化，这个方法对归一化的必要性没有比较直观的解释。下面是我的早期记录。

　　下面把数据映射为各个方向上的坐标称为编码（即降维），再反映射回来称为解码。

　　定义矩阵$X\in R^{m\times n}$为数据集，由$m$个$n$维行向量数据$\{x_1,x_2,...,x_m\}$组成。首先求得数据均值与方差（以下是按元素进行的运算）：

$\displaystyle\overline{x}=\sum\limits_{i=1}^{m}\frac{x_i}{m}$

$\displaystyle s=\frac{1}{m-1} \sum\limits_{i=1}^{m}\left(x_i-\overline{x}\right)^2 $

　　然后规范化样本：

$\displaystyle x_i^\ast = \frac{x_i-\overline{x}}{\sqrt{s}}, i=1,2,\ldots,m$

　　获得规范化后的数据集$X^*$。随后PCA将把$X^*$编码成矩阵$Y^*\in R^{m\times k},k<n$，通过右乘列向量单位正交的方向矩阵$D\in R^{n\times k}$得到编码后的矩阵$Y^*$：

$Y^* = X^*\cdot D$

　　而再解码回来是再乘上$D^T$：

$\hat{X^*} = Y^*\cdot D^{T}$

　　因为矩阵$D$并不可逆，所以编码是不可逆的，也就是说解码并不能完全恢复回$X^*$。为了最小化精度损失，容易想到，优化$D$来最小化$\hat{X^*}$与$X^*$之差的Frobenius范数：

$ \begin{aligned} &D = \displaystyle\text{arg}\min\limits_D{||\hat{X^*} - X^*||^2_F}\\ &D = \text{arg}\min\limits_D{|| X^*DD^T - X^*||^2_F}\\ &D = \text{arg}\min\limits_D{\text{Tr}[(X^*DD^T - X^*)^T(X^*DD^T - X^*)]}\\ &D = \text{arg}\min\limits_D{\text{Tr}[DD^TX^{*T}X^*DD^T-DD^TX^{*T}X^*-X^{*T}X^*DD^T+X^{*T}X^*]}\\ \end{aligned} $

　　去掉与$D$不相关的$X^{*T}X^*$项，再由迹的性质，矩阵乘法循环右移迹的值不变：

 $ \begin{aligned} &D = \text{arg}\min\limits_D{\text{Tr}[X^*DD^TDD^TX^{*T}-X^*DD^TX^{*T}-D^TX^{*T}X^*D]}\\ &D = \text{arg}\min\limits_D{\text{Tr}[X^*DD^TX^{*T}-X^*DD^TX^{*T}-D^TX^{*T}X^*D]}\\ &D = \text{arg}\max\limits_D{\text{Tr}[D^TX^{*T}X^*D]}\\ \end{aligned} $

　　容易发现当$D$的列向量取$X^{*T}X^*$的前$k$大特征值对应的特征向量时，有这个迹的值最大，等于前$k$大特征值之和，且此时压缩的精度损失最少。

性质

　　经过上面的推导，我们可以发现，$X^{*T}X^*$实际上就是$X$行向量每个元素之间的协方差矩阵$\Sigma_x$（没有除以$m-1$）。定义$\lambda_1,\lambda_2,...,\lambda_n$为$\Sigma_x$从大到小排列的特征值。我们可以推出下面两个性质。

　　性质一、编码得到的矩阵$Y^*\in R^{m\times k}$的行向量各个元素对应的协方差矩阵$\displaystyle\Sigma_y = \text{diag}(\lambda_1,\lambda_2,...,\lambda_k)$ 。

　　性质二、如果不进行压缩，即编码矩阵形如$D\in R^{n\times n}$，从而$Y^*\in R^{m\times n}$，那么$Y^*$的行向量各个元素方差之和等于$X^*$（$X$也一样）的行向量各个元素方差之和，也就是：

$\displaystyle\sum\limits_{i=1}^n \Sigma_{y_{ii}}= \sum\limits_{i=1}^n \Sigma_{x_{ii}} = \sum\limits_{i=1}^n (X^{*T}X^*)_{ii}$

　　首先证明性质一。

　　由之前的定义，$X^*$是$X$经过行标准化后的矩阵，所以对$X^*$的任意列进行求和有

$\displaystyle\sum\limits_{i=1}^mX^*_{ij}=0,\;\;j=1,2,...,n$

　　而对于$Y^*=X^*D$，对$Y^*$的任意列进行求和有

$ \begin{aligned} \sum\limits_{i=1}^m Y^*_{ij} &= \sum\limits_{i=1}^mX^*_{i:}D_{:j} = \sum\limits_{i=1}^m\sum\limits_{k=1}^nX^*_{ik}D_{kj}\\ &= \sum\limits_{k=1}^nD_{kj}\sum\limits_{i=1}^mX^*_{ik}= 0,\;\;j=1,2,...,k \end{aligned} $

　　每列求和都为0，说明，$Y^*$行向量每个元素的均值都为0，无需再计算行向量的均值进行标准化。因此，可以直接计算$Y^*$的协方差矩阵：

$ \begin{aligned} \Sigma_y = Y^{*T}Y^* = D^TX^{*T}X^*D=\text{diag}(\lambda_1,\lambda_2,...,\lambda_k) \end{aligned} $

　　有性质一以后，性质二就容易了：

　　如果不进行压缩，由性质一可得，$Y^*$行向量各元素方差之和就是$X^{*T}X^*$的特征值之和，矩阵特征值之和又等于矩阵的迹，而$X^{*T}X^*$的迹就是$X^*$行向量各个元素的方差之和。所以得到性质二。 

　　另外，映射矩阵$D$的列向量在人脸识别中就是所谓的特征脸（PCA+SVM实现人脸识别）。

代码实现

numpy实现

　　实验代码：

import matplotlib.pyplot as plt
import pylab                #原本在jupyter里才能显示的图片，可以用窗口显示
import numpy  as np 

def PCA_img(img,n):
    t = np.mean(img,0)                                          #以行向量为压缩单元，求行平均
    s = np.std(img,0)                                           #求行标准差
    img=(img-t)/s                                               #每行减去平均值进行规范化
    v,vectors = np.linalg.eig(np.dot(img.T,img)/(len(img)-1))   #生成规范化矩阵后，转置乘自身再除以（行数-1），即为行向量各个元素的协方差矩阵，求协方差矩阵的特征值与特征向量
    indv = v.argsort()                                          #将特征值排序，从小到大
    vectors = vectors[:,indv]                                   #将特征向量按特征值大小排序
    vectors = vectors[:,-n:]                                    #选择前n大特征值对应的特征向量，作为映射矩阵。在人脸识别中就是特征脸
    img = np.dot(np.dot(img,vectors),vectors.T)                 #规范化矩阵先乘映射矩阵再乘映射矩阵的转置，也就是先压缩再解压缩
    img = np.multiply(img,s)+t                                  #把减掉的均值加回去变为原矩阵
    return img
        
img=plt.imread("aaa.jpg") #读取图片  

img.flags["WRITEABLE"] = True
img = img.astype(float) 
n = 50
img[:,:,0] = PCA_img(img[:,:,0],n)
img[:,:,1] = PCA_img(img[:,:,1],n)
img[:,:,2] = PCA_img(img[:,:,2],n)
fig = plt.figure()           
ax = fig.add_subplot(111)    
ax.imshow(img/255)  
print(img)           
pylab.show()              

　　原图是1200×1920的图：

　　压缩为1200×50后再解压回去：

sklearn包PCA

　　sklearn中计算PCA比直接使用numpy矩阵计算快得多，这是因为它里面实现了比较高效的随机算法，虽然每次算出来的值不一样，但是相差不大。代码如下：

import sklearn.decomposition as dc 

pca = dc.PCA(n_components = 10)       #规定PCA要取前k大特征值对应的特征向量的个数，即要映射到几维，PCA默认压缩行，比如100*100压缩到100*10
pca.fit(A)                            #将要用于压缩的矩阵传入，用来“学习”压缩矩阵
D = pca.components_                   #获取“学习”到的编码矩阵（压缩矩阵）D。原矩阵A经过AD^T会得到压缩后的矩阵（这里的D是上面推算的D^T）
C = pca.transform(B)                  #传入新的矩阵B，用“学习”到的编码矩阵进行压缩，即BD^T。在人脸特征提取中，同一个人的不同人脸照片，用同一个压缩矩阵压缩后的坐标向量通常相似，因此可以用于人脸识别的降维