線性降維：主成分分析原理及仿真

時間 2020-04-30

標籤線性成分分析原理仿真欄目應用數學简体版

原文原文鏈接

鶯嘴啄花紅溜，燕尾點波綠皺。
指冷玉笙寒，吹徹小梅春透。
依舊，依舊，人與綠楊俱瘦。
——《如夢令·春景》秦觀

更多精彩內容請關注微信公衆號「優化與算法」html

一、背景

隨着信息技術的發展，數據量呈現爆照式增加，高維海量數據給傳統的數據處理方法帶來了嚴峻的挑戰，所以，開發高效的數據處理技術是很是必要的。數據降維是解決維度災難的一種有效方法，之因此對數據進行降維是由於：算法

在原始的高維空間中，包含有冗餘信息以及噪音信息，在實際應用例如圖像識別中形成了偏差，下降了準確率；而經過降維,咱們但願減小冗餘信息所形成的偏差，提升識別（或其餘應用）的精度。
但願經過降維算法來尋找數據內部的本質結構特徵。

主成分分析（Principal Component Analysis，PCA）是一種常見的線性降維方法，普遍應用於圖像處理、人臉識別、數據壓縮、信號去噪等領域。微信

二、算法思想與原理

設原數據大小爲 \(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的這種思想。
假設某地的房子面積和總價在理想狀況下符合以下正比關係：spa

表1. 理想房子面積與總價關係

這些數據在二維平面中能夠表示爲：3d

圖1. 理想房子面積與總價關係

然而在實際中並不是如此，因爲各類緣由總會產生一些波動，有可能會是如表2所示的關係。code

表2. 實際房子面積與總價關係

一樣的，這些數據在二維平面中能夠表示爲：htm

圖2. 實際房子面積與總價關係

若把表一中的數據看成一個矩陣來看，就變成：blog

表3. 表1中數據的矩陣形式

從表1中的數據來看，在理想狀況下，房子面積和總價符合正比關係，也就是說房子面積和總價具備很大的相關性，若把它當成一個 \(2\times 5\) 的矩陣的話，它的秩爲1，這說明用表1中的10個數據來描述房子面積和總價的這種關係顯得有些冗餘，咱們只須要一行數據（\(5\)個）就能夠描述清楚這種關係。在實際狀況中，表2中的兩行數據雖然不知足徹底的正比關係，可是它們也很是接近正比關係，所以，他們之間的相關性也是很是大的。那麼咱們如何將 \(\bf{X}\) 中的數據用一行數（也就是從二維降到一維）來表示呢？由前面相關性分析可知，應該是要進行「去相關性」操做，也就是要使兩行數據經過一個變換，使其相關性變小。

數學裏一般用方差和協方差來描述數據的相關性，方差定義以下：

\[Var({\bf{x}}) = \frac{1}{N}\sum\limits_{i = 1}^N {{{\left( {{x_i} - E({\bf{x}})} \right)}^2}} \]

其中 \(E({\bf{x}}) = \frac{1}{N}\sum\limits_{i = 1}^N {\left( {{x_i}} \right)}\)爲向量 \(\bf{x}\) 的指望。

很明顯，方差是描述一個向量（維度）內數據的分散程度（相關性）的度量，若若是數據是像表3同樣的矩陣形式呢？那麼就用協方差來描述數據之間的相關性：

\[Cov( \bf{x},\bf{y}) \ =\ \frac{1}{N}\sum ^{N}_{i=1}( x_{i} \ -E(\bf{x}))( y_{i} -E(\bf{y})) =E[( \bf{x}-E(\bf{x}))( \bf{y}-E(\bf{y}))\]

數量高於2個的 \(N\) 維數據一般用協方差矩陣來表示:

\[ \bf{C}\ =\begin{pmatrix} cov(\bf{x},\bf{x}) & cov( \bf{x},\bf{y})\\ cov( \bf{y},\bf{x}) & cov( \bf{y},\bf{y}) \end{pmatrix}\]

假設數據已經通過中心化處理，數據矩陣 \(\bf{X}\) 的協方差矩陣就能夠寫成：

\[ \bf{C_{\bf{X}}}\ =\begin{pmatrix} cov(\bf{X}_{1},\bf{X}_{1}) & cov(\bf{X}_{1},\bf{X}_{2})\\ cov( \bf{X}_{2},\bf{X}_{1}) & cov( \bf{X}_{2},\bf{X}_{2}) \end{pmatrix}\]

其中，\(\bf{X}_{1}\) 和 \(\bf{X}_{2}\) 分別爲 \(\bf{X}\) 的第一行和第二行。
在這個例子中，協方差具體表示什麼意義呢？

圖3. 協方差的意義

如圖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}}}\) 變成對角矩陣就好了。

圖4. 如何選擇投影方向？

咱們發現 \(\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{C}}_{\bf{X}}} = {{\bf{P}}^{\rm T}}{\bf{\Lambda P}} \]

固然也能夠經過奇異值分解來解決問題：

\[{{\bf{C}}_{\bf{X}}} = {\bf{U\Sigma V}} \]

通過奇異值分解後獲得的 \(\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\) 維空間：

\[{\bf{Y}} = {\bf{P}}_K^{\rm T}{\bf{X}} \]

咱們能夠看到最後結果，那條投影線跟咱們直觀感覺的同樣：

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)'))

奇異值分解圖像壓縮：

圖5. SVD重構的圖像與原圖對比

PCA圖像壓縮：

圖6. PCA重構的圖像與原圖對比

從圖5和圖6基本看不出來PCA和SVD壓縮圖像有多大區別，其特徵值貢獻比差別很大確不能說明什麼問題，圖像重建質量可使用峯值信噪比（PSNR）來評價，感興趣的能夠去試一試。