白话论文：A Tutorial on Principal Component Analysis
最近自己也在阅读一些论文，发现大部分文章都没有那么好懂，阅读起来总是需要很多思考和推理的时间。所以想把自己读过的文章都讲一讲，记录一下自己思考和推导的过程。希望有一天也能讲到自己写的文章- -
这次的文章是 A Tutorial on Principal Component Analysis（https://arxiv.org/abs/1404.1100），文章讲解了主成分分析（PCA）的动机与方法。如标题所说，这篇文章严格来说并不是一篇论文，更像一篇 PCA 的教程。
本文使用的符号
感觉这篇文章较为不好的一点是，作者使用的符号在各个小节会改来改去，矩阵的大小一会儿 $m \times n$ 一会儿 $n \times m$。这里规定一下这篇博客使用的符号。
$X$ 是一个 $m \times n$ 的矩阵，表示我们测量的数据。
测量一共进行了 $m$ 次，每次都测量了 $n$ 个变量，所以我们用 $x_i^T$ 表示 $X$ 里的第 $i$ 行，即第 $i$ 次测量的结果。
为什么需要 PCA
现实中的各种数据内容复杂，种类繁多，然而它们背后隐藏的规律大部分都是较为简单的。现实中的数据，除了真正有用的部分以外，还可能有较多的冗余和一定的噪声。如何将繁杂的数据化繁为简，发掘数据背后简单的规律，是所有科学研究人员都会遇到的问题。
文章举了一个简单的例子说明 PCA 的动机。假如我们要研究一个挂在弹簧上的小球的运动，根据物理定律我们知道，挂在单个弹簧上的小球在理想状态下应该做的是一维的直线运动，我们只要记录各个时刻小球沿着这条直线的 x 坐标，就能表达小球的运动。
然而，假如我们不知道这个物理规律，我们可能会“保险起见”（万一小球的运动由很多因素决定呢？），设置三个摄像机，分别记录小球在每个摄像机视角下的 x 坐标和 y 坐标。这样，原本只需要 1 个变量表示的现象，用了 6 个变量，这就造成了数据的冗余；另外，小球在现实中的运动并不是一个理想状态，各种因素都可能对小球的运动造成些许的影响，这就给数据带来了噪声。
PCA 能做的，就是帮我们对数据进行重新表示，重新表示后的数据不仅精简，还能在一定程度上减小噪声，为之后的数据分析带来方便。
PCA 作出的假设
PCA 作出的第一个假设就是线性性：重新表示后的数据，每个变量是重新表示前的变量的线性组合。这样，PCA 的目标就是要求出 $XP = Y$ 中的矩阵 $P$，使得 $X$ 在经过 $P$ 所蕴含的线性变换变成 $Y$ 之后，能够变得“更好”。
我们还可以从投影的角度解释这个线性变换。将 $P$ 的第 $i$ 列用 $p_i$ 表示，不难发现，数据重新表示后，第 $i$ 次测量的第 $j$ 个变量是由 $x_i^Tp_j$ 得到的，这个式子可以看作将 $x_i$ 投影到 $p_j$ 上，并把长度乘以 $|p_j|$（如果 $|p_j| = 1$ 那就真的是投影了）。所以，PCA 的目标也可以看成：在原来的数据空间中选出若干个“更好”的方向，将原来的数据投影到这些方向上。
这个较强的假设虽然方便了问题的处理，但仍然存在着一定的问题。如果我们要研究的是小球的圆周运动，显然仅使用一个变量 $\theta$ 来表示圆周运动的相位是最合适的。然而，$\theta$ 并不是小球的 x 坐标和 y 坐标的线性组合，这就是 PCA 没法完成的事情。当然啦，如果我们对要研究的问题事先有一定的了解，我们可以利用核函数把原来的变量进行一定的变换（例如把 xy 坐标变成极坐标），然后再使用 PCA 效果就会更好。
既然我们希望要选出“更好”的方向，那么什么样的方向才是“好”的呢？PCA 作出的第二个假设是：如果投影后的数据方差越大，那么这个方向越重要。
这个假设貌似来源于对噪声的一种度量方法——信噪比（signal-to-noise ratio，SNR），$\text{SNR} = \frac{\sigma_s^2}{\sigma_n^2}$，其中 $\sigma_s^2$ 是信号的方差，$\sigma_n^2$ 是噪声的方差。在一定量的噪声下，自然是信号的方差越大，数据越准确。例如，下面的示意图就画出了信号的方向和噪声的方向。可以看出，信号的方向上方差较大，而噪声的方向上方差较小。PCA 就是要把信号的方向选出来。
这个假设也可以从数据的冗余角度来解释。例如下图，最左边的数据冗余度较低（因为我们很难通过 $r_1$ 的值去推测 $r_2$ 的值），在较多的方向上数据的方差都较大；而最右边的数据冗余度较高（通过 $r_1$ 的值较容易推测出 $r_2$ 的值），数据只在“左上-右下”方向上显现出较大的方差。
为了减小数据的冗余，PCA 还希望选出来的方向之间相关度较低。说到相关度，果然容易想到协方差。假设我们的测量结果中，每一列的均值都是 0（如果不是 0，我们就求一下这一列的均值，然后把这一列都减去均值即可），显然有 $$\sum_{i=1}^m x_i = 0$$ 我们可以通过样本，算出 $n$ 个变量之间的协方差矩阵 $$C_X = \frac{1}{n}\sum_{i=1}^m x_i x_i^T - \frac{1}{n^2}(\sum_{i=1}^m x_i)(\sum_{j=1}^m x_j^T) = \frac{1}{n}X^TX$$ 由于投影后数据每一列均值是 0 的特性不会改变，容易发现投影后的协方差矩阵为 $$C_Y = \frac{1}{n}Y^TY$$ 既然我们希望不同方向之间的相关度较低，那么 $C_Y$ 最好情况下应该是一个对角矩阵，对角线上的元素表示的就是数据投影到各个方向之后的方差。
为了方便计算，PCA 还作出了第三个假设：选出的各个方向互相正交。也就是说，$P$ 是一个正交矩阵。这是一个非常强的假设，也让 PCA 在很多情况下无法使用（例如下图，很显然两个重要的方向并不正交）。
但是这个假设却极大地方便了我们的计算，让我们能够直接通过线性代数的方法推出 $P$ 的值。例如，我们现在可以把 $C_Y$ 表示成 $C_X$ 和 $P$ 的形式：$$C_Y = \frac{1}{n}Y^TY = \frac{1}{n}(XP)^T(XP) \\ = \frac{1}{n}(XP)^T(XP) = P^T(\frac{1}{n}X^TX)P = P^TC_XP$$
通过特征值分解（EVD）求解 PCA
接下来我们就要利用刚才的三个假设求出 $P$ 的值，也就是求出用来重新表示数据的方向。我们首先利用的方法是特征值分解（EVD）。
我们知道，$C_X = X^TX$ 是一个实对称矩阵，而实对称矩阵是一定能够正交对角化的，即 $X^TX = EDE^T$，其中 $E$ 是一个正交矩阵，每一列是 $X^TX$ 的一个特征向量，$D$ 是一个对角矩阵，对角线上的元素是与 $E$ 顺序对应的特征值，不妨令 $d_{1, 1} \ge d_{2, 2} \ge \dots \ge d_{n, n}$。把这个式子代入 $C_Y$ 得 $$C_Y = P^TC_XP = P^TEDE^TP$$ 此时，只要我们选择 $P = E$，由于它们是正交矩阵，我们有 $P^{-1} = E^{-1} = E$，那么 $$C_Y = D$$ 这就是我们想要的结果，通过选择一个合适的 $P$，把 $C_Y$ 变成一个对角矩阵。
概括地说，$X$ 的主成分就是 $X^TX$ 的特征向量，而且特征向量对应的特征值越大，这个主成分也就越重要（因为数据在这个方向上的投影方差越大）。
通过奇异值分解（SVD）求解 PCA
PCA 的求解还可以通过奇异值分解（SVD）的方式进行。这里先简单介绍一下 SVD。
SVD 是说，任何一个矩阵 $X$ 都能表示成 $XV = U\Sigma$ 的形式。其中 $V$ 是 $n \times n$ 的矩阵（记每一列为 $v_i$），$U$ 是 $m \times m$ 的矩阵（记每一列为 $u_i$），$\Sigma$ 是 $m \times n$ 的矩阵（记它的元素为 $\sigma_{i, j}$），只有 $i = j$ 时 $\sigma_{i, j}$ 可能非 0，记为 $\sigma_i$。
这些矩阵满足这样的条件：$v_i$ 是 $X^TX$ 的特征向量；$\sigma_i = \sqrt{\lambda_i}$，其中 $\lambda_i$ 是 $v_i$ 对应的特征值；$u_i = \frac{1}{\sigma_i}Xv_i$。另外要求 $V$ 是正交矩阵，那么 $$u_i^Tu_j = \frac{1}{\sigma_i\sigma_j}v_i^T(X^TXv_j) = \frac{\sigma_j}{\sigma_i}v_i^Tv_j = \begin{cases} 1 & i = j \\ 0 & i \ne j \end{cases}$$ 所以 $U$ 也是正交矩阵。
既然 $V$ 和 $U$ 都是正交矩阵，那么 $XV = U\Sigma$ 就可改为 $X = U\Sigma V^T$。把这个式子代入 $C_Y$ 有 $$C_Y = P^TC_XP \\ = P^TV\Sigma^TU^TU\Sigma V^TP = P^TVDV^TP$$ 其中 $d_{i, i} = \sigma_i^2 = \lambda_i$。如同上一节一样选择 $P = V$，我们就有 $$C_Y = D$$ 这也是我们想要的结果。$X$ 的主成分仍然是 $X^TX$ 的特征向量，而且特征向量对应的奇异值越大，这个主成分也就越重要。
