Principal Components Analysis

　　主成分分析（Principal Components Analysis, PCA）是一种常见的数据分析方式，常用于高维数据的降维，可用于提取数据的主要特征分量，其目标是找到一组低维的基来表示原数据。一般从理论上解释 PCA 原理有两个方向：

1. 一是从最大方差理论（最大可分型）
2. 二是从最小平方误差理论（最近重构性）。

1. PCA 最大方差理论

　　我们想将高维空间的数据映射到一个低维空间，映射的结果考虑是要尽量使得数据分散得开，因为如有数据重叠，那么可以看成信息丢失了。数值的分散程度可以用方差或协方差来衡量，于是目标也就是要最大化投影方差，最大化方差尽量使得投影后的数据之间具有较大的区分度。

　　当然也可以从熵的角度进行理解，熵越大所含信息越多。这样考虑的一个依据是来自信号处理领域的信噪比，我们认为信号具有较大方差，噪声具有较小方差，信号与噪声之比称为信噪比。信噪比越大意味着数据的质量越好，反之，信噪比越小意味着数据的质量越差。如下图拿二维数据为例：

1.1 一维向量视角（拉格朗日乘子法）

　　对于给定的一组数据点 $\{v_1, v_2, \ldots, v_n \}$，其中所有向量均为列向量，中心化（零均值化）后的表示为

　　其中：

　　我们知道，向量内积在几何上表示为第一个向量投影到第二个向量上的长度，因此向量 $x_i$ 在单位方向向量 $w$ 上的投影坐标可以表示为：

　　所以我们的目标是寻找一个一维基 $w$，使得所有数据变换为这个基上的坐标表示后，方差值最大。

　　计算方差之前，我们先计算投影后坐标矩阵的均值：

　　投影后的均值为 $0$，所以在计算方差时就不用考虑均值了，这也是为什么 PCA 要做中心化！那么方差计算如下：

　　仔细一看，其实 $\Sigma = \frac{1}{n} \sum_{i=1}^{n} \boldsymbol{x}_i \boldsymbol{x}_i^{\mathrm{T}} $ 就是样本协方差矩阵。另外，由于 $w$ 是单位方向向量，即有 $\omega^{T} \omega=1$。因此我们要求解一个最大化问题，可表示为

　　引入拉格朗日乘子，并对 $w$ 求导令其等于 0，便可以推出：

　　此时方差可以写成:

　　熟悉线性代数的读者马上就会发现，原来，$x$ 投影后的方差就是协方差矩阵的特征值。我们要找到最大的方差也就是协方差矩阵最大的特征值，最佳投影方向就是最大特征值所对应的特征向量。次佳投影方向位于最佳投影方向的正交空间中，是第二大特征值对应的特征向量，以此类推。

1.2 多维矩阵视角（矩阵对角化）

　　设原始 $n\times m$ 维的数据矩阵 $V_{n\times m}$，对应的均值矩阵为 $U_{n\times m}=\frac{1}{n}V$，中心化后的矩阵记为 $X_{n \times m}=V-U$，则样本的协方差矩阵为 $C$：

　　可见 $C$ 是一个对称矩阵，其对角线分别对应各个变量的方差，而第 $i$ 行 $j$ 列和 $j$ 行 $i$ 列元素相同，表示 $i$ 和 $j$ 两个变量的协方差。

　　再设 $P$ 是一组基按行组成的矩阵，设 $Y=PX$，则 $Y$ 为 $X$ 对 $P$ 做基变换后的数据。设 $Y$ 的协方差矩阵为 $D$，我们推导一下 数据变换前后的协方差 $C$ 与 $D$ 的关系：

　　有了这个公式后，我们再总结一下 PCA 的优化目标有两个：

1. 将为后同一维度的方差最大
2. 不同维度之间的相关性为零

　　根据线性代数可知：

1. 同一元素的协方差就表示该元素的方差
2. 不同元素之间的协方差表示它们的相关性

　　这两个优化目标偶可以用协方差矩阵来表示：

　　样本降维后希望得到的最理想的协方差矩阵（包含最多信息），即优化目标就是令类似上面的 $C_Y$ 矩阵的 $D$ 矩阵的对角线元素之和最大，即 $\max \operatorname{tr}(D)$。当然，我们要求 $P$ 作为基是正交的，所以：

　　则最终的优化目标可以表示为：

　　根据拉格朗日乘子法可以得到：

　　其中迹（tr）表示对角线元素的和，这里用到其重要公式：

　　对拉格朗日乘子法结果求导并取零点有：

　　当 $\frac{\partial f}{\partial P} = 0$ 时，得到：

　　即优化目标只要对协方差矩阵做特征值分解即可。

总结一下 PCA 求解步骤：

1. 对样本数据进行中心化处理。
2. 求样本协方差矩阵：$\Sigma = \frac{1}{n} \sum_{i=1}^{n} \boldsymbol{x}_i \boldsymbol{x}_i^{\mathrm{T}} $。
3. 对协方差矩阵进行特征值分解，将特征值从大到小排列。
   • 通过奇异值分解（SVD 求取 $\Sigma$ 的 特征向量（eigenvectors）
     • $\left(U, S, V^{T}\right)=SVD(\Sigma)$

4. 从 $U$ 中取出前 $d$ 列左奇异特征向量 $U_\text{sub}= \{\omega_{1}, \omega_{2}, \dots, \omega_{d}\}$，通过以下映射将 $n$ 维样本映射到 $d$ 维

　　新的 $x_i^\prime$ 的第 $d$ 维就是 $x_i$ 在第 $d$ 个主成分 $w_d$ 方向上的投影，通过选取最大的 $d$ 个特征值对应的特征向量，我们将方差较小的特征（噪声）抛弃，使得每个 $n$ 维列向量 $x_i$ 被映射为 $d$ 维列向量 $x_i^\prime$，定义降维后的信息占比为

　　降维后的维度 $d^\prime$ 可以通过如下几种方式得到：

1. 由用户根据经验事先指定
2. 通过在 $d^\prime$ 不同的地位空间中对 $k$ 近邻分类器（或其他小开销的学习器）进行交叉验证来选取较好的 $d^\prime$ 值
3. 还可以从重构的角度设置一个重构阈值，例如 $t=95 \%$，然后通过选取使得下式成立的最小 $d^\prime$ 值

　　PCA 通过线性变换将向量投影到低维空间。对向量进行投影就是对向量左乘一个矩阵，得到结果向量，所以其决策函数为： \(y=W x\)

　　其中 $W=U_\text{sub}^T$。

　　先找一个主成分，然后再找第二个一直往后找。在找第二个时，还要加上一个新的约束，与第一个主成分相互垂直，不然第二次还是找到第一个主成分。

2. PCA 最小平方误差理论

　　为什么说 PCA 是线性降维方法呢？可以看下其决策函数就可以看到，就是一个线性计算，那么我们就可以将其转换一下思路，其目标也是求解一个线性函数使得对应直线能够更好地拟合样本点集合。如果我们从这个角度定义PCA的目标，那么问题就会转化为一个回归问题。

　　数据集中每个点 $x_k$ 到 $d$ 维超平面 $D$ 的距离为

　　其中 $\widetilde{\boldsymbol{x}_{k}}$ 表示 $x_k$ 在超平面 $D$ 上的投影向量。如果该超平面由 $d$ 个标准正交基构成，根据线性代数理论可以由这组基线性表示

　　其中 $w_i^T x_k$ 表示 $x_k$ 在 $w _i$ 方向上投影的长度。因此，$\widetilde{\boldsymbol{x}_k}$ 实际上就是 $x_k$ 在 $W$ 这组标准正交基下的坐标。而 PCA 要优化的目标为

　　由向量内积的性质，我们知道 $\boldsymbol{x}_k^{\mathrm{T}} \widetilde{\boldsymbol{x}_k}=\widetilde{\boldsymbol{x}_k}^{\mathrm{T}} \boldsymbol{x}_k$，于是将上式中的每一个距离展开

　　其中第一项 $x_k^Tx_k$ 与选取的 $W$ 无关，是常数，则有：

　　注意到，其中 $\omega_i^{\mathrm{T}} \boldsymbol{x}_k$ 和 $\omega_j^{\mathrm{T}} \boldsymbol{x}_k$ 表示投影长度，都是数字。且当 $i \neq j$ 时，$\omega_i^{\mathrm{T}} \omega_j=0$，因此有:

　　注意到，$\sum_{i=1}^{d} \omega_{i}^{\mathrm{T}} \boldsymbol{x}_k \boldsymbol{x}_k^{\mathrm{T}} \boldsymbol{\omega}_i$ 实际上就是矩阵 $\boldsymbol{W}^{\mathrm{T}} \boldsymbol{x}_k \boldsymbol{x}_k^{\mathrm{T}} \boldsymbol{W}$ 的迹（对角线元素之和），于是可以距离公式化简

　　因此式之前式子可以写成

　　根据矩阵乘法的性质，因此优化问题可以转化为 $\arg \max _ {W} \sum_{k=1}^{n} \operatorname{tr}\left(\boldsymbol{W}^{\mathrm{T}} \boldsymbol{x}_k \boldsymbol{x}_k^{\mathrm{T}} \boldsymbol{W}\right)$，这等价于求解带约束的优化问题：

　　如果我们对 $W$ 中的 $d$ 个基 $\omega_{1}, \omega_{2}, \ldots, \omega_{d}$ 依次求解，就会发现和最大方差理论的方法完全等价。比如当 $d=1$ 时，我们实际求解的问题是

　　最佳直线 $w$ 与最大方差法求解的最佳投影方向一致，即协方差矩阵的最大特征值所对应的特征向量，差别仅是协方差矩阵 $\Sigma$ 的一个倍数，以及常数 $\sum_{k=1}^{n} \boldsymbol{x}_k^{\mathrm{T}} \boldsymbol{x}_k$ 偏差，但这并不影响我们对最大值的优化。

3. 核化 PCA

　　传统的 PCA 只能做线性降维，对于非线性的情况效果就不尽人意。那么是不是可以像 SVM 那样，我们先把数据映射到高维空间，在高维空间就能用 PCA 线性降维了！

　　TODO 待重新整理，周志华图书中说的不清楚，可参考附录 3，说的比较清楚。

　　假设 $z_i = \phi\left(\boldsymbol{x}_{i}\right)$，$i=1,2, \ldots, m$， 是样本 $x_i$ 在高维特征空间中由 $w$ 确定的超平面上的像，因为高维特征空间可以用 PCA 线性降维，则具体就是求解如下的式子：

　　做一下变换求解 $\boldsymbol{\omega}$:

　　其中 $\boldsymbol{\alpha}_i=\frac{1}{\lambda} \boldsymbol{z}_i^{\mathrm{T}} \boldsymbol{\omega}$，一般情况下我们不知道 $\phi$ 的具体形式，于是引入核函数：

　　如果我们记 $\mathbf{K}$ 为 $\kappa$ 对应的核矩阵，$\mathbf{K}_{ij}=\kappa\left(\boldsymbol{x}_i, \boldsymbol{x}_j\right)$，$\mathbf{A}=\left(\boldsymbol{\alpha}_1 ; \boldsymbol{\alpha}_2 ; \ldots ; \boldsymbol{\alpha}_m\right)$，则上面式子可以化简为：

　　这里就跟传统的 PCA 一样了，只剩下特征值分解问题了，取 $\mathbf{K}$ 最大的 $d^\prime$ 个特征值对应的特征向量即可。

　　对于新样本 $x$，其投影后的第 $j$ $\left(j=1,2, \ldots, d^{\prime}\right)$ 维坐标为：

　　其中 $\boldsymbol{\alpha}_i$ 已经规范化，$\boldsymbol{\alpha}_i^j$ 是 $\boldsymbol{\alpha}_i$ 的第 $j$ 个分量。我们可以看到为获得投影后的坐标，KPCA 需要做所有样本的加和，计算开销比较大。

4. 应用

4.1 选择合适的维度

　　看特征值累加的比例，比如我们去 95% 来选择维度，可以自己这么做：

pca = PCA()
pca.fit(X_train)
cumsum = np.cumsum(pca.explained_variance_ratio_)
d = np.argmax(cumsum >= 0.95) + 1

　　在 Sklearn 里初始化时就能选定：

pca = PCA(n_components=0.95)
X_reduced = pca.fit_transform(X_train)
pca.n_components_
np.sum(pca.explained_variance_ratio_)

4.2 核化 PCA

from sklearn.decomposition import KernelPCA

rbf_pca = KernelPCA(n_components = 2, kernel="rbf", gamma=0.04)
X_reduced = rbf_pca.fit_transform(X)

选择合适的核并调整参数：

from sklearn.model_selection import GridSearchCV
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import Pipeline

clf = Pipeline([
        ("kpca", KernelPCA(n_components=2)),
        ("log_reg", LogisticRegression(solver="liblinear"))
    ])

param_grid = [{
        "kpca__gamma": np.linspace(0.03, 0.05, 10),
        "kpca__kernel": ["rbf", "sigmoid"]
    }]

grid_search = GridSearchCV(clf, param_grid, cv=3)
grid_search.fit(X, y)

print(grid_search.best_params_)

4.3 如何选择用什么类型的 PCA

1. Vanilla PCA 是默认的 PCA 模式，一般只在数据能够全量放进内存的情况下使用。
2. 增量 PCA 则能解决数据太大不能放进内存的情况，但是速度会相对稍逊一筹；并且对于一些在线任务，增量 PCA 能接受随时加入的新样本，会更加合适一些。
3. 随机化 PCA 当你期望大剂量的降低数据维度时，且数据能放进内存时可做首选，此时比 Vanilla PCA 会更快。
4. 而核化 PCA 主要针对的是非线性的数据集。

References

1. Linear Discriminant Analysis – Bit by Bit
2. 《机器学习》周志华
3. Kernel Principal Components Analysis
4. Feature-Engineering中文版第6章：降维：用 PCA 压缩数据集
5. 降维——PCA（非常详细）
6. PCA算法原理：为什么用协方差矩阵