In [1]:
# setup
from __future__ import division, print_function, unicode_literals

import numpy as np
import os

np.random.seed(42)

%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['axes.labelsize'] = 14
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12

终于到了激动人心的维度灾难！

机器学习时涉及到几千甚至数百万维的特征的训练实例，这会导致训练过程非常慢，而且很难找到一个很好地解，甚至发生维数灾难。

我们可以将一个极大维度维度问题降低其特征维度将一个棘手的问题转变成一个可以较为容易解决的问题。

降维肯定会丢失一些信息，因此即使加快了训练的速度，但也会让你的系统表现的稍微差一点。

降维对数据可视化也非常重要，二维和三维是可以画出的，我们可以通过视觉直观发现一些重要的信息，比如聚类。

我们将讨论两种主要的降维方法：主成分分析 PCA 和核主成分分析 Kernel PCA 和局部线性嵌入 LLE。

维数灾难

  • 高纬度的物体的大多数点都分布在边界处
  • 高纬度的数据集有很大风险分布非常稀疏
  • 维度越高，过拟合风险越大

这些事增大维度以后我们将会面临的事情，统称 维数灾难 。

理论上说，解决方案可以是增加训练集的大小，不幸的是，需要的训练实例呈指数增长。

降维

降维有两个主要方法：

  • 投影
  • 流形学习

投影 Projection

训练实例并不是所有维度上均匀分布，许多特征几乎是常数，而其他特征则高度相关，结果所有训练实例实际上位于高维空间的低维子空间内。

但在很多情况下，子空间可能会扭曲和转动，比如瑞士滚动卷数据集。

流形学习

瑞士卷是一个二维流形的例子，二维流行是一中二维形状，它可以在更高维空间中弯曲或扭曲。

一个 d 维流形是类似于 d 维超平面的 n 维空间的一部分。

瑞士卷中， d = 2 , n = 3 。它是 2d 平面在 三维 中卷曲。

我们可以根据 流形假设 manifold assumption 对流形进行建模，做到降维的效果。这种假设经常在实践中被证实。

In [5]:
from sklearn.datasets import make_swiss_roll
from mpl_toolkits.mplot3d import proj3d
# 瑞士卷
X, t = make_swiss_roll(n_samples=1000, noise=0.2, random_state=42)
# 大小
axes = [-11.5, 14, -2, 23, -12, 15]

fig = plt.figure(figsize=(6, 5))
ax = fig.add_subplot(111, projection='3d')
# 瑞士卷显示
ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=t, cmap=plt.cm.hot)
# 视角
ax.view_init(10, -70)
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])

plt.show()
In [6]:
plt.figure(figsize=(11, 4))

plt.subplot(121)
plt.scatter(X[:, 0], X[:, 1], c=t, cmap=plt.cm.hot)
plt.axis(axes[:4])
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$x_2$", fontsize=18, rotation=0)
plt.grid(True)

plt.subplot(122)
plt.scatter(t, X[:, 1], c=t, cmap=plt.cm.hot)
plt.axis([4, 15, axes[2], axes[3]])
plt.xlabel("$z_1$", fontsize=18)
plt.grid(True)

plt.show()

PCA

主成分分析 Principal Component Analysis ：首先周到接近数据集分布的超平面，然后将所有数据都投影到这个超平面上。

保留（最大）方差

假设我们的二维数据集要向一维进发，有三个轴可供选择，应该选择在投影后有最大方差的轴，因为他可能比其他投影损失更少的信息。

选择这个轴使得原始数据集投影到该轴上的均方距离最小，这就是 PCA 背后的思想。

主成分 Principle Components

PCA 寻找训练集中会的最大方差的轴。

以一条直线为例，我们选择一条轴，然后选择与其正交的第二个轴，计算最大的残差。如果是更高维，我们可以找到第三个轴，第四个轴...定义第 i 个轴的单位矢量为第 i 个主成分 PC 。第一个 PC 是 c1 ，第二个 PC 是 c2 ...

那么如何找到训练集的主成分 PC 呢？有一种称为奇异值分解 (SVD) 的标准矩阵分解技术，可以将训练集矩阵 X 分解成三个矩阵 $U\cdot \sum \cdot V^{T}$ ，其中的 $V^{T}$ 包含我们想要的所有主成分。

$$V^{T}=\begin{pmatrix} |& |& & |\\ c_{1}& c_{2}& ...&c_{n} \\ |& |& & | \end{pmatrix}$$

In [84]:
np.random.seed(4)
m = 60
w1, w2 = 0.1, 0.3
noise = 0.1

angles = np.random.rand(m) * 3 * np.pi / 2 - 0.5
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)
In [9]:
# 数据中心化
X_centered = X - X.mean(axis=0)
U, s, Vt = np.linalg.svd(X_centered)
c1 = Vt.T[:, 0]  # Vt[0]
c2 = Vt.T[:, 1]  # Vt[1]

投影到 d 维空间

我们确认主成分以后，就可以通过将数据集投影到前 d 个主成分构成的超平面上，从而将数据集的维数降至 d 维。选择这个超平面可以确保投影将保留尽可能多的方差。

训练集投影到 d 维空间：

$$X_{d-proj}=X\cdot W_{d}$$

其中 w_{d} 表示前 d 个主成分的矩阵，即 V^{T} 的前 d 列组成的矩阵。

In [28]:
# 将训练集投影到由前两个主成分定义的超平面上：
W2 = Vt.T[:, :2]
X2D = X_centered.dot(W2)
# 这个套路可以将任何一个数据集降维并尽可能地保留原始数据集了。

使用 sklearn

sklearn 的 PCA 类使用 SVD 分解来实现，就像上面的做法，下面应用 PCA 将数据集的维度降至两维。

ps: 它自动进行了数据中心化处理

In [30]:
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
X2D = pca.fit_transform(X)
# PCA 转化以后，使用 components_ 访问每一个主成分

方差解释率 Explained Variance Ratio

另一个非常有用的信息是每个主成分的方差解释率，可以通过 explained_varianceratio 获得，它表示位于每个主成分轴上的数据方差的比例。

比如下面是一个三维数据集前两个分量的方差解释率。

In [39]:
pca.explained_variance_ratio_
Out[39]:
array([0.84248607, 0.14631839])

选择正确的维度

通常我们倾向于加起来的方差解释率能够占比足够的维度的数量，而不是任意选择要降低到的维度数量。下面计算出保留训练集方差95%所需的最小维度数。

In [47]:
# 配置代理
import socket
import socks
SOCKS5_PROXY_HOST = '127.0.0.1' 
SOCKS5_PROXY_PORT = 1086
default_socket = socket.socket
socks.set_default_proxy(socks.SOCKS5, SOCKS5_PROXY_HOST, SOCKS5_PROXY_PORT) 
socket.socket = socks.socksocket
# 在 mnist 上训练随机森林，然后画出每个像素的重要性
# 获取 mnist
from six.moves import urllib
from sklearn.datasets import fetch_mldata
try:
    mnist = fetch_mldata('MNIST original')
except urllib.error.HTTPError as ex:
    print("Could not download MNIST data from mldata.org, trying alternative...")
    from scipy.io import loadmat
    mnist_alternative_url = "https://github.com/amplab/datascience-sp14/raw/master/lab7/mldata/mnist-original.mat"
    mnist_path = "./mnist-original.mat"
    response = urllib.request.urlopen(mnist_alternative_url)
    with open(mnist_path, "wb") as f:
        content = response.read()
        f.write(content)
    mnist_raw = loadmat(mnist_path)
    mnist = {
        "data": mnist_raw["data"].T,
        "target": mnist_raw["label"][0],
        "COL_NAMES": ["label", "data"],
        "DESCR": "mldata.org dataset: mnist-original",
    }
    print("Success!")

Could not download MNIST data from mldata.org, trying alternative...
Success!
In [88]:
from sklearn.model_selection import train_test_split

X = mnist["data"]
y = mnist["target"]

X_train, X_test, y_train, y_test = train_test_split(X, y)
In [89]:
pca=PCA()
pca.fit(X_train)
cumsum=np.cumsum(pca.explained_variance_ratio_)
d=np.argmax(cumsum>=0.95)+1
d
Out[89]:
154

当然，我们也可以设置 n_components = d 来设置 PCA ，有一个更好的选择，不指定你想要保留的主成分个数，而是将 n_components 设置为 0.0 到 1.0 之间的浮点数，表明希望保留的方差比率。

In [90]:
pca = PCA(n_components=0.95)
X_reduced = pca.fit_transform(X_train)
pca.n_components_
Out[90]:
154

另一种方法是用画出方差解释率的函数，曲线中有一个肘部，方差解释器停止快速增长。

In [91]:
pca = PCA(n_components=400)
pca.fit(X_train)
Out[91]:
PCA(copy=True, iterated_power='auto', n_components=400, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
In [107]:
cumsum = np.cumsum(pca.explained_variance_ratio_)
plt.plot(range(400), cumsum, 'k-')
plt.axis([4, 400, 0, 1])
plt.plot(154, cumsum[154], 'ko')
plt.plot([154 ,154], [0, cumsum[154]], 'k--')
plt.plot([0, 400], [cumsum[154], cumsum[154]], 'k--')
plt.xlabel("$Dimensions$", fontsize=18)
plt.ylabel("$Explained Variance$", fontsize=18)
Out[107]:
Text(0,0.5,'$Explained Variance$')

PCA 压缩

我们将 MNIST 数据集保留 95% 后，784个特征剩下150个，保留大部分方差的同时，数据集改不到原来的 20% 。这是一个合理的压缩比例。

通过应用 PCA 投影的你变化，也可以将缩小的数据集压缩回 784 维。当然这并不会返回最原始的数据，因为投影丢失了一些信息，这些信息在 5% 的方差内。

原始数据和重构数据之间的均方距离被称为重构误差 reconstruction error 。

In [113]:
pca = PCA(n_components=154)
X_reduced = pca.fit_transform(X_train)
X_recovered = pca.inverse_transform(X_reduced)
In [114]:
def plot_digits(instances, images_per_row=5, **options):
    size = 28
    images_per_row = min(len(instances), images_per_row)
    images = [instance.reshape(size,size) for instance in instances]
    n_rows = (len(instances) - 1) // images_per_row + 1
    row_images = []
    n_empty = n_rows * images_per_row - len(instances)
    images.append(np.zeros((size, size * n_empty)))
    for row in range(n_rows):
        rimages = images[row * images_per_row : (row + 1) * images_per_row]
        row_images.append(np.concatenate(rimages, axis=1))
    image = np.concatenate(row_images, axis=0)
    plt.imshow(image, cmap = matplotlib.cm.binary, **options)
    plt.axis("off")
In [115]:
plt.figure(figsize=(7, 4))
plt.subplot(121)
# 每隔 2100 张选出一张照片
plot_digits(X_train[::2100])
plt.title("Original", fontsize=16)
plt.subplot(122)
plot_digits(X_recovered[::2100])
plt.title("Compressed", fontsize=16)
Out[115]:
Text(0.5,1,'Compressed')

PCA 逆变换，回到原来的数据维度：

$$X_{recovered}=X_{d-proj}\cdot W_{d}^{T}$$

增量 PCA (Incremental PCA)

之前我们实现的方式是处理整个训练集以便 SVD 算法运行。幸运的是，我们已经开发了增量 PCA (IPCA) 算法：我们可以将训练集分批，并以此只对一个批量使用 IPCA 算法。这对大型训练集非常有用，或者在线学习。

In [123]:
from sklearn.decomposition import IncrementalPCA

n_batches = 100
inc_pca = IncrementalPCA(n_components=154)
for X_batch in np.array_split(X_train, n_batches):
    print('.', end='')
    # 最小批次调用，而不是整体调用 fit
    inc_pca.partial_fit(X_batch)
X_reduced = inc_pca.transform(X_train)

....................................................................................................
In [124]:
X_recovered_inc_pca = inc_pca.inverse_transform(X_reduced)
plt.figure(figsize=(7, 4))
plt.subplot(121)
plot_digits(X_train[::2100])
plt.subplot(122)
plot_digits(X_recovered_inc_pca[::2100])
plt.tight_layout
Out[124]:
<function matplotlib.pyplot.tight_layout(pad=1.08, h_pad=None, w_pad=None, rect=None)>
In [126]:
# np.memmap 允许我们操作存储在磁盘上的二进制文件中的大型数组。
m, n = X_train.shape
X_mm = np.memmap('mnist-original.mat', dtype='float32', mode='write', shape=(m, n))
batch_size = m // n_batches
inc_pca = IncrementalPCA(n_components=154, batch_size=batch_size)
inc_pca.fit(X_mm)

/Users/ronnie/.pyenv/versions/3.6.4rc1/envs/chatbot/lib/python3.6/site-packages/sklearn/decomposition/incremental_pca.py:267: RuntimeWarning: invalid value encountered in true_divide
  explained_variance_ratio = S ** 2 / np.sum(col_var * n_total_samples)
Out[126]:
IncrementalPCA(batch_size=525, copy=True, n_components=154, whiten=False)

随机 PCA (Randomized PCA)

sklearn 提供了另一种执行 PCA 的选择，称为随机 PCA 。这是一种随机算法，可以快速找到前 d 个主成分的近似值，它的复杂度是 O(m.d^2)+O(d^3) ，当 d 远小于 n 时，它比之前的算法快得多。

In [127]:
rnd_pca = PCA(n_components=154, svd_solver='randomized')
X_reduced = rnd_pca.fit_transform(X_train)

核 PCA (Kernel PCA)

核技巧，将实例隐式映射到高维空间的技术，让 svm 可以应用于非线性分类和回归。

同样的技巧可以应用于 PCA ，从而可以执行复杂的非线性投影来降低维度。这就是所谓的核 PCA 。它通常能够很好地保留投影后的簇，甚至可以展开分布近似于扭曲流形的数据集。

In [128]:
# 这段代码使用 sklearn 的 KernelPCA 类来执行带有 RBF 核的 kPCA 。
X, t = make_swiss_roll(n_samples=1000, noise=0.2, random_state=42)
from sklearn.decomposition import KernelPCA

rbf_pca = KernelPCA(n_components=2, kernel='rbf', gamma=0.04)
X_reduced = rbf_pca.fit_transform(X)
In [129]:
# 我们使用了三种核将瑞士卷降到2维
# 线性核
lin_pca = KernelPCA(n_components=2, kernel='linear', fit_inverse_transform=True)
# rbf 核
rbf_pca = KernelPCA(n_components = 2, kernel="rbf", gamma=0.0433, fit_inverse_transform=True)
# sigmoid 核 (Logistic 核)
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$")):
    # 使用 pca 降维后的输入
    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)

plt.show()

选择一种核并调整超参数

由于 kPCA 是无监督学习，因此没有明显的性能指标可以帮助我们选择最佳的核方法和超参数值。

但是降维通常是监督学习任务的准备步骤，我们可以简单实用网格搜索来选择可以让该任务最佳表现的核方法和超参数。

下面我们创建一个两步流水线，首先实用 kPCA 将维度降至 2 维，然后应用 Logistic 回归进行分类。然后使用 Grid SearchCV 为 kPCA 找到最佳的核和 gamma 值，以便在最后获得最佳的分类准确性。

In [131]:
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()),
])

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_)

{'kpca__gamma': 0.043333333333333335, 'kpca__kernel': 'rbf'}

/Users/ronnie/.pyenv/versions/3.6.4rc1/envs/chatbot/lib/python3.6/site-packages/sklearn/model_selection/_search.py:737: DeprecationWarning: The default of the `iid` parameter will change from True to False in version 0.22 and will be removed in 0.24. This will change numeric results when test-set sizes are unequal.
  DeprecationWarning)

bestparams 中可以查看模型最好的核和超参数

例一种完全为非监督的方法，是选择最低重建误差的核和超参数。

我们在对原始数据进行 kPCA 操作的时候，相当于将原始数据映射到无穷维，然后投影到 2D 平面上，此时我们反向的重构点将位于特征空间（无限维）中，而不是原始空间。我们无法找出重建点。

我们可以在原始空间中找到一个贴近重建点的点，叫重建前图像。

这需要设置 fit_inverse_transform = True 。

In [132]:
rbf_pca = KernelPCA(n_components=2, kernel='rbf', gamma=0.0433,
                   fit_inverse_transform=True)
X_reduced = rbf_pca.fit_transform(X)
# 回到从前的数据
X_preimage = rbf_pca.inverse_transform(X_reduced)
from sklearn.metrics import mean_squared_error
# 计算从前的数据和原始数据之间的方差
mean_squared_error(X, X_preimage)
Out[132]:
32.78630879576613

LLE

局部线性嵌入 Locally Linear Embedding 是另一种非常有效的非线性降维 NLDR 方法。

这是一种流形学习技术，不依赖前面的投影。

LLE 首先测量每个训练实例与其最近邻之间的线性关系，然后寻找能最好地保留这些局部关系的训练集的低维表示。

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

下面使用 sklearn 的 LocallyLinearEmbedding 类来展开瑞士卷。得到的二维数据集将展示。

In [133]:
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, random_state=42)
# 经过 LLE 后，将得到展开的二维矩阵
X_reduced = lle.fit_transform(X)
In [134]:
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)

plt.show()

$$\widehat{w}=argmin_{w}\sum _{i=1}^{m}(x^{(i)}-\sum _{j=1}^{m}w_{i,j}x^{(j)})^{2}$$

$$subject\quad to\left\{\begin{matrix} w_{i,j}=0\\ \sum_{j=1}^{m}w_{i,j}=1 \end{matrix}\right.$$

上面是 LLE 的第一步，寻找 k 近邻，如果 xj 不是 xi 的 k 近邻，那么距离权重设置为 0。

经过这一步，权重矩阵 w 包含对实例的线性关系。现在需要将训练实例投影到一个 d 维空间去。同时尽可能保留这些局部关系。

如果 Z^{i} 是 x^{i} 在这个 d 维空间的图像，那么我们想要 z 和 \sum wx 之间的平方距离尽可能小。

$$\widehat{z}=argmin_{z}\sum _{i=1}^{m}(z^{(i)}-\sum _{j=1}^{m}\widehat{w}_{i,j}z^{(j)})^{2}$$

上面是算法的第二步，和保持实例固定不变，找最佳权重相反；我们保持权重不变，在低维空间中找到实例空间的最佳位置。 z 是包含所有 z^{i} 的矩阵。

时间复杂度上，查找 k 近邻为 O(mlog(m)nlog(k)) 。优化权重为 O(mnk^3) 。建立低维表示是 O(dm^2) 。

最后一项 m^2 导致这个算法在处理大数据集的时候表现较差。

其他降维

In [135]:
# mds
from sklearn.manifold import MDS

mds = MDS(n_components=2, random_state=42)
X_reduced_mds = mds.fit_transform(X)
In [136]:
# isomap
from sklearn.manifold import Isomap

isomap = Isomap(n_components=2)
X_reduced_isomap = isomap.fit_transform(X)
In [137]:
# t-sne
from sklearn.manifold import TSNE

tsne = TSNE(n_components=2, random_state=42)
X_reduced_tsne = tsne.fit_transform(X)
In [ ]:
# lda
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

lda = LinearDiscriminantAnalysis(n_components=2)
X_mnist = mnist["data"]
y_mnist = mnist["target"]
lda.fit(X_mnist, y_mnist)
X_reduced_lda = lda.transform(X_mnist)
In [138]:
titles = ["MDS", "Isomap", "t-SNE"]

plt.figure(figsize=(11,4))

for subplot, title, X_reduced in zip((131, 132, 133), titles,
                                     (X_reduced_mds, X_reduced_isomap, X_reduced_tsne)):
    plt.subplot(subplot)
    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)

plt.show()