张小贤的平凡之路---机器学习之PCA降维---周志华《机器学习》

1. 维数灾难

许多学习算法都涉及距离计算，而高维空间会给距离计算带来很大的麻烦，例如当维数很高时甚至连计算内积都不再容易。事实上，在高维情形下出现的数据样本稀疏、距离计算困难等问题，是所有机器学习方法共同面临的严重障碍，被称为“维数灾难”（curse of dimensionality）。缓解维数灾难的一个重要途径是降维（dimension reduction），亦称“维数简约”：通过某种数学变换将原始高维属性空间转变为一个低维“子空间”（subspace），在这个子空间中样本密度大幅提高，距离计算也变得更为容易。

降维、

123

2. 主成分分析（principal component analisys, PCA）

主成分分析是最常用的一种降维方法，也是一种常用的无监督学习方法，这一方法利用正交变换把由线性相关变量表示的观测数据转换为少数几个由线性无关变量表示的数据，线性无关的变量称为主成分。主成分的个数通常小于原始变量的个数，所以主成分分析属于降维方法。主成分分析主要用于发现数据中的基本结构，即数据中变量之间的关系，是数据分析的有力工具，也用于其他机器学习方法的前处理。

统计分析中，数据的变量之间可能存在相关性，以致增加了分析的难度。于是，考虑由少数不相关的变量来代替相关的变量，用来表示数据，并且要求能够保留数据中的大部分信息。

如何表示数据之间的相关性呢？可以采用相似性矩阵：

2.1 相似性矩阵

若两个随机变量$X$和$Y$相互独立，则$E[(X-E(X))(Y-E(Y))]=0$，因而若上述数学期望不为零，则$X$和$Y$必不是相互独立的，亦即它们之间存在着一定的关系。
协方差与方差之间有如下关系：
$D(X+Y)=D(X)+D(Y)+2Cov(X，Y)$
$D(X-Y)=D(X)+D(Y)-2Cov(X，Y)$
协方差与期望值有如下关系：
$Cov(X，Y)=E(XY)-E(X)E(Y)$
协方差的性质：
（1）$Cov(X，Y)=Cov(Y，X)$；
（2）$Cov(aX，bY)=abCov(X，Y)$，（$a$，$b$是常数）；
（3）$Cov(X1+X2，Y)=Cov(X1，Y)+Cov(X2，Y)$。
由协方差定义，可以看出$Cov(X，X)=D(X)$，$Cov(Y，Y)=D(Y)$。
协方差作为描述$X$和$Y$相关程度的量，在同一物理量纲之下有一定的作用，但同样的两个量采用不同的量纲使它们的协方差在数值上表现出很大的差异。为此引入如下概念：
定义${\rho }_{XY}=\frac{Cov(X,Y)}{\sqrt{D(X)}\sqrt{D(Y)}}$称为随机变量$X$和$Y$的$(Pearson)$相关系数。
定义
若${\rho }_{XY}=0$，则称X与Y不线性相关。
即${\rho }_{XY}=0$的充分必要条件是$Cov(X，Y)=0$，亦即不相关和协方差为零是等价的。

由相关系数的表达式可以看出，两个变量之间的相关系数可以通过对数据进行规范化，使得各自的数据方差为$1$，然后采用该两个变量之间的协方差矩阵来表示。

2.2 协方差矩阵

主成分分析需要回答的问题是如何从多个变量中综合为少数几个代表性变量，既能够代表原始变量的绝大多数信息，又互不相关，并且在新的综合变量基础上，可以进一步的统计分析。考虑下图：

该图是一个二维的图像，若需要降维，将二维的图像降为一维的，采用哪种度量会最大程度地保留原来的信息呢？很明显的是采用$x$轴，因为该两点在$x$轴上的投影之间的差距能够完全反应原来两个数据之间的差距；为什么不选择$y$轴呢？因为在$y$轴上的投影无法反应两个数据点之间的差距。而这种差距怎么形容呢？可以采用方差来形容，即该两点在$x$轴投影的方差很大，而在$y$轴上的方差很小。这个例子非常巧合，巧合的是两点的投影刚刚好和原来的坐标系重合了。

如上图所示，这时候如何投影呢？应当将各数据点投影到 $y-x=0$ 这条直线上，因为在该直线上的投影能够反映数据之间的差距，这时候投影的方向和原来的坐标轴完全不重合了。

那么问题又来了，是怎么寻找到这条 $y-x=0$ 的投影直线呢？因为实际上可以投影的直线有成千上万条，为什么说 $y-x=0$ 是我们需要寻找的那条呢？

2.3 投影与内积

我们先应该分析投影的数学意义，投影实际上就是求两个向量的内积。以二维向量为例，

内积运算将两个向量映射为一个实数。其计算方式非常容易理解，但是其意义并不明显。下面我们分析内积的几何意义。假设$A$和$B$是两个$n$维向量，我们知道n维向量可以等价表示为$n$维空间中的一条从原点发射的有向线段，为了简单起见我们假设$A$和$B$均为二维向量，则，。则在二维平面上$A$和$B$可以用两条发自原点的有向线段表示，现在我们从$A$点向$B$所在直线引一条垂线。我们知道垂线与B的交点叫做$A$在$B$上的投影，再设$A$与$B$的夹角是$a$，则投影的矢量长度为$|A|cos(a)$，其中$|A|$是向量$A$的模，也就是$A$线段的标量长度。注意这里我们专门区分了矢量长度和标量长度，标量长度总是大于等于$0$，值就是线段的长度；而矢量长度可能为负，其绝对值是线段长度，而符号取决于其方向与标准方向相同或相反。现在事情似乎是有点眉目了：$A$与$B$的内积等于$A$到$B$的投影长度乘以$B$的模。再进一步，如果我们假设$B$的模为$1$，即让，那么就变成了：也就是说，设向量$B$的模为$1$，则$A$与$B$的内积值等于$A$向$B$所在直线投影的矢量长度！这就是内积的一种几何解释，也是我们得到的第一个重要结论。

也即是说，投影依然能够反映原始数据之间的差异。

那么投影方向怎么找呢？

在主成分分析中，首先对给定数据进行规范化，使得数据每一变量的平均值为$0$，方差为$1$。之后对数据进行正交变换，原来由线性相关变量表示的数据，通过正交变换变成由若干个线性无关的新变量表示的数据。新变量是可能的正交变换中变量的方差和（信息保存）最大的，方差表示在新变量上信息的大小。将新变量依次称为第一主成分、第二主成分等。这就是主成分分析的基本思想。通过主成分分析，可以利用主成分近似地表示原始数据，这可理解为发现数据的基本结构；也可以把数据由少数主成分表示，这可理解为对数据降维。

上面的陈述很好地解释了这个方向应该怎么找，即：寻找方差最大的那个方向。因此需要对协方差矩阵进行求解特征值等操作。

选择方差最大的那个方向投影是因为方差最大可以很好地保留/映射原数据的分布。协方差最大的方向，表示该维度上的变量相关性越强，因此可以采用一个新变量来表示原来两个向量，而这个新向量就是投影的方向。

2.3.1 方差最大的解释

假设有两个变量$x_1$和$x_2$，三个样本点$A$、$B$、$C$，样本分布在由$x_1$和$x_2$轴组成的坐标系中，如下图所示。对坐标系进行旋转变换，得到新的坐标轴$y_1$，表示新的变量$y_1$。样本点$A$、$B$、$C$在$y_1$轴上投影，得到$y_1$轴的坐标值$A^{\prime }$、$B^{\prime }$、$C^{\prime }$。坐标值的平方和$O{A}^{\prime 2}+O{B}^{\prime 2}+O{C}^{\prime 2}$表示样本在变量$y_1$上的方差和。主成分分析旨在选取正交变换中方差最大的变量，作为第一主成分，也就是旋转变换中坐标值的平方和最大的轴。注意到旋转变换中样本点到原点的距离的平方和$O{A}^2+O{B}^2+O{C}^2$保持不变，根据勾股定理，坐标值的平方和$O{A}^{\prime 2}+O{B}^{\prime 2}+O{C}^{\prime 2}$最大等价于样本点到$y_1$轴的距离的平方和$A{A}^{\prime 2}+B{B}^{\prime 2}+C{C}^{\prime 2}$最小。所以，等价地，主成分分析在旋转变换中选取离样本点的距离平方和最小的轴，作为第一主成分，第二主成分在保证与已选坐标值正交的条件下，类似地进行。

这段陈述解释了为什么第二主成分需要和第一主成分正交的原因，也可以通过特征值求解的概念来理解。

2.3.2 从数学的角度理解降维

由上一小节的方差最大的解释可以很好地理解：实际上降维要求原始空间中样本之间的距离在低维空间中得以继续保持，即得到“多维缩放”（Multiple Dimensinal Scaling，简称MDS）。

假定$m$个样本在原始空间的距离矩阵为$D\in {\mathbb{R}}^{m\times m}$，其第$i$行第$j$列的元素$dis{t}_{ij}$为样本$x_i$到$x_j$的距离。我们的目标是获得样本在$d^{\prime }$维空间的表示$Z\in {\mathbb{R}}^{d^{\prime }\times m},d^{\prime }⩽d$，且任意两个样本在$d^{\prime }$维空间中的欧式距离等于原始空间中的距离，即$z_i-z_j=dis{t}_{ij}$.

令$B=Z^{T}Z\in {\mathbb{R}}^{m\times m}$，其中$B$为降维后样本的内积矩阵，$b_{ij}=z_i^T{z}_j$，有
$$dis{t}_{ij}^2=||z_i|{|}^2+||z_j|{|}^2-2z_i^T{z}_j=b_{ii}+b_{jj}-2b_{ij}$$
为便于讨论，令降维后的样本$Z$被中心化，即$\sum _{i=1}^m{z}_i=0$。显然，矩阵$B$的行与列之和均为零，即$\sum _{i=1}^m{b}_{ij}=\sum _{j=1}^m{b}_{ij}=0$.易知：
$$\sum _{i=1}^{m}dis{t}_{ij}^2=tr(\text{B})+m{b}_{jj},$$
$$\sum _{j=1}^{m}dis{t}_{ij}^2=tr(\text{B})+m{b}_{ii},$$
$$\sum _{i=1}^m\sum _{j=1}^{m}dis{t}_{ij}^2=2m\times tr(\text{B}),$$
其中，$tr(\bullet )$表示矩阵的迹(trace)，$tr(\text{B})=\sum _{i=1}^m||z_i|{|}^2$，令
$$dis{t}_{i.}^2=\frac{1}{m}\sum _{j=1}^{m}dis{t}_{ij}^2,$$
$$dis{t}_{.j}^2=\frac{1}{m}\sum _{i=1}^{m}dis{t}_{ij}^2,$$
$$dis{t}_{..}^2=\frac{1}{m^2}\sum _{i=1}^m\sum _{j=1}^{m}dis{t}_{ij}^2,$$
可得：
$$b_{ij}=-\frac{1}{2}(dis{t}_{ij}^2-dis{t}_{i.}^2-dis{t}_{.j}^2+dis{t}_{..}^2)$$
由此即可通过降维前后保持不变的距离矩阵$D$求取内积矩阵$B$.

对矩阵$B$做特征值分解， $B=V\Lambda V^T$，其中$\Lambda =diag({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_{d\ast })$
$$Z={\lambda }_{\ast }^{1/2}{V}_{\ast }^T\in d\ast \times m$$

2.4 PCA步骤

实际上$PCA$是一种通过线性变换进行的降维方法，投影实际上是一种线性映射，即通过向量之间的内积。类似地，采用非线性变换进行降维的方法有核化线性降维等。

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

1）将原始数据按列组成$n$行$m$列矩阵$X$

2）将$X$的每一行（代表一个属性字段）进行零均值化，即减去这一行的均值

3）求出协方差矩阵

4）求出协方差矩阵的特征值及对应的特征向量

5）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前$k$行组成矩阵$P$

6）即为降维到$k$维后的数据

根据上面对PCA的数学原理的解释，我们可以了解到一些PCA的能力和限制。PCA本质上是将方差最大的方向作为主要特征，并且在各个正交方向上将数据“离相关”，也就是让它们在不同正交方向上没有相关性。

因此，PCA也存在一些限制，例如它可以很好的解除线性相关，但是对于高阶相关性就没有办法了，对于存在高阶相关性的数据，可以考虑$Kernel$ PCA，通过$Kernel$函数将非线性相关转为线性相关，关于这点就不展开讨论了。另外，PCA假设数据各主特征是分布在正交方向上，如果在非正交方向上存在几个方差较大的方向，PCA的效果就大打折扣了。最后需要说明的是，PCA是一种无参数技术，也就是说面对同样的数据，如果不考虑清洗，谁来做结果都一样，没有主观参数的介入，所以PCA便于通用实现，但是本身无法个性化的优化。

2.5 编程实现

2.5.1 基础PCA

[Row Col]=size(X);
covX=cov(X);                 %求样本的协方差矩阵(散步矩阵除以(n-1)即为协方差矩阵)
[V D]=eigs(covX);                        %求协方差矩阵的特征值D和特征向量V
meanX=mean(X);                        样本均值m
% 所有样本X减去样本均值m，再乘以协方差矩阵（散步矩阵）的特征向量V，
% 即为样本的主成份SCORE
tempX= repmat(meanX,Row,1);
SCORE2=(X-tempX)*V                  %主成份：SCORE
pcaData2=SCORE2(:,1:k)

2.5.2 快速PCA（fast-PCA）

%% 使用快速PCA算法实现的方法
[pcaData3 COEFF3] = fastPCA(X, k )
%   其中fastPCA函数的代码实现如下：
function [pcaA V] = fastPCA( A, k )
% 快速PCA
% 输入：A --- 样本矩阵，每行为一个样本
%      k --- 降维至 k 维
% 输出：pcaA --- 降维后的 k 维样本特征向量组成的矩阵，每行一个样本，
% 列数 k 为降维后的样本特征维数
%      V --- 主成分向量
[r c] = size(A);  % 样本均值
meanVec = mean(A);   % 计算协方差矩阵的转置 covMatT
Z = (A-repmat(meanVec, r, 1));
covMatT = Z * Z';   % 计算 covMatT 的前 k 个本征值和本征向量
[V D] = eigs(covMatT, k);   % 得到协方差矩阵 (covMatT)' 的本征向量
V = Z' * V;   % 本征向量归一化为单位本征向量
for i=1:k
    V(:,i)=V(:,i)/norm(V(:,i));
end     % 线性变换（投影）降维至 k 维
pcaA = Z * V;   % 保存变换矩阵 V 和变换原点 meanVec

Ref1
Ref2
Ref3