Cpt5: 数据降维

我们发现即使我们对数据进行了预处理，使用这些干净的数据直接进行分析和建模可能仍然存在一些问题。机器学习发展至今，需要处理的数据集也跟随问题的复杂性不断增加，数据集的维度也在不断增加。

这些高维数据集可能包含许多特征，但其中一些特征可能是冗余的，或者与目标变量没有强相关性。这些冗余或无关的特征可能会导致模型过拟合，降低模型的泛化能力。

因此，我们也有相对应的解决方案。在这一章中，我们来看看有关数据降维，你需要了解的基本概念和方法。

什么是数据降维？

首先我们需要认识到，我们为什么这么需要数据降维。

假设我们的数据集只有两个维度，我们从一个单位正方形中选择一个随机的点，那么这个点距离正方形边界距离小于0.001的概率是 <0.4%。

哦哦但是如果我们考虑一个10000维的数据集，我们在这个超立方体中选择一个随机的点，那么这个点距离边界距离小于0.001的概率是 >99.99999%。

同样如果还是这两个数据集，我们在里面随机选择两个点，那么对于2D数据集来说，平均距离是0.52，而对于10000维的数据集来说，平均距离是408.25。边界距离的数据告诉我们，在高维空间中，数据点更可能位于边界附近，而不是中心区域。而平均距离的数据告诉我们，在高维空间中，数据点之间的距离变得非常大，这可能导致模型难以捕捉到数据之间的关系。

因此，我们首先想到的第一对策就是，能不能把这些高维数据“精简化”，看看有没有方法丢弃一些我们不需要的特征，或者把一些特征合并成一个新的特征，这样我们就可以在一个更低维的空间中进行分析和建模了。

我们常说的 **数据降维 (Dimensionality Reduction)**就是这样的一个过程，它是指将高维数据映射到一个更低维的空间中，同时尽可能保留原始数据的结构和信息。这样一来，我们使用这些数据训练出来的模型就更加简单，泛化能力也更强了。

投影

工程师想出来的直接对策之一就是 投影(Projection)。投影的思想简单粗暴，就是把高维的数据想办法投影到一个更低维的空间中去。

为什么我们能够这样做，同时不需要考虑夸张的信息丢失？在绝大多数的数据集中，特征之间是存在相关性的。有些特征是冗余的，而有些特征之间存在强相关性。那么我们可以选择合并这些相关特征，同时丢弃一些冗余的特征，这样我们就可以在一个更简单的空间中进行分析和建模了。

如图所示，加入我们有一个三维数据集，而我们想把它投影到一个二维的空间中去，那么我们就可以选择一个平面，把这个三维数据集投影到这个平面上去。这样我们就得到了一个新的二维数据集。

但是投影不是万金油，它有一个很大的问题，就是投影比较吃数据集的分布情况。

下图是一个瑞士卷形状的数据集，你可以很清楚地在3D图像中看出来。mika狂喜 但是当我们把它投影到一个二维的空间中去的时候，我们就完全看不出来这个数据集的结构了：

直接丢失 $x_3$ 的信息会导致Squash，被压扁了。不是很好。

因此，我们需要一些更高级的降维方法，来帮助我们更好地保留数据集的结构和信息。其中一个解决方案就是 Manifold Learning：

流形学习

流形学习(Manifold Learning) 就能够很好的解决瑞士卷问题。

首先我们要知道什么是流形。流形(Manifold) 是一个数学概念，它是一个局部类似于欧几里得空间的拓扑空间。简单来说，流形就是一个在局部看起来像平坦空间的曲面或者空间。比如说，地球表面就是一个流形，在局部看起来像平坦的二维空间，但实际上它是一个三维的球面。

虽然我们的数据集可能存在于一个高维空间中，但它实际上可能存在于一个更低维的流形上。对于瑞士卷，我们可以将它看作是一个二维平面卷起来，形成三维结构。因此，我们可以使用流形学习的方法来找到这个二维流形，并将数据点映射到这个流形上去。这样我们就能够在一个更低维的空间中进行分析和建模了，同时也能够保留数据集的结构和信息。

经过流形学习得出的低维展开一般是是简单可分的，但是经过流形学习降维的数据可能还不如原始数据更难以分离。比如下图还是瑞士卷形状的数据集，但是分布不一：一个顺着流形面方向分布，而另一个垂直于流形面方向分布：

所以我们还需要一些其他的降维方法来帮助我们更好地保留数据集的结构和信息。

主成分分析 (PCA)

主成分分析的内容在之前的笔记中也有所提及，有兴趣可以来看一看

但是鉴于我觉得当时写的比较乱七八糟的，我们在这里重新整理一下吧。

主成分分析，全称Principal Component Analysis，简称PCA，可以说是目前最流行的数据降维方法了。

还记得刚才我们说的投影吗？PCA做的事情就是在所有可能的投影中，找到一个最优的投影，这个投影能够最大程度地保留数据集的结构和信息。因此，我们要找到那个能够使得数据点之间的距离最大化的超平面。

方差

为了选择最佳超平面，我们需要最大限度地保留数据集的方差。我们可以将这个过程理解为在高维空间中寻找一个“最有信息量”的投影方向。

就像你要为一根法棍拍一张照片。如果你从侧面拍摄，你会看到法棍的长度和宽度，这样你就能够获得更多的信息。相反，如果你从顶部拍摄，你只能看到法棍的宽度，而无法看到它的长度，这样你就会丢失一些信息。

在数据集中，方差能够衡量数据在某个方向上的分布广度，也可以理解成“信息量”。因此，我们希望找到一个投影方向，使得数据在这个方向上的方差最大化，这样我们就能够保留更多的信息了。方差大的方向意味着数据点在这个方向上分布更广，法棍就看起来更长；反之，方差小的方向意味着数据点在这个方向上分布更窄，法棍就看起来更不像法棍。

这就引出了主成分的概念。

主成分

主成分(Principal Component) 是指在数据集中方差最大的方向，也可以理解成数据集中最有信息量的方向。第一个主成分是方差最大的方向，第二个主成分是与第一个主成分正交且方差次大的方向，依此类推。

你能看出我们尝试将一个长条形的蓝色的数据集投影到一维轴线上。数据在c1方向上的方差最大，因此c1是第一个主成分。数据在c2方向上的方差次大，因此c2是第二个主成分。如果这时候我们选择按照c1方向保留信息，那么我们就尽可能保留了最多的原始数据信息。

那么我们肯定不能靠人用眼睛来判断主成分方向。我们需要告诉计算机如何使用数学工具计算出任何数据集的主成分。

奇异值分解

PCA有传统的旧方法来算出主成分的方法，这需要我们对数据集计算协方差矩阵，然后对这个矩阵进行特征值分解，最后根据特征值的大小来选择主成分的方向。但是对于现代的数据集来说，这个方法可能会非常慢，因为计算协方差矩阵和进行特征值分解的计算复杂度都很高。那么有没有方法去找到主成分的方向，而不需要计算协方差矩阵和进行特征值分解呢？

隆重介绍，奇异值分解(Singular Value Decomposition, SVD)。与其计算协方差矩阵并进行特征值分解，不如直接对数据矩阵进行奇异值分解。

SVD的核心逻辑是：对于一个给定的数据矩阵 $X$，其中：

$$X \in \mathbb{R}^{n \times d}$$

我们说 $n$ 是数据点的数量，$d$ 是特征的数量。那么我们可以对这个矩阵进行奇异值分解，得到三个矩阵 $U$、$\Sigma$ 和 $V^T$，使得：

$$X = U \Sigma V^T$$

其中：

• $U \in \mathbb{R}^{n \times n}$，列向量是左奇异向量，列向量两两正交，且都是单位向量。这对于数据来说意味着，每个数据点在不同方向上的投影是相互独立的，并且每个投影的长度都是1，这样我们就能够更好地理解数据在不同方向上的分布情况了。

• $\Sigma \in \mathbb{R}^{n \times d}$ 是准对角矩阵，只有主对角线可能有非零元素，其他地方一定是0。我们叫这些非零元素为奇异值，它们按照从大到小的顺序排列。

  很多人不理解奇异值究竟代表什么含义，其实它们代表了数据在不同方向上的重要程度。每个奇异值 $\sigma_i$ 就表示了数据在第i个主成分方向上的方差大小，也就是数据在这个方向上的分布广度。

  我们会将奇异值按照从大到小的顺序排列，这样我们就能够知道哪些主成分方向是最重要的，哪些主成分方向是次要的了。通常来说，我们会选择前k个奇异值对应的主成分方向来进行降维，这样我们就能够保留数据中最重要的信息了。

• $V^T \in \mathbb{R}^{d \times d}$，列向量是右奇异向量，同样也是正交矩阵。在这个矩阵中， $V$ 的列向量对应于数据集的特征方向，也就是我们之前提到的主成分方向。

数学逻辑就是上面这样，那么直觉上怎么理解？我们把 $X$ 看作是一个线性变换。我们输入一个 $R^d$ 的向量，经过这个线性变换之后，我们得到一个 $R^n$ 的向量。而SVD告诉我们这个变换可以拆成三步，也就是 $$X = U \Sigma V^T$$。

想象一下我们把一团数据标准化，然后尝试进行SVD。首先：

1. 由于这个矩阵代表的是原始空间的主成分方向，$V^T$代表我们在输入空间内做一个旋转，尝试将坐标轴与数据集的主方向对齐。

2. 接下来我们要做的是把这些数据在这些方向上“拉开”，让它们的分布反映原始数据的方差。$\Sigma$表示我们在这个新的坐标系中进行缩放，而每个奇异值 $$\sigma_i$$ 就表示我们在第i个主成分方向上的缩放程度。

3. 最后，我们已经在主成分空间中完成了投影和拉伸，现在我们要做的就是把这些数据“搬回”原始空间，让它们的位置和方向能够正确表示。$U$就是这样一个搬运工，表示我们要在输出空间内再做一次旋转，来将数据点映射到最终的空间中去。

投影到主成分空间

一旦我们找到并且选择了前k个主成分方向，我们就可以将数据点投影到这些主成分方向上去，得到一个新的低维数据集。

在数学上，我们通过下面的变换来投影数据点：

$$x = [x_1, x_2, ..., x_d], x \in \mathbb{R}^d$$

$$\downarrow xW, W \in \mathbb{R}^{d \times k}$$

$$z = [z_1, z_2, ..., z_k], z \in \mathbb{R}^k$$

何意味？我们有一个原始数据点 $x$，通过乘以一个权重矩阵 $W$，我们得到了一个新的数据点 $z$，这个新的数据点就是在主成分空间中的表示了。这个权重矩阵 $W$ 的列向量就是我们之前选择的前k个主成分方向。

代码实践

使用python的这些库，我们可以很轻松的使用封装好的PCA类来进行数据降维。还记得之前的葡萄酒数据集吗？我们来看看如何对这个数据集进行降维：

手动挡

首先，我们来做标准化和数据分割：

1
2
3
4
5
6
7
8
9
10
11
12
13
14

import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data', header=None)
df_wine.columns = ['Class', 'Alcohol', 'Malic acid', 'Ash', 'Alcalinity of ash',    'Magnesium', 'Total phenols', 'Flavanoids', 'Nonflavanoid phenols', 'Proanthocyanins', 'Color intensity', 'Hue', 'OD280/OD315 of diluted wines', 'Proline']

# iloc[:, 1:]表示选择所有行，除了第一列以外的所有列作为特征；iloc[:, 0]表示选择所有行，第一列作为标签
X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, stratify=y, random_state=0)
sc = StandardScaler()
X_train_std = sc.fit_transform(X_train)
X_test_std = sc.fit.transform(X_test)

随后，我们来直接对标准化后的训练数据进行奇异值分解：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

import numpy as np

# 直接对标准化后的训练数据进行奇异值分解
U, s, Vt = np.linalg.svd(X_train_std)
print(Vt.shape) # (13, 13)
print(Vt[:2])


# 选择前两个主成分方向
print(X_train[0])
W = Vt.T[:, :2] # .T表示转置
print(W)
X_train_pca = X_train_std.dot(W) 
# 这里我们直接使用了矩阵乘法来进行投影，X_train_std是一个(n_samples, n_features)的矩阵，W是一个(n_features, n_components)的矩阵，结果是一个(n_samples, n_components)的矩阵，也就是我们降维后的数据集了。

print(X_train_pca[0])

经过奇异值分解后，我们就可以得到一个新的数据集，这个数据集在主成分空间中表示了原始数据集的结构和信息了。

自动挡

说得好，但是有没有更简单的方法？

我们直接用传奇伟大scikit-learn库中的PCA类来进行数据降维，几行就ok了：

1
2
3

from sklearn.decomposition import PCA
pca = PCA(n_components=2)
X_train_pca = pca.fit_transform(X_train_std)

完事了。你可以使用matplotlib来画出前后两个数据集的分布情况，你会发现结果是一样的。

方差比例

方差比例告诉我们每个主成分方向上保留了原始数据集的方差的比例，解释了每个主成分方向上保留了多少信息。

你可以通过explained_variance_ratio_属性来查看：

1

print(pca.explained_variance_ratio_)

就好比说这个输出：

1
2

array([0.39809123, 0.12190877, 0.09455811, 0.07313287, 0.06475628, 0.05660977, 0.04934809, 0.03530219, 0.02545516, 0.02144817, 0.01807942, 0.01608624,
       0.01597406])

就说明有39%的方差被第一个主成分方向保留，12%的方差被第二个主成分方向保留，依此类推。

通过这些数据，你就能够知道如果我们选择前k个主成分方向来进行降维，我们就能够保留原始数据集的多少信息。通常来说，我们会选择前k个主成分方向，使得它们的方差比例之和达到一个预设的阈值，比如说95%：

1
2
3
4
5

pca = PCA()
pca.fit(X_train_std)
cumsum = np.cumsum(pca.explained_variance_ratio_)
k = np.argmax(cumsum >= 0.95) + 1
print(k)

除了这种写法，你也可以直接在PCA类中设置n_components=0.95：

1
2
3

pca = PCA(n_components=0.95)
X_train_pca = pca.fit_transform(X_train_std)
print(X_train_pca.shape)

数据压缩

PCA不仅可以用来进行数据降维，还可以用来进行数据压缩。比如我们在MNIST上运用PCA压缩数据，我们就可以用更少的特征来表达原始数据，同时能够尽可能保留信息。

原始的MNIST数据集有784个特征，每一个的分辨率都是28x28。使用PCA后，我们可以将这个数据集降维到一个更低的维度，比如说50维，这样我们就能够用更少的特征来表达原始数据了。

1
2
3
4
5
6
7
8

from sklearn.datasets import fetch_openml
from sklearn.decomposition import PCA
mnist = fetch_openml('mnist_784', version=1)
X = mnist.data.astype('float64')
y = mnist.target.astype('int64')
pca = PCA(n_components=50)
X_reduced = pca.fit_transform(X)
print(X_reduced.shape) # (70000, 50)

更好的是，我们还可以通过逆变换来将降维后的数据重新映射回原始空间中去，这样我们就能够得到一个近似于原始数据的重建数据了：

1
2

X_approx = pca.inverse_transform(X_reduced)
print(X_approx.shape) # (70000, 784)

其他PCA的变体

随机PCA

随机PCA (Randomized PCA) 是PCA的一种变体，它使用了一个随机化的算法来计算PCA。相比于传统的SVD算法，随机PCA在处理大规模数据集时效率更高，能够更快的找到前d个主成分方向。

如果我们将 svd_solver 的超参数设置为 'randomized'，scikit-learn会使用一个随机化的算法来计算PCA，这个算法在处理大规模数据集时效率更高，能够更快的找到前d个主成分方向。相比于SVD来说，当k比d小很多时，随机化算法的计算复杂度更低，因此在处理大规模数据集时更为高效。

1
2

pca = PCA(n_components=2, svd_solver='randomized')
X_train_pca = pca.fit_transform(X_train_std)

默认来讲，svd_solver 的参数默认设置为 'auto'，这意味着scikit-learn会根据数据集的大小和维度来自动选择使用SVD还是随机化算法。当数据集较小或者维度较低时，scikit-learn会选择使用SVD；当数据集较大或者维度较高时，scikit-learn会选择使用随机化算法。

具体来说，当 $d > 500 \text{AND} k < 0.8d$ 的时候，scikit-learn会选择使用随机化算法来计算PCA；否则，scikit-learn会选择使用SVD来计算PCA。

如果你想强制要求scikit-learn使用SVD来计算PCA，你可以将 svd_solver 的参数设置为 'full'：

1
2

pythonpca = PCA(n_components=2, svd_solver='full')
X_train_pca = pca.fit_transform(X_train_std)

增量PCA

随机PCA在计算过程中，要求我们将整个数据集完整地加载到内存中去，这对于一些非常大的数据集来说可能会导致内存不足的问题。

而 增量PCA (Incremental PCA) 则可以解决这个问题，它允许我们在数据集无法完全加载到内存中的情况下进行PCA计算。这样一来，我们不仅能够处理更大规模的数据集，甚至还能处理流式数据。

比如说，我们可以把MNIST数据集分成100个小批次，每次只加载一个小批次的数据到内存中去进行PCA计算：

1
2
3
4
5
6
7

from sklearn.decomposition import IncrementalPCA
n_batches = 100
inc_pca = IncrementalPCA(n_components=87)
for X_batch in np.array_split(X, n_batches):
    inc_pca.partial_fit(X_batch)
X_reduced = inc_pca.transform(X)
print(X_reduced.shape) # (70000, 87)

局部线性嵌入

讲完PCA，我们来简单看一下另外一个流行的非线性降维方法。

局部线性嵌入 (Locally Linear Embedding, LLE) 是一种非线性降维方法，它通过保持数据点在局部邻域内的线性关系来进行降维。

LLE的核心思想是：每个数据点都可以通过它的k个最近邻来进行线性重构，然后我们希望在低维空间中保持这些线性重构关系不变。随后，LLE会在低维空间中找到一个新的表示，使得这些线性重构关系得到最好的保留。这种方法能够很好地展开流形结构的数据集，尤其是当数据集没有太多噪声的时候。

还记得之前的瑞士卷数据集吗？我们可以选择去使用流形学习展开，但是如果我们使用LLE来进行降维，我们就能够更好地保留数据集的结构和信息。

LLE的数学原理

对于每一个训练样本 $\mathbf{x}^{(i)}$，我们首先要告诉LLE要寻找它的k个最近邻，例如 n_neighbors=10。

随后，LLE会尝试将 $\mathbf{x}^{(i)}$ 重建为一个它的k个最近邻的线性组合。要做到这一点，LLE会找到一组权重 $w_{i,j}$，使得 $\mathbf{x}^{(i)}$ 与

$$\sum_{j=1}^m w_{i,j} \mathbf{x}^{(j)}$$

的距离尽可能最小化。并且在这个过程中，我们假定当 $\mathbf{x}^{(j)}$ 不是 $\mathbf{x}^{(i)}$ 的k个最近邻时，权重 $w_{i,j}$ 为0。

所以总的来说，LLE要找到一组权重矩阵 $W$，使得每个数据点 $\mathbf{x}^{(i)}$ 都能够通过它的k个最近邻的线性组合来进行重建：

$$\hat{W} = \arg\min_W \sum_{i=1}^m \left| \mathbf{x}^{(i)} - \sum_{j=1}^m w_{i,j} \mathbf{x}^{(j)} \right|^2$$

在这个过程中，我们还需要满足刚才提到的约束条件。这样重建过程能够确保一定依赖于邻居，并且权重具有归一性，不会导致重建结果的缩放问题。

经过这一步的计算，包含权重 $\hat{w}_{i,j}$ 的权重矩阵 $\hat{W}$会将这些局部线性关系编码到全局训练数据集中去。然后，我们会将训练数据映射到一个k维的空间中去 ($k < d$)，同时尽可能保留这些局部线性关系。我们通过下面的优化问题来实现这个目标：

如果 $\mathbf{z}^{(i)}$ 表示 $\mathbf{x}^{(i)}$ 在k维空间中的表示，那么我们希望找到一组新的表示 $\mathbf{z}^{(i)}$，使得

$$ \sum_{i=1}^m \hat{w}_{i,j} \mathbf{z}^{(j)} $$

的距离尽可能最小化。这意味着原始空间中的线性关系要在低维空间中尽可能得到保留。

其他降维方法

除了PCA和LLE之外，还有很多其他的降维方法，有些也能够在scikit-learn中找到对应的实现。比如下面这些：

• 随机投影 (Random Projection)：通过使用一个随机矩阵来将数据投影到一个更低维的空间中去。这个方法的计算效率非常高，适用于处理大规模数据集。

• 多维度缩放 (Multidimensional Scaling, MDS)：通过保持数据点之间的距离关系来进行降维。这个方法适用于处理非线性结构的数据集。

• Isomap: 通过保持数据点之间的地理距离来进行降维。这个方法适用于处理具有流形结构的数据集。

• t分布随机邻域嵌入 (t-Distributed Stochastic Neighbor Embedding, t-SNE)：通过保持数据点之间的概率分布关系来进行降维。这个方法适用于处理高维数据集，尤其是当数据集具有复杂结构的时候。

• 线性判别分析 (Linear Discriminant Analysis, LDA)：通过最大化类间距离和最小化类内距离来进行降维。这个方法适用于处理有标签的数据集，尤其是当数据集具有多个类别的时候。

这些方法各有优缺点，适用于不同类型的数据集和不同的降维需求。选择合适的降维方法需要根据具体的数据集和任务来进行判断。

练习

学完理论，现在可以来一起看看这些练习，巩固一下我们学到的内容吧：

1. 对数据集特征进行降维的主要目的是什么？我们对特征进行降维，有没有什么潜在的风险？

2. 什么是 Curse of Dimensionality？它对数据分析和建模有什么影响？

3. 一旦数据集的维度被降维，有没有什么方法将这些降维后的数据重新映射回原始空间中去？如果有，这个过程是如何实现的？

4. 如果我们的数据集是非线性的，那么PCA对于这个数据集来说还会有效吗？如果不有效，我们应该使用什么样的降维方法来处理非线性数据集？

5. 如果我们对一个1000维的数据集执行PCA，并且我们将可解释方差比例设置为0.95，那么我们最终得到的降维后的数据集的维度会是多少？这个维度是如何确定的？

6. 在什么情况下我们会使用原版PCA，增量PCA和随机PCA？

7. 我们如何评判一个降维算法的性能？

8. 链式使用两个降维方法是不是一个好主意？

9. 我们载入MNIST数据集，并将其分成训练集和测试集。前60000个样本作为训练集，后10000个样本作为测试集。

我们在这些数据集上训练一个随机森林分类器，并记录下训练时长，并评判模型的性能。

然后，我们使用PCA，选择0.95作为可解释方差比例来对数据进行降维，然后我们在这些数据上训练一个新的随机森林分类器，并记录下训练时长，并评判模型的性能。

那我问你：

• 使用PCA后的数据训练模型的时长和使用原始数据训练模型的时长哪个更长？为什么？

• 使用PCA后的数据训练模型的性能如何？与使用原始数据训练的模型相比，是否有提升或下降？