機器學習筆記(14)——sklearn降維方法舉例(RandomProjection,TSVD,t-SNE)

機器學習筆記(14)——sklearn降維方法舉例(RandomProjection,TSVD,t-SNE)
目錄

sklearn降維方法舉例

以datasets.digits資料為例匯入相關包

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import time
from sklearn.datasets import load_digits

大樣本資料的視覺化是一個相對比較麻煩的事情,

一般情況下我們都要用到降維的方法先處理特徵。我們找一個例子來看看,可以怎麼做,比如我們資料集取經典的『手寫數字集』

Each datapoint is a 8×8 image of a digit.

Classes->10

Samples per class->180

Samples total->1797

Dimensionality->64

Features integers->0-16

匯入資料,0,1,2,3,4,5,6這7個數字的手寫影象

digits = load_digits(n_class=7)#integer
X = digits.data
y = digits.target
n_samples,n_features = X.shape

輸出影象

from matplotlib import offsetbox
def plot_embedding(X, title=None):
x_min, x_max = np.min(X, 0), np.max(X, 0)
X = (X - x_min) / (x_max - x_min)#正則化
print X
plt.figure(figsize=(10, 10))
ax = plt.subplot(111)
for i in range(X.shape[0]):
plt.text(X[i, 0], X[i, 1], str(digits.target[i]),
color=plt.cm.Set1(y[i] / 10.),
fontdict={'weight': 'bold', 'size': 12})
#列印彩色字型
if hasattr(offsetbox, 'AnnotationBbox'):
# only print thumbnails with matplotlib > 1.0
shown_images = np.array([[1., 1.]])  # just something big
for i in range(digits.data.shape[0]):
dist = np.sum((X[i] - shown_images) ** 2, 1)
if np.min(dist) < 4e-3:
# don't show points that are too close
continue
shown_images = np.r_[shown_images, [X[i]]]
imagebox = offsetbox.AnnotationBbox(
offsetbox.OffsetImage(digits.images[i], cmap=plt.cm.gray_r),
X[i])
ax.add_artist(imagebox)#輸出圖上輸出圖片
plt.xticks([]), plt.yticks([])
if title is not None:
plt.title(title)

1.隨機投射(Random-Projection)

首先我們來看一個問題, 如果你手頭有一組資料X∈Rn” role=”presentation”>X∈RnX∈RnX \in R^n, 它的維數太高, 從而不得不進行降維至Rk” role=”presentation”>RkRkR^k, 你會怎麼辦?

相信不少人反射性地求它的x′x” role=”presentation”>x′xx′xx’x,然後轉化為延最大方差展開的問題。 當然這個方法需要求解固有值和固有向量。對於測試用的小資料一切似乎都很好, 但如果我給你的是50G大小csv格式的Twitter發言呢?如果需要對全部單詞降維, 恐怕半個月吧。 因此我們需要個更加快速的方法。

那麼現在再來想, 既然要快, 那我們就用個最簡單的方法: 隨機在高維空間裡選幾個單位向量ei” role=”presentation”>eieie_i, 但注意, 這裡我們僅僅要求是單位向量, 並沒有要求它們之間必須正交&lt;ei,ej&gt;=0″ role=”presentation”><ei,ej>=0<ei,ej>=0=0, 因此你可以隨便選。 最後, 我們把高維資料投影到選擇的這一組基上。也就完成了降維。

隨便選的軸如何能保證降維效果? 但它是有數學依據的, 叫做Johnson-Lindenstrauss Lemma。這個定理保證了任何降維方法的精度上下限。

Johnson-Lindenstrauss Lemma

假設我們有資料x∈Rn” role=”presentation”>x∈Rnx∈Rnx \in R^n, 而我們通過一種方法f(x)” role=”presentation”>f(x)f(x)f(x)將其降維成y∈Rk” role=”presentation”>y∈Rky∈Rky \in R^k, 那麼, 將為前後任意兩點a,b” role=”presentation”>a,ba,ba,b之間的距離有不等式保證:

(1−ϵ)‖a−b‖2≤‖f(a)−f(b)‖2≤(1 ϵ)‖a−b‖2″ role=”presentation”>(1−ϵ)∥a−b∥2≤∥f(a)−f(b)∥2≤(1 ϵ)∥a−b∥2(1−ϵ)‖a−b‖2≤‖f(a)−f(b)‖2≤(1 ϵ)‖a−b‖2(1-\epsilon )\| a-b \| ^2 \le \| f(a)-f(b) \| ^2 \le (1 \epsilon ) \| a-b \|^2

對於隨機對映來說, 只要注意到下面幾點就行了:

1.不定式中的精度僅僅受制於維數、資料大小等因素, 與將為方法無關。

2.在維數差不是很大時, 總可以保證一個相對較高的精度, 不論用什麼方法。

3.到這裡一切就很明顯了, 既然精度有上界, 那我們也就不必擔心軸的選取,那麼最簡單的方法自然就是隨機挑選了, 這也就產生的Random Projection這一方法。

Random Projection

簡單來說,Random Projection流程就是

選擇影射矩陣R∈RK×N” role=”presentation”>R∈RK×NR∈RK×NR \in \mathcal{R}^{K \times N}。

用隨機數填充影射矩陣。 可以選擇uniform或者gaussian。

歸一化每一個新的軸(影射矩陣中的每一行)。

對資料降維y=RX” role=”presentation”>y=RXy=RXy=R\mathcal{X}。

上個小節裡面的JL-lemma保證的降維的效果不會太差。

從直觀上來看看。

首先假設我們有資料X={x|fi(θ) σ}” role=”presentation”>X={x|fi(θ) σ}X={x|fi(θ) σ}\mathcal{X}=\{ x | f^{i}(\theta) \sigma \}, 其中θ” role=”presentation”>θθ\theta是一組引數, σ” role=”presentation”>σσ\sigma則是高斯噪聲。 回想PCA方法, 我們很自然的認為所有的特徵都是正交的, 我們既沒有考慮它是不是有多箇中心, 也沒有考慮是不是有特殊結構, 然而對於實際資料, 很多時候並不是這樣。 比如我們把取f(θ)=Si(θi)N(A,σi)∈R3″ role=”presentation”>f(θ)=Si(θi)N(A,σi)∈R3f(θ)=Si(θi)N(A,σi)∈R3f(\theta)=S^{i}(\theta _i)N(A, \sigma^i) \in R^3, 其中Si∈SO(3)” role=”presentation”>Si∈SO(3)Si∈SO(3)S^{i} \in SO(3), 這樣我們得到的資料可能會像×” role=”presentation”>××\times或者∗” role=”presentation”>∗∗*一樣, 顯然用PCA並不能得到好的結果。

在這種情況下, 考慮在非常高維度的空間, 可以想象random projection這種撞大運的方法總會選到一些超平面能夠接近資料的特徵, 同時也能夠想象如果目標維數過低, 那麼命中的概率就會大大降低。

實現:

from sklearn import random_projection
rp = random_projection.SparseRandomProjection(n_components=2, random_state=42)
#隨機對映演算法
start_time = time.time()
X_projected = rp.fit_transform(X)#隨機投射
plot_embedding(X_projected, "Random Projection of the digits (time: %.3fs)" % (time.time() - start_time))
plt.show()

這裡寫圖片描述

時間很短但效果一般

官方描述:

def sparse_random_matrix(n_components, n_features, density='auto',
random_state=None):
"""Generalized Achlioptas random sparse matrix for random projection
Setting density to 1 / 3 will yield the original matrix by Dimitris
Achlioptas while setting a lower value will yield the generalization
by Ping Li et al.
If we note :math:`s = 1 / density`, the components of the random matrix are
drawn from:
- -sqrt(s) / sqrt(n_components)   with probability 1 / 2s
-  0                              with probability 1 - 1 / s
-  sqrt(s) / sqrt(n_components)   with probability 1 / 2s
Read more in the :ref:`User Guide <sparse_random_matrix>`.
Parameters
----------
n_components : int,
Dimensionality of the target projection space.
n_features : int,
Dimensionality of the original source space.
density : float in range ]0, 1] or 'auto', optional (default='auto')
Ratio of non-zero component in the random projection matrix.
If density = 'auto', the value is set to the minimum density
as recommended by Ping Li et al.: 1 / sqrt(n_features).
Use density = 1 / 3.0 if you want to reproduce the results from
Achlioptas, 2001.
random_state : integer, RandomState instance or None (default=None)
Control the pseudo random number generator used to generate the
matrix at fit time.
Returns
-------
components: numpy array or CSR matrix with shape [n_components, n_features]
The generated Gaussian random matrix.
"""
_check_input_size(n_components, n_features)
density = _check_density(density, n_features)
rng = check_random_state(random_state)
if density == 1:
# skip index generation if totally dense
components = rng.binomial(1, 0.5, (n_components, n_features)) * 2 - 1
return 1 / np.sqrt(n_components) * components
else:
# Generate location of non zero elements
indices = []
offset = 0
indptr = [offset]
for i in xrange(n_components):
# find the indices of the non-zero components for row i
n_nonzero_i = rng.binomial(n_features, density)
indices_i = sample_without_replacement(n_features, n_nonzero_i,
random_state=rng)
indices.append(indices_i)
offset  = n_nonzero_i
indptr.append(offset)
indices = np.concatenate(indices)
# Among non zero components the probability of the sign is 50%/50%
data = rng.binomial(1, 0.5, size=np.size(indices)) * 2 - 1
# build the CSR structure by concatenating the rows
components = sp.csr_matrix((data, indices, indptr),
shape=(n_components, n_features))
return np.sqrt(1 / density) / np.sqrt(n_components) * components

2.TSVD(截斷奇異值分解, TruncatedSVD,PCA的一種實現)

截斷奇異值分解(Truncated singular value decomposition,TSVD)是一種矩陣因式分解(factorization)技術,將矩陣 M 分解成 U , Σ 和 V 。它與PCA很像,只是SVD分解是在資料矩陣上進行,而PCA是在資料的協方差矩陣上進行。通常,SVD用於發現矩陣的主成份。對於病態矩陣,目前主要的處理辦法有預調節矩陣方法、區域分解法、正則化方法等,截斷奇異值分解技術TSVD就是一種正則化方法,它犧牲部分精度換去解的穩定性,使得結果具有更高的泛化能力。

對於原始資料矩陣A(N*M) ,N代表樣本個數,M代表維度,對其進行SVD分解:

A=UΣVT,Σ=diag(δ1,δ2,…,δn)=[δ100..00δ20..000δ3..0…………000..δn]” role=”presentation”>A=UΣVT,Σ=diag(δ1,δ2,…,δn)=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢δ100..00δ20..000δ3..0…………000..δn⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥A=UΣVT,Σ=diag(δ1,δ2,…,δn)=[δ100..00δ20..000δ3..0…………000..δn]

A = U\Sigma V^T, \Sigma=diag(\delta_1,\delta_2,…,\delta_n)=\left[
\begin{matrix}
\delta_1&0&0&.&.&0\\
0&\delta_2&0&.&.&0\\
0&0&\delta_3&.&.&0\\
.&.&.&.&.&.\\
.&.&.&.&.&.\\
0&0&0&.&.&\delta_n
\end{matrix}
\right]

上式中的δ” role=”presentation”>δδ\delta就是資料的奇異值,且δ1″ role=”presentation”>δ1δ1\delta_1>δ2″ role=”presentation”>δ2δ2\delta_2>δ3″ role=”presentation”>δ3δ3\delta_3…,通常如果A非常病態,delta的後面就越趨向於0,δ1δn” role=”presentation”>δ1δnδ1δn\frac{\delta_1}{\delta_n}就是資料的病態程度,越大說明病態程度越高,無用的特徵越多,通常會擷取前p個最大的奇異值,相應的U擷取前p列,V擷取前p列,這樣A依然是N*M的矩陣,用這樣的計算出來的A代替原始的A,就會比原始的A更穩定。

TSVD與一般SVD不同的是它可以產生一個指定維度的分解矩陣。例如,有一個 n×n 矩陣,通過SVD分解後仍然是一個 n×n 矩陣,而TSVD可以生成指定維度的矩陣。這樣就可以實現降維了。

演示一下TSVD工作原理

In[1]:
import numpy as np
from scipy.linalg import svd
D = np.array([[1, 2], [1, 3], [1, 4]])
D
Out[1]:
array([[1, 2],
[1, 3],
[1, 4]])
In[2]:
U, S, V = svd(D, full_matrices=False)
#我們可以根據SVD的定義,用 UU , SS 和 VV 還原矩陣 DD 
U.shape, S.shape, V.shape
Out[2]:
((3, 2), (2, 1), (2, 2))
In[3]:
np.dot(U.dot(np.diag(S)), V)#還原
Out[3]:
array([[ 1.,  2.],
[ 1.,  3.],
[ 1.,  4.]])

TruncatedSVD返回的矩陣是 U 和 S 的點積。如果我們想模擬TSVD,我們就去掉最新奇異值和對於 U 的列向量。例如,我們想要一個主成份,可以這樣:

In[4]:
new_S = S[0]
new_U = U[:, 0]
new_U.dot(new_S)
Out[4]:
array([-2.20719466, -3.16170819, -4.11622173])

TruncatedSVD有個“陷阱”。隨著隨機數生成器狀態的變化,TruncatedSVD連續地擬合會改變輸出的符合。為了避免這個問題,建議只用TruncatedSVD擬合一次,然後用其他變換。這正是管線命令的另一個用處。

實現:

from sklearn import decomposition
X_pca = decomposition.TruncatedSVD(n_components=2).fit_transform(X)
start_time = time.time()
plot_embedding(X_pca,"Principal Components projection of the digits (time: %.3fs)" % (time.time() - start_time))
plt.show()

這裡寫圖片描述

效果稍微好一些

3.t-SNE(t-分佈鄰域嵌入演算法)

流形學習方法(Manifold Learning),簡稱流形學習,可以將流形學習方法分為線性的和非線性的兩種,線性的流形學習方法如我們熟知的主成份分析(PCA),非線性的流形學習方法如等距對映(Isomap)、拉普拉斯特徵對映(Laplacian eigenmaps,LE)、區域性線性嵌入(Locally-linear embedding,LLE)。

這裡寫圖片描述
t-SNE詳細介紹:http://lvdmaaten.github.io/tsne/

from sklearn import manifold
#降維
tsne = manifold.TSNE(n_components=2, init='pca', random_state=0)
start_time = time.time()
X_tsne = tsne.fit_transform(X)
#繪圖
plot_embedding(X_tsne,
"t-SNE embedding of the digits (time: %.3fs)" % (time.time() - start_time))
plt.show()
#這個非線性變換降維過後,僅僅2維的特徵,就可以將原始資料的不同類別,在平面上很好地劃分開。
#不過t-SNE也有它的缺點,一般說來,相對於線性變換的降維,它需要更多的計算時間。

這裡寫圖片描述