sklearn 降维方法概述

降维方法

• 现实中的许多数据都是稀疏的(sparse)，高维数据处理的时间和空间复杂度都十分大，因此需要对数据进行降维
• 对数据进行降维，会在一定程度上降低数据的精度，同时也会增加机器学习模型处理流程的复杂度。

主要的降维方法

映射(Projection)

• 现实中的许多数据的特征都是相关的，或者特征为常数，可以利用映射的方法将高维数据映射到低维

流行学习(Manifold Learning)

• 流行学习依靠流行假设：即自然界的大多数高维数据集与一个低维的流行很接近
• 通常伴随流程假设的还有另外一个假设：高维数据的建模任务在转换为低维空间之后，任务的难度会下降很多。但是现实的情况往往是：的确是降低了模型秋季的难度，但是数据更加难以划分

PCA(Principal Component Analysis)

• 主成分分析是目前最常用的降维方法，它是首先找到数据的超平面，然后将数据投影到这个超平面上
• 在进行超平面投影的时候，需要注意的是：要尽可能地保留数据的方差特征
• 主成分首先找到训练集上方差最大的投影方向，然后再找与第一个方向正交的第二个主方向，数据沿着它的投影方向投影之后的方差是第二大的；以此类推，可以找出数据集的若干个主成分
• 我们可以使用SVD方法，很方便地找到数据集的主成分

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np


angle = np.pi / 5
stretch = 5
m = 200

np.random.seed(3)
X = np.random.randn(m, 2) / 10
X = X.dot(np.array([[stretch, 0],[0, 1]])) # stretch
X = X.dot([[np.cos(angle), np.sin(angle)], [-np.sin(angle), np.cos(angle)]]) # rotate

u1 = np.array([np.cos(angle), np.sin(angle)])
u2 = np.array([np.cos(angle - 2 * np.pi/6), np.sin(angle - 2 * np.pi/6)])
u3 = np.array([np.cos(angle - np.pi/2), np.sin(angle - np.pi/2)])

X_proj1 = X.dot(u1.reshape(-1, 1))
X_proj2 = X.dot(u2.reshape(-1, 1))
X_proj3 = X.dot(u3.reshape(-1, 1))

plt.figure(figsize=(6,3))
plt.subplot2grid((3,2), (0, 0), rowspan=3)
plt.plot([-1.4, 1.4], [-1.4*u1[1]/u1[0], 1.4*u1[1]/u1[0]], "k-", linewidth=1)
plt.plot([-1.4, 1.4], [-1.4*u2[1]/u2[0], 1.4*u2[1]/u2[0]], "k--", linewidth=1)
plt.plot([-1.4, 1.4], [-1.4*u3[1]/u3[0], 1.4*u3[1]/u3[0]], "k:", linewidth=2)
plt.plot(X[:, 0], X[:, 1], "bo", alpha=0.5)
plt.axis([-1.4, 1.4, -1.4, 1.4])
plt.arrow(0, 0, u1[0], u1[1], head_width=0.1, linewidth=5, length_includes_head=True, head_length=0.1, fc='k', ec='k')
plt.arrow(0, 0, u3[0], u3[1], head_width=0.1, linewidth=5, length_includes_head=True, head_length=0.1, fc='k', ec='k')
plt.text(u1[0] + 0.1, u1[1] - 0.05, r"$\mathbf{c_1}$", fontsize=22)
plt.text(u3[0] + 0.1, u3[1], r"$\mathbf{c_2}$", fontsize=22)
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$x_2$", fontsize=18, rotation=0)
plt.grid(True)

plt.subplot2grid((3,2), (0, 1))
plt.plot([-2, 2], [0, 0], "k-", linewidth=1)
plt.plot(X_proj1[:, 0], np.zeros(m), "bo", alpha=0.3)
plt.gca().get_yaxis().set_ticks([])
plt.gca().get_xaxis().set_ticklabels([])
plt.axis([-2, 2, -1, 1])
plt.grid(True)

plt.subplot2grid((3,2), (1, 1))
plt.plot([-2, 2], [0, 0], "k--", linewidth=1)
plt.plot(X_proj2[:, 0], np.zeros(m), "bo", alpha=0.3)
plt.gca().get_yaxis().set_ticks([])
plt.gca().get_xaxis().set_ticklabels([])
plt.axis([-2, 2, -1, 1])
plt.grid(True)

plt.subplot2grid((3,2), (2, 1))
plt.plot([-2, 2], [0, 0], "k:", linewidth=2)
plt.plot(X_proj3[:, 0], np.zeros(m), "bo", alpha=0.3)
plt.gca().get_yaxis().set_ticks([])
plt.axis([-2, 2, -1, 1])
plt.xlabel("$z_1$", fontsize=18)
plt.grid(True)

plt.show()

X_centered = X - X.mean( axis=0 )
U, s, V = np.linalg.svd( X_centered )
c1 = V.T[:,0]
c2 = V.T[:,1]
print( c1, c2 )

[-0.79644131 -0.60471583] [-0.60471583  0.79644131]

• 可以通过找到数据的前d个主成分，将高维数据转换为d维数据，只需要将原始数据与几个主成分相乘即可

W2 = V.T[:,:2]
X2D = X_centered.dot( W2 )
plt.figure(figsize=(3,3))
plt.plot( np.zeros((X2D.shape[0], 1)), X2D[:,0],'ro', linewidth=4, alpha=0.3 )
plt.plot( X2D[:,1], np.zeros((X2D.shape[0], 1)), 'bo', linewidth=4, alpha=0.3 )
plt.axis("equal")
plt.show()

• 直接使用sklearn也可以PCA求解
• PCA模块可以输出每个主成分的方差权重(可以理解为重要性权重)
• PCA中，需要选择正确的主成分个数，可以通过计算前k个主成分是否到某个百分比来确定，也可以直接在PCA传参中确定主成分需要占得比例

from sklearn.decomposition import PCA
pca = PCA( n_components=2 )
X2D = pca.fit_transform( X )
print( pca.explained_variance_ratio_ )
cumsum = np.cumsum( pca.explained_variance_ratio_ )
d = np.argmax( cumsum >= 0.95 ) + 1 
print( d )
pca = PCA( n_components=0.95 )
X_reduced = pca.fit_transform( X )
print( X_reduced.shape )

[ 0.95369864  0.04630136]
1
(200, 1)

from sklearn.datasets import fetch_mldata
mnist = fetch_mldata('MNIST original')
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)
pca = PCA()
pca.fit(X_train)
cumsum = np.cumsum(pca.explained_variance_ratio_)
d = np.argmax(cumsum >= 0.95) + 1
print( d )

154

import matplotlib
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")
    return

pca = PCA(n_components = 154)
X_reduced = pca.fit_transform(X_train)
X_recovered = pca.inverse_transform(X_reduced)
plt.figure(figsize=(6, 3))
plt.subplot(121)
plot_digits(X_train[::2100])
plt.title("Original", fontsize=10)
plt.subplot(122)
plot_digits(X_recovered[::2100])
plt.title("Compressed", fontsize=10)
plt.show()

增量式PCA

• PCA的一个缺点就是它在运行的时候需要利用所有的训练集。针对这个缺点，开发出增量式PCA，将训练集分成很多个batch，每次使用训练集中的一部分
• 也可以使用numpy中的memmap类，可以对储存在二进制文件中的大型矩阵进行操作，只有当用到矩阵的该部分时，才会将其加载到内存

Randomized PCA

• cklearn也有一种随机PCA的方法，他的训练时间复杂度是$O(m \cdot d^2)+O(d^3)$，而不是之前介绍的PCA的$O(m \cdot n^2)+O(d^3)$，因此当原始数据维数很高的时候，这种方法可以大大加快PCA的速度

from sklearn.decomposition import IncrementalPCA
from time import clock
n_batches = 100
inc_pca = IncrementalPCA( n_components=154 )
start = clock()
for X_batch in np.array_split( X_train, n_batches ):
    inc_pca.partial_fit( X_batch )
finish = clock()
print( "time (s) : ", (finish-start) )
X_train_reduced = inc_pca.transform( X_train )

time (s) :  37.81862586160648

rnd_pca = PCA(n_components=154, svd_solver="randomized")
start = clock()
X_reduced = rnd_pca.fit_transform(X_train)
finish = clock()
print( "time (s) ： ", (finish-start) )

time (s) ：  9.015620531600973

Kernel PCA

• kernel技巧就是：利用一些数学方面的技巧，将数据映射到高维空间，使得SVM可以处理这种非线性的分类或者回归任务。同样的技巧也可以被用于PCA
• kPCA可以保护一些投影后的样本簇，同时将数据展开

from sklearn.decomposition import KernelPCA
from sklearn.datasets import make_swiss_roll

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

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=(9, 3))
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.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()

• kernel选择和超参数调节
• 可以采用网格搜索的方法，搜索当前模型超参数的最优解

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'}

from sklearn.metrics import mean_squared_error
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 )
print( mean_squared_error(X, X_preimage) )

32.7863087958

LLE(Locally Linear Embedding )

• LLE不依赖之前所采用的投影方法，它是找出与它最近邻的一些样本呢，然后尝试用低维特征去表示这些关系，同时使得局部的样本关系可以使得被保存。它属于流行学习方法，主要思想是用近邻的点去线性表征该样本，在低维空间中，这种线性关系也可以基本保证不变。
• LLE算法的主要流程
  
  • 寻找每个样本点的k个最近邻点
  • 由每个样本点 的近邻点计算出该样本点的局部重建权值矩阵
  • 由该样本点的局部重建权值矩阵和其近邻点计算出该样本点的输出值
• LLE可以对一些噪声不大的难以分割的数据进行很好的展开，但是他也有局限性：数据集不能是闭合流形，不能是稀疏的数据集，不能是分布不均匀的数据集等等。
• 参考文献：https://ask.hellobi.com/blog/guodongwei1991/10912

from sklearn.manifold import LocallyLinearEmbedding

lle = LocallyLinearEmbedding(n_components=2, n_neighbors=10, random_state=42)
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)

plt.show()

其他的一些降维方法

• MDS(Multidimensional Scaling )
• Isomap
• t-SNE(t-Distributed Stochastic Neighbor Embedding)
• LDA(Linear Discriminant Analysis)