模式识别系列｜特征降维(1)主成分分析PCA

1-PCA概述

主成分分析是一种常用的降维方法，它不使用标签信息，通过将原始坐标空间的数据（$d\times 1$）投影到新的正交空间（$k\times 1$）中实现数据降维，所谓的主成分就是指数据在新空间的基的方向。PCA以方差作为信息损失衡量的标准，使得数据降维过程中信息损失最小，即降维后数据的方差要尽量大。PCA首先找到所有数据方差最大的方向，并将其作为新的坐标空间的第一个轴的方向，然后在这个方向的垂直超平面上寻找第二个投影后方差最大的方向，并作为新坐标空间第二个轴的方向，以此类推，直到找到需要的k个方向，也就是K个主成分，显然这k个新的基方向是两两垂直的。PCA的主要过程可以用“扭动坐标轴，保留K个轴”来形容。

为什么要以方差最大为依据呢？降维是为了数据更好地表示与计算，显然我们不希望降维后的数据成了一坨，使得原本分界明显的数据掺和在一起。例如，将数据投影到一维坐标系中，显然绿色的投影更好一些，因为其分散程度大，也就是方差更大。

对n个d维数据构成的数据集$X$（$d\times n$），要降维到k维的PCA过程如下：

1、将样本去均值得到新的$X$，即$X$的每一列减去$X$的按行求和取平均

2、计算协方差矩阵$S=\frac{1}{n}X{X}^T$

3、对协方差矩阵进行特征值特征向量分解

4、将特征值从大到小排序，将前K个特征值对应的特征向量组合成变换阵U（$d\times k$）

5、将样本降维$Y=U^{T}X$（$k\times n$）

2-理论推导

在推导之前，我们首先要明确三件事：一个是向量的内积与投影，第二个是基变换，第三个是协方差矩阵。

2.1-向量的内积与投影：

假设两个向量$a=(a_1,a_2,\cdots ,a_n{)}^T,b=(b_1,b_2,\cdots ,b_n{)}^T$，其内积表示为
$$a^{T}b=a_1{b}_1+a_2{b}_2+\cdots +a_n{b}_n=|a||b|\cos \theta$$其中$\theta$为两向量夹角。如图所示二维空间两个向量，向量a到b的投影应为b上由原点到绿线截断的部分，根据几何知识，投影长度为$|a|\cos \theta$，于是 当向量b模长为1时，向量a与向量b的内积就是向量a到向量b的投影的矢量长度。

2.2-基的表示与变换：

我们平时表示的诸如$(3,2{)}^T,(1,0{)}^T$等坐标都默认使用$(1,0{)}^T,(0,1{)}^T$的基，也就是说对平面上一个矢量，例如$(3,2{)}^T$实际上应该为$\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}3\\ 2\end{array}\right]$，其中左边表示基，右侧表示在这组基下的表示。

现在我们想将这个$(3,2{)}^T$在新的基下进行表示，其实只需要将原始基投影到新的基上即可，即点在新的基下的表示等于原始基到新基的转换乘以原始基乘以点在原始基下的表示。

例如，将原始坐标系中的点用蓝色线表示的基（45度）进行表示，蓝色线所示单位向量为$(\frac{\sqrt{2}}{2},\frac{\sqrt{2}}{2}{)}^T$和$(-\frac{\sqrt{2}}{2},\frac{\sqrt{2}}{2}{)}^T$。由2.1节可知，当b模长为1时，a到b的投影长度就是a与b的内积。因此，原始基到新的基的变换为${\left[\begin{array}{cc}\frac{\sqrt{2}}{2}& -\frac{\sqrt{2}}{2}\\ \frac{\sqrt{2}}{2}& \frac{\sqrt{2}}{2}\end{array}\right]}^T\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$，因此，原始坐标系下$(3,2{)}^T$点在新基的表示为
$${\left[\begin{array}{cc}\frac{\sqrt{2}}{2}& -\frac{\sqrt{2}}{2}\\ \frac{\sqrt{2}}{2}& \frac{\sqrt{2}}{2}\end{array}\right]}^T\left[\begin{array}{c}3\\ 2\end{array}\right]=\left[\begin{array}{c}\frac{5}{\sqrt{2}}\\ -\frac{1}{\sqrt{2}}\end{array}\right]$$综上，要得到原始空间一点在新空间的表示，只需要知道新空间各轴在原空间的单位矢量表示，然后将其组成矩阵U，U的转置与原空间点坐标的乘积即为新空间的表示。如例子中的$(\frac{\sqrt{2}}{2},\frac{\sqrt{2}}{2}{)}^T$和$(-\frac{\sqrt{2}}{2},\frac{\sqrt{2}}{2}{)}^T$。

注：左乘一个矩阵，既可以以基变换的角度去看，也可以以矢量运动的角度去看，这与在特征值与特征向量中的理解并不矛盾。

2.3-协方差矩阵：

协方差矩阵的每个元素是各个向量元素之间的协方差，因此它是对称矩阵，对角线元素为各维度的方差，其计算公式为
$$S=\frac{1}{n}\sum _{i=1}^n(x_i-\bar{x})(x_i-\bar{x}{)}^T=\left[\begin{array}{cc}\frac{1}{n}\sum _{i=1}^n{a}_i^2& \frac{1}{n}\sum _{i=1}^n{a}_i{b}_i\\ \frac{1}{n}\sum _{i=1}^n{a}_i{b}_i& \frac{1}{n}\sum _{i=1}^n{b}_i^2\end{array}\right]$$
其中$\bar{x}=\frac{1}{n}\sum _{i=1}^n{x}_i$为所有样本的均值，等式最右侧只表示了两个维度，意在说明协方差矩阵对角线元素为样本个维度的方差，其余元素为协方差，可以表示两个维度线性相关的程度，当两个向量垂直时，其协方差为0。

2.4-PCA推导

主成分分析的基本思想是以方差最大化为标准，寻找$k$个（小于d）正交轴，使得数据在其上投影得到新的数据。主要步骤是首先寻找一个样本投影后方差最大的方向，然后在该方向的垂直的超平面内寻找下一个方差最大的方向，直到找到所有$k$个方向。

我们首先寻找第一个方向$u_1$，它能使投影后的数据的方差最大。投影后，样本点$x_i$变为$u_1^T{x}_i$，计算投影后的方差
1 n ∑ i = 1 n { u 1 T x i − u 1 T x ˉ } 2 = 1 n ∑ i = 1 n u 1 T ( x i − x ˉ ) ( x i − x ˉ ) T u 1 = u 1 T S u 1 \frac{1}{n}\sum\limits _{i=1}^{n}\{u_1^Tx_i-u_1^T\bar{x}\}^2=\frac{1}{n}\sum\limits _{i=1}^{n}u_1^T(x_i-\bar{x})(x_i-\bar{x})^Tu_1=u_1^TSu_1 n1​i=1∑n​{u1T​xi​−u1T​xˉ}2=n1​i=1∑n​u1T​(xi​−xˉ)(xi​−xˉ)Tu1​=u1T​Su1​
为了找到最优方向，限制投影向量的模长为1，即$u_1^T{u}_1=1$。求解该限制条件下二次函数最优解，使用拉格朗日乘子法（关于拉格朗日乘子法的理解，参考另一博客），使拉格朗日乘数对$u_1$的导数为0，有
$$\frac{\partial \left(u_1^{T}S{u}_1+{\lambda }_1(1-u_1^T{u}_1)\right)}{\partial u_1}=0⟹S{u}_1={\lambda }_1{u}_1$$上式显然是协方差矩阵的特征值定义式（关于特征值特征向量的理解，可以参考另一博客），这是我们求解$u_1^{T}S{u}_1$最大时得到的条件，等式左右两边同时乘以$u_1^T$（由限制条件知$u_1^T=u_1^{-1}$）可得
$$u_1^{T}S{u}_1={\lambda }_1$$即最大方差就是协方差矩阵$S$的最大特征值。至此我们已经找到了第一个投影轴，即协方差矩阵最大特征值对应的特征向量。下面我们要找剩余$k-1$个投影轴，并且和第一个投影轴两两正交。接下来我们找第二个投影方向$u_2$，它能使投影后数据的方差最大，同理也满足上式，如果对投影轴不加限制，我们得到的仍然是最大特征值对应的特征向量。投影轴需要两两正交，又因为协方差矩阵的特征向量两两正交，因此第二个投影方向为协方差矩阵的第二大特征值对应的特征向量，以此类推。

注意：PCA寻找k个投影轴是寻找以原数据分别向第一个、第二个、第$n$个投影轴投影后的方差最大，用投影轴两两正交进行限制，保证不会得到$k$个相同的向量。

3-几何理解

PCA的过程可以用“旋转坐标系后去掉某些轴”来形容，旋转依据最大方差找到的投影轴来旋转，剔除是指去掉不需要的特征向量方向。

4-计算过程

对n个d维数据构成的数据集$X$（$d\times n$），要降维到k维的PCA过程如下：

1、将样本去均值得到新的$X$，即$X$的每一列减去$X$的按行求和取平均

2、计算协方差矩阵$S=\frac{1}{n}X{X}^T$

3、对协方差矩阵进行特征值特征向量分解

4、将特征值从大到小排序，将前K个特征值对应的特征向量组合成变换阵U（$d\times k$）

5、将样本降维$Y=U^{T}X$（$k\times n$）

如果不是任务规定$k$值大小，而自己又不知道$k$取多少合适，可以规定降维后方差损失不超过一个阈值，如信息损失不能超过0.01，根据如下公式确定$k$值。
$$1-\frac{\sum _{i=1}^k{\lambda }_i}{\sum _{i=1}^n{\lambda }_i}<0.01$$

4.1-样本数小于特征数时的计算

通常情况下，样本数n是大于样本维度d的，如果样本数n小于样本维度d，则上述计算过程应做相应改变，因为此时样本矩阵$X$秩最大为n，协方差矩阵$S=\frac{1}{n}X{X}^T$是一个$d\times d$的矩阵，显然不是满秩的，求解出的特征向量有d-n+1位为0。

此时，PCA过程应做如下修改：

2、计算协方差矩阵$S^{\prime }=\frac{1}{n}{X}^{T}X$

3、进行特征值特征向量分解得到${\lambda }_i$和$v_i$（虽然协方差矩阵变了，特征向量变了，但特征值没变）

4、计算得到原始需求的特征向量$u_i=\frac{1}{(n{\lambda }_i{)}^2}X{v}_i$，根据特征值排序，得到变换阵U

上述过程其实是利用了$X{X}^T$与$X^{T}X$特征值与特征向量的关系，二者具有相同的特征值。

4.2-matlab代码

[d, n] = size(X);
X = X - mean(X, 2);  %去均值
S = 1/n .* (X * X'); %协方差
[Ve, Va] = eig(S);   %特征值分解
[~,I] = sort(sum(Va), 'descend'); %特征值排序
k=2;                 %定义将到的维数
U = zeros(d, k);
for i = 1:k
   U(:,i) = Ve(:, I(i)) ; %取前k大特征值对应特征向量组成变换阵
end
Y = U' * X;          %投影

5-实例

我们构造一个二维的例子如下，将其降维到一维上。

样本生成代码如下：

for i = 1:10
   X(1,i) = i + 2*rand;
   X(2,i) = 3*i + 5 + 3*rand;
end

去均值后的样本、协方差矩阵、特征值、特征向量以及变换后的数据如下图所示：

最终实现的降维效果如下图所示，其中右图为一维表示。

同样地，我们将MNIST数据集（28*28）降维到3维和2维进行可视化如图所示：

由效果图可知，PCA对如此高维度进行降维还是有些吃力，但是还是可以有效地将数据进行聚类。

参考

[1] 通俗的主成分分析理解
[2] 国科大模式识别课程ppt(密码：uews)