今天看了用主成分分析簡化數(shù)據(jù)，就順便用MNIST數(shù)據(jù)集做了下實驗，想直觀地看一下效果，并通過完成這個小demo深入理解下原理。

我發(fā)現(xiàn)“是什么、能做什么、怎么用、效果是什么、原理是什么、優(yōu)缺點是什么”這樣的思路能讓我更好地接受一個新知識，之所以把原理放在效果后面，是因為我比較喜歡先看看它的作用，可視化意義之后能提起我對一個知識的興趣，加深對它意義的理解，后面看數(shù)學(xué)原理會容易，所以整篇文章就以這樣的思路組織整理。

主成分分析是什么

主成分分析（Principal Component Analysis，PCA），一種降維方法，在PCA中，數(shù)據(jù)從原來的坐標系轉(zhuǎn)換到了新的坐標系，新坐標系由數(shù)據(jù)本身決定，在新坐標系中，第一個坐標軸選擇的是原始數(shù)據(jù)中方差最大的方向，第二個坐標軸選擇的是和第一個坐標軸正交且具有最大方差的方向。該過程一直重復(fù)，重復(fù)次數(shù)為原始數(shù)據(jù)中特征的數(shù)目。我們會發(fā)現(xiàn)，大部分方差都包含在最前面的幾個新坐標軸中。因此，我們可以忽略余下的坐標軸，即對數(shù)據(jù)進行了降維處理。

初看這段話感覺是抽象的。方差大意味著什么？方差是衡量源數(shù)據(jù)和期望值相差的度量值，方差越大，數(shù)據(jù)差別越大。選擇方差最大的方向，就是選擇數(shù)據(jù)差別最大的方向。重復(fù)特征數(shù)目次，就是說找第一個特征（第一維）方差最大的方向（即覆蓋數(shù)據(jù)點最多的一條直線），做第一個軸，正交且最大方差方向做第二個軸，在此基礎(chǔ)上再看第二個特征（第二維），找方差最大方向做第一個軸，正交且最大方差方向做第二個軸，依次類推。這樣執(zhí)行后會發(fā)現(xiàn)前幾個坐標軸已經(jīng)差不多囊括所有大差異了，剩下的就不要了，所以實現(xiàn)了降維。

上面從理論上講了主成分分析和它是如何一步一步實現(xiàn)降維的，有一個感性認識。

主成分分析能做什么

降維，在多個指標中只取重要的幾個指標，能使復(fù)雜問題簡單化，就像說話說重點一樣。

主成分分析怎么用

要做的事就是使用tensorflow里的MNIST數(shù)據(jù)集，取前100張圖片中所有的手寫數(shù)字7圖片，對他們進行主成分分析，輸出經(jīng)過降維反變換回去的圖片，對比差異，看看降維后的效果。

• 引入MNIST數(shù)據(jù)集、numpy和PIL的Image

import tensorflow.examples.tutorials.mnist.input_data as input_data
import numpy as np
from PIL import Image

• 獲得MNIST數(shù)據(jù)集的所有圖片和標簽

mnist = input_data.read_data_sets("MNIST_data/", one_hot=False)
imgs = mnist.train.images
labels = mnist.train.labels

這里可以看看imgs和labels的type和shape，對于一個python初學(xué)者來說總是想搞清楚各個變量的類型和長相。

print(type(imgs))             # <type 'numpy.ndarray'>
print(type(labels))           # <type 'numpy.ndarray'>
print(imgs.shape)             # (55000, 784)
print(labels.shape)           # (55000,)

• 取前1000張圖片里的100個數(shù)字7

origin_7_imgs = []
for i in range(1000):
      if labels[i] == 7 and len(origin_7_imgs) < 100:
          origin_7_imgs.append(imgs[i])

看看shape

print(np.array(origin_7_imgs).shape)   # (100, 784)

• 把10張圖片排成2x5的表格
  由于一張圖片是一個784維的一維數(shù)組，變成我們想看的圖片就需要把它reshape成28x28的二維數(shù)組，然后再用Image里的方法，把它拼成一張2x5的大圖。
• 由于tensorflow中MNIST都是灰度圖（L），所以shape是（55000，784），每張圖的dtype是float32，如果是彩色圖（RGB），shape可能是（55000，784，3），圖的dtype是uint8，從array轉(zhuǎn)到Image需要用下面的方法：

  def array_to_img(array):
      array=array*255
      new_img=Image.fromarray(array.astype(np.uint8))
      return new_img
  

• 拼圖

 def comb_imgs(origin_imgs, col, row, each_width, each_height, new_type):
     new_img = Image.new(new_type, (col* each_width, row* each_height)) 
     for i in range(len(origin_imgs)):
         each_img = array_to_img(np.array(origin_imgs[i]).reshape(each_width, each_width))
         # 第二個參數(shù)為每次粘貼起始點的橫縱坐標。在本例中，分別為（0，0）（28，0）（28*2，0）依次類推，第二行是（0，28）（28，28），（28*2，28）類推
         new_img.paste(each_img, ((i % col) * each_width, (i / col) * each_width)) 
     return new_img
 ```
- 效果圖
```Python
 ten_origin_7_imgs=comb_imgs(origin_7_imgs, 10, 10, 28, 28, 'L')
 ten_origin_7_imgs.show()

數(shù)字7原圖

• 實現(xiàn)主成分分析算法（詳細代碼解析在文章后面的原理部分）

def pca(data_mat, top_n_feat=99999999):
    """ 
    主成分分析：  
    輸入：矩陣data_mat ，其中該矩陣中存儲訓(xùn)練數(shù)據(jù)，每一行為一條訓(xùn)練數(shù)據(jù)  
         保留前n個特征top_n_feat，默認全保留
    返回：降維后的數(shù)據(jù)集和原始數(shù)據(jù)被重構(gòu)后的矩陣（即降維后反變換回矩陣）
    """  

    # 獲取數(shù)據(jù)條數(shù)和每條的維數(shù) 
    num_data,dim = data_mat.shape  
    print(num_data)  # 100
    print(dim)   # 784

    # 數(shù)據(jù)中心化，即指變量減去它的均值
    mean_vals = data_mat.mean(axis=0)  #shape:(784,)
    mean_removed = data_mat - mean_vals # shape:(100, 784)

    # 計算協(xié)方差矩陣（Find covariance matrix）
    cov_mat = np.cov(mean_removed, rowvar=0) # shape：(784, 784)

    # 計算特征值(Find eigenvalues and eigenvectors)
    eig_vals, eig_vects = linalg.eig(mat(cov_mat)) # 計算特征值和特征向量，shape分別為（784，）和(784, 784)

    eig_val_index = argsort(eig_vals)  # 對特征值進行從小到大排序，argsort返回的是索引，即下標

    eig_val_index = eig_val_index[:-(top_n_feat + 1) : -1] # 最大的前top_n_feat個特征的索引
    # 取前top_n_feat個特征后重構(gòu)的特征向量矩陣reorganize eig vects, 
    # shape為(784, top_n_feat)，top_n_feat最大為特征總數(shù)
    reg_eig_vects = eig_vects[:, eig_val_index] 
    
    # 將數(shù)據(jù)轉(zhuǎn)到新空間
    low_d_data_mat = mean_removed * reg_eig_vects # shape: (100, top_n_feat), top_n_feat最大為特征總數(shù)
    recon_mat = (low_d_data_mat * reg_eig_vects.T) + mean_vals # 根據(jù)前幾個特征向量重構(gòu)回去的矩陣，shape:(100, 784)
    
    return low_d_data_mat, recon_mat

• 調(diào)用PCA進行降維

low_d_feat_for_7_imgs, recon_mat_for_7_imgs = pca(np.array(origin_7_imgs), 1) # 只取最重要的1個特征
print(low_d_feat_for_7_imgs.shape) # (100, 1)
print(recon_mat_for_7_imgs.shape) # (100, 784)

• 看降維后只用1個特征向量重構(gòu)的效果圖

low_d_img = comb_imgs(recon_mat_for_7_imgs, 10, 10, 28, 28, 'L')
low_d_img.show()

數(shù)字7降維后的圖

主成分分析效果是什么

降維前后對比圖

不難發(fā)現(xiàn)降維后數(shù)字7長得規(guī)則多了，或許降維后再用tensorflow入門教程的softmax進行分類accuracy會更高。

主成分析的原理是什么

前面轉(zhuǎn)坐標軸從理論上考慮，這里主要從數(shù)學(xué)的角度考慮。
第一個主成分是數(shù)據(jù)差異最大（方差最大）的方向，第二個主成分是數(shù)據(jù)差異次大且與第一個主成分正交的方向。通過數(shù)據(jù)集的協(xié)方差矩陣及其特征值分析，就能求得這些主成分的值。

統(tǒng)計學(xué)中的幾個概念

• 平均值
  這個最為熟悉最不容易忘記，描述樣本集合的中間點。

• 標準差
  描述樣本集合中各個點到平均值的距離。

• 方差
  標準差的平方。

  
  
  方差

• 協(xié)方差
  方差是用來描述一維數(shù)據(jù)的，協(xié)方差用來描述二維數(shù)據(jù)，用來描述兩個隨機變量之間的關(guān)系，如果是正值，則說明兩變量正相關(guān)，負值，說明負相關(guān)，0，說明不相關(guān)，即相互獨立。

  
  
  協(xié)方差
  
  

  從公式可以看出協(xié)方差的一些性質(zhì)：
  1、cov(X, X) = var(X)
  2、cov(X,Y) = cov(Y, X)

• 協(xié)方差矩陣
  協(xié)方差可以描述二維數(shù)據(jù)，但是對于多維數(shù)據(jù)來說，我們只能兩個點兩個點地計算多次協(xié)方差，一個n維數(shù)據(jù)，我們需要計算C(n, 2)=A(n,2)/2=n!/((n-2)!*2)個協(xié)方差，自然就需要用矩陣來組織這些數(shù)據(jù)。所以協(xié)方差矩陣的定義為

  
  
  協(xié)方差矩陣
  
  

  比如數(shù)據(jù)集有三個維度，X，Y，Z，則協(xié)方差矩陣為

  
  
  三維協(xié)方差矩陣
  
  可見，矩陣的對角線為方差，由于cov(X,Y) = cov(Y, X)，所以是一個對稱矩陣。

  注意，協(xié)方差矩陣計算的是不同維度之間的協(xié)方差，不是不同樣本之間的協(xié)方差。

結(jié)合代碼分析原理

目的就是找出差異最大的方向，也就是影響最大的幾個特征，數(shù)學(xué)上通過協(xié)方差矩陣來找差異最大的特征，排序，最后找到降維后的特征矩陣。

  # 數(shù)據(jù)中心化，即指變量減去它的均值
  mean_vals = data_mat.mean(axis=0)  #shape:(784,)
  mean_removed = data_mat - mean_vals # shape:(100, 784)
  # 計算協(xié)方差矩陣（Find covariance matrix）
  cov_mat = np.cov(mean_removed, rowvar=0)

協(xié)方差矩陣需要計算平均值，上面強調(diào)了計算的是不同維度的協(xié)方差，數(shù)據(jù)每行是一個樣本，每列是一個維度，因此計算的是列的平均值，即axis=0，因此shape為（784，）。使用np的cov函數(shù)計算協(xié)方差矩陣，api入下：

numpy.cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None, aweights=None)[source]
詳細的API請點這里

**rowvar代表是否轉(zhuǎn)置。在API里，默認rowvar是True，也就是行是variable，列是observation，我們這里列是observation，行是variable。**

 eig_vals, eig_vects = linalg.eig(mat(cov_mat)) # 計算特征值和特征向量

• mat(cov_mat)：將輸入轉(zhuǎn)成矩陣。和matrix，asmatrix不同，如果輸入已經(jīng)是舉著嗯或者ndarray，它不會制作副本。相當于matrix(data, copy=False)
  詳細API請點這里
• linalg.eig(a)：計算特征值和特征向量
  詳細API請點這里

矩陣乘法對應(yīng)了一個變換，在這個變換的過程中，原向量主要發(fā)生旋轉(zhuǎn)、伸縮的變化。如果矩陣對某一個向量或某些向量只發(fā)生伸縮變換，不對這些向量產(chǎn)生旋轉(zhuǎn)的效果，那么這些向量就稱為這個矩陣的特征向量，伸縮的比例就是特征值。

eig_val_index = argsort(eig_vals)  # 對特征值進行從小到大排序，argsort返回的是索引，即下標

numpy.argsort(a, axis=-1, kind='quicksort', order=None)
詳細API請點這里

eig_val_index = eig_val_index[:-(top_n_feat + 1) : -1] # 最大的前top_n_feat個特征的索引
reg_eig_vects = eig_vects[:, eig_val_index]

這里有一個語法問題，[::]代表切片，[開始：結(jié)束：步長]，負號代表從后往前每隔步長個取一次，比如有一個array[1, 2, 3, 4, 5]，取[:-4:-2]，0是第一個，-1是最后一個(在這里是5的下標)，從最后一個往前數(shù)，一直數(shù)到-4（在這里是2的下標），每兩個取1個數(shù)，最后得到的array是[5, 3]。

[:, eig_val_index]代表第一維不變，第二維要eig_val_index個，所以它的shape是（784，top_n_feat）

 # 將數(shù)據(jù)轉(zhuǎn)到新空間
  low_d_data_mat = mean_removed * reg_eig_vects # shape: (100, top_n_feat), top_n_feat最大為特征總數(shù)
  recon_mat = (low_d_data_mat * reg_eig_vects.T) + mean_vals # 根據(jù)前幾個特征向量重構(gòu)回去的矩陣，shape:(100, 784)

一個shape是（100，784）的矩陣，乘以一個shape是（784，top_n_feat）的矩陣，最后得到降維的矩陣（100， top_n_feat）

recon_mat再將矩陣變回（100，784），得到降維后再重構(gòu)的矩陣。

主成分分析的優(yōu)缺點是什么

• 優(yōu)點：降低數(shù)據(jù)的復(fù)雜性，識別最重要的特征
• 缺點：不一定需要，且可能損失有用信息
  適用數(shù)據(jù)類型：數(shù)值型數(shù)據(jù)