机器学习算法之主成分分析(Principle Component Analysis, PCA)

基变换的矩阵表示

众所周知，在平面直角坐标系上，一个数据点或者说一个向量的坐标实际就是这个向量与两个坐标轴的单位向量做点积的结果，比如有一个向量 $\mathbf{x}$ 和 两个坐标轴基向量 $\mathbf{i}$ 和 $\mathbf{j}$:
$$\mathbf{x}=\left(\begin{array}{c}1\\ 2\end{array}\right),\mathbf{i}=\left(\begin{array}{c}1\\ 0\end{array}\right),\mathbf{j}=\left(\begin{array}{c}0\\ 1\end{array}\right)$$, $\mathbf{x}$ 的横坐标和纵坐标分别是：
$$\begin{array}{rl}{x}_1& ={\mathbf{x}}^T\mathbf{i}=1\\ x_2& ={\mathbf{x}}^T\mathbf{j}=2\end{array}$$
当然在平面上我们也可以选择另外两个线性无关的向量作为基向量比如取这么一组基向量 $\mathbf{i}$ 和 $\mathbf{j}$：
$$\mathbf{i}=\left(\begin{array}{c}\frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\end{array}\right),\mathbf{j}=\left(\begin{array}{c}-\frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\end{array}\right)$$
$\mathbf{x}$ 的在新基下的横坐标和纵坐标分别是：
$$\begin{array}{rl}{x}_1^{\prime }& ={\mathbf{x}}^T\mathbf{i}=\frac{5}{\sqrt{2}}\\ x_2^{\prime }& ={\mathbf{x}}^T\mathbf{j}=\frac{1}{\sqrt{2}}\end{array}$$
求新坐标的过程可以用矩阵运算来表示 ：
$$\left(\begin{array}{c}{\mathbf{i}}^T\\ {\mathbf{j}}^T\end{array}\right)\mathbf{x}=\left(\begin{array}{c}{x}_1^{\prime }\\ x_2^{\prime }\end{array}\right)$$即$$\left(\begin{array}{cc}\frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\\ -\frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\end{array}\right)\left(\begin{array}{c}1\\ 2\end{array}\right)=\left(\begin{array}{c}\frac{3}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\end{array}\right)$$
左边矩阵第一行是第一个基向量，第二行是第二个基向量，乘以原向量，结果就是新基下的坐标。假如我们有很多数据点$(1,2),(2,3),(3,4)$，则它们在新基下的坐标就是：
$$\left(\begin{array}{cc}\frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\\ -\frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\end{array}\right)\left(\begin{array}{ccc}1& 2& 3\\ 2& 3& 4\end{array}\right)=\left(\begin{array}{ccc}\frac{3}{\sqrt{2}}& \frac{5}{\sqrt{2}}& \frac{7}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\end{array}\right)$$
写成通用的形式就是：
$$\left(\begin{array}{c}{\mathbf{p}}_1^{\mathbf{T}}\\ {\mathbf{p}}_2^{\mathbf{T}}\\ ⋮\\ {\mathbf{p}}_{\mathbf{k}}^{\mathbf{T}}\end{array}\right)\left(\begin{array}{cccc}{\mathbf{a}}_1& {\mathbf{a}}_2& \dots & {\mathbf{a}}_{\mathbf{m}}\end{array}\right)=\left(\begin{array}{cccc}{\mathbf{p}}_1^{\mathbf{T}}{\mathbf{a}}_1& {\mathbf{p}}_1^{\mathbf{T}}{\mathbf{a}}_2& \dots & {\mathbf{p}}_1^{\mathbf{T}}{\mathbf{a}}_{\mathbf{m}}\\ {\mathbf{p}}_2^{\mathbf{T}}{\mathbf{a}}_1& {\mathbf{p}}_2^{\mathbf{T}}{\mathbf{a}}_2& \dots & {\mathbf{p}}_2^{\mathbf{T}}{\mathbf{a}}_{\mathbf{m}}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{p}}_{\mathbf{k}}^{\mathbf{T}}{\mathbf{a}}_1& {\mathbf{p}}_{\mathbf{k}}^{\mathbf{T}}{\mathbf{a}}_2& \dots & {\mathbf{p}}_{\mathbf{k}}^{\mathbf{T}}{\mathbf{a}}_{\mathbf{m}}\end{array}\right)$$
这里的 ${\mathbf{p}}_i$ 和 ${\mathbf{a}}_i$ 都是列向量。${\mathbf{a}}_i$ 表示第 $i$ 个样本数据看到这里我们就明白了:

一个矩阵乘以一个列向量可以解释为把该列向量变换到以矩阵的行向量为基向量的空间中去。

两个矩阵相乘就可以看作是把右边矩阵的每一个列向量变换到以左边矩阵的所有行向量组成的基向量的空间中去。

最大可分性

选择不同的基可以对同样一组数据给出不同的表示，如果基的数量少于向量本身的维数，则可以达到降维的效果。

那么如果我们想把一组 N 维的数据降到 K维，我们应该如何选择最优的 K 个基呢？或者说我们该如何选择 K 个基，才能尽可能多地保留原来数据的信息呢？

一种直观的看法是：希望投影后的投影值尽可能分散，因为如果重叠就会有样本消失。当然这个也可以从熵的角度进行理解，投影值越分散，投影值的方差就越大，熵就越大所含信息越多。

总体方差和样本方差

数值的分散可以用 方差来表示，方差分为总体方差和样本方差。在数据是一维的情况下，设有一个一维随机变量 $A$, 则总体方差是：
$$Var(A)=\frac{1}{m}\sum _{i=1}^m(a_i-\mu {)}^2$$
而为了保证样本方差是总体方差的无偏估计，样本方差的定义中，是除以 $m-1$ 而非样本数 $m$:
$$S(A)=\frac{1}{m-1}\sum _{i=1}^m(a_i-\mu {)}^2$$
我们要降维的数据属于样本，故应该用样本方差来衡量数据的分散程度。为了简便计算，我们对 数值先做归一化，使它们的均值为 $0$ :
$$S(A)=\frac{1}{m-1}\sum _{i=1}^m{a}_i^2\text{\hspace{1em}(1)}$$
我们的问题就转化为寻找一个基向量，使得样本点投影到该基向量方向上的投影值尽可能分散，即投影值的方差最大。

协方差

在一维的情况下，我们用方差来衡量随机变量自身的分散程度，如果有多个随机变量，我们用协方差来衡量它们之间的相关性。， A 和 B 两个随机变量的协方差是：
$$Cov(\mathbf{A},\mathbf{B})=\frac{1}{m-1}\sum _{i=1}^m(a_i-{\mu }_a)(b_i-{\mu }_b)$$
如果协方差为0，则表示这两个一维随机变量是线性无关的。
为了简便，我们同样先对数据做零均值化：
$$Cov(\mathbf{A},\mathbf{B})=\frac{1}{m-1}\sum _{i=1}^m{a}_i{b}_i\text{\hspace{1em}(2)}$$
协方差为零表示两个随机变量线性无关。 为了使保留下来的信息尽可能多，不冗余，我们选择的 K 个基向量应该是线性无关的，即相互正交。一般我们会把基向量归一化成单位向量。我们的目标就是即把一组N维数据降到K维，需要找到K个N维正交单位向量，使原始数据变换到这组基向量空间后，各个随机变量之间的协方差为0，而自身的变量方差最大。 这里所说的变量就是数据的特征，数据从N维降到K维就是从N个数据特征降到K个数据特征的过程。

协方差矩阵

设有两个特征，即两个随机变量 $\mathbf{A}$ 和 $\mathbf{B}$ , 有 $m$ 个样本，把它们用矩阵表示，每列是一个样本，每行是一个特征：
$$X=\left(\begin{array}{c}{\mathbf{a}}^T\\ {\mathbf{b}}^T\end{array}\right)=\left(\begin{array}{cccc}{a}_1& a_2& \dots & a_m\\ b_1& b_2& \dots & b_m\end{array}\right)$$
其中 $\mathbf{a}$ 和 $\mathbf{b}$ 分别是：
$$\mathbf{a}=\left(\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_m\end{array}\right),\mathbf{b}=\left(\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_m\end{array}\right)$$
相当于我们现在有一个随机向量 $\mathbf{v}$：
$$\mathbf{v}=\left(\begin{array}{c}A\\ B\end{array}\right)$$
然后有m个服从$\mathbf{v}$的分布的样本, 每列一个样本组成了 $X$。在多维空间下，我们用协方差矩阵来衡量随机向量的各个维度对应的随机变量之间的协方差。当然了一维随机变量与自身的协方差就是方差即$Cov(A,A)=S(A)$。协方差矩阵公式如下：
$$Cov(\mathbf{v})=\left(\begin{array}{cc}Cov(A,A)& Cov(A,B)\\ Cov(B,A)& Cov(B,B)\end{array}\right)$$
考虑$(1)$ 和 $(2)$,有：
$$\begin{array}{rl}Cov(\mathbf{v})& =\left(\begin{array}{cc}\frac{1}{m-1}\sum _{i=1}^m{a}_i^2& \frac{1}{m-1}\sum _{i=1}^m{a}_i{b}_i\\ \frac{1}{m-1}\sum _{i=1}^m{b}_i{a}_i& \frac{1}{m-1}\sum _{i=1}^m{b}_i^2\end{array}\right)\\ & =\frac{1}{m-1}\left(\begin{array}{cc}{\mathbf{a}}^{\mathbf{T}}\mathbf{a}& {\mathbf{a}}^{\mathbf{T}}\mathbf{b}\\ {\mathbf{b}}^{\mathbf{T}}\mathbf{a}& {\mathbf{b}}^{\mathbf{T}}\mathbf{b}\end{array}\right)\\ & =\frac{1}{m-1}\left(\begin{array}{c}{\mathbf{a}}^T\\ {\mathbf{b}}^T\end{array}\right)\left(\begin{array}{cc}\mathbf{a}& \mathbf{b}\end{array}\right)\\ & =\frac{1}{m-1}X{X}^T\end{array}$$
至此，我们得到了原始数据矩阵和随机向量$\mathbf{v}$的协方差矩阵的关系式：
$$Cov(\mathbf{v})=\frac{1}{m-1}X{X}^T$$
目前的随机向量 $\mathbf{v}$ 是2维的，当它是$n$ 维的时候，显然也有：
设有一个 $n$ 维随机向量$\mathbf{v}$, 它每一维都是一个一维随机变量，有 $m$ 个服从 $\mathbf{v}$ 分布的 $n$ 维向量按照每列一个 $n$ 维向量共 $m$ 列组成原始数据矩阵 $X$，则随机向量 $v$ 的协方差矩阵 $C$ 和数据矩阵 $X$ 的关系满足 :
$$C=\frac{1}{m-1}X{X}^T\text{\hspace{1em}(3)}$$

矩阵对角化

现在我们的目标就是使得基变换后降维的随机变量的协方差矩阵的对角线元素尽可能的大，而对角线以外的元素为0，这样就能保证在各个投影轴上的投影值尽可能分散，而各个投影轴直接是线性无关即相互正交，使得降维后的各个新特征的信息不冗余。设我们的基变换矩阵是 $P$, 原始数据矩阵 $X$ 变换后的新数据矩阵是 $Y=PX$ , 变换后新随机向量 $\mathbf{v}\prime$ 的协方差矩阵是 $D$, 我们有：
$$\begin{array}{rl}D& =\frac{1}{m-1}Y{Y}^T\\ & =\frac{1}{m-1}(PX)(PX{)}^T\\ & =\frac{1}{m-1}PX{X}^T{P}^T\\ & =P(\frac{1}{m-1}X{X}^T)P^T\\ & =PC{P}^T\end{array}$$
到这里我们发现，我们要找的基变换矩阵 $P$ 是把原协方差矩阵 $C$ 对角化并且使得 $PC{P}^T$ 的对角线元素从从上往下是从大到小排列的矩阵，如果我们要把原来的数据从 $n$ 维降到 $k$ 维，那么只要取 $P$ 的前 $k$ 行作为降维的基变换矩阵即可达到我们的目标。
那么我们该如何找到矩阵 $P$ 呢？协方差矩阵 $C$ 是一个实对称矩阵， 实对称矩阵有以下性质：

1. 实对称矩阵的不同特征值对应的特征向量相互正交
2. 设特征向量 $\lambda$ 的重数为 $r$ ，必定存在 $r$ 个线性无关的特征向量对应于 $\lambda$，因此可以将这 $r$ 个线性无关的特征向量单位正交化。

注意我们的协方差矩阵 $C$ 是大小是 $n\times n$ ，对于这样的实对称矩阵一定可以对角化，我们一定可以找到 $n$ 个正交单位向量 ${\mathbf{e}}_1,{\mathbf{e}}_2,\dots ,{\mathbf{e}}_{\mathbf{n}}$, 把它们按列组成矩阵 $E=({\mathbf{e}}_1,{\mathbf{e}}_2,\dots ,{\mathbf{e}}_{\mathbf{n}})$,则我们有以下结论：
$$E^{T}CE=\Lambda =\left(\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & \ddots & \\ & & & {\lambda }_n\end{array}\right)$$
其中 $\Lambda$ 是对角阵，对角线上的元素是 $C$ 的特征值。现在我们可以发现，$\Lambda$ 实际就是 $D$, 其对角线元素即 $C$ 的每个特征值就是投影后在各个轴上的方差：
$$P=E^T\text{\hspace{1em}(4)}$$
$P$ 是协方差矩阵 $C$ 的特征向量单位化后按行排列出的矩阵，其中每一行都是 $C$ 的一个特征向量。如果设 $P$ 按照 $\Lambda$ 中特征值从大到小，将特征向量从上到下排列，则用 P 的前 K 行组成的矩阵乘以原始数据矩阵 $X$，就得到了我们需要的降维后的数据矩阵 $Y$。

拉格朗日乘子法

我们也可以用拉格朗日乘子法来推导找到 $P$ ，设原数据矩阵 $$X=({\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_{\mathbf{n}})$$
样本点 ${\mathbf{x}}_i$ 在某个单位正交基 ${\mathbf{w}}_{\mathbf{j}}$ 方向上的投影值是 ${\mathbf{x}}_{\mathbf{i}}^{\mathbf{T}}{\mathbf{w}}_{\mathbf{j}}$, 我们有原始数据在 ${\mathbf{w}}_{\mathbf{j}}$ 方向上的投影值的方差：
$$\begin{array}{rl}{D}_j& =\frac{1}{m-1}\sum _{i=1}^m({\mathbf{x}}_{\mathbf{i}}^{\mathbf{T}}{\mathbf{w}}_{\mathbf{j}}{)}^2\\ & =\frac{1}{m-1}\sum _{i=1}^m({{\mathbf{x}}_{\mathbf{i}}^{\mathbf{T}}{\mathbf{w}}_{\mathbf{j}})}^T({\mathbf{x}}_{\mathbf{i}}^{\mathbf{T}}{\mathbf{w}}_{\mathbf{j}})\\ & =\frac{1}{m-1}\sum _{i=1}^m{\mathbf{w}}_{\mathbf{j}}^{\mathbf{T}}{\mathbf{x}}_{\mathbf{i}}{\mathbf{x}}_{\mathbf{i}}^{\mathbf{T}}{\mathbf{w}}_{\mathbf{j}}\\ & ={\mathbf{w}}_{\mathbf{j}}^{\mathbf{T}}(\frac{1}{m-1}\sum _{i=1}^m{\mathbf{x}}_{\mathbf{i}}{\mathbf{x}}_{\mathbf{i}}^{\mathbf{T}}){\mathbf{w}}_{\mathbf{j}}\\ & ={\mathbf{w}}_{\mathbf{j}}^{\mathbf{T}}(\frac{1}{m-1}X{X}^T){\mathbf{w}}_{\mathbf{j}}\\ & ={\mathbf{w}}_{\mathbf{j}}^{\mathbf{T}}C{\mathbf{w}}_{\mathbf{j}}\end{array}$$
现在我们的目标和条件是：
{ m a x { w j T C w j } w j T w j = 1 \begin{cases} max \{\bm{w_j^T}C\bm{w_j}\} \\ \bm{w_j^Tw_j=1} \end{cases} {max{wjT​Cwj​}wjT​wj​=1​
用拉格朗日函数表示：
$$L({\mathbf{w}}_{\mathbf{j}})={\mathbf{w}}_{\mathbf{j}}^{\mathbf{T}}C{\mathbf{w}}_{\mathbf{j}}+\lambda (1-{\mathbf{w}}_{\mathbf{j}}^{\mathbf{T}}{\mathbf{w}}_{\mathbf{j}})$$
对 ${\mathbf{w}}_{\mathbf{j}}$ 求导， 令其为 $0$ ：
$$C{\mathbf{w}}_{\mathbf{j}}=\lambda {\mathbf{w}}_{\mathbf{j}}$$
这里表明，${\mathbf{w}}_{\mathbf{j}}$是 $C$ 的特征向量，$\lambda$ 是对应的特征值。而且有：
$$D_j={\mathbf{w}}_{\mathbf{j}}^{\mathbf{T}}C{\mathbf{w}}_{\mathbf{j}}={\mathbf{w}}_{\mathbf{j}}^{\mathbf{T}}\lambda {\mathbf{w}}_{\mathbf{j}}=\lambda {\mathbf{w}}_{\mathbf{j}}^{\mathbf{T}}{\mathbf{w}}_{\mathbf{j}}=\lambda$$
我们发现投影后的最大方差就是基向量 ${\mathbf{w}}_{\mathbf{j}}$ 对应的特征值。我们要找的 $k$ 个投影用的基向量就是 $C$ 的最大的 $k$ 个特征值对应的单位特征向量。

总结PCA的步骤

PCA的求解步骤如下：

1.把 $m$ 个 $n$ 维向量表示的样本按列组成原始数据矩阵 $X$

2.把 $X$ 的每个特征，即每一行进行零均值化，即每个元素减去这行的均值

3.求出协方差矩阵 $C=\frac{1}{m-1}X{X}^T$ 的特征值和对应的单位特征向量

4.把 $C$ 的单位特征向量按照特征值从大到小，按行组成矩阵，取前 $k$ 行组成基变换矩阵 $P$

5.降维后的 $k$ 行 $m$ 列数据矩阵就是 $Y=PX$

PCA的特点

1. 缓解维度灾难：PCA降维后，使得样本的采样密度增大了。
2. 降噪：当数据收到噪声影响时，最小的特征值对应的特征向量往往与噪声相关，降维的过程中将它们舍弃，可以起到降噪的效果。
3. 过拟合：由于PCA舍去了一些信息，而这些信息可能是分类或者回归需要的关键信息，只是在训练集上没有重要表现，导致在训练集上过拟合。
4. 特征独立：PCA不仅使数据降到低维，还使得降维之后的数据特征相互独立。

细节

零均值化

当使用PCA对训练集进行降维时，也需要对验证集和测试集做同样的PCA。而且对验证集、测试集零均值化用的均值必须是训练集的均值，因为验证集和测试集是不可观测的。

PCA和SVD的关系

PCA本质是特征值分解，而特征值和特征向量只有方阵才有，而对与任意矩阵都可以做奇异值分解。

PCA

对于一个实对称方阵 $A$ ，总可以对角化:
$$A=P\Lambda P^{-1}$$
如果 $P$ 是正交矩阵，有 $P^T=P^{-1}$, 进而有：
$$A=P\Lambda P^T$$

SVD

一个矩阵$A_{n,m}$ 可以 进行奇异值分解：
$$A_{n,m}=U_{n,n}{\Sigma }_{n,m}{V^T}_{m,m}\approx U_{n,k}{\Sigma }_{k,k}{V^T}_{k,m}$$
这里的 $U_{n,n}$ 和 $V_{m,m}$ 都是正交矩阵,。 满足 $U_{n,n}{U^T}_{n,n}=I_n$, $V_{m,m}{V^T}_{m,m}=I_m$。后面的约等于是因为 $\Sigma$ 有多个奇异值，对角线上的奇异值按从大到小从上往下排列，那么后面很多其实很小接近于 $0$, 不妨近似取前 $k$ 个较大的。
$$\begin{array}{rl}{A}^{T}A& =(U\Sigma V^T{)}^T(U\Sigma V^T)=V{\Sigma }^2{V}^T\\ A{A}^T& =(U\Sigma V^T)(U\Sigma V^T{)}^T=U{\Sigma }^2{U}^T\end{array}$$可见， $V$ 的每一列是 $A^{T}A$ 的特征向量，$U$ 的每一列是 $A{A}^T$ 的特征向量。${\Sigma }^2$ 的对角线元素是 $A^{T}A$ 和 $A{A}^T$ 的特征值。记它们的特征值矩阵维 $\Lambda$:
$$\Lambda ={\Sigma }^2$$ 可见 $A^{T}A$ 和 $A{A}^T$ 的特征值是 $A$ 的对应奇异值的平方。
在PCA中，我们需要对协方差矩阵 $C=\frac{1}{m-1}X{X}^T$ 进行特征值分解，如果令 $$A=\frac{1}{\sqrt{m-1}}X$$, 则 $C=A{A}^T$, 对 $A$ 进行奇异值分解可以等价于对 $C$ 进行特征值分解。
可以用 SVD 替代来求解 PCA：

1. 当样本维度很高时，协方差矩阵计算太慢；
2. 方阵特征值分解计算效率不高；
3. SVD 出了特征值分解这种求解方式外，还有更高效更准球的迭代求解方式，避免了 $X{X}^T$ 的计算；
4. 其实 PCA 与 SVD 的右奇异向量的压缩效果相同。

Reference

1. 【机器学习】降维——PCA（非常详细）, 阿泽, 知乎

2. 从SVD到PCA——奇妙的数学游戏 ，Lee的白板报，开源中国

3. Statistics-and-PCA， MIT