机器学习之主成分分析PCA
PCA(主成分分析)是一种常见的数据降维算法，其核心思想是找到一个维数更低的仿射集，然后将所有的数据点在其上做投影，以得到新的在更低维空间中的点作为新的数据。那么问题来了，如何选择这样的仿射集，以降维的同时不至于损失过多的信息呢？一般来说有两种思路：
- 最近重构性：样本点到该仿射集的距离要尽量小；
- 最大可分性：样本点到该放射集的投影要尽可能分开。
下面我们构建数学模型，并且说明，两种思路其实是一回事。这里的数学推导是本人自己作为练习详细给出的，因为大多数机器学习的书基本上是含糊其辞地混过去了一些细节，让人怀疑作者本身有没有搞明白。
1.数学推导：
首先强调一下，所有推导中出现的向量若没有指出其为行向量则均是列向量。现在的待求解问题是：
问题：已知数据集 $D=\lbrace x_{i}\in\mathbb{R}^{M}\rbrace_{i=1}^{N}$, 现在我们需要找到一个$m$维仿射集($m<M$)，使得$D$中的点到该仿射集距离之平方和最小。
我们知道，任何$m$维仿射集$H$可以表示为如下形式：
\begin{equation}H=\lbrace{Wu+b}\mid u\in\mathbb{R}^{m}\rbrace\end{equation}
其中矩阵$W\in\mathcal{M}_{M\times m}$满足：$W^{T}W=Id$, $b\in\mathbb{M}$. 任何仿射集只不过是线性空间平移得到的，就如（1）中的仿射集合H是由线性空间$L_{W}\triangleq\lbrace Wu\mid u\in\mathbb{R}^{m}\rbrace$平移得到，这时如果我们将矩阵W写成一组列向量的形式：$W=[w_{1},...,w_{m}],$ 则我们知道$w_{i}$,$i=1,...,m$是$L_{W}$的一组标准正交基。
一点$x\in\mathbb{R}^{M}$到超平面$H$的投影：
\begin{equation}P_{H}(x)=b+WW^{T}(x-b)\end{equation}
其实我们不妨设投影为$P_{H}(x)=b+Wu$, 则必然有：$x-P_{H}(x)\perp L_{W}$, 也就是：$W^{T}(x-P_{H}(x))=0,$于是可以得到$u=W^{T}(x-b)$.这时候很容易算出到$H$之距离之平方：
\begin{equation}d^{2}_{H}(x)=\Vert x-P_{H}(x)\Vert^{2}=\Vert (\text{Id}-WW^{T})(x-b)\Vert^{2}=(x-b)^{T}(\text{Id}-WW^{T})(x-b)\end{equation}
所有这个时候所有的点到达超平面距离之和的平方为：
\begin{equation}S(W,b)=\sum_{i=1}^{N}(x_{i}-b)^{T}(\text{Id}-WW^{T})(x_{i}-b)\end{equation}
于是我们只需要求优化问题：
\begin{split}\text{min }& \ S(W,b)\newline \text{s.t}& \ W^{T}W=\text{Id} \end{split}
的某个特定解就够了。注意到上述优化问题的极值一定存在，因为约束条件构成一个欧氏空间中的紧集，并且注意约束条件其实是有多达$m\times m$个，即：$w_{i}^{T}w_{j}=\delta_{ij}$, $i,j=1,...,m$. 所以我们这个时候利用Lagrange乘子法可知道极小值一定满足：
存在常数$a_{ij}$使得:
\begin{equation}\text{d}S=\sum_{i,j}a_{ij}\text{d}(w_{i}^{T}w_{j}).\end{equation}
首先右边完全没有关于db的项，所以左边必然也没有任何\text{d}b有关的项，也就是：
\begin{split}\text{d}_{b}S=2(Nb-\sum_{i=1}^{N}x_{i})^{T}(\text{Id}-WW^{T})db=0,\end{split}
我们很自然得到条件:
\begin{equation}(Nb-\sum_{i=1}^{N}x_{i})^{T}(\text{Id}-WW^{T})=0,\end{equation}
当然也等价于: $$(\text{Id}-WW^{T})(Nb-\sum_{i=1}^{N}x_{i})=0$$
注意到$(\text{Id}-WW^{T})y$其实就是y的垂直于$L_{W}$之分量，对于任意的$y\in\mathbb{R}^{M}$, 所以$(\text{Id}-WW^{T})y=0$的解集自然是线性空间$L_{W}$自身，于是(6)的解集为：
\begin{equation}{b=\frac{1}{N}(\sum_{i=1}^{N}x_{i})+Wu, u\in\mathbb{R}^{m}}\end{equation}
这时为简单起见，我们不妨现在直接取：
\begin{equation}b=\frac{1}{N}\sum_{i=1}^{N}x_{i}\end{equation}
由以上条件可以知道，最优解对应的超平面必然过数据集的重心。下面我们继续计算$d_{W}S$.容易计算：
$$S(w,b)=\sum_{i=1}^{N}\Vert x_{i}-b\Vert^{2}-\sum_{i=1}^{N}(x_{i}-b)^{T}WW^{T}(x_{i}-b)$$
这个时候我们寻求将上式表示为矩阵形式。首先可以定义一个矩阵$X$,其第$i$行为行向量$(x_{i}-b)^{T}$。则我们很容易知道：
\begin{equation}\begin{split}S(w,b)=&\sum_{i=1}^{N}\Vert x_{i}-b\Vert^{2}-\sum_{i=1}^{N}tr[XWW^{T}X^{T}]\newline =&\sum_{i=1}^{N}\Vert x_{i}-b\Vert^{2}-\sum_{i=1}^{N}tr[W^{T}X^{T}XW]\end{split}\end{equation}
于是由上式很容易求出：
$$\text{d}_{W}S=-2tr(W^{T}X^{T}X\cdot \text{d}W).$$
而注意到（5）的右边：
\begin{equation}\sum_{i,j}a_{ij}\text{d}(w_{i}^{T}w_{j})=tr(A^{T}\cdot \text{d}(W^{T}W)),\end{equation}
其中我们定义矩阵：$A\triangleq (a_{ij})$, 这时容易计算得到：
\begin{equation}\sum_{i,j}a_{ij}\text{d}(w_{i}^{T}w_{j})=tr((A+A^{T})W^{T}\cdot dW),\end{equation}
所以我们立即得到：
\begin{equation}W^{T}X^{T}X=-\frac{1}{2}(A+A^{T})W^{T},\end{equation}
我们令$B=-\frac{1}{2}(A+A^{T})$, 则$B$是对称$m\times m$矩阵，其必然可以对角化，也就是存在$m\times m$正交阵$\Sigma$使得$B=\Sigma\text{diag}(\lambda_{1},...,\lambda_{m})\Sigma^{T}$, 这时上式等价于：
\begin{equation}W^{T}X^{T}X=\Sigma\text{diag}(\lambda_{1},...,\lambda_{m})\Sigma^{T} W^{T},\end{equation}
两边取转置$T$立马得到：
\begin{equation}X^{T}X(W\Sigma)=(W\Sigma) \text{diag}(\lambda_{1},...,\lambda_{m})\end{equation}
这时不妨设$W\Sigma=[v_{1},...,v_{m}]$, 由上式我们知道所有的列向量$v_{i}$只不过是$X^{T}X$的特征向量，且$X^{T}Xv_{i}=\lambda_{i}v_{i}$, 这时求:
$$S(w,b)=\sum_{i=1}^{N}\Vert x_{i}-b\Vert^{2}-\sum_{i=1}^{m}\lambda_{i},$$
为使其极小我们需要使得$\lambda_{1},...,\lambda_{m}$正好是$X^{T}X$最大的$m$个特征值，再求得相应的模长为1的特征向量$v_{1},...,v_{m}$使其相互正交，然后我们不妨令(14)中$\Sigma=\text{Id}$, 则$(W,b)=([v_{1},...,v_{m}],\frac{1}{N}\sum_{i=1}^{N}x_{i})$就是最优化问题的一个解。
算法
输入：数据集$D=\lbrace x_{i}\in\mathbb{R}^{M}\rbrace_{i=1}^{N}$，正整数m<M
输出：$M\times m$矩阵$W$
Step1：去中心化：$x_{i}=x_{i}-\frac{1}{N}\sum_{i=1}^{N}x_{i}$；
Step2：令$X=(x_{1},...,x_{N})^{T}$, 计算矩阵$X^{T}X$；
Step3：计算$X^{T}X$的特征值，并选择出$m$个最大的特征值$\lambda_{1},...,\lambda_{m}$, 得到相应的特征向量$w_{1},...,w_{m}$, 使其相互正交，模长为1，
输出$W=[w_{1},...,w_{m}]$
分析：
现在由上述PCA算法，我们求出了$W,b$自然也就已经得到了仿射集$H=\lbrace Wu+b\mid u\in\mathbb{R}^{m}\rbrace$，注意到这时候$x_{i}$到$H$的投影：
$$x_{i}^{\prime}=b+WW^{T}(x_{i}-b)$$
结合$b=\frac{1}{N}\sum_{i=1}^{N}x_{i}$很容易得到：
$$\frac{1}{N}\sum_{i=1}^{N}x^{\prime}_{i}=b$$
并且：
$$\Vert x^{\prime}_{i}-b\Vert^{N}=(x_{i}-b)^{T}WW^{T}(x_{i}-b)$$
结合第一节中我们已经得到的公式$S(w,b)=\sum_{i=1}^{N}\Vert x_{i}-b\Vert^{2}-\sum_{i=1}^{N}(x_{i}-b)^{T}WW^{T}(x_{i}-b)$我们可以得到：
\begin{equation}\sum_{i=1}^{N}\Vert x_{i}-b\Vert^{2}=\sum_{i=1}^{N}\Vert x^{\prime}_{i}-b\Vert^{2}+S(W,b)\end{equation}
该公式中第一项表示的是数据集点的分开程度，而右边第一项是到$H$的投影点的分开程度，第二项正好是距离平方和，由于左边是常数，距离平方和最小的时候，自然投影点之间越分开，于是我们论证了开头所说的最近重构性和最大可分性的等价性。
