PCA (Principle Component Analysis)，主成分分析，是用的最广泛的数据预处理方法之一。同时PCA是一个almost always works的经典方法，基本上几乎所有的数据用PCA处理过以后再做后续的分类，往往会得到精度提升，不像很多其他方法只适用于某些数据而不适用于另一些数据。为什么PCA这么有效呢？我想主要是因为数据中普遍多多少少存在噪声，用PCA抽取主成分以后，噪声被滤除很多，使得分类变得更容易。另外，一些模型本身对数据空间有一些bias，用PCA做变换之后也会有影响，但这里的影响比较不确定。
说回PCA本身，PCA的想法是把数据从原空间投影到一个低维空间，使得这些低维空间上的投影（1）Variance最大，或者（2）Reconstruction error最小。把数据从高维向低维投影，这本身便是一个去噪的过程。因为维度减小，必有信息损失，我们希望去掉噪声的信息，保留有用的信息。这里提到的两个约束条件，variance其实代表了信息量，设想投影到某一个维度上variance很小，表示数据在这个维度上基本上没怎么变化，这就说明这个维度不是一个很有信息量的维度。而reconstruction error则更直观了，它就表示损失的信息量。降维投影和最优化约束二者综合起来，达到去掉无用信息，尽可能保留有用信息的目的。可以证明最大variance和最小reconstruction error这两者是等价的。
PCA也有可能失败，设想另一个极端，如果数据中噪声非常大，而有用信息非常微弱，则这里大variance的维度是噪声的维度，降维时反而是有用信息被去掉了。所幸的是，我们所遇到的绝大部分数据中噪声都是相对小的，要不然分析数据也就没有希望了。
这里介绍最大variance约束下PCA的推导。
一个向量(mathbf{x}inmathbb{R}^D)在一个方向(mathbf{w})上的投影是(left(mathbf{x}^ opfrac{mathbf{w}}{||mathbf{w}||} ight)frac{mathbf{w}}{||mathbf{w}||})，其中(frac{mathbf{w}}{||mathbf{w}||})是单位长方向向量，(mathbf{x}^ opfrac{mathbf{w}}{||mathbf{w}||})是(mathbf{x})的投影在这个方向上的长度。若(mathbf{w})为单位向量，投影可简化为((mathbf{x}^ opmathbf{w})mathbf{w}=(mathbf{w}mathbf{w}^ op)mathbf{x}). 进一步，如果投影到一个低维空间，而非一条线（一维空间），设该空间的一组单位正交基为(mathbf{W}=(mathbf{w}_1, ..., mathbf{w}_K)^ op)，其中(forall k, mathbf{w}_k^ opmathbf{w}_k=1)，而(forall k eq j, mathbf{w}_k^ opmathbf{w}_j=0)。则(mathbf{x})在该空间的投影(hat{mathbf{x}}=sum_k alpha_k mathbf{w}_k)满足(mathbf{x}=hat{mathbf{x}}+mathbf{x}_perp)，(mathbf{x}_perp)为垂直分量，因此对任意k，有
$$mathbf{w}_k^ opmathbf{x}=sum_j alpha_j mathbf{w}_k^ opmathbf{w}_j + mathbf{w}_k^ opmathbf{x}_perp=alpha_k$$
从而
$$hat{mathbf{x}} = sum_k (mathbf{w}_k^ opmathbf{x})mathbf{w}_k = left(sum_k mathbf{w}_kmathbf{w}_k^ op ight)mathbf{x}=mathbf{M}mathbf{x}$$
其中(mathbf{M}=sum_k mathbf{w}_kmathbf{w}_k^ op).
设给定数据集(mathbf{X}=(mathbf{x}_1, ..., mathbf{x}_n)^ op)，将每一个数据点投影到K维空间，得到的投影数据集(hat{mathbf{X}}=(hat{mathbf{x}}_1, ..., hat{mathbf{x}}_n)^ op)的variance为
$$mathrm{Var}(hat{mathbf{X}})=frac{1}{n}sum_i (hat{mathbf{x}}_i - ar{hat{mathbf{x}}})^ op(hat{mathbf{x}}_i - ar{hat{mathbf{x}}})$$
其中(ar{hat{mathbf{x}}})为投影的均值，
$$ar{hat{mathbf{x}}}=frac{1}{n}sum_i hat{mathbf{x}_i} = mathbf{M}left(frac{1}{n}sum_imathbf{x}_i ight)=mathbf{M}ar{mathbf{x}}$$
(ar{mathbf{x}})为原数据的均值。从而
$$mathrm{Var}(hat{mathbf{X}})=frac{1}{n}sum_i (mathbf{x}_i - ar{mathbf{x}})^ opmathbf{M}^ opmathbf{M}(mathbf{x}_i - ar{mathbf{x}})$$
利用单位正交基的性质，矩阵(mathbf{M}^ opmathbf{M})可以化为
$$mathbf{M}^ opmathbf{M}=left(sum_k mathbf{w}_kmathbf{w}_k^ op ight)^ opleft(sum_j mathbf{w}_jmathbf{w}_j^ op ight)=sum_ksum_j mathbf{w_k}left(mathbf{w}_k^ opmathbf{w}_j ight)mathbf{w}_j^ op=sum_k mathbf{w}_kmathbf{w}_k^ op$$
这其中用到了(mathbf{w}_k^ opmathbf{w}_j)仅在(k=j)时为1而其他情况为0的性质。从而(mathrm{Var}(hat{mathbf{X}}))可进一步化为
$$egin{align*}mathrm{Var}(hat{mathbf{X}})&=frac{1}{n}sum_i (mathbf{x}_i - ar{mathbf{x}})^ opleft(sum_k mathbf{w}_kmathbf{w}_k^ op ight)(mathbf{x}_i - ar{mathbf{x}})\&=frac{1}{n}sum_i sum_k mathbf{w}_k^ op(mathbf{x}_i - ar{mathbf{x}})(mathbf{x}_i - ar{mathbf{x}})^ opmathbf{w}_k\&=sum_k mathbf{w}_k^ opleft(frac{1}{n}sum_i (mathbf{x}_i - ar{mathbf{x}})(mathbf{x}_i - ar{mathbf{x}})^ op ight)mathbf{w}_k\&= sum_k mathbf{w}_k mathrm{Cov}(mathbf{x})mathbf{w}_kend{align*}$$
记(mathbf{V}=mathrm{Cov}(mathbf{x}))为原数据集的协方差矩阵，我们可以得到PCA的优化问题为
$$egin{align*}max_mathbf{W} & sum_k mathbf{w}_k^ opmathbf{V}mathbf{w}_k \ mathrm{s.t.} & forall k, mathbf{w}_k^ opmathbf{w}_k=1 \ & forall k eq j, mathbf{w}_k^ opmathbf{w}_j=0end{align*}$$
这个优化问题可以用拉格朗日乘子方法来解，考虑第一组约束(forall k, mathbf{w}_k^ opmathbf{w}_k=1)，引入拉格朗日乘子后目标函数变为
$$L=sum_k mathbf{w}_k^ opmathbf{V}mathbf{w}_k - sum_k lambda_k left(mathbf{w}_k^ opmathbf{w}_k - 1 ight)$$
其取最优解时满足
$$egin{align*}frac{partial L}{partial mathbf{w}_k} &=2mathbf{V}mathbf{w}_k - 2lambda_k mathbf{w}=0 \ frac{partial L}{partial lambda_k}&=mathbf{w}_k^ opmathbf{w}_k-1=0 end{align*}$$
其中第一式即
$$mathbf{V}mathbf{w}_k=lambda_k mathbf{w}_k$$
这正是特征值（eigen value）和特征向量（eigen vector）的定义，说明取最优解时，所有(mathbf{w}_k)均为(mathbf{V})的特征向量，而(lambda_k)为对应的特征值。上式两边同时左乘(mathbf{w}_k^ op)，得到
$$mathbf{w}_k^ opmathbf{V}mathbf{w}_k = lambda_k mathbf{w}_k^ opmathbf{w}_k$$
从而目标函数变为
$$L=sum_k lambda_k mathbf{w}_k^ opmathbf{w}_k - sum_k lambda_k left(mathbf{w}_k^ opmathbf{w}_k - 1 ight) = sum_k lambda_k$$
要使该目标函数最大，则需使每个(lambda_k)都尽量大。又考虑到约束条件(forall k eq j, mathbf{w}_k^ opmathbf{w}_j=0)，说明最优解时各特征向量两两垂直，从而对应的特征值需两两不同（若某个特征值为m重特征值则至多重复m次），要使(sum_k lambda_k)最大，只需取(mathbf{V})最大的K个特征值以及对应的单位特征向量即可。
因此，PCA优化问题的解为原数据协方差矩阵的最大K个特征值对应的特征向量。每个特征值表示了对应那个维度的variance大小，这一点可以从(mathrm{Var}(hat{mathbf{X}})=sum_k mathbf{w}_k^ opmathbf{V}mathbf{w}_k=sum_k lambda_k)中很容易看出来，选取最大的K个特征值，也就是选取了最大variance的K个维度。
优化问题求解完毕后，PCA可以用来降维。将K个单位特征向量记为矩阵形式(mathbf{W}=(mathbf{w}_1, ..., mathbf{w}_K)^ opinmathbb{R}^{K imes D})，任意原空间的数据向量(mathbf{x})的投影为(hat{mathbf{x}}=mathbf{W}mathbf{x}inmathbb{R}^{K})，在原空间的表示为(mathbf{x}'=mathbf{W}^ opmathbf{W}mathbf{x}inmathbb{R}^D)，此即为其reconstruction。写成矩阵形式后，用程序实现将非常简洁方便。
由于PCA学到的每一个低维空间基向量(mathbf{w}_kinmathbb{R}^D)在原空间中，从而也可以可视化。每一个基向量都是一个有代表性的模式组成部分，应用于图像数据时，大特征值对应的特征向量往往是低频平滑的图像，而特征值越小特征向量对应的越是高频变化快的图像。用于人脸图像时，这些基向量还有一个名字叫eigen faces。注意这些模式和用kmeans或者mixture of Gaussians学出来的模式有本质不同，在mixture model中，每一个模式就是一个有代表性的数据点，而这里PCA的所有component线性组合起来才构成一个数据模式。
最后说一说实现。PCA用任何代数工具包都很容易实现，以MATLAB为例，对任何数据矩阵(mathbf{X})（每行一个数据），首先求协方差矩阵((mathbf{X}-ar{mathbf{X}})^ op(mathbf{X}-ar{mathbf{X}}))，再用eigs函数求最大的K个特征值及对应特征向量即可。此处计算协方差矩阵无需除以n，因为这个因子的作用只是减小特征值，并不影响特征向量。代码非常简单，三行：
xmean = mean(X, 1); X = X - ones(size(X, 1), 1) * xmean; [W, ~] = eigs(double(X'*X), [], K);
Matlab算出的W为DxK矩阵，是上面所说W的转置。使用时要降维，只需做矩阵乘法XW。用SVD分解也可以求得W，但经过试验发现用eigs直接求会更快且更省内存一些。另外减掉均值后再做投影应该往往更好一些。