主成分分析(PCA)直观理解与数学推导

近期在完成信息论的作业，发现网上的资料大多是直观解释，对其中的数学原理介绍甚少，并且只介绍了向量降维，而没有介绍向量重构的问题（重构指的是：根据降维后的低维向量来恢复原始向量），因此在这里做一个总结，配合另一篇博客看效果可能会更好参考博客。

简介

PCA降维就是要把$m$维空间的n个样本点$x_i,i=1\cdots n$映射成$l$维的低维空间中的向量$y_i,i=1\cdots n$，其中$l\ll m$。这个映射不是随意的，而是要确保利用低维的向量$y_i$重构出的$x_i^{{}^{\prime }}$与原向量$x_i$的误差最小化。这个过程可以描述成
$$y_{l\times 1}=W_{l\times m}{x}_{m\times 1}$$
于是PCA的任务就是找到这样一个W矩阵。

算法数学推导

我们考虑一个二维空间中的数据压缩成一维的问题：
给定五个样本点$x_i,i=1\cdots 5$，每个样本点都是一个二维的列向量，把它们拼成一个矩阵$X=(x_1,x_2,x_3,x_4,x_5)=\left[\begin{array}{ccccc}1& 1& 2& 4& 2\\ 1& 3& 3& 4& 4\end{array}\right]$。
这个例子和文章开头提到的那篇博客中的例子是一致的，可以参考那篇博客的图片。
我们知道每个向量是通过一组基和在这组基下的坐标来描述的，例如$x_1=[1,1{]}^T$是该向量在单位正交基${\xi }_1=[0,1{]}^T,{\xi }_2=[1,0{]}^T$下的坐标表示，坐标实际上就是向量向各个基上的投影值。同样的，如果我们要将一个二维向量压缩成一维向量，那么只需要找到一条直线，直线的单位方向向量$w$作为基，然后用向量向这条直线的投影值就可以描述压缩后的一维向量。

一条直线可以用它经过的点$\mu$和单位方向向量$w$来描述，即$x=\mu +\alpha w$，（这里我们用的$\mu$是上述五个样本点的均值$[2,3{]}^T$）这里的$\alpha$就可以理解为坐标，$w$是基。那么我们要寻找的二维向量$x_i$经过压缩后的一维向量其实就是${\alpha }_i$，现在需要确定$w$，这样才能求出二维向量向直线的投影值。正如前一节所述，这个$w$不是任意的，而是应该确保重构后的误差最小化，它可以描述成如下的优化问题：
$$\underset{w}{\min }f(w)=\frac{1}{2}\sum _{i=1}^n{\Vert \left(\mu +{\alpha }_{i}w\right)-x_i\Vert }_2^2$$
先求出${\alpha }_i$的取值，
$$\frac{\partial f}{\partial {\alpha }_i}=w^T(\mu +{\alpha }_{i}w-x_i)=0$$
由于$w$是单位向量，即$w^{T}w=1$，由上式可得：
$${\alpha }_i=w^T(x_i-\mu )$$
这个公式就给出了${\alpha }_i$的求解方法。下面继续推导确定$w$的过程，定义散布矩阵如下：
$$S=\sum _{i=1}^n(x_i-\mu )(x_i-\mu {)}^T$$
对于上面的例子，我们把$X$中的每个样本$x_i$都减去均值$\mu$得到一个新的矩阵记为$\tilde{X}=\left[\begin{array}{ccccc}-1& -1& 0& 2& 0\\ -2& 0& 0& 1& 1\end{array}\right]$，那么上面的散步矩阵其实可以简单地记为
$$S=\tilde{X}{\tilde{X}}^T$$
说明它是一个对称矩阵。
将上面求出的${\alpha }_i$的表达式代入到$f(w)$中，得到：
$$\begin{array}{rl}f(\mathbf{w})& =\frac{1}{2}\sum _{i=1}^n{\Vert {\alpha }_i\mathbf{w}-\left({\mathbf{x}}_i-\mathbf{\mu }\right)\Vert }_2^2\\ & =\frac{1}{2}\left(\sum _{i=1}^n{\alpha }_i^2\Vert \mathbf{w}{\Vert }_2^2-2\sum _{i=1}^n{\alpha }_i{\mathbf{w}}^{\mathrm{T}}\left({\mathbf{x}}_i-\mathbf{\mu }\right)+\sum _{i=1}^n{\Vert {\mathbf{x}}_i-\mathbf{\mu }\Vert }_2^2\right)\\ & =-\frac{1}{2}\sum _{i=1}^n{\alpha }_i^2+\frac{1}{2}\sum _{i=1}^n{\Vert {\mathbf{x}}_i-\mathbf{\mu }\Vert }_2^2\\ & =-\frac{1}{2}\sum _{i=1}^n{\left[{\mathbf{w}}^{\mathrm{T}}\left({\mathbf{x}}_i-\mathbf{\mu }\right)\right]}^2+\frac{1}{2}\sum _{i=1}^n{\Vert {\mathbf{x}}_i-\mathbf{\mu }\Vert }_2^2\\ & =-\frac{1}{2}\sum _{i=1}^n{\mathbf{w}}^{\mathrm{T}}\left({\mathbf{x}}_i-\mathbf{\mu }\right){\left({\mathbf{x}}_i-\mathbf{\mu }\right)}^{\mathrm{T}}\mathbf{w}+\frac{1}{2}\sum _{i=1}^n{\Vert {\mathbf{x}}_i-\mathbf{\mu }\Vert }_2^2\\ & =-\frac{1}{2}{\mathbf{w}}^{\mathrm{T}}\left(\sum _{i=1}^n\left({\mathbf{x}}_i-\mathbf{\mu }\right){\left({\mathbf{x}}_i-\mathbf{\mu }\right)}^{\mathrm{T}}\right)\mathbf{w}+\frac{1}{2}\sum _{i=1}^n{\Vert {\mathbf{x}}_i-\mathbf{\mu }\Vert }_2^2\\ & =-\frac{1}{2}{\mathbf{w}}^{\mathrm{T}}S\mathbf{w}+\frac{1}{2}\sum _{i=1}^n{\Vert {\mathbf{x}}_i-\mathbf{\mu }\Vert }_2^2\end{array}$$
上式的第二项与$w$无关，因此要极小化$f(w)$，只要使第一项极小化，于是优化问题转化为
$$\begin{array}{rl}\min & -\frac{1}{2}{w}^{T}Sw\\ s.t.& w^{T}w=1\end{array}$$
这个优化问题可以用拉格朗日乘子法求解，令拉格朗日函数为
$$L(w,\lambda )=-\frac{1}{2}{w}^{T}Sw+\frac{\lambda }{2}(w^{T}w-1)$$
令
$$\frac{\partial L}{\partial w}=-Sw+\lambda w=0$$
从而$$Sw=\lambda w$$
到这里，结果已经逐渐清晰了，我们要求的$w$正是矩阵$S$的特征向量。稍作变形：
$$w^{T}Sw=\lambda w^{T}w=\lambda$$
我们要最小化$-\frac{1}{2}{w}^{T}Sw$，就是要最大化$w^{T}Sw$，则$w$应该是$S$的最大特征值${\lambda }_{\max }$对应的特征向量。

算法总结

至此，我们可以总结一下二维向量压缩成一维的PCA的方法：
(1)求矩阵$$S=\sum _{i=1}^n(x_i-\mu )(x_i-\mu {)}^T=\tilde{X}{\tilde{X}}^T$$
(2)求$S$最大的特征值对应的特征向量，即$w$
(3)求${\alpha }_i=w^T(x_i-\mu )$

于是$X=(x_1,x_2,x_3,x_4,x_5)$经过压缩之后得到的结果就是$Y=({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5)$。

投影到方向向量$w$所对应的直线之后，$w$成了唯一的一个基，于是一维空间中的样本$x_i^{{}^{\prime }}$可以由基向量$w$表示：
$$x_i^{{}^{\prime }}=\mu +{\alpha }_{i}w$$
在原来的2维空间中，我们用基的系数来表示样本$x_i$，而在1维空间中，同样以基$w$的系数${\alpha }_i$来表示一维向量，它被称为主成分。

将二维向量压缩成一维向量${\alpha }_i$有时候只是为了减少传输时的数据量，一维向量是无法直接使用的，需要根据一维向量重构出原来的二维向量。如何重构呢？其实上面关于的$x_i^{{}^{\prime }}$的公式已经给出了答案。

二维样本PCA降维的例子

还是上面的例子，我们按照这个流程走一遍：
$$\mu =[2,3{]}^T$$
$$\tilde{X}=X-\mu =\left[\begin{array}{ccccc}-1& -1& 0& 2& 0\\ -2& 0& 0& 1& 1\end{array}\right]$$
$$S=\tilde{X}{\tilde{X}}^T=\left[\begin{array}{ccccc}-1& -1& 0& 2& 0\\ -2& 0& 0& 1& 1\end{array}\right]\left[\begin{array}{cc}-1& -2\\ -1& 0\\ 0& 0\\ 2& 1\\ 0& 1\end{array}\right]=\left[\begin{array}{cc}6& 4\\ 4& 6\end{array}\right]$$
S的最大特征值为10，对应特征向量$w=[\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}{]}^T$
于是降维后的表示为：
$$Y=({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5)=w^T\tilde{X}=[\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}]\left[\begin{array}{ccccc}-1& -1& 0& 2& 0\\ -2& 0& 0& 1& 1\end{array}\right]=[-3/\sqrt{2},-1/\sqrt{2},0,3/\sqrt{2},-1/\sqrt{2}]$$

要重构第一个样本$x_1{}^{\prime }$，方法是：
$$x_1^{{}^{\prime }}=\mu +{\alpha }_{1}w=[2,3{]}^T+\frac{-3}{\sqrt{2}}[\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}{]}^T=[\frac{1}{2},\frac{3}{2}{]}^T$$
当然，与原本的$x_1=[1,1{]}^T$还是有一些误差的。

更高维的情况

前面介绍的是二维降为一维的情况，更一般地，对于$m$维向量$x_i$如果要降维为$l$维的$y_i$，算法也是类似的，不加证明地给出以下步骤：
(1)求矩阵$$S=\sum _{i=1}^n(x_i-\mu )(x_i-\mu {)}^T=\tilde{X}{\tilde{X}}^T$$
(2)求$S$的所有特征值，从大到小排列，选取前$l$个特征值所对应的特征向量，即$W=(w_1,w_2,\cdots ,w_l)$。
(3)求各个样本点$x_i$对应于基$w_1,w_2,\cdots ,w_l$的系数，即主成分，${\alpha }_{i,k}=w_k^T(x_i-\mu ),k=1\cdots l$，得到低维的表示$y_i=({\alpha }_{i,1},{\alpha }_{i,2},\cdots {\alpha }_{i,l}{)}^T$
这个过程可以写成矩阵的形式：
$$Y=W^T\tilde{X}$$
$Y=(y_1,y_2,\cdots ,y_n)$是压缩后的样本点组成的矩阵。
(4) 原向量的重构：
$$x_i^{{}^{\prime }}=\mu +\sum _{k=1}^L{\alpha }_{i,k}{w}_k$$
写成矩阵的形式为：
$$X^{{}^{\prime }}=\mu +WY$$

PCA用于图像处理

原图：

PCA处理后：(选40个特征向量)