降維方法

• 現實中的許多資料都是稀疏的(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⋅d2)+O(d3)，而不是之前介紹的PCA的O(m⋅n2)+O(d3)，因此當原始資料維數很高的時候，這種方法可以大大加快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可以對一些噪聲不大的難以分割的資料進行很好的展開，但是他也有侷限性：資料集不能是閉合流形，不能是稀疏的資料集，不能是分佈不均勻的資料集等等。

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)