Book Notes on Pattern Recognition and Machine Learning (PRML)
PCA 主成分分析，被广泛应用于数据降维，特征提取等领域。我们有两种广泛使用的 PCA 定义，从它们可以引出相同的 PCA 算法。PCA 的第一种定义将数据正交投影在低维线性空间，使得投影数据的 variance 最大化的投影方式。PCA 的第二种定义是最小化原数据和投影数据间均方距离的线性投影方式。接下来我们分别来讨论这两种定义
我们假设观测数据集为 \(\{x_1,...,x_n\}\) ，其中每个观测数据 \(x_i\) 的维度为 D。我们的目标是将数据 \(x_i\) 降维到 M 维，使得投影后的数据方差最大。在这里我们先假设降维维度 M 是已知的，考虑最简单的降为一维的情况 (M=1)
我们用一个 D 维向量 \(u\) 来表示原空间的方向。因为我们只关心其方向，我们令向量 \(u\) 为单位向量 \((u^Tu=1)\) 。每个数据向量都会被投影为一个标量，即 \(u^Tx_i\) 。我们写出投影数据的方差
\[\frac{1}{n}\sum_{i=1}^n \{u^Tx_i - u^T\overline{x}\}^2, \quad where \ \ \overline{x} = \frac{1}{n}\sum_{i=1}^n x_i\]我们将协方差矩阵提取出来，上式可以写成矩阵向量的形式
\[u^TSu, \qquad where\ \ S = \frac{1}{n}\sum_{i=1}^n (x_i-\overline{x})(x_i-\overline{x})^T\]根据 PCA 定义，我们想要使得投影后的数据方差最大，此时我们的目标问题即为
\[max_{u}\ u^TSu, \quad s.t.\ u^Tu=1\]通过引入拉格朗日乘子 \(\lambda\) 来消除约束，解决这个问题
\[L = u^TSu + \lambda (1-u^Tu)\]我们对 \(u\) 求导并令其为零可以发现
\[\frac{\partial L}{\partial u} = 0 \Rightarrow Su=\lambda u\]因为我们引入的 \(\lambda\) 是一个常数，这意味着 \(u\) 是协方差矩阵 S 的一个 eigenvector，我们在等式两边同时乘 \(u^T\)，结合单位向量的假设，可以发现
\[u^TSu = \lambda\]为使目标方差 \(u^TSu\) 最大，我们需要 \(\lambda\) 取最大值。而 \(\lambda\) 是协方差矩阵 S 的一个 eigenvalue，这意味着我们需要选取协方差矩阵 S 中最大的 eigenvalue 对应的 eigenvector 作为 \(u\) ，此时这个 eigenvector 也被称为第一主成分。
至此我们就得到了 PCA 的建立及求解方式。对应非一维的数据降维，我们对应的选取前 M 大 eigenvalue 对应的 eigenvector 组合起来作为投影矩阵 \(u\)
该过程中主要的计算量消耗在 eigenvector 的计算上，对于尺寸为 \(D\times D\) 的矩阵，计算的时间复杂度为 \(O(D^3)\) 。因为我们只关心前 M 大的 eigenvector，我们可以用 power method 简化，将复杂度降为 \(O(MD^2)\)
Minimum error formulation 的核心思想是通过旋转坐标系以最小化投影误差。这里同样假设观测数据集为 \(\{x_1,...,x_n\}\) ，其中每个观测数据 \(x_i\) 的维度为 D 。我们首先引入一组 D 维的完全正交基向量 \(\{u_1,...,u_D\}\) ，我们可以根据这组基向量来表示出每一个观测数据，即
\[x_n = \sum_{i=1}^D \alpha_{ni}u_i\]我们可以认为原来的观测数据 \(x_n = \{x_{n1},...,x_{nD}\}\) 经过旋转得到新的坐标，即这里的 \(\{\alpha_{n1},...,\alpha_{nD}\}\) ，这意味着 \(\alpha_{ni} = x_n^Tu_i\) ，我们代入上式可得
\[x_n = \sum_{i=1}^D (x_n^Tu_i)u_i\]因为我们需要将数据降到 M 维，我们用前 M 个基向量来表示这个 M 维的线性子空间。此时观测数据的估计值可以写成
\[\tilde{x}_n = \sum_{i=1}^M z_{ni}u_i + \sum_{i=M+1}^D b_iu_i\]第一项就是我们希望得到的维数为 M 的子空间。第二项数据中包含较多冗余信息的部分，对于这些冗余的维度，我们对每个维度用一个常数 \(b_i\) 来约束。根据 PCA 第二种定义，我们要最小化原数据和投影数据间均方距离，即
\[min_{u,z,b}\ J = \frac{1}{n}\sum_{i=1}^n ||x_i-\tilde{x}_i||^2 = \frac{1}{n}\sum_{i=1}^n ||x_i-\sum_{j=1}^M z_{ij}u_j - \sum_{j=M+1}^D b_ju_j||^2\]我们首先对 \(z_{ij}\) 求导并令其为零，结合基向量的正交性，可以得到
\[z_{ij} = x_i^Tu_j,\quad where\ \ j=1,...,M\]接着我们对 \(b_{j}\) 求导并令其为零，结合基向量的正交性，可以得到
\[b_j = (\frac{1}{n}\sum_{i=1}^n x_i)^Tu_j = \overline{x}^Tu_j,\quad where\ \ j=M+1,...,D\]现在我们将 \(z_{ij}\) 和 \(b_{j}\) 的形式代入优化目标中，可以得到
\[J = \frac{1}{n}\sum_{i=1}^n\sum_{j=M+1}^{D}\{(x_i-\overline{x})^Tu_j\}u_j = \frac{1}{n}\sum_{i=1}^n\sum_{j=M+1}^{D}(x_i^Tu_j-\overline{x}^Tu_j)^2 = \sum_{i=M+1}^D u_i^TSu_i\]这里的 S 同样是协方差矩阵，我们得到的形式等价于第一种定义。我们可以得到映射后的数据 \(\hat{x}\) 根据公式
\[\hat{x} = \sum_{i=1}^M (x^Tu_i)u_i\]在高维数据中，有时维数会高于样本数，即 \(D>N\) 。在这种情况下，如果我们用之前提到的方法计算 eigenvector，时间复杂度 \(O(D^3)\) 是难以接受的。因此需要进一步的修改优化。
我们用 X 来表示输入数据集，是一个 \((N\times D)\) 的矩阵，第 n 行表示的是第 n 个归一化的数据 \((x_n-\overline{x})^T\) ，此时的协方差矩阵可以写出 \(S=N^{-1}X^TX\) ，此时的 eigenvector 方程为
\[\frac{1}{N}X^TXu = \lambda u\]我们在等式两边同时乘 X ，可以得到
\[\frac{1}{N}XX^T(Xu) = \lambda (Xu) \Rightarrow \frac{1}{N}XX^Tv = \lambda v\]通过这种方式，我们将问题转化为求 \(XX^T\) 矩阵的 eigenvector，时间复杂度从原问题的 \(O(D^3)\) 变为 \(O(N^3)\)。在求出 eigenvector v 之后，我们需要还原到待求的原 eigenvector u，通过在等式两边再同时乘 \(X^T\)
\[\frac{1}{N}X^TX(X^Tv) = \lambda (X^Tv)\]可以发现我们又回到了原问题的协方差矩阵，此时 \(X^Tv\) 是对应的 eigenvector。最后我们进行一下归一化，限制 \(\vert\vert u \vert\vert=1\) ，得到最后的投影矩阵 u
\[u_i = \frac{X^Tv_i}{\sqrt{N\lambda_i}}\]总结来说，在高维数据的 PCA中，我们首先计算 \(XX^T\) 的 eigenvector \(\{v_i\}\) ，再用上式去还原出原协方差矩阵的 eigenvector \(\{u_i\}\) 作为最终结果
结合之前学过的 Kernel 方法，我们可以通过 kernel 的方式来引入非线性变换。在上面的标准 PCA 中，我们对原数据进行矩阵相乘，将其从 D 维空间线性变换到 M 维空间，即 \(x \rightarrow u^Tx\)。现在我们考虑用一个变换函数 \(\phi(x)\) 实现映射变换，即 \(x \rightarrow \phi(x)\)。对于某个数据 \(x_i\) 而言，它是一个 D 维向量。而变换后的 \(\phi(x_i)\) 则变成一个 M 维向量。
我们假设映射后的数据同样均值为零，即 \(\sum_i \phi(x_i) = 0\) ，那么我们可以写出 MxM 特征空间的协方差矩阵
\[C = \frac{1}{n}\sum_{i=1}^n\phi(x_i)\phi(x_i)^T\]我们用 \(\{v_1,...,v_M\}\) 来表示协方差矩阵 C 的特征向量。这些特征向量也可以写为特征空间数据的线性组合，即
\[v_i = \sum_{j=1}^n \alpha_{ij}\phi(x_j)\]代入到矩阵 C 的特征值方程中，对于第 k 个特征值 \(\lambda_k\) 可以得到，
\[(\frac{1}{n}\sum_{i=1}^n\phi(x_i)\phi(x_i)^T) (\sum_{j=1}^n \alpha_{kj}\phi(x_j)) = \lambda_k \sum_{i=1}^n \alpha_{ki}\phi(x_i)\]为了凑出 kernel 的形式，我们在等式两边同时乘 \(\phi(x_s)^T\) ，此时上式就可以写成矩阵的形式
\[K^2\alpha_k = \lambda_knK\alpha_k\]这里的 K 表示 kernel 矩阵，\(K_{ij} = K(x_i,x_j) = \phi(x_i)^T\phi(x_j)\) 。\(\alpha_k\) 则表示一个 N 维列向量。我们可以去掉上式等号两边的 K，此时得到的解与原式的不同在于零特征值对应的特征向量，而这种特征向量的差别对主成分的分解不会产生影响。
那么现在的问题就和标准 PCA 一样了，我们希望得到矩阵 K 的特征向量，即
\[K \alpha_k = \lambda_kn\alpha_k\]此外，我们还需要保证特征向量 \(\alpha_i\) 的归一化，这一点可以通过归一化特征空间特征向量来满足，即
\[v_i^Tv_i = \alpha_i^TK\alpha_i = \lambda_i n \alpha_i^T\alpha_i = 1\]在确定系数 \(\alpha_i\) 后，就可以得到原问题的特征向量 \(v_i\) ，进而实现变换。原数据 x 映射在特征空间第 m 维的值同样可以表示为 kernel 的形式
\[y_m(x) = \phi(x)^Tv_m = \sum_{i=1}^n \alpha_{mi}\phi(x)^T\phi(x_i) = \sum_{i=1}^n\alpha_{mi}k(x,x_i)\]至此我们就完成了 kernel PCA 的建模。最开始我们假设映射后的数据同样均值为零，即 \(\sum_i \phi(x_i) = 0\) 。但这个假设在实际情况中通常并不成立，所以我们需要人为地进行处理
\[\tilde{\phi}(x_k) = \phi(x_k) - \frac{1}{n}\sum_{i=1}^n\phi(x_i)\]此时的 kernel 矩阵也会发生变化
\[\tilde{K}_{ij} = \tilde{\phi}(x_i)^T\tilde{\phi}(x_j) \Rightarrow \tilde{K} = K-1_nK-K1_n+1_nk1_n\]这里的 \(1_n\) 为 nxn 的矩阵，其中每个元素都是 1/n 。这种变化对我们的建模并没有影响。特别地，如果取线性 kernel \(k(x_i,x_j)=x_i^Tx_j\) ，那么就可以得到标准的 PCA 模型。
最后是 kernel PCA 的缺点。第一个缺点是它需要找到 NxN 尺寸的 kernel 矩阵 K 的特征向量，而在传统 PCA 中，我们只需要找的是 DxD 尺寸的协方差矩阵。当数据量大的时候就需要用到近似。第二个缺点是它必须保留所有的特征向量用于变换，而在传统的 PCA 中，我们只需要保留前 M 大的特征向量即可。
PCA 也可以表示为概率隐变量模型的极大似然解。传统的 PCA 可以看成是将 D 维的数据投影到 M 维的线性子空间。而 probabilistic PCA 则可以看成一种从隐空间到数据空间的映射。probabilistic PCA 是线性高斯框架下的应用，这意味着所有的边缘分布，条件分布都是高斯分布。我们从生成模型的角度来引入 probabilistic PCA
我们假设观测数据 x 为 D 维向量，由隐变量 z 控制，z 是 M 维向量。我们首先假设隐变量 z 先验分布为标准高斯分布，即 \(z \sim N(0,I)\) ，之后我们可以证明，z 的高斯分布参数对建模没有影响。
同样，条件概率分布也是高斯分布，我们假设 \(x\vert z \sim N(Wz+u,\sigma^2I)\)
这里 x 可以看成是 z 的线性变换，由矩阵 W (\(D \times M\)) 控制。也就是说，观测数据可以表示为
这里 \(\epsilon\) 表示噪声，服从均值为零，方差为 \(\sigma^2I\) 的高斯分布，对应条件分布中的方差项。
在表示出隐变量 z 的先验分布和条件分布后，我们可以写出观测数据 x 的分布，即
\[p(x) = \int p(x\vert z)p(z)dz \Rightarrow x \sim N(u,WW^T+\sigma^2I)\]观测数据的分布由参数 \(u,W,\sigma\) 控制。通过观察我们可以发现一个现象，如果我们将 \(W\) 替换为 \(W_0 = WR\) ，这里 R 是一个正交矩阵。此时出现在方差项的 \(W_0W_0^T\) 可以写成
\[W_0W_0^T = WRR^TW = W^TW\]我们对 W 的变换并不会影响观测数据 x 的分布，也就是说有一系列 W 都会产生相同的预测分布。这里的矩阵 R 可以看成是隐变量空间的旋转，R 的存在表明隐变量空间具有旋转不变性的特征。
回到之前对隐变量 z 先验分布的假设。如果我们不约束其高斯分布的参数，记 \(z \sim N(m,\Sigma)\) ，此时可以写出 x 的分布 \(x\sim N(Wm+u,\sigma^2I+W\Sigma W^T)\)。我们令 \(\hat{u} = Wm+u, \hat{W}=W\Sigma^{1/2}\) 可以得到和上面标准正态分布假设一样的结果。说明 z 的具体分布形式对建模没有影响。
接下来我们考虑用极大似然来确定模型参数，
\[\ln p(X|u,W,\sigma) = -\frac{ND}{2}\ln (2\pi)-\frac{N}{2}\ln |C|-\frac{1}{2}\sum_{i=1}^N(x_i-u)^TC^{-1}(x_i-u)\]这里的 C 就是观测数据 x 分布的方差，\(C=WW^T+\sigma^2I\)
我们直接写出通过求导确定的参数解析解
\[u = \frac{1}{N}\sum_{i=1}^Nx_i\] \[W = U_M(L_M-\sigma^2I)^{1/2}R\]这里 \(U_M\) 是一个 DxM 的矩阵，每一列都是协方差矩阵 S 的一个特征向量。\(L_M\) 是一个 MxM 的对角矩阵，每个对角值是 \(U_M\) 中特征向量对应的特征值。R 是一个任意的 MxM 正交矩阵。可以证明，只有当 \(U_M\) 中的 M 列是 S 的特征值前 M 大对应的特征向量时，可以取得似然方程的最大值（上式的其他解都对应的是鞍点）
\[\sigma^2 = \frac{1}{D-M}\sum_{i=M+1}^D \lambda_i\]接下来我们对模型进行分析。考虑观测数据 x 的协方差矩阵 C，假设预测分布在某方向上的方差由单位向量 v 决定，即 \(v^TCv\) 。如果我们首先假设 v 垂直与主成分空间，是剩余特征向量的线性组合，此时我们根据 \(v^TU=0\) 可以得到 \(v^TCv = \sigma^2\) 。如果我们再假设 v 就是某个被保留在 U 中的特征向量，此时我们根据 \(v=u_i\) 可以得到 \(v^TCv = \lambda_i\) 。这说明我们的模型完整地保留了数据在主成分方向上的方差，而用 \(\sigma^2\) 去近似估计数据在其余方向的方差。
在最开始我们提到 probabilistic PCA 则可以看成一种从隐空间 \(\{z\}\) 到数据空间 \(\{x\}\) 的映射，那么我们也可以利用贝叶斯定理来反向映射。
和之前的假设一致，我们有
\[z \sim N(0,I),\ x\vert z \sim N(Wz+u,\sigma^2I) \Rightarrow z\vert x \sim N(M^{-1}W^T(x-u),\sigma^{-2}M)\]这里 \(M = W^TW+\sigma^2I\) ，是一个 MxM 的矩阵。代入极大似然的解 \(u = \overline{x}\) ，我们写出后验分布的均值
\[E[z|x] = M^{-1}W^T(x-\overline{x})\]通过观察可以发现，此时的方差 \(\sigma^{-2}M\) 与 观测数据 x 是相互独立的。我们希望方差尽可能地小，取 \(\sigma^2 \rightarrow 0\) ，此时后验分布的均值可以化简为
\[E[z|x] = (W^TW)^{-1}W^T(x-\overline{x})\]这表示从数据空间投影到隐空间，我们又回到了传统的 PCA