今天看了用主成分分析簡化資料，就順便用MNIST資料集做了下實驗，想直觀地看一下效果，並通過完成這個小demo深入理解下原理。

我發現“是什麼、能做什麼、怎麼用、效果是什麼、原理是什麼、優缺點是什麼”這樣的思路能讓我更好地接受一個新知識，之所以把原理放在效果後面，是因為我比較喜歡先看看它的作用，視覺化意義之後能提起我對一個知識的興趣，加深對它意義的理解，後面看數學原理會容易，所以整篇文章就以這樣的思路組織整理。

主成分分析是什麼

主成分分析（Principal Component Analysis，PCA），一種降維方法，在PCA中，資料從原來的座標系轉換到了新的座標系，新座標系由資料本身決定，在新座標系中，第一個座標軸選擇的是原始資料中方差最大的方向，第二個座標軸選擇的是和第一個座標軸正交且具有最大方差的方向。該過程一直重複，重複次數為原始資料中特徵的數目。我們會發現，大部分方差都包含在最前面的幾個新座標軸中。因此，我們可以忽略餘下的座標軸，即對資料進行了降維處理。

初看這段話感覺是抽象的。方差大意味著什麼？方差是衡量源資料和期望值相差的度量值，方差越大，資料差別越大。選擇方差最大的方向，就是選擇資料差別最大的方向。重複特徵數目次，就是說找第一個特徵（第一維）方差最大的方向（即覆蓋資料點最多的一條直線），做第一個軸，正交且最大方差方向做第二個軸，在此基礎上再看第二個特徵（第二維），找方差最大方向做第一個軸，正交且最大方差方向做第二個軸，依次類推。這樣執行後會發現前幾個座標軸已經差不多囊括所有大差異了，剩下的就不要了，所以實現了降維。

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

主成分分析能做什麼

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

主成分分析怎麼用

要做的事就是使用tensorflow裡的MNIST資料集，取前100張圖片中所有的手寫數字7圖片，對他們進行主成分分析，輸出經過降維反變換回去的圖片，對比差異，看看降維後的效果。

引入MNIST資料集、numpy和PIL的Image

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

獲得MNIST資料集的所有圖片和標籤

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

這裡可以看看imgs和labels的type和shape，對於一個python初學者來說總是想搞清楚各個變數的型別和長相。

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

取前1000張圖片裡的100個數字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維的一維陣列，變成我們想看的圖片就需要把它reshape成28x28的二維陣列，然後再用Image裡的方法，把它拼成一張2x5的大圖。

由於tensorflow中MNIST都是灰度圖（L），所以shape是（55000，784），每張圖的dtype是float32，如果是彩色圖（RGB），shape可能是（55000，784，3），圖的dtype是uint8，從array轉到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)) # 第二個引數為每次貼上起始點的橫縱座標。在本例中，分別為（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)) returnnew_img

效果圖

ten_origin_7_imgs=comb_imgs(origin_7_imgs, 10, 10, 28, 28, 'L') ten_origin_7_imgs.show()

實現主成分分析演算法（詳細程式碼解析在文章後面的原理部分）

def pca(data_mat, top_n_feat=99999999): """ 主成分分析： 輸入：矩陣data_mat ，其中該矩陣中儲存訓練資料，每一行為一條訓練資料 保留前n個特徵top_n_feat，預設全保留 返回：降維後的資料集和原始資料被重構後的矩陣（即降維後反變換回矩陣）

""" # 獲取資料條數和每條的維數 num_data,dim = data_mat.shape print(num_data) # 100 print(dim) # 784 # 資料中心化，即指變數減去它的均值 mean_vals = data_mat.mean(axis=0) #shape:(784,) mean_removed = data_mat - mean_vals # shape:(100, 784) # 計算協方差矩陣（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個特徵後重構的特徵向量矩陣reorganize eig vects, # shape為(784, top_n_feat)，top_n_feat最大為特徵總數 reg_eig_vects = eig_vects[:, eig_val_index] # 將資料轉到新空間 low_d_data_mat = mean_removed * reg_eig_vects # shape: (100, top_n_feat), top_n_feat最大為特徵總數 recon_mat = (low_d_data_mat * reg_eig_vects.T) + mean_vals # 根據前幾個特徵向量重構回去的矩陣，shape:(100, 784)

return low_d_data_mat, recon_mat

呼叫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個特徵向量重構的效果圖

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

主成分分析效果是什麼

不難發現降維後數字7長得規則多了，或許降維後再用tensorflow入門教程的softmax進行分類accuracy會更高。

主成分分析的原理是什麼

前面轉座標軸從理論上考慮，這裡主要從數學的角度考慮。

第一個主成分是資料差異最大（方差最大）的方向，第二個主成分是資料差異次大且與第一個主成分正交的方向。通過資料集的協方差矩陣及其特徵值分析，就能求得這些主成分的值。

統計學中的幾個概念

平均值

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

標準差

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

方差

標準差的平方。

協方差

方差是用來描述一維資料的，協方差用來描述二維資料，用來描述兩個隨機變數之間的關係，如果是正值，則說明兩變數正相關，負值，說明負相關，0，說明不相關，即相互獨立。

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

協方差矩陣

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

比如資料集有三個維度，X，Y，Z，則協方差矩陣為

可見，矩陣的對角線為方差，由於cov(X,Y) = cov(Y, X)，所以是一個對稱矩陣。

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

結合程式碼分析原理

目的就是找出差異最大的方向，也就是影響最大的幾個特徵，數學上通過協方差矩陣來找差異最大的特徵，排序，最後找到降維後的特徵矩陣。

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

協方差矩陣需要計算平均值，上面強調了計算的是不同維度的協方差，資料每行是一個樣本，每列是一個維度，因此計算的是列的平均值，即axis=0，因此shape為（784，）。使用np的cov函式計算協方差矩陣，api入下：

numpy.cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None, aweights=None)[source]

**rowvar代表是否轉置。在API裡，預設rowvar是True，也就是行是variable，列是observation，我們這裡列是observation，行是variable。

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

• mat(cov_mat)：將輸入轉成矩陣。和matrix，asmatrix不同，如果輸入已經是舉著嗯或者ndarray，它不會製作副本。相當於matrix(data, copy=False) 詳細API請點這裡（https://docs.scipy.org/doc/numpy/reference/generated/numpy.mat.html）
• linalg.eig(a)：計算特徵值和特徵向量 詳細API請點這裡（https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html）

矩陣乘法對應了一個變換，在這個變換的過程中，原向量主要發生旋轉、伸縮的變化。如果矩陣對某一個向量或某些向量只發生伸縮變換，不對這些向量產生旋轉的效果，那麼這些向量就稱為這個矩陣的特徵向量，伸縮的比例就是特徵值。

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

numpy.argsort(a, axis=-1, kind='quicksort', order=None) 詳細API請點這裡（https://docs.scipy.org/doc/numpy/reference/generated/numpy.argsort.html）

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

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

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

# 將資料轉到新空間 low_d_data_mat = mean_removed * reg_eig_vects # shape: (100, top_n_feat), top_n_feat最大為特徵總數 recon_mat = (low_d_data_mat * reg_eig_vects.T) + mean_vals # 根據前幾個特徵向量重構回去的矩陣，shape:(100, 784)

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

recon_mat再將矩陣變回（100，784），得到降維後再重構的矩陣。

主成分分析的優缺點是什麼

優點：降低資料的複雜性，識別最重要的特徵

缺點：不一定需要，且可能損失有用資訊

適用資料型別：數值型資料