機(jī)器學(xué)習(xí)基礎(chǔ)與實(shí)踐:PCA的目的

寫在前面：本來(lái)想?yún)⒖糚RML來(lái)寫，但是發(fā)現(xiàn)里面涉及到比較多的數(shù)學(xué)知識(shí)，寫出來(lái)可能不好理解，我決定還是用最通俗的方法解釋PCA，并舉一個(gè)實(shí)例一步步計(jì)算，然后再進(jìn)行數(shù)學(xué)推導(dǎo)，最后再介紹一些變種以及相應(yīng)的程序。

正文：

在數(shù)據(jù)處理中，經(jīng)常會(huì)遇到特征維度比樣本數(shù)量多得多的情況，如果拿到實(shí)際工程中去跑，效果不一定好。一是因?yàn)槿哂嗟奶卣鲿?huì)帶來(lái)一些噪音，影響計(jì)算的結(jié)果；二是因?yàn)闊o(wú)關(guān)的特征會(huì)加大計(jì)算量，耗費(fèi)時(shí)間和資源。所以我們通常會(huì)對(duì)數(shù)據(jù)重新變換一下，再跑模型。數(shù)據(jù)變換的目的不僅僅是降維，還可以消除特征之間的相關(guān)性，并發(fā)現(xiàn)一些潛在的特征變量。

一、PCA的目的

PCA是一種在盡可能減少信息損失的情況下找到某種方式降低數(shù)據(jù)的維度的方法。通常來(lái)說(shuō)，我們期望得到的結(jié)果，是把原始數(shù)據(jù)的特征空間（n個(gè)d維樣本）投影到一個(gè)小一點(diǎn)的子空間里去，并盡可能表達(dá)的很好（就是說(shuō)損失信息最少）。常見(jiàn)的應(yīng)用在于模式識(shí)別中，我們可以通過(guò)減少特征空間的維度，抽取子空間的數(shù)據(jù)來(lái)最好的表達(dá)我們的數(shù)據(jù)，從而減少參數(shù)估計(jì)的誤差。注意，主成分分析通常會(huì)得到協(xié)方差矩陣和相關(guān)矩陣。這些矩陣可以通過(guò)原始數(shù)據(jù)計(jì)算出來(lái)。協(xié)方差矩陣包含平方和與向量積的和。相關(guān)矩陣與協(xié)方差矩陣類似，但是第一個(gè)變量，也就是第一列，是標(biāo)準(zhǔn)化后的數(shù)據(jù)。如果變量之間的方差很大，或者變量的量綱不統(tǒng)一，我們必須先標(biāo)準(zhǔn)化再進(jìn)行主成分分析。

二、PCA VS MDA

提到PCA，可能有些人會(huì)想到MDA（Multiple Discriminate Analysis,多元判別分析法），這兩者都是線性變換，而且很相似。只不過(guò)在PCA中，我們是找到一個(gè)成分（方向）來(lái)把我們的數(shù)據(jù)最大化方差，而在MDA中，我們的目標(biāo)是最大化不同類別之間的差異（比如說(shuō)，在模式識(shí)別問(wèn)題中，我們的數(shù)據(jù)包含多個(gè)類別，與兩個(gè)主成分的PCA相比，這就忽略了類別標(biāo)簽）。

換句話說(shuō)，通過(guò)PCA，我們把整個(gè)數(shù)據(jù)集（不含類別標(biāo)簽）投射到一個(gè)不同的子空間中，在MDA中，我們?cè)噲D決定一個(gè)合適的子空間來(lái)區(qū)分不同類別。再換種方式說(shuō)，PCA是找到數(shù)據(jù)傳播最廣的時(shí)候的最大方差的軸axis，MDA是最大化類別與類別之間的區(qū)別。

上文我們提到了子空間，那么怎么樣去尋找“好的”子空間呢？

假設(shè)我們的目標(biāo)是減少d維的數(shù)據(jù)集，將其投影到k維的子空間上（看k

下文中我們會(huì)計(jì)算數(shù)據(jù)中的特征向量（主成分），然后計(jì)算散布矩陣（scatter_matrix）中（也可以從協(xié)方差矩陣中計(jì)算）。每個(gè)特征向量與特征值相關(guān)，即特征向量的“長(zhǎng)度”或“大小”。如果發(fā)現(xiàn)每個(gè)特征值都很小，那就可以說(shuō)明我們的原始數(shù)據(jù)就已經(jīng)是一個(gè)“好的”空間了。但是，如果有些特征值比其他值要大得多，我們只需要關(guān)注那些特別大的特征值，因?yàn)檫@些值包含了數(shù)據(jù)分布情況的絕大部分信息。反之，那些接近于0的特征值包含的信息幾乎沒(méi)有，在新的特征空間里我們可以忽略不計(jì)。

三、PCA的過(guò)程

通常來(lái)說(shuō)有以下六步：

1.去掉數(shù)據(jù)的類別特征（label），將去掉后的d維數(shù)據(jù)作為樣本

2.計(jì)算d維的均值向量（即所有數(shù)據(jù)的每一維向量的均值）

3.計(jì)算所有數(shù)據(jù)的散布矩陣（或者協(xié)方差矩陣）

4.計(jì)算特征值（e1,e2,...,ed）以及相應(yīng)的特征向量（lambda1,lambda2,...,lambdad）

5.按照特征值的大小對(duì)特征向量降序排序，選擇前k個(gè)最大的特征向量，組成d*k維的矩陣W（其中每一列代表一個(gè)特征向量）

6.運(yùn)用d*K的特征向量矩陣W將樣本數(shù)據(jù)變換成新的子空間。（用數(shù)學(xué)式子表達(dá)就是，其中x是d*1維的向量，代表一個(gè)樣本，y是K*1維的在新的子空間里的向量）

1.數(shù)據(jù)準(zhǔn)備----生成三維樣本向量

首先隨機(jī)生成40*3維的數(shù)據(jù)，符合多元高斯分布。假設(shè)數(shù)據(jù)被分為兩類，其中一半類別為w1，另一半類別為w2

#coding:utf-8 import numpy as np np.random.seed(4294967295) mu_vec1 = np.array([0,0,0]) cov_mat1 = np.array([[1,0,0],[0,1,0],[0,0,1]]) class1_sample = np.random.multivariate_normal(mu_vec1, cov_mat1, 20).T assert class1_sample.shape == (3,20)#檢驗(yàn)數(shù)據(jù)的維度是否為3*20，若不為3*20，則拋出異常 mu_vec2 = np.array([1,1,1]) cov_mat2 = np.array([[1,0,0],[0,1,0],[0,0,1]]) class2_sample = np.random.multivariate_normal(mu_vec2, cov_mat2, 20).T assert class1_sample.shape == (3,20)#檢驗(yàn)數(shù)據(jù)的維度是否為3*20，若不為3*20，則拋出異常

運(yùn)行這段代碼后，我們就生成了包含兩個(gè)類別的樣本數(shù)據(jù)，其中每一列都是一個(gè)三維的向量，所有數(shù)據(jù)是這樣的矩陣：

結(jié)果：

2.作圖查看原始數(shù)據(jù)的分布

from matplotlib import pyplot as plt from mpl_toolkits.mplot3d import Axes3D from mpl_toolkits.mplot3d import proj3d fig = plt.figure(figsize=(8,8)) ax = fig.add_subplot(111, projection='3d') plt.rcParams['legend.fontsize'] = 10 ax.plot(class1_sample[0,:], class1_sample[1,:], class1_sample[2,:], 'o', markersize=8, color='blue', alpha=0.5, label='class1') ax.plot(class2_sample[0,:], class2_sample[1,:], class2_sample[2,:], '^', markersize=8, alpha=0.5, color='red', label='class2') plt.title('Samples for class 1 and class 2') ax.legend(loc='upper right') plt.show()

結(jié)果：

3.去掉數(shù)據(jù)的類別特征

1 all_samples = np.concatenate((class1_sample, class2_sample), axis=1) 2 assert all_samples.shape == (3,40)#檢驗(yàn)數(shù)據(jù)的維度是否為3*20，若不為3*20，則拋出異常

4.計(jì)算d維向量均值

mean_x = np.mean(all_samples[0,:]) mean_y = np.mean(all_samples[1,:]) mean_z = np.mean(all_samples[2,:]) mean_vector = np.array([[mean_x],[mean_y],[mean_z]]) print('Mean Vector: ', mean_vector)

結(jié)果：

1 print('Mean Vector: ', mean_vector) 2 Mean Vector:, 3 array([[ 0.68047077], 4 [ 0.52975093], 5 [ 0.43787182]]))

5.計(jì)算散步矩陣或者協(xié)方差矩陣

a.計(jì)算散步矩陣

散布矩陣公式：

其中m是向量的均值：（第4步已經(jīng)算出來(lái)是mean_vector）

1 scatter_matrix = np.zeros((3,3)) 2 for i in range(all_samples.shape[1]): 3 scatter_matrix += (all_samples[:,i].reshape(3,1) - mean_vector).dot((all_samples[:,i].reshape(3,1) - mean_vector).T) 4 print('Scatter Matrix: ', scatter_matrix)

結(jié)果：

1 print('Scatter Matrix: ', scatter_matrix) 2 ('Scatter Matrix:, 3 array([[ 46.81069724, 13.95578062, 27.08660175], 4 [ 13.95578062, 48.28401947, 11.32856266], 5 [ 27.08660175, 11.32856266, 50.51724488]]))

b.計(jì)算協(xié)方差矩陣

如果不計(jì)算散布矩陣的話，也可以用python里內(nèi)置的numpy.cov()函數(shù)直接計(jì)算協(xié)方差矩陣。因?yàn)樯⒉骄仃嚭蛥f(xié)方差矩陣非常類似，散布矩陣乘以（1/N-1）就是協(xié)方差，所以他們的特征空間是完全等價(jià)的（特征向量相同，特征值用一個(gè)常數(shù)（1/N-1，這里是1/39）等價(jià)縮放了）。協(xié)方差矩陣如下所示：

1 cov_mat = np.cov([all_samples[0,:],all_samples[1,:],all_samples[2,:]]) 2 print('Covariance Matrix: ', cov_mat)

結(jié)果：

1 >>> print('Covariance Matrix: ', cov_mat) 2 Covariance Matrix:, 3 array([[ 1.20027429, 0.35784053, 0.69452825], 4 [ 0.35784053, 1.23805178, 0.29047597], 5 [ 0.69452825, 0.29047597, 1.29531397]]))

6.計(jì)算相應(yīng)的特征向量和特征值

# 通過(guò)散布矩陣計(jì)算特征值和特征向量 eig_val_sc, eig_vec_sc = np.linalg.eig(scatter_matrix) # 通過(guò)協(xié)方差矩陣計(jì)算特征值和特征向量 eig_val_cov, eig_vec_cov = np.linalg.eig(cov_mat) for i in range(len(eig_val_sc)): eigvec_sc = eig_vec_sc[:,i].reshape(1,3).T eigvec_cov = eig_vec_cov[:,i].reshape(1,3).T assert eigvec_sc.all() == eigvec_cov.all() print('Eigenvector {}: {}'.format(i+1, eigvec_sc)) print('Eigenvalue {} from scatter matrix: {}'.format(i+1, eig_val_sc[i])) print('Eigenvalue {} from covariance matrix: {}'.format(i+1, eig_val_cov[i])) print('Scaling factor: ', eig_val_sc[i]/eig_val_cov[i]) print(40 * '-')

結(jié)果：

Eigenvector 1: [[-0.84190486] [-0.39978877] [-0.36244329]] Eigenvalue 1 from scatter matrix: 55.398855957302445 Eigenvalue 1 from covariance matrix: 1.4204834860846791 Scaling factor: 39.0 ---------------------------------------- Eigenvector 2: [[-0.44565232] [ 0.13637858] [ 0.88475697]] Eigenvalue 2 from scatter matrix: 32.42754801292286 Eigenvalue 2 from covariance matrix: 0.8314755900749456 Scaling factor: 39.0 ---------------------------------------- Eigenvector 3: [[ 0.30428639] [-0.90640489] [ 0.29298458]] Eigenvalue 3 from scatter matrix: 34.65493432806495 Eigenvalue 3 from covariance matrix: 0.8885880596939733 Scaling factor: 39.0 ----------------------------------------

其實(shí)從上面的結(jié)果就可以發(fā)現(xiàn)，通過(guò)散布矩陣和協(xié)方差矩陣計(jì)算的特征空間相同，協(xié)方差矩陣的特征值*39 = 散布矩陣的特征值

當(dāng)然，我們也可以快速驗(yàn)證一下特征值-特征向量的計(jì)算是否正確，是不是滿足方程（其中為協(xié)方差矩陣，v為特征向量，lambda為特征值）

1 for i in range(len(eig_val_sc)): 2 eigv = eig_vec_sc[:,i].reshape(1,3).T 3 np.testing.assert_array_almost_equal(scatter_matrix.dot(eigv), eig_val_sc[i] * eigv,decimal=6, err_msg='', verbose=True)

得出結(jié)果未返回異常，證明計(jì)算正確

注：np.testing.assert_array_almost_equal計(jì)算得出的結(jié)果不一樣會(huì)返回一下結(jié)果：

>>> np.testing.assert_array_almost_equal([1.0,2.33333,np.nan], ... [1.0,2.33339,np.nan], decimal=5) ... : AssertionError: Arrays are not almost equal (mismatch 50.0%) x: array([ 1. , 2.33333, NaN]) y: array([ 1. , 2.33339, NaN])

可視化特征向量

from matplotlib import pyplot as plt from mpl_toolkits.mplot3d import Axes3D from mpl_toolkits.mplot3d import proj3d from matplotlib.patches import FancyArrowPatch 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) fig = plt.figure(figsize=(7,7)) ax = fig.add_subplot(111, projection='3d') ax.plot(all_samples[0,:], all_samples[1,:], all_samples[2,:], 'o', markersize=8, color='green', alpha=0.2) ax.plot([mean_x], [mean_y], [mean_z], 'o', markersize=10, color='red', alpha=0.5) for v in eig_vec_sc.T: a = Arrow3D([mean_x, v[0]], [mean_y, v[1]], [mean_z, v[2]], mutation_scale=20, lw=3, arrowstyle="-|>", color="r") ax.add_artist(a) ax.set_xlabel('x_values') ax.set_ylabel('y_values') ax.set_zlabel('z_values') plt.title('Eigenvectors') plt.show()

結(jié)果：

7.根據(jù)特征值對(duì)特征向量降序排列

我們的目標(biāo)是減少特征空間的維度，即通過(guò)PCA方法將特征空間投影到一個(gè)小一點(diǎn)的子空間里，其中特征向量將會(huì)構(gòu)成新的特征空間的軸。然而，特征向量只會(huì)決定軸的方向，他們的單位長(zhǎng)度都為1，可以用代碼檢驗(yàn)一下：

1 for ev in eig_vec_sc: 2 numpy.testing.assert_array_almost_equal(1.0, np.linalg.norm(ev))

因此，對(duì)于低維的子空間來(lái)說(shuō)，決定丟掉哪個(gè)特征向量，就必須參考特征向量相應(yīng)的特征值。通俗來(lái)說(shuō)，如果一個(gè)特征向量的特征值特別小，那它所包含的數(shù)據(jù)分布的信息也很少，那么這個(gè)特征向量就可以忽略不計(jì)了。常用的方法是根據(jù)特征值對(duì)特征向量進(jìn)行降序排列，選出前k個(gè)特征向量

# 生成（特征向量，特征值）元祖 eig_pairs = [(np.abs(eig_val_sc[i]), eig_vec_sc[:,i]) for i in range(len(eig_val_sc))] #對(duì)（特征向量，特征值）元祖按照降序排列 eig_pairs.sort(key=lambda x: x[0], reverse=True) #輸出值 for i in eig_pairs: print(i[0])

結(jié)果：

1 84.5729942896 2 39.811391232 3 21.2275760682

8.選出前k個(gè)特征值最大的特征向量

本文的例子是想把三維的空間降維成二維空間，現(xiàn)在我們把前兩個(gè)最大特征值的特征向量組合起來(lái)，生成d*k維的特征向量矩陣W

1 matrix_w = np.hstack((eig_pairs[0][1].reshape(3,1), eig_pairs[1][1].reshape(3,1))) 2 print('Matrix W: ', matrix_w)

結(jié)果：

1 >>> print('Matrix W: ', matrix_w) 2 Matrix W:, 3 array([[-0.62497663, 0.2126888 ], 4 [-0.44135959, -0.88989795], 5 [-0.643899 , 0.40354071]]))

9.將樣本轉(zhuǎn)化為新的特征空間

最后一步，把2*3維的特征向量矩陣W帶到公式中，將樣本數(shù)據(jù)轉(zhuǎn)化為新的特征空間

matrix_w = np.hstack((eig_pairs[0][1].reshape(3,1), eig_pairs[1][1].reshape(3,1))) print('Matrix W: ', matrix_w) transformed = matrix_w.T.dot(all_samples) assert transformed.shape == (2,40), "The matrix is not 2x40 dimensional." plt.plot(transformed[0,0:20], transformed[1,0:20], 'o', markersize=7, color='blue', alpha=0.5, label='class1') plt.plot(transformed[0,20:40], transformed[1,20:40], '^', markersize=7, color='red', alpha=0.5, label='class2') plt.xlim([-4,4]) plt.ylim([-4,4]) plt.xlabel('x_values') plt.ylabel('y_values') plt.legend() plt.title('Transformed samples with class labels') plt.show()

到這一步，PCA的過(guò)程就結(jié)束了。其實(shí)python里有已經(jīng)寫好的模塊，可以直接拿來(lái)用，但是我覺(jué)得不管什么模塊，都要懂得它的原理是什么。matplotlib有matplotlib.mlab.PCA()，sklearn也有專門一個(gè)模塊Dimensionality reduction專門講PCA，包括傳統(tǒng)的PCA，以及增量PCA，核PCA等等，除了PCA以外，還有ZCA白化等等，在圖像處理中也經(jīng)常會(huì)用到。