机器学习--sklearn之主成分分析(PCA)

PCA原理

假设在$n$维空间内有$m$个点{$x^{(1)},x^{(2)},\dots \dots ,x^{(m)}$ }， 为了降低维度，对于每一个$n$维向量的样本点$x^{(i)}$，我们希望找到一个对应$x^{(i)}$的$l$维编码向量$c^{(i)}$。
假设编码函数为$c=f_{(x)}$，通过这个编码函数将原始的$m$维数据降维至$l$维。
假设解码函数为$g_{(c)}=g_{(f_{(x)})}\approx x$，通过将上一步得到的编码向量解码至$n$维并与原始向量作比较得到误差。
用矩阵形式表示函数，则$g_{(c)}=Dc$，其中$D$是一个$n\times l$的矩阵。
但是我们发现，如果我们按比例缩小所有点对应的编码向量$c$中的第$i$个位置的元素$c_i$，那么只要等比例放大$D_{:,i}$，则解码后的结果不变，如：
假设$n=3,m=4,l=2$，即样本点有4个，它们的维度为3，经过编码后，它们的维度被降低到了2维。
$$C_1=\left[\begin{array}{cccc}1& 2& 3& 4\\ 5& 6& 7& 8\end{array}\right]$$
$$D_1=\left[\begin{array}{cc}1& 2\\ 4& 5\\ 7& 8\end{array}\right]$$
则通过$g_{(c)}=Dc$解码可得
$$g_{(c)1}=\left[\begin{array}{cccc}11& 14& 17& 20\\ 29& 28& 47& 56\\ 47& 62& 77& 92\end{array}\right]$$
而在按0.5的比例缩小所有样本点的初始向量的第一个维度后，
$$C_2=\left[\begin{array}{cccc}0.5& 1& 1.5& 2\\ 5& 6& 7& 8\end{array}\right]$$
此时只要将$D_1$的第一列放大两倍至$D_2$：
$$D_2=\left[\begin{array}{cc}2& 2\\ 8& 5\\ 14& 8\end{array}\right]$$
即可得：
$$g_{(c)2}=\left[\begin{array}{cccc}11& 14& 17& 20\\ 29& 28& 47& 56\\ 47& 62& 77& 92\end{array}\right]=g_{(c)1}$$
因此，每个样本点编码后的向量不唯一。
为了使问题具有唯一解，我们限制$D$中所有列向量都有单位范数。
为了使编码问题更加简单，$PCA$限制$D$的列向量彼此正交。
为了为每个样本点寻找一个最优编码向量$c^{\ast }$，我们使用$L2$范数的平方来衡量原始向量$x$和重构向量$g_{(c^{\ast })}$的距离。
即对于一个样本点来说，有：
$$c^{\ast }=\underset{c}{\mathrm{argmin}}||x-g_{(c)}|{|}_2^2$$
$$=(x-g_{(c)}{)}^T(x-g_{(c)})$$
$$=x^{T}x-x^T{g}_{(c)}-g_{(c)}^{T}x+g_{(c)}^T{g}_{(c)}$$
由于$x^T{g}_{(c)}$和$g_{(c)}^{T}x$都是标量，故二者相等。且由于第一项$x^{T}x$不依赖于$c$，可以忽略，所以：
$$c^{\ast }=\underset{c}{\mathrm{argmin}}-2x^T{g}_{(c)}+g_{(c)}^T{g}_{(c)}$$
用矩阵形式表示$g_{(c)}$后，得到：
$$c^{\ast }=\underset{c}{\mathrm{argmin}}-2x^{T}Dc+c^T{D}^{T}Dc$$
由于$D^{T}D=I_l$，
$$c^{\ast }=\underset{c}{\mathrm{argmin}}-2x^{T}Dc+c^{T}c$$
将$-2x^{T}Dc+c^T{D}^{T}Dc$对$c$求导得：
$${▽}_c(-2x^{T}Dc+c^{T}c)=0$$
$$-2D^{T}x+2c=0$$
$$c=D^{T}x$$
故编码操作被化简为一个简单的矩阵-向量乘法操作。即：
$$f_{(x)}=D^{T}x$$
进一步，PCA解码操作为
$$r_{(x)}=g_{(D^{T}x)}=D{D}^{T}x$$
接下来，我们要挑选编码矩阵$D$。
因为要用相同的矩阵D对所有点进行编码，我们不能再孤立地看待每个点。反之，我们必须最小化所有维数和所有点上的误差矩阵的$Frobenius$范数，即：
$$D^{\ast }=\underset{D}{\mathrm{argmin}}\sqrt{\sum _{i,j}(x_j^{(i)}-r(x^{(i)}{)}_j{)}^2}\text{subject to }D{D}^T=I_l$$
首先，我们考虑$l=1$的情况，即将$n$维原始向量编码为1维向量。此时$D$为一个单一向量$d$。则问题被简化为：
$$d^{\ast }=\underset{d}{\mathrm{argmin}}\sum _i||x^{(i)}-d{d}^T{x}^{(i)}|{|}_2^2\text{subject to }||d|{|}_2=1$$
通常，标量$d^T{x}^{(i)}$应放在向量$d$的左边，且考虑到标量的转置和其自身相等，上式可以写作：
$$d^{\ast }=\underset{d}{\mathrm{argmin}}\sum _i||x^{(i)}-x^{(i)T}dd|{|}_2^2\text{subject to }||d|{|}_2=1$$
此时，将表示各点的$m$个向量上下堆叠成一个$m\times n$的矩阵，记为$X$，其中，$X_{i,:}=x^{(i)T}$。所以这个问题可以被重新表述为：
$$d^{\ast }=\underset{d}{\mathrm{argmin}}||X-Xd{d}^T|{|}_F^2\text{subject to }||d|{|}_2=1$$
此处由$x^{(i)}-x^{(i)T}dd$变为$X-Xd{d}^T$是合理的，例如，$n=3,m=3$，设：
$$x^{(1)}=\left[\begin{array}{c}1\\ 1\\ 1\end{array}\right]x^{(2)}=\left[\begin{array}{c}2\\ 2\\ 2\end{array}\right]x^{(3)}=\left[\begin{array}{c}3\\ 3\\ 3\end{array}\right]$$
$$d=\left[\begin{array}{c}5\\ 5\\ 5\end{array}\right]$$
则：
$$X=\left[\begin{array}{ccc}1& 1& 1\\ 2& 2& 2\\ 3& 3& 3\end{array}\right]$$
$$x^{(1)}-x^{(1)T}dd=\left[\begin{array}{c}1\\ 1\\ 1\end{array}\right]-\left[\begin{array}{ccc}1& 1& 1\end{array}\right]\left[\begin{array}{c}5\\ 5\\ 5\end{array}\right]\left[\begin{array}{c}5\\ 5\\ 5\end{array}\right]=\left[\begin{array}{c}-74\\ -74\\ -74\end{array}\right]$$
同理可得：
$$x^{(2)}-x^{(2)T}dd=\left[\begin{array}{c}-149\\ -149\\ -149\end{array}\right]x^{(3)}-x^{(3)T}dd=\left[\begin{array}{c}-224\\ -224\\ -224\end{array}\right]$$
而
$$X-Xd{d}^T=\left[\begin{array}{ccc}1& 1& 1\\ 2& 2& 2\\ 3& 3& 3\end{array}\right]-\left[\begin{array}{c}15\\ 30\\ 45\end{array}\right]\left[\begin{array}{ccc}5& 5& 5\end{array}\right]=\left[\begin{array}{ccc}-74& -74& -74\\ -149& -149& -149\\ -224& -224& -224\end{array}\right]$$
故可以用$X-Xd{d}^T$代替$x^{(i)}-x^{(i)T}dd$。
暂时不考虑约束，此$Frobenius$范数可做以下简化：
$$\underset{d}{\mathrm{argmin}}||X-Xd{d}^T|{|}_F^2$$
$$=\underset{d}{\mathrm{argmin}}Tr((X-Xd{d}^T{)}^T(X-Xd{d}^T))$$
$$=\underset{d}{\mathrm{argmin}}Tr(X^{T}X)-Tr(X^{T}Xd{d}^T)-Tr(d{d}^T{X}^{T}X)+Tr(d{d}^T{X}^{T}Xd{d}^T)$$
$$=\underset{d}{\mathrm{argmin}}-2Tr(X^{T}Xd{d}^T)+Tr(X^{T}Xd{d}^{T}d{d}^T)$$
这些化简利用的性质有，
1）与$d$无关的项不影响$argmin$。
2）循环改变迹运算中相乘矩阵的顺序不影响结果。
此时，我们再来考虑约束条件
$$\underset{d}{\mathrm{argmin}}-2Tr(X^{T}Xd{d}^T)+Tr(X^{T}Xd{d}^{T}d{d}^T)\text{subject to }{d}^{T}d=1$$
$$=\underset{d}{\mathrm{argmin}}-2Tr(X^{T}Xd{d}^T)+Tr(X^{T}Xd{d}^T)\text{subject to }{d}^{T}d=1$$
$$=\underset{d}{\mathrm{argmin}}-Tr(X^{T}Xd{d}^T)\text{subject to }{d}^{T}d=1$$
$$=\underset{d}{\mathrm{argmax}}Tr(X^{T}Xd{d}^T)\text{subject to }{d}^{T}d=1$$
$$=\underset{d}{\mathrm{argmax}}Tr(d^T{X}^{T}Xd)\text{subject to }{d}^{T}d=1$$
故这个优化问题可以通过特征分解来求解。具体来讲，最优的$d$是$X^{T}X$最大特征值对应的特征向量。
以上推导特定于$l=1$的情况，仅得到了一个主成分。更一般地，当我们需要得到主成分的基时，矩阵$D$由前$l$个最大的特征值对应的特征向量组成。

矩阵分解

在求特征向量的时候，要对$X^{T}X$经行矩阵分解操作，下面列举两种方法：

特征值分解矩阵

如果一个向量v是矩阵A的特征向量，将一定可以表示成下面的形式：

其中，λ是特征向量v对应的特征值，一个矩阵的一组特征向量是一组正交向量。
对于矩阵A，有一组特征向量v，将这组向量进行正交化单位化，就能得到一组正交单位向量。特征值分解，就是将矩阵A分解为如下式：

其中，Q是矩阵A的特征向量组成的矩阵，$\Sigma$则是一个对角阵，对角线上的元素就是特征值。

奇异值分解矩阵（SVD）

奇异值分解是一个能适用于任意矩阵的一种分解的方法，对于任意矩阵A总是存在一个奇异值分解：

假设A是一个$m\ast n$的矩阵，那么得到的U是一个$m\ast m$的方阵，U里面的正交向量被称为左奇异向量。 $\Sigma$是一个$m\ast n$的矩阵， $\Sigma$除了对角线其它元素都为0，对角线上的元素称为奇异值。$V^T$是$V$的转置矩阵，是一个$n\ast n$的矩阵，它里面的正交向量被称为右奇异值向量。而且一般来讲，我们会将 $\Sigma$上的值按从大到小的顺序排列。

SVD分解矩阵A的步骤：
(1) 求$A{A}^T$的特征值和特征向量，用单位化的特征向量构成$U$。
(2) 求$A^{T}A$的特征值和特征向量，用单位化的特征向量构成$V$。
(3) 将$A{A}^T$或者$A^{T}A$的特征值求平方根，然后构成$\Sigma$。

PCA应用（sklearn）

一个简单的小例子

用PCA给鸢尾花数据集中的四个特征进行降维操作。

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

iris = load_iris()
y = iris.target
X = iris.data

#调用PCA
pca = PCA(n_components=2)           #实例化
pca = pca.fit(X)                    #拟合模型
X_dr = pca.transform(X)             #获取新矩阵

#要展示三种分类的分布，需要对三种鸢尾花分别绘图
colors = ['red', 'black', 'orange']
names = iris.target_names
 
plt.figure()
for i in [0, 1, 2]:
    plt.scatter(X_dr[y == i, 0]  #布尔索引切片
                ,X_dr[y == i, 1]
                ,alpha=.7        #指画出的图像的透明度
                ,c=colors[i]
                ,label=names[i]
               )
plt.legend()#图例
plt.title('PCA of IRIS dataset')
plt.show()

#属性explained_variance_，查看降维后每个新特征向量上所带的信息量大小（可解释性方差的大小）
pca.explained_variance_  # 方差从大到小排列，第一个最大，依次减小

#属性explained_variance_ratio，查看降维后每个新特征向量所占的信息量占原始数据总信息量的百分比
#又叫做可解释方差贡献率
pca.explained_variance_ratio_

函数介绍：

pca = PCA(n_components, svd_solver, random_state)：
n_components的可选参数：

• 数字表示要降到几维。
• ‘mle’表示让PCA用大似然估计自选降到几维。
• 输入[0,1]之间的浮点数，并且让参数svd_solver ==‘full’，如PCA(n_components=0.97, svd_solver=“full”)表示希望降维后的总解释性方差占比大于n_components 指定的百分比。
• 当参数n_components中不填写任何值，则默认返回min(X.shape)个特征，一般来说，样本量都会大于特征数目，所以什么都不填就相当于转换了新特征空间，但没有减少特征的个数。一般来说，不会使用这种输入方式。但我们却可以使用这种输入方式来画出累计可解释方差贡献率曲线，以此选择好的n_components的整数取值。

svd_solver的可选参数：

• “auto”：基于X.shape和n_components的默认策略来选择分解器：如果输入数据的尺寸大于500x500且要提取的特征数小于数据最小维度min(X.shape)的80%，就启用效率更高的“randomized”方法。否则,精确完整的SVD将被计算，截断将会在矩阵被分解完成后有选择地发生。
• “full”：从scipy.linalg.svd中调用标准的LAPACK分解器来生成精确完整的SVD，适合数据量比较适中、计算时间充足的情况，生成的精确完整的SVD的结构为:
  
• “arpack”：从scipy.sparse.linalg.svds调用ARPACK分解器来运行截断奇异值分解(SVD truncated)，分解时就将特征数量降到n_components中输入的数值k，可以加快运算速度，适合特征矩阵很大的时候,但一般用于特征矩阵为稀疏矩阵的情况，此过程包含一定的随机性。截断后的SVD分解出的结构为：
  
• “randomized”：通过Halko等人的随机方法进行随机SVD。在"full"方法中，分解器会根据原始数据和输入的n_components值去计算和寻找符合需求的新特征向量，但是在"randomized"方法中，分解器会先生成多个随机向量，然后一一去检测这些随机向量中是否有任何一个符合我们的分解需求，如果符合，就保留这个随机向量，并基于这个随机向量来构建后续的向量空间。这个方法已经被Halko等人证明，比"full"模式下计算快很多，并且还能够保证模型运行效果。适合特征矩阵巨大，计算量庞大的情况。

而参数random_state在参数svd_solver的值为"arpack" or “randomized"的时候生效，可以控制这两种SVD模式中的随机模式。通常我们就选用"auto”，所以一般用不到这个参数。

PCA中的SVD

PCA和SVD涉及了大量的矩阵计算，两者都是运算量很大的模型，但其实，SVD有一种惊人的数学性质，即是它可以不计算协方差矩阵，直接找出一个新特征向量组成的n维空间，而这个n维空间就是奇异值分解后的右矩阵$V^T$。
右奇异矩阵$V^T$有着如下性质：

$k$就是n_components，是我们降维后希望得到的维度。若X为(m,n)的特征矩阵，就是结构为(n,n)的矩阵，取这个矩阵的前k行（进行切片），即将V转换为结构为(k,n)的矩阵。而$V_{(k,n)}^T$与原特征矩阵X相乘，即可得到降维后的特征矩阵X_dr。也就是说，奇异值分解可以不计算协方差矩阵等等结构复杂计算冗长的矩阵，就直接求出新特征空间和降维后的特征矩阵。
在sklearn中，矩阵$U$和$\Sigma$虽然会被计算出来，但完全不会被用到，也无法调取查看或者使用，因此我们可以认为，$U$和$\Sigma$在PCA.fit()过后就被遗弃了。奇异值分解追求的仅仅是$V$，只要有了$V$，就可以计算出降维后的特征矩阵。在transform过程之后，fit中奇异值分解的结果除了$V(k,n)$以外，就会被舍弃，而$V(k,n)$会被保存在属性components_ 当中，可以调用查看：

PCA(2).fit(X).components_

可以用这个函数看最终降维得到的新的特征空间的样子，比如对人脸图像降维后，得到的特征空间应该为五官或轮廓等。

PCA数据降噪

另外，使用接口inverse_transform可以将降维后的数据还原回原始数据（不会完全还原）：

X_inverse = pca.inverse_transform(X_dr)

基于此，我们可以使用PCA算法对数据进行降噪处理。

导入没有噪声的数据集

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

digits = load_digits()
digits.data.shape  # (1797, 64)
set(digits.target.tolist())  # 查看target有哪几个数
digits.images.shape  # (1797, 8, 8)

画出原始数据

def plot_digits(data):
    #data的结构必须是（m,n），并且n要能够被分成（8,8）这样的结构
    fig, axes = plt.subplots(4,10,figsize=(10,4)
                            ,subplot_kw = {"xticks":[],"yticks":[]}
                            )
    for i, ax in enumerate(axes.flat):
        ax.imshow(data[i].reshape(8,8),cmap="binary")
        
plot_digits(digits.data)

得到图像为：

在数据集上加噪声

rng = np.random.RandomState(42)
 
#在指定的数据集中，随机抽取服从正态分布的数据
#两个参数，分别是指定的数据集，和从这个数据集中抽取出来的正太分布的方差
noisy = rng.normal(digits.data,2)
plot_digits(noisy)

得到图像：

所以加了噪声之后，数字变得不是很容易识别了。

PCA降噪

pca = PCA(0.5,svd_solver='full').fit(noisy)
X_dr = pca.transform(noisy)
X_dr.shape  # (1797, 6)

降噪后的图像

without_noise = pca.inverse_transform(X_dr)
plot_digits(without_noise)

得到图像：

总结

PCA参数列表

PCA属性列表

PCA接口列表