主成分分析降維（MNIST資料集）
今天看了用主成分分析簡化資料，就順便用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/" 這裡可以看看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的大圖。
def array_to_img(array):
array=array*255
new_img=Image.fromarray(array.astype(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))
return new_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]
詳細的API請點這裡
rowvar代表是否轉置。在API裡，預設rowvar是True，也就是行是variable，列是observation，我們這裡列是observation，行是variable。
eig_vals, eig_vects = linalg.eig(mat(cov_mat)) # 計算特徵值和特徵向量矩陣乘法對應了一個變換，在這個變換的過程中，原向量主要發生旋轉、伸縮的變化。如果矩陣對某一個向量或某些向量只發生伸縮變換，不對這些向量產生旋轉的效果，那麼這些向量就稱為這個矩陣的特徵向量，伸縮的比例就是特徵值。
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]這裡有一個語法問題，[::]代表切片，[開始：結束：步長]，負號代表從後往前每隔步長個取一次，比如有一個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），得到降維後再重構的矩陣。
主成分分析的優缺點是什麼
- 優點：降低資料的複雜性，識別最重要的特徵
- 缺點：不一定需要，且可能損失有用資訊
適用資料型別：數值型資料