機器學習基礎與實踐（三）----數據降維之PCA

時間 2019-11-14

標籤機器學習基礎實踐數據 pca 简体版

原文原文鏈接

寫在前面：原本這篇應該是上週四更新，可是上週四寫了一篇深度學習的反向傳播法的過程，就推遲更新了。原本想參考PRML來寫，可是發現裏面涉及到比較多的數學知識，寫出來可能很差理解，我決定仍是用最通俗的方法解釋PCA，並舉一個實例一步步計算，而後再進行數學推導，最後再介紹一些變種以及相應的程序。（數學推導及變種下次再寫好了）html

正文：python

　　在數據處理中，常常會遇到特徵維度比樣本數量多得多的狀況，若是拿到實際工程中去跑，效果不必定好。一是由於冗餘的特徵會帶來一些噪音，影響計算的結果；二是由於無關的特徵會加大計算量，耗費時間和資源。因此咱們一般會對數據從新變換一下，再跑模型。數據變換的目的不只僅是降維，還能夠消除特徵之間的相關性，並發現一些潛在的特徵變量。spring

　　1、PCA的目的併發

　　PCA是一種在儘量減小信息損失的狀況下找到某種方式下降數據的維度的方法。一般來講，咱們指望獲得的結果，是把原始數據的特徵空間（n個d維樣本）投影到一個小一點的子空間裏去，並儘量表達的很好（就是說損失信息最少）。常見的應用在於模式識別中，咱們能夠經過減小特徵空間的維度，抽取子空間的數據來最好的表達咱們的數據，從而減小參數估計的偏差。注意，主成分分析一般會獲得協方差矩陣和相關矩陣。這些矩陣能夠經過原始數據計算出來。協方差矩陣包含平方和與向量積的和。相關矩陣與協方差矩陣相似，可是第一個變量，也就是第一列，是標準化後的數據。若是變量之間的方差很大，或者變量的量綱不統一，咱們必須先標準化再進行主成分分析。dom

　　2、PCA VS MDA函數

　　提到PCA，可能有些人會想到MDA（Multiple Discriminate Analysis,多元判別分析法），這二者都是線性變換，並且很類似。只不過在PCA中，咱們是找到一個成分（方向）來把咱們的數據最大化方差，而在MDA中，咱們的目標是最大化不一樣類別之間的差別（好比說，在模式識別問題中，咱們的數據包含多個類別，與兩個主成分的PCA相比，這就忽略了類別標籤）。學習

換句話說，經過PCA，咱們把整個數據集（不含類別標籤）投射到一個不一樣的子空間中，在MDA中，咱們試圖決定一個合適的子空間來區分不一樣類別。再換種方式說，PCA是找到數據傳播最廣的時候的最大方差的軸axis，MDA是最大化類別與類別之間的區別。spa

　　上文咱們提到了子空間，那麼怎麼樣去尋找「好的」子空間呢？3d

　　假設咱們的目標是減小d維的數據集，將其投影到k維的子空間上（看k<d）。因此，咱們如何來肯定k呢？如何知道咱們選擇的特徵空間可以很好的表達原始數據呢？code

　　下文中咱們會計算數據中的特徵向量（主成分），而後計算散佈矩陣（scatter_matrix）中（也能夠從協方差矩陣中計算）。每一個特徵向量與特徵值相關，即特徵向量的「長度」或「大小」。若是發現每一個特徵值都很小，那就能夠說明咱們的原始數據就已是一個「好的」空間了。可是，若是有些特徵值比其餘值要大得多，咱們只須要關注那些特別大的特徵值，由於這些值包含了數據分佈狀況的絕大部分信息。反之，那些接近於0的特徵值包含的信息幾乎沒有，在新的特徵空間裏咱們能夠忽略不計。

　　3、PCA的過程

　　一般來講有如下六步：

　　1.數據準備----生成三維樣本向量

　　首先隨機生成40*3維的數據，符合多元高斯分佈。假設數據被分爲兩類，其中一半類別爲w1，另外一半類別爲w2

 1 #coding:utf-8
 2 import numpy as np
 3 
 4 np.random.seed(4294967295) 
 5 
 6 mu_vec1 = np.array([0,0,0])
 7 cov_mat1 = np.array([[1,0,0],[0,1,0],[0,0,1]])
 8 class1_sample = np.random.multivariate_normal(mu_vec1, cov_mat1, 20).T
 9 assert class1_sample.shape == (3,20)#檢驗數據的維度是否爲3*20，若不爲3*20，則拋出異常
10 
11 mu_vec2 = np.array([1,1,1])
12 cov_mat2 = np.array([[1,0,0],[0,1,0],[0,0,1]])
13 class2_sample = np.random.multivariate_normal(mu_vec2, cov_mat2, 20).T
14 assert class1_sample.shape == (3,20)#檢驗數據的維度是否爲3*20，若不爲3*20，則拋出異常

　　運行這段代碼後，咱們就生成了包含兩個類別的樣本數據，其中每一列都是一個三維的向量，全部數據是這樣的矩陣：

　　結果：

　　2.做圖查看原始數據的分佈

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

　　結果：

　　3.去掉數據的類別特徵

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

　　4.計算d維向量均值

1 mean_x = np.mean(all_samples[0,:])
2 mean_y = np.mean(all_samples[1,:])
3 mean_z = np.mean(all_samples[2,:])
4 
5 mean_vector = np.array([[mean_x],[mean_y],[mean_z]])
6 
7 print('Mean Vector:\n', mean_vector)

　　結果：

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

　　5.計算散步矩陣或者協方差矩陣

　　a.計算散步矩陣

　　散佈矩陣公式：

　　其中m是向量的均值：（第4步已經算出來是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:\n', scatter_matrix)

　　結果：

1  print('Scatter Matrix:\n', 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.計算協方差矩陣

　　若是不計算散佈矩陣的話，也能夠用python裏內置的numpy.cov()函數直接計算協方差矩陣。由於散步矩陣和協方差矩陣很是相似，散佈矩陣乘以（1/N-1）就是協方差，因此他們的特徵空間是徹底等價的（特徵向量相同，特徵值用一個常數（1/N-1，這裏是1/39）等價縮放了）。協方差矩陣以下所示：

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

　　結果：

1 >>> print('Covariance Matrix:\n', 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.計算相應的特徵向量和特徵值

 1 # 經過散佈矩陣計算特徵值和特徵向量
 2 eig_val_sc, eig_vec_sc = np.linalg.eig(scatter_matrix)
 3 
 4 # 經過協方差矩陣計算特徵值和特徵向量
 5 eig_val_cov, eig_vec_cov = np.linalg.eig(cov_mat)
 6 
 7 for i in range(len(eig_val_sc)):
 8     eigvec_sc = eig_vec_sc[:,i].reshape(1,3).T
 9     eigvec_cov = eig_vec_cov[:,i].reshape(1,3).T
10     assert eigvec_sc.all() == eigvec_cov.all()
11 
12 print('Eigenvector {}: \n{}'.format(i+1, eigvec_sc))
13 print('Eigenvalue {} from scatter matrix: {}'.format(i+1, eig_val_sc[i]))
14 print('Eigenvalue {} from covariance matrix: {}'.format(i+1, eig_val_cov[i]))
15 print('Scaling factor: ', eig_val_sc[i]/eig_val_cov[i])
16 print(40 * '-')

　　結果：

 1 Eigenvector 1:
 2     [[-0.84190486]
 3      [-0.39978877]
 4      [-0.36244329]]
 5     Eigenvalue 1 from scatter matrix: 55.398855957302445
 6     Eigenvalue 1 from covariance matrix: 1.4204834860846791
 7     Scaling factor:  39.0
 8     ----------------------------------------
 9     Eigenvector 2:
10     [[-0.44565232]
11      [ 0.13637858]
12      [ 0.88475697]]
13     Eigenvalue 2 from scatter matrix: 32.42754801292286
14     Eigenvalue 2 from covariance matrix: 0.8314755900749456
15     Scaling factor:  39.0
16     ----------------------------------------
17     Eigenvector 3:
18     [[ 0.30428639]
19      [-0.90640489]
20      [ 0.29298458]]
21     Eigenvalue 3 from scatter matrix: 34.65493432806495
22     Eigenvalue 3 from covariance matrix: 0.8885880596939733
23     Scaling factor:  39.0
24     ----------------------------------------

　　其實從上面的結果就能夠發現，經過散佈矩陣和協方差矩陣計算的特徵空間相同，協方差矩陣的特徵值*39 = 散佈矩陣的特徵值

　　固然，咱們也能夠快速驗證一下特徵值-特徵向量的計算是否正確，是否是知足方程（其中爲協方差矩陣，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)

　　得出結果未返回異常，證實計算正確

　　注：np.testing.assert_array_almost_equal計算得出的結果不同會返回一下結果：

 1 >>> np.testing.assert_array_almost_equal([1.0,2.33333,np.nan],
 2 ...                                      [1.0,2.33339,np.nan], decimal=5)
 3 ...
 4 <type 'exceptions.AssertionError'>:
 5 AssertionError:
 6 Arrays are not almost equal
 7 
 8 (mismatch 50.0%)
 9  x: array([ 1.     ,  2.33333,      NaN])
10  y: array([ 1.     ,  2.33339,      NaN])

　　可視化特徵向量

 1 from matplotlib import pyplot as plt
 2 from mpl_toolkits.mplot3d import Axes3D
 3 from mpl_toolkits.mplot3d import proj3d
 4 from matplotlib.patches import FancyArrowPatch
 5 
 6 
 7 class Arrow3D(FancyArrowPatch):
 8     def __init__(self, xs, ys, zs, *args, **kwargs):
 9         FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
10         self._verts3d = xs, ys, zs
11 
12     def draw(self, renderer):
13         xs3d, ys3d, zs3d = self._verts3d
14         xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
15         self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
16         FancyArrowPatch.draw(self, renderer)
17 
18 fig = plt.figure(figsize=(7,7))
19 ax = fig.add_subplot(111, projection='3d')
20 
21 ax.plot(all_samples[0,:], all_samples[1,:], all_samples[2,:], 'o', markersize=8, color='green', alpha=0.2)
22 ax.plot([mean_x], [mean_y], [mean_z], 'o', markersize=10, color='red', alpha=0.5)
23 for v in eig_vec_sc.T:
24     a = Arrow3D([mean_x, v[0]], [mean_y, v[1]], [mean_z, v[2]], mutation_scale=20, lw=3, arrowstyle="-|>", color="r")
25     ax.add_artist(a)
26 ax.set_xlabel('x_values')
27 ax.set_ylabel('y_values')
28 ax.set_zlabel('z_values')
29 
30 plt.title('Eigenvectors')
31 
32 plt.show()

　　結果：

　　7.根據特徵值對特徵向量降序排列

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

　　所以，對於低維的子空間來講，決定丟掉哪一個特徵向量，就必須參考特徵向量相應的特徵值。通俗來講，若是一個特徵向量的特徵值特別小，那它所包含的數據分佈的信息也不多，那麼這個特徵向量就能夠忽略不計了。經常使用的方法是根據特徵值對特徵向量進行降序排列，選出前k個特徵向量

1 # 生成（特徵向量，特徵值）元祖
2 eig_pairs = [(np.abs(eig_val_sc[i]), eig_vec_sc[:,i]) for i in range(len(eig_val_sc))]
3 
4 #對（特徵向量，特徵值）元祖按照降序排列
5 eig_pairs.sort(key=lambda x: x[0], reverse=True)
6 
7 #輸出值
8 for i in eig_pairs:
9     print(i[0])

　　結果：

 1 84.5729942896
 2 39.811391232     
 3 21.2275760682

　　8.選出前k個特徵值最大的特徵向量

　　本文的例子是想把三維的空間降維成二維空間，如今咱們把前兩個最大特徵值的特徵向量組合起來，生成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:\n', matrix_w)

　　結果：

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

　　9.將樣本轉化爲新的特徵空間

　　最後一步，把2*3維的特徵向量矩陣W帶到公式中，將樣本數據轉化爲新的特徵空間

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

　　結果：

　　到這一步，PCA的過程就結束了。其實python裏有已經寫好的模塊，能夠直接拿來用，可是我以爲無論什麼模塊，都要懂得它的原理是什麼。matplotlib有matplotlib.mlab.PCA()，sklearn也有專門一個模塊Dimensionality reduction專門講PCA，包括傳統的PCA，也就是我上文寫的，以及增量PCA，核PCA等等，除了PCA之外，還有ZCA白化等等，在圖像處理中也常常會用到，內容太多，下次再寫。

　　最後推薦一個博客，動態展現了PCA的過程：http://setosa.io/ev/principal-component-analysis/ 寫的也很清楚，能夠看一下；再推薦一個維基百科的，講的真的是詳細啊https://en.wikipedia.org/wiki/Principal_component_analysis

------------------------------------本博客全部內容以學習、研究和分享爲主，如需轉載，請聯繫本人，標明做者和出處，而且是非商業用途，謝謝！--------------------------------