降维实现

2019-05-26 本文已影响0人单调不减

1、维度诅咒

许多机器学习问题的每个训练实例涉及数千甚至数百万的特征。正如我们将要看到的，这不仅使训练变得非常缓慢，而且还可能使得找到一个好的解决方案变得更加困难。这个问题通常被称为维度的诅咒。

在高维空间中的一个麻烦在于，高维超立方体中的大多数点非常靠近边界。例如，如果选择单位平方（1×1平方）内的随机点，它只有0.4％的可能性位于距离边界小于0.001的位置（换句话说，随机点在任何维度上都不太可能“极端”)，但是在10000维单元超立方体（1×1×⋯×1立方体，具有1万1个）中，这个概率大于99.999999％。

在高维空间中的另一个麻烦在于，大多数训练实例可能彼此相距很远。如果在单位平方中随机选取两个点，这两个点之间的距离平均约为0.52。如果在单位立方体中选取两个随机点，则平均距离大约为0.66。但是在1000000维超立方体中随机选取两个点呢？平均距离大概是408.25（大约1,000,000 / 6）！这一事实意味着高维数据集存在非常稀疏的风险：大多数训练实例可能彼此相距很远。当然，这也意味着新实例可能远离任何训练实例，使得预测比较低维度更不可靠，因为它们将基于更大的推断。简而言之，训练集的维度越多，过拟合的风险就越大。

理论上，维度诅咒的一个解决方案是增加训练集的大小以达到足够密度的训练实例。不幸的是，在实践中，达到给定密度所需的训练实例的数量随着维度的数量呈指数增长。对于100个特征（远少于MNIST问题），我们需要比宇宙中可观察的原子更多的训练实例，才能使训练实例平均距离在0.1之内（假设它们在所有维度上均匀分布），这是无法忍受的。

2、降维的用处

在现实世界的问题中，通常可以减少特征的数量，将棘手的问题变成易处理的问题。

例如，考虑MNIST图像（在第3章中介绍）：图像边框上的像素几乎总是白色，因此可以从训练集中完全丢弃这些像素而不会丢失太多信息。
此外，两个相邻像素通常高度相关：如果将它们合并为单个像素（例如，通过取两个像素强度的平均值），将不会丢失太多信息。

降维确实会丢失一些信息，但它会加快训练速度，在某些情况下，减少训练数据的维数还可以过滤掉一些噪音和不必要的细节，从而得到更好的性能。

当然，降维在实际运用中最重要的一个用处是有助于对数据进行可视化。将维数减少到两（或三）维后我们可以在图上绘制高维训练集，通过观察来发现数据中潜在的模式，从而获得一些重要的见解。

3、降维的主要方法

降维的两类主要方法分别是投影和流形学习。

3.1 投影

在大多数现实生活的问题中，训练实例并不是在所有维度上均匀分布的。许多特征几乎是常数，而其他特征则高度相关（如前面讨论的 MNIST）。结果，所有训练实例实际上位于（或接近）高维空间的低维子空间内。例如：

此时我们就可以通过把高维空间中的数据投影到一个低维的超平面上来实现降维。但是，投影并不总是降维的最佳方法，在很多情况下，子空间可能会扭曲和转动，比如：

面对这种瑞士卷数据集，我们需要采用流行学习。

3.2 流形学习

许多降维算法通过对训练实例所在的流形进行建模从而达到降维目的，这叫做流形学习。它依赖于流形猜想，也称为流形假设，它认为大多数现实世界的高维数据集大都靠近一个更低维的流形。这种假设经常在实践中被证实。

回到 MNIST 数据集：所有手写数字图像都有一些相似之处。它们由连线组成，边界是白色的，大多是在图片中中间的等等。如果随机生成图像，只有一小部分看起来像手写数字。换句话说，如果尝试创建数字图像，那么自由度远低于生成任何随便一个图像时的自由度。这些约束往往会将数据集压缩到较低维流形中。

流形假设通常包含着另一个隐含的假设：你在的上的工作（例如分类或回归）如果在流形的较低维空间中表示，那么它们会变得更简单。比如在三维空间中（图左），分类边界会相当复杂，但在二维展开的流形空间中（图右），分类边界是一条简单的直线：

但是这个假设并不总是成立。例如在原始三维空间中，决策边界位于x1 = 5（图左）。这个决策边界看起来非常简单，但在展开的流形中却变得更复杂了（图右）：

简而言之，如果在训练模型之前降低训练集的维数，那训练速度肯定会加快，但并不总是会得出更好的训练效果；这一切都取决于数据集。

4、PCA

主成分分析（Principal Component Analysis，PCA）是目前为止最流行的降维算法。它首先找到接近数据集分布的超平面，然后将所有的数据都投影到这个超平面上。

关于PCA寻找超平面的主要思想有两种解释：

选择这个超平面使得将原始数据集投影到该超平面上的均方距离最小。这一点上一节的投影已经说过，训练实例并不是在所有维度上均匀分布的，许多特征几乎是常数，而其他特征则高度相关，我们在寻找超平面的过程中，自然希望损失的是几乎是常数的那部分特征的信息（因为根据香农信息论的思想，最不出乎意料意味着包含最少的信息量），而“几乎是常数”或“变化不大”就体现在各样本点到超平面的距离。因此我们希望原始数据集投影到该超平面上的均方距离最小。
选择保持最大方差的超平面。这种解释和上面的解释其实是一个硬币的两面。

如何找到训练集的主要组成部分呢？幸运的是，有一种称为奇异值分解（SVD）的标准矩阵分解技术可以将训练集矩阵分解为 $M=U\Sigma V^T$ ，其中 $V^T$ 包含我们想要的所有主成分。（ $V^T$ 具体含义见https://www.jianshu.com/p/8c7dac32620f）

接下来我们来实际操作一下，首先生成三维数据集，然后使用NumPy的svd()函数获取训练集的所有主成分，并提取前两个PC：

#随机生成3维数据集
np.random.seed(4)
m = 60
w1, w2 = 0.1, 0.3
noise = 0.1
#angles为60个-0.5到3pi/2-0.5之间的随机数
angles = np.random.rand(m) * 3 * np.pi / 2 - 0.5
#X为60*3的随机矩阵
X = np.empty((m, 3))

X[:, 0] = np.cos(angles) + np.sin(angles)/2 + noise * np.random.randn(m) / 2
X[:, 1] = np.sin(angles) * 0.7 + noise * np.random.randn(m) / 2
X[:, 2] = X[:, 0] * w1 + X[:, 1] * w2 + noise * np.random.randn(m)

#PCA using SVD decomposition
X_centered = X-X.mean(axis=0)
U,s,Vt = np.linalg.svd(X_centered)
c1 = Vt.T[:,0]
c2 = Vt.T[:,1]
s

out:
array([6.77645005, 2.82403671, 0.78116597])

注意这里的SVD分解得到的s(sigma)并不是矩阵的形式，而是 $\Sigma$ 对角元素组成的数组。

m,n = X.shape

S = np.zeros(X_centered.shape)
#将上面得到的数组s作为对角元素构成对角矩阵
S[:n,:n] = np.diag(s)
#numpy的allclose方法，比较两个array是不是每一元素都相等，默认在1e-05的误差范围内
np.allclose(X_centered,U.dot(S).dot(Vt))

out:
True

上述过程检验了SVD分解的正确性。

接下来我们取前两个主成分，并将训练集投影到由前两个主成分定义的平面上：

#提取前两个主成分
W2 = Vt.T[:,:2]
#将原数据集投影到上面主成分定义的平面上
X2D_using_svd  = X_centered.dot(W2)
X2D_using_svd 

out:
array([[-1.26203346, -0.42067648],
       [ 0.08001485,  0.35272239],
       [-1.17545763, -0.36085729],
       [-0.89305601,  0.30862856],
       [-0.73016287,  0.25404049]])

可以看到，我们已经得到了原三位数据在二维平面上的投影，实现了降维。

下面我们使用Scikit-Learn的PCA类来运行SVD从而实现PCA（请注意，它会自动处理数据居中）：

from sklearn.decomposition import PCA

pca = PCA(n_components = 2)
X2D = pca.fit_transform(X)
X2D

out:
array([[ 1.26203346,  0.42067648],
       [-0.08001485, -0.35272239],
       [ 1.17545763,  0.36085729],
       [ 0.89305601, -0.30862856],
       [ 0.73016287, -0.25404049]])

np.allclose(X2D,-X2D_using_svd)

out:
True

我们发现，使用PCA类降维得到的结果刚好是我们之前结果的相反数，这其实并不奇怪。主成分的方向不是稳定的：如果稍微扰乱训练集并再次运行PCA，一些新PC可能指向与原始PC相反的方向。但是，它们通常仍然位于相同的轴上。也就是说，我们所得到的新的降维结果之所以是之前结果的相反数，只是因为求出的主成分方向与之前相反，因此两者所表示的在标准基下的坐标其实是相同的。

下面我们验证这一点（PCA对象对其计算的主成分提供的访问(access)）：

pca.components_

out:
array([[-0.93636116, -0.29854881, -0.18465208],
       [ 0.34027485, -0.90119108, -0.2684542 ]])

Vt[:2]

out:
array([[ 0.93636116,  0.29854881,  0.18465208],
       [-0.34027485,  0.90119108,  0.2684542 ]])

接下来我们由投影得到的二位数据还原到三维空间：

X3D_inv = pca.inverse_transform(X2D)
np.allclose(X3D_inv,X)

out:
False

可以看到在投影步骤中有一些信息丢失，所以恢复的3D点并不完全等于原始的3D点。

PCA中一个非常有用的信息是每个主成分的解释方差比，可通过explain_varianceratio变量获得。它表示沿着每个主成分的轴的数据集方差的比例：

pca.explained_variance_ratio_
1 - pca.explained_variance_ratio_.sum()

out:
array([0.84248607, 0.14631839])
0.011195535570688975

这告诉我们84.2％的数据集的方差位于第一个轴上，14.6％位于第二个轴。对于第三个轴，留下不到1.2％，因此可以合理地假设它可能携带很少的信息。

以下是如何使用SVD方法计算解释的方差比（回想一下s是矩阵S的对角线）：

np.square(s) / np.square(s).sum()

out:
array([0.84248607, 0.14631839, 0.01119554])

接下来的部分就是数据可视化了：

首先是一个用于绘制3D箭头的实用程序类：

from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d import proj3d

class Arrow3D(FancyArrowPatch):
    def __init__(self, xs, ys, zs, *args, **kwargs):
        FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
        self._verts3d = xs, ys, zs

    def draw(self, renderer):
        xs3d, ys3d, zs3d = self._verts3d
        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
        self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
        FancyArrowPatch.draw(self, renderer)

将平面表示为x1和x2的函数：

axes = [-1.8, 1.8, -1.3, 1.3, -1.0, 1.0]

x1s = np.linspace(axes[0], axes[1], 10)
x2s = np.linspace(axes[2], axes[3], 10)
x1, x2 = np.meshgrid(x1s, x2s)

C = pca.components_
R = C.T.dot(C)
z = (R[0, 2] * x1 + R[1, 2] * x2) / (1 - R[2, 2])

在该平面上绘制三维数据集，平面和投影：

from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(6, 3.8))
ax = fig.add_subplot(111, projection='3d')

X3D_above = X[X[:, 2] > X3D_inv[:, 2]]
X3D_below = X[X[:, 2] <= X3D_inv[:, 2]]

ax.plot(X3D_below[:, 0], X3D_below[:, 1], X3D_below[:, 2], "bo", alpha=0.5)

ax.plot_surface(x1, x2, z, alpha=0.2, color="k")
np.linalg.norm(C, axis=0)
ax.add_artist(Arrow3D([0, C[0, 0]],[0, C[0, 1]],[0, C[0, 2]], mutation_scale=15, lw=1, arrowstyle="-|>", color="k"))
ax.add_artist(Arrow3D([0, C[1, 0]],[0, C[1, 1]],[0, C[1, 2]], mutation_scale=15, lw=1, arrowstyle="-|>", color="k"))
ax.plot([0], [0], [0], "k.")

for i in range(m):
    if X[i, 2] > X3D_inv[i, 2]:
        ax.plot([X[i][0], X3D_inv[i][0]], [X[i][1], X3D_inv[i][1]], [X[i][2], X3D_inv[i][2]], "k-")
    else:
        ax.plot([X[i][0], X3D_inv[i][0]], [X[i][1], X3D_inv[i][1]], [X[i][2], X3D_inv[i][2]], "k-", color="#505050")
    
ax.plot(X3D_inv[:, 0], X3D_inv[:, 1], X3D_inv[:, 2], "k+")
ax.plot(X3D_inv[:, 0], X3D_inv[:, 1], X3D_inv[:, 2], "k.")
ax.plot(X3D_above[:, 0], X3D_above[:, 1], X3D_above[:, 2], "bo")
ax.set_xlabel("$x_1$", fontsize=18)
ax.set_ylabel("$x_2$", fontsize=18)
ax.set_zlabel("$x_3$", fontsize=18)
ax.set_xlim(axes[0:2])
ax.set_ylim(axes[2:4])
ax.set_zlim(axes[4:6])

save_fig("dataset_3d_plot")
plt.show()

投影后的二维平面如下：


fig = plt.figure()
ax = fig.add_subplot(111, aspect='equal')

ax.plot(X2D[:, 0], X2D[:, 1], "k+")
ax.plot(X2D[:, 0], X2D[:, 1], "k.")
ax.plot([0], [0], "ko")
ax.arrow(0, 0, 0, 1, head_width=0.05, length_includes_head=True, head_length=0.1, fc='k', ec='k')
ax.arrow(0, 0, 1, 0, head_width=0.05, length_includes_head=True, head_length=0.1, fc='k', ec='k')
ax.set_xlabel("$z_1$", fontsize=18)
ax.set_ylabel("$z_2$", fontsize=18, rotation=0)
ax.axis([-1.5, 1.3, -1.2, 1.2])
ax.grid(True)
save_fig("dataset_2d_plot")

关于PCA还有一个重要问题，那就是主成分数目的选择。通常倾向于选择方差部分加起来足够大的维数（例如，95％），而不是随意选择要减少的维数。当然，除非正在减少数据可视化的维度——在这种情况下，通常希望将维度降低到2或3。

可以将n_components设置为介于0.0和1.0之间的浮点数，而不是指定要保留的主成分的数量，这表示希望保留的方差比：

pca = PCA(n_components=0.95)
X_reduced = pca.fit_transform(X)

5、Kernel PCA

核技巧可以应用于PCA，使得我们有可能执行复杂的非线性投影以减少维数。这称为核PCA（kPCA）。核PCA通常擅长在投影后保留实例簇，或者有时甚至展开靠近扭曲流形的数据集。

下面我们使用瑞士卷数据进行测试，实验Linear Kernel、RBF Kernel和Logistic Kernel的降维效果：

from sklearn.datasets import make_swiss_roll
X, t = make_swiss_roll(n_samples=1000, noise=0.2, random_state=42)

from sklearn.decomposition import KernelPCA

lin_pca = KernelPCA(n_components = 2, kernel="linear", fit_inverse_transform=True)
rbf_pca = KernelPCA(n_components = 2, kernel="rbf", gamma=0.0433, fit_inverse_transform=True)
sig_pca = KernelPCA(n_components = 2, kernel="sigmoid", gamma=0.001, coef0=1, fit_inverse_transform=True)

y = t > 6.9

plt.figure(figsize=(11, 4))
for subplot, pca, title in ((131, lin_pca, "Linear kernel"), (132, rbf_pca, "RBF kernel, $\gamma=0.04$"), (133, sig_pca, "Sigmoid kernel, $\gamma=10^{-3}, r=1$")):
    X_reduced = pca.fit_transform(X)
    if subplot == 132:
        X_reduced_rbf = X_reduced
    
    plt.subplot(subplot)
    #plt.plot(X_reduced[y, 0], X_reduced[y, 1], "gs")
    #plt.plot(X_reduced[~y, 0], X_reduced[~y, 1], "y^")
    plt.title(title, fontsize=14)
    plt.scatter(X_reduced[:, 0], X_reduced[:, 1], c=t, cmap=plt.cm.hot)
    plt.xlabel("$z_1$", fontsize=18)
    if subplot == 131:
        plt.ylabel("$z_2$", fontsize=18, rotation=0)
    plt.grid(True)

save_fig("kernel_pca_plot")
plt.show()

由于kPCA是一种无监督学习算法，没有明显的性能指标来帮助我们选择最佳内核和超参数值。但是，降维通常是监督学习任务（例如，分类）的准备步骤，因此可以简单地使用网格搜索来选择使该任务具有最佳性能的内核和超参数。

6、LLE

局部线性嵌入（Locally Linear Embedding，LLE）是另一种非常强大的非线性降维（NLDR）技术。它是一种流形学习技术，不依赖于先前算法的预测。简而言之，LLE的工作原理是:

首先测量每个训练实例如何与其最近邻居(c.n.)线性相关，
然后寻找训练集的低维表示，其中这些局部关系得到最好的保留。

这使得它特别擅长展开扭曲的流形，特别是在没有太多噪音的情况下。

下面我们看一下LLE在瑞士卷数据集上的表现：

X, t = make_swiss_roll(n_samples=1000, noise=0.2, random_state=41)

from sklearn.manifold import LocallyLinearEmbedding

lle = LocallyLinearEmbedding(n_components=2,n_neighbors=10)
X_reduced = lle.fit_transform(X)

plt.title("Unrolled swiss roll using LLE", fontsize=14)
plt.scatter(X_reduced[:, 0], X_reduced[:, 1], c=t, cmap=plt.cm.hot)
plt.xlabel("$z_1$", fontsize=18)
plt.ylabel("$z_2$", fontsize=18)
plt.axis([-0.065, 0.055, -0.1, 0.12])
plt.grid(True)

save_fig("lle_unrolling_plot")
plt.show()