主成分分析（Principal Component Analysis，PCA）

维度灾

在学习PCA之前我们需要明白为什么要使用它？或者说为什么需要进行降维？首先我们需要了解一下什么是维度诅咒？

维度诅咒是指当数据集的维度增加时，数据将变得极其稀疏，使得对于聚类、离群点分析等关心的关于两点的距离、密度等信息变得失去意义，而且子空间的组合也会呈指数增长。

如果我们在一个单位平面（1*1的正方形）中随机选择一个点，它距离边界的距离小于0.001的概率只有大约0.4%；但是在一个10000维的单位超立方体中，这个概率大于99.999999%，在高维超立方体中大多数的点都非常接近边界。

这就说明了当数据集的维度很高时，数据会变得非常的稀疏，我们无法去想象这么一个空间，也无法使用现有的方法在这么高的维度上处理问题。而且维度越高，发生过拟合的风险就越大。而在实际中我们使用的数据集往往具有很多的特征，属于一个高维的数据集，而且总存在着一些问题，我们无法方便的进行可视化和使用它。

总的来说进行降维的动机概括如下：
• 消除噪声和冗余的特征
• 加速后续的训练算法
• 方便数据的可视化，便于从中了解到最重要的特征
• 用于数据压缩，节省存储空间

但是进行降维后，由于我们舍弃了一部分特征，必然会有其他问题的出现，主要的弊端如下所示：
• 丢失部分信息，可能使后面训练算法的性能降低，但这是无法避免的
• 可能是计算密集型的
• 为机器学习流水线增添了复杂度
• 转换后的特征可解释性下降

但权衡来说利大于弊所以需要进行降维，而其中最常用的是投影和流形学习，今天看其中一种方法就是主成分分析（PCA）。

数学支撑

这部分内容可见 PCA的推导和解释 ，文中讲了向量的内积、向量的投影、基、协方差、SVD等相关的数学基础知识，告诉我们如何从公式推导理解PCA的算法过程。

主成分分析（Principal Component Analysis，PCA）

PCA属于投影的一种，那么什么是投影？在很多问题中，我们得到的训练数据往往是高维数据集，同时数据不是在每一个维度均匀分布的，很多的特征之间有着一定的关联。因此可以说高维的数据空间近似于受一个低得多的低维子空间所影响，所以我们可以将高维的数据投影到一个和它相关联的低维的子空间中，尽可能的保留数据之间的信息，又方便我们使用。

比如在如下的三维空间中有很多数据，他们都紧靠着一个超平面，如果将其投影到这个超平面上，数据之间的某些关联并不会损失太多，比如两点之间的距离。

根据这样的想法将其进行投影，得到如下所示的二维的新的数据集，从中可以看出效果还是很好的。

但是将超平面稍微转动一下，我们就可以得到无限多的超平面，使得三维的数据集可以投影到二维平面上，如何选择最优的那个超平面，成为了一个新的问题。仔细回想一下之前我们是不是也遇到过类似的问题呢？没错就是在支持向量机中，由于感知机得到的分类线（以二维情况为例）是随机选择的一条，所以在支持向量机中根据间隔最大化的原则，在保证分类正确的前提下，尽可能使得数据点离分类边界远一些，也就是使间隔大一些，这样得到的模型效果好，泛化能力强。

在PCA中我们同样是寻找那个最优的超平面，但要使得数据点距离超平面足够的近，满足最近重构性，；又要使得在超平面上的投影尽可能分开，满足最大可分性。根据在SVM中的思想，我们需要找到一个标准来衡量我们超平面的选择，根据PCA对超平面的要求，我们可以从两个角度思考这个问题

1. 最近重构性
   在线性回归中我们希望在二维平面的数据集中，找到一条直线可以很好的拟合数据，同时又要防止过拟合。那里通常使用的衡量标准是均方误差，可以使用最小二乘法求得闭式解，也可以通过梯度下降法求得结果。它衡量是预测值和真实值之间的差距，在PCA中没有这个需求，但同样可以借鉴这种思想，我们可以使用点到直线的距离来判断是否满足最近重构性。

   
   

2. 最大可分性
   要使投影后的图满足最大可分性，即使得投影到直线上的点距离就要尽量的远，也可以说数据在某一个轴上的投影方差最大。

   如下图所示，在二维平面上我们有蓝色的数据点分布，将其标准化后我们的数据分布如下所示。假设我们选择了三条直线去做投影，结果如右下图所示，显然方向向量为C1的投影效果最好，因为它保留了最多的点的信息，同时投影方差最大；C2方向的次之，如果我们主成分只选择一个的话，就选择C1。

   
   

下面我们看一下PCA是如何求解的，以一个小的实验为例。假设我们研究理想弹簧的物理运动，小球的质量为m，弹簧的质量很小，摩擦力很小。将弹簧拉伸离平衡位置一段距离后放开，则弹簧会沿着平衡点位置来回的运动，我们希望可以找到它运动的规律，用一个数学表达式进行表达。

如果我们对于这个物理现象有足够的认知，就会知道可以将其表达为关于时间的函数，这时就只含有一个变量X。但是在正常情况下，我们无法知道哪几个维度{x,y,z}组成的三维空间是最优的，所以任取三个方向{a,b,c}放置三个摄像机来观察弹簧的运动，三个摄像机之间的放置角度也许并没有90°。

现在我们的问题 Q1是：如何通过摄像机获取的数据集来得到关于x的方程式？ 在真实的情况中，由于空气阻力、摩擦力以及摄像机的问题，我们得到的数据往往是含有噪声的，我们需要从被污染的数据中找到想要的规律。而PCA的目的是计算最有意义的基来重新表达一个有噪声的数据集。希望这个新的基能够滤除噪声，揭示隐藏的结构。所以在这个例子中，我们希望可以通过PCA证实：沿着X轴的方向是最重要的维度。

我们将每个时间点的样本作为数据集中一个独立的数据，每一个数据都是多维的（记录了动量、位置等）。其中某一时刻，摄像机记录了球的位置$（xA,yA）$，则一个样本就可以使用六维列向量表示$$\mathbf{X}=\left(\begin{array}{c}{x}_A\\ y_A\\ x_B\\ y_B\\ x_c\\ y_C\end{array}\right)$$

这样每一个X都是一个m维的向量，其中m表示测量的类型数，它们都位于一个由某些标准正交基所组成的多维空间中，所有的X都可以表示成这些基的线性组合，所以现在问题就转换成了Q2：如何找到这些标准正交基？

假如我们只看摄像机A，通常使用（1，0）和（0，1）来作为标准正交基，而不是其他的基，因为它反映了我们数据采集的方法。则位置（2，2）就可以表示成表示向上2个单位，向左2个单位。

在线性代数中，（1，0）和（0，1）可以被表示为单独的行向量，将其组合可以得到2*2的矩阵Ⅰ。将其推广到m维，同样可以得到m * m的矩阵Ⅰ如下所示，将其记为B

其中$bi$是含有m个分量的正交基向量，我们数据集中的每一个数据都可以表示为它们的线性组合。但是这是最普通的基，为了更好的表达我们的数据，我们的问题转换为Q3：能否通过PCA找到一组比原始基更好的基来表达X呢？

设X是一个$m\ast n$的矩阵（n = 6,m = 72000）,它记录了原始的数据，Y是另一个$m\ast n$的矩阵，由线性变换P得到$PX=Y$其中$pi$表示P的行，$xi$表示X的列，$yi$表示Y的列。

其中关于P我们可以认为它是X到Y的线性转换矩阵、在几何意义上，P表示将X变成Y的旋转或拉伸的表示或认为P的每一行都是一组表达X列的新的基向量。因此可以得到如下的表达式

Y可以表示成列向量的形式

其中$y_i$的每一个系数都可以表示成$x_i$与p中对应行的点积，也可以说$y_i$是$x_i$在P上的投影，因此P的行组成了一组表示X的列的新的基向量。根据前面的分析知道，P的行向量就是X的主成分，找到P就解决了问题，那么问题就又转换成了Q4：如何找到这个P，而且它是表达X最优的选择?

为了得到最优的P，就必须知道我们希望Y表现出什么特征，前面假设的条件只有线性一个，为了得到上面问题的答案，还需多加一些假设。

影响表达能力的三个因素

在此之前我们看一下在线性假设条件下，影响表达能力的三个因素：噪声、旋转和冗余。

1. 噪声和旋转

为了提取到数据中有用的东西，数据集中的噪声自然越小越好，但是噪声的衡量不是绝对的，无法说多大尺度就是噪声。通常采用信噪比（SNR）来衡量数据，因为在信号处理中认为信号具有较大的方差，噪声具有较小的方差，信噪比越大越好。因此SNR可以由如下的等式表示$$SNR=\frac{{\sigma }_{sigmal}^2}{{\sigma }_{noise}^2}$$

我们将从摄像机A中得到的数据绘制成图，因为理想情况下弹簧的运动是直线的，任何偏离直线的运动产生的数据就必然是噪声。而信号和噪声的方差分别由图中指示的两个方向上的数据产生，它们的长度之比就是信噪比的大小。当SNR>>1时，数据表示成一条很瘦的直线，当SNR =1 时，表现为一个圆，其他更糟的情况应该也可以想象出来。

将方差和SNR关于直线角度的关系绘制成图，如下所示，从中可以看出方差最大的方向不是X轴的方向，也不是Y轴的方向，而是如上图长的直线所示的方向，因此也证实了原始基并不是表达数据集最好的选择。我们感兴趣的特征的方向是最大方差和最大信噪比的方向。

将角度是p时方差最大，我们需要将原始基旋转到这个角度，即可以最好拟合数据的直线的平行方向上，这时就可以得到关于弹簧问题的运动情况。但是如何找到这个最优的旋转角度p，以及如何将其推广到任意维的空间中？

2. 冗余
在弹簧的例子中，我们使用了三个摄像机来获取它运动的数据，但是是否有必要这样进行记录？事实上它们会记录大量重复的信息，故是没必要的。如下图所示，假如有如下三种数据的分布，它反映了r1和r2这两种任意测量类型之间的一系列可能的图。

（a）r1和r2完全不相关，它们是静态独立的。
（b）r1和r2存在一定的相关性
（c）中r1和r2是高度相关的，这时只需记录一个变量 r 就可以很好的表达数据，因为另一个变量的数据可以通过它的最佳拟合直线和 r 推导出。

PCA的求解

讲了这么多，我们如何来求解PCA呢？主要有两种方法。

1. 使用特征向量的协方差求解
   假设$X$是$m\ast n$的矩阵，m表示测量类型数，n表示样本数，我们的目标就是找到最优的P，使得$Y=PX$,得到它的协方差矩阵$$C_Y=\frac{1}{n-1}Y{Y}^T$$

   则P的行就是X的主成分。首先将Cy写成关于P的表达式$$C_Y=\frac{1}{n-1}Y{Y}^T=\frac{1}{n-1}(PX)(PX{)}^T=\frac{1}{n-1}PX{X}^T{p}^T=\frac{1}{n-1}P(X{X}^T)P^T$$ $$C_Y=\frac{1}{n-1}PA{P}^T$$

   其中$A$是一个对称阵，且$A\equiv X{X}^T$根据 奇异值分解中关于对称阵的分解，我们有$A=ED{E}^T$，其中D是特征值组成的对角阵，E为特征值对应的特征向量所组成矩阵。$A有$$r\le m$个标准正交特征向量（r是矩阵的秩），当A退化或所有数据都占据一个子空间$r\le m$时，r小于m，为保持正交性约束，可以通过选择$(m-r)$额外的标准正交向量(这些额外的向量不影响最终的解，因为与这些方向相关的方差为零)来“填充”矩阵E。

   假设P的每一行是$X{X}^T$的特征向量。有$P\equiv E^T$成立，将其带入$C_Y$的表达式得到$$A=P^{T}DP$$由$$P^{-1}=P^T$$来计算Cy$$C_Y=\frac{1}{n-1}PA{P}^T=\frac{1}{n-1}P(P^{T}DP)P^T=\frac{1}{n-1}(P{P}^T)D(P{P}^T)=\frac{1}{n-1}(P{P}^{-1})D(P{P}^{-1})$$ $$C_Y=\frac{1}{n-1}D$$

   显然将$C_Y$对角化就得到了P。总结一下我们可以得到：（1） $X$的主成分是$X{X}^T$的特征向量或者P的行（2）$C_Y$的第i个对角线值是$X$沿$Pi$方向的方差。

2. 使用SVD求解PCA
   关于SVD可见奇异值分解，根据公式求解即可。下面我们来看一下这两种方法有没有什么关联，假设X依然是一个$m\ast n$的矩阵，定义一个Y为$n\ast m$的新矩阵，其中Y的每一列均值为零$$Y\equiv \frac{1}{n-1}{X}^T$$

   接下来计算$Y^{T}Y$ $$Y^{T}Y=(\frac{1}{n-1}{X}^T)(\frac{1}{n-1}{X}^T)=\frac{1}{n-1}{X}^T{}^T{X}^T=\frac{1}{n-1}X{X}^T$$ $$Y^{T}Y=C_X$$

   这样就得到了关于X的协方差矩阵，可见两种方法是密切相关的。

PCA的算法步骤：

设有$m$条$n$维数据。

1）将原始数据按列组成$n$行$m$列矩阵$X$
2）将$X$的每一行（代表一个属性字段）进行零均值化，即减去这一行的均值
3）求出协方差矩阵$C=\frac{1}{m}XXT$
4）求出协方差矩阵的特征值及对应的特征向量
5）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前$k$行组成矩阵$P$
6）$Y=PX$即为降维到$k$维后的数据

对PCA的分析

1. 为什么通过PCA等到的K个主成分就可以很好的表示数据呢？
   通常，人们会发现前K个主成分具有较大方差，然后急剧下降。我们可以得出结论，我们关注的特征只存在于在前k维。

2. PCA有什么优缺点呢？
   它是一个无参数的方法，这是它的优点，同样也是它的缺点。因为它无需用户在算法运行前设置参数，减少了工作量，但同时如果用户预先知道了要处理的问题的一些特征，那么将其结合到任务的分析中是有一定意义的。

3. 上面我们所讨论的PCA是在线性的假设条件下进行的，如果是非线性的变换呢？这时就要使用在SVM中提到的核技巧，这种方法的PCA称为核PCA。常见的核变换包括傅里叶变换和高斯变换。这个过程是参数化的，因为用户在根据自己的先验知识选择核时，需要指定相关的参数，但它也是更优的。

参考：

周志华 《机器学习》第十章-降维与度量学习

https://blog.csdn.net/baimafujinji/article/details/79376378

https://www.cnblogs.com/hadoop2015/p/7419087.html

http://blog.codinglabs.org/articles/pca-tutorial.html