鶯嘴啄花紅溜，燕尾點波綠皺。
指冷玉笙寒，吹徹小梅春透。
依舊，依舊，人與綠楊俱瘦。
——《如夢令·春景》 秦觀
更多精彩內容請關注微信公眾號 “優化與算法”
1、背景
隨着信息技術的發展，數據量呈現爆照式增長，高維海量數據給傳統的數據處理方法帶來了嚴峻的挑戰，因此，開發高效的數據處理技術是非常必要的。數據降維是解決維度災難的一種有效方法，之所以對數據進行降維是因為：
- 在原始的高維空間中，包含有冗余信息以及噪音信息，在實際應用例如圖像識別中造成了誤差，降低了准確率；而通過降維,我們希望減少冗余信息所造成的誤差，提高識別（或其他應用）的精度。
- 希望通過降維算法來尋找數據內部的本質結構特征。
主成分分析（Principal Component Analysis，PCA）是一種常見的線性降維方法，廣泛應用於圖像處理、人臉識別、數據壓縮、信號去噪等領域。
2、算法思想與原理
設原數據大小為 \(N\times M\) ，經過PCA降維后的數據大小為 \(N \times K\)，其中 \(K<M\)。
PCA的中心思想是：將高維數據投影到具有最大方差的低維空間中。
2.1 從不同角度理解PCA
PCA原理可以從不同角度來理解：
- 從空間變換的角度來看，這里我們把 \(N\) 稱為數據維度，\(M\) 稱為原子個數，因此，PCA也可以直白的理解為將 \(N\) 維空間中相關性較大的 \(M\) 個原子組成的空間壓縮成相關性最小的 \(K\) 個原子組成的空間。
- 從分類的角度來看，這里我們把 \(N\) 稱為樣本數量，\(M\) 稱為樣本類數，因此，PCA也可以普通的理解為將 \(M\) 個樣本類別重新整理成具有最大區分度的 \(K\) 個類別。
- 從信號表示的角度來看，這里我們把 \(N\) 稱為單個信號長度，\(M\) 稱為信號數量，因此，PCA也可以粗略的理解為將 \(M\) 個很相似的信號通過適當的投影（或者說變換\映射）用 \(K\) 個信號來簡單的表示原來的。就比如說一只貓的叫聲為“喵喵喵~”，另一只貓的叫聲也為“喵喵喵~”，雖然他們的叫聲會存在一定的差異，但是對於人類來說幾乎是一樣的，這兩只貓發出的聲音信號可以直接用一個信號來表示。
- 從特征提取的角度來說，這里我們把 \(N\) 稱為數據維度，\(M\) 稱為特征數量，因此，PCA也可以近似的理解為將數據集中 \(M\) 個特征用 \(K\) 個特征來表示，其中去掉的 \(M-K\) 個特征是無關緊要的。
2.2 PCA直觀思想
為了方便理解，下面以一個具體實例來說明PCA的這種思想。
假設某地的房子面積和總價在理想情況下符合如下正比關系：
|
這些數據在二維平面中可以表示為：
|
然而在實際中並非如此，由於各種原因總會產生一些波動，有可能會是如表2所示的關系。
|
同樣的，這些數據在二維平面中可以表示為：
|
若把表一中的數據當作一個矩陣來看，就變成：
|
從表1中的數據來看，在理想情況下，房子面積和總價符合正比關系，也就是說房子面積和總價具有很大的相關性，若把它當成一個 \(2\times 5\) 的矩陣的話，它的秩為1，這說明用表1中的10個數據來描述房子面積和總價的這種關系顯得有些冗余，我們只需要一行數據（\(5\)個）就可以描述清楚這種關系。在實際情況中，表2中的兩行數據雖然不滿足完全的正比關系，但是它們也非常接近正比關系，因此，他們之間的相關性也是非常大的。那么我們如何將 \(\bf{X}\) 中的數據用一行數（也就是從二維降到一維）來表示呢？由前面相關性分析可知，應該是要進行“去相關性”操作，也就是要使兩行數據通過一個變換，使其相關性變小。
數學里通常用方差和協方差來描述數據的相關性，方差定義如下：
其中 \(E({\bf{x}}) = \frac{1}{N}\sum\limits_{i = 1}^N {\left( {{x_i}} \right)}\)為向量 \(\bf{x}\) 的期望。
很明顯，方差是描述一個向量（維度）內數據的分散程度（相關性）的度量，若如果數據是像表3一樣的矩陣形式呢？那么就用協方差來描述數據之間的相關性：
數量高於2個的 \(N\) 維數據通常用協方差矩陣來表示:
假設數據已經經過中心化處理，數據矩陣 \(\bf{X}\) 的協方差矩陣就可以寫成：
其中，\(\bf{X}_{1}\) 和 \(\bf{X}_{2}\) 分別為 \(\bf{X}\) 的第一行和第二行。
在這個例子中，協方差具體表示什么意義呢？
|
如圖3所示，將各個數據點分別對橫軸和縱軸做投影，其中綠色線與橫軸的交點為數據在橫軸上的投影，藍線與縱軸的交點為數據在縱軸上的投影。那么橫軸上交點的分散程度（方差）就對應了 \(\bf{C}\) 中的 \(cov(\bf{x},\bf{x})\)，而縱軸上的交點的分散程度就對應了 \(\bf{C}\) 中的 \(cov(\bf{y},\bf{y})\)。\(cov( \bf{y},\bf{x}\) 則表示 \(\bf{x}\) 和 \(\bf{y}\) 之間的相關程度。
那么由PCA的中心思想：將高維數據投影到具有最大方差的低維空間中，我們就可以知道只要將數據投影到方差最大的低維（\(K\) 維）空間中就行，那么如何找到方差最大的低維空間呢？如圖4所示，數據不僅可以投影到橫軸和縱軸，還可以360度旋轉，投影到任何直線（1維空間）上，對於圖4中的 \(pc1\) 方向和 \(pc2\) 方向，可以直觀的看出來肯定要投影到 \(pc1\) 方向要好一點，因為投影之后的數據點更分散。
數據投影到一個維度越分散在一定程度上就說明兩類數據的相關性越小，這就指導我們一條思路：只需要將 \(\bf{C_{\bf{X}}}\) 中的兩個協方差 \(cov(\bf{X}_{1},\bf{X}_{2})\) 和 \(cov(\bf{X}_{2},\bf{X}_{1})\) 變成0就可以了，這就意味着將 \(\bf{C_{\bf{X}}}\) 變成對角矩陣就行了。
|
我們發現 \(\bf{C_{\bf{X}}}\) 是一個對稱矩陣，自然就會想到特征值分解。
\(\bf{C_{\bf{X}}}\) 是 \(\bf{M}\) 階對稱矩陣，則必有正交矩陣 \(\bf{P}\)使得：
\({{\bf{P}}^{ - 1}}{{\bf{C}}_{\bf{X}}}{\bf{P}} = {{\bf{P}}^{\rm T}}{{\bf{C}}_{\bf{X}}}{\bf{P}} = {\bf{\Lambda }}\)
其中 \(\bf{\Lambda }\) 是以 \(\bf{C_{\bf{X}}}\) 的 \(\bf{M}\) 個特征值為對角元
因此我們只要找到這個正交矩陣 \(\bf{P}\) 和它對應的特征值，然后取前 \(K\) 個最大的特征值及其對應的特征向量就可以了，這實際上就是進行特征值分解：
當然也可以通過奇異值分解來解決問題：
經過奇異值分解后得到的 \(\bf{\Sigma}\) 中最大的 \(K\) 個奇異值及其對應的 \(\bf{U}\) 中的 \(K\) 個列向量組成的矩陣就是變換矩陣。如果用右奇異矩陣 \(\bf{U}\) 的話就是對 \(\bf{X}\) 的列進行降維了。
求協方差矩陣 \({{\bf{C}}_{\bf{X}}} = {{\bf{X}}}{\bf{X}^{\rm T}}\) (這里實際上是散度矩陣，是協方差矩陣的 \(N-1\) 倍)的復雜度為 \(O\left( {M{N^2}} \right)\)，當數據維度過大時，還是很耗時的，此時用奇異值分解有一定的優勢，比如有一些SVD的實現算法可以先不求出協方差矩陣也能求出右奇異矩陣 \(\bf{V}\)，也就是說，PCA算法可以不用做特征分解而是通過SVD來完成，這個方法在樣本量很大的時候很有效。
至於為什么協方差矩陣經過特征值（奇異值）分解后前 \(K\) 個最大特征值對應的分量就是方差最大的 \(K\) 維空間呢？這里給出一個簡單的理解方法，數學證明就不貼出來了。\(K\) 個特征值與\(K\) 個特征向量構成了一個 \(K\) 維空間，而特征值正是衡量每一維空間所占的分量的，特征值越大，說明這一維度所包含的信息量越大，權重就越大。前 \(K\) 個最大特征值組成的空間當然是權重最大的，從而矩陣 \(\bf{X}\) 投影到這 \(K\) 維空間中足夠分散，也就是說方差足夠大。
2.3 PCA的理論推導
這里引用一下
這里是引用主成分分析（PCA）原理總結這篇博客的證明。
2.4 PCA實例求解
還是以某地房子面積與總價的問題作為例子來總結PCA算法步驟。
步驟1. 數據去中心化：
步驟2：計算散度矩陣(協方差矩陣)
步驟3：對協方差矩陣進行特征值（奇異值分解）：
步驟4：取最大的 \(K\) 個特征值與特征向量構成一個 變換矩陣 \(\bf{P}_K\)（這里 \(K=1\)）：
步驟5：將原始數據投影到 \(K\) 維空間：
我們可以看到最后結果，那條投影線跟我們直觀感受的一樣：
3. PCA的應用
PCA在特征工程、計算機視覺、人臉識別、去噪等領域有很廣泛的應用。這里僅僅測試一下PCA在圖像壓縮上的效果。matlab代碼如下：
function [P,r] = fun_pca(X,K)
% 輸入：
% X：經過中心化處理的數據
% K: 需要降低到的維度
% 輸出：
% P: 投影矩陣
% r:前K個奇異值占所有奇異值的比
Cov = X*X';
[u,s] = svd(Cov) ;
P = u(:,1:K) ;
s = diag(s) ;
r = sum(abs(s(1:K)))/sum(abs(s)) ;
end
PCA for image compression
clear
img = imread('xinhuan.jpg');
img = im2double(img);
K=50;
% 保留主成分個數
img_matrix = [img(:,:,1),img(:,:,2),img(:,:,3)] ;
% 將RGB三通道圖像合並成一個矩陣
img_matrix_mean = img_matrix - mean(img_matrix) ;
% 求均值
Proj = fun_pca(img_matrix_mean,K) ;
% 用PCA算法獲得投影矩陣P
img_K = Proj'*img_matrix_mean ;
% 壓縮到K維后的圖像
%% 圖像重建
img_rec = Proj*img_K ;
% 圖像解壓縮
img_rec = img_rec + mean(img_matrix) ;
% 重建后的圖像加上均值
c_n = size(img_matrix,2)/3 ;
img_rec = cat(3,img_rec(:,1:c_n),img_rec(:,c_n+1:2*c_n),img_rec(:,2*c_n+1:end)) ;
figure(1)
subplot(1,2,1)
imshow(img),title('原圖')
subplot(1,2,2)
imshow(img_rec),title(strcat('壓縮后的圖像(', '奇異值數量：10)'))
SVD for image compression
clear
img = imread('xinhuan.jpg');
img = im2double(img);
K = 10 ;
for i = 1:3
[u,s,v] = svd(img(:,:,i)) ;
s = diag(s) ;
r = sum(abs(s(1:K)))/sum(abs(s)) ;
s = diag(s(1:K)) ;
W(:,:,i) = u(:,1:K)*s*v(:,1:K)' ;
end
figure(1)
subplot(2,2,1)
imshow(img),title('原圖')
subplot(2,2,2)
imshow(P),title(strcat('壓縮后的圖像(', '奇異值數量：10)'))
奇異值分解圖像壓縮：
|
PCA圖像壓縮：
|
從圖5和圖6基本看不出來PCA和SVD壓縮圖像有多大區別，其特征值貢獻比差異很大確不能說明什么問題，圖像重建質量可以使用峰值信噪比（PSNR）來評價，感興趣的可以去試一試。
4. 討論
PCA的主要優點有：
1）僅僅需要以方差衡量信息量，不受數據集以外的因素影響。
2）各主成分之間正交，可消除原始數據成分間的相互影響的因素。
3）計算方法簡單，主要運算是特征值分解，易於實現。
PCA算法的主要缺點有：
1）主成分各個特征維度的含義具有一定的模糊性，不如原始樣本特征的解釋性強。
2）方差小的非主成分也可能含有對樣本差異的重要信息，因降維丟棄可能對后續數據處理有影響。
更多精彩內容請關注訂閱號優化與算法和加入QQ討論群1032493483獲取更多資料
