ICA快速算法原理和程序

本文格式為Word版，下載可任意編輯——ICA快速算法原理和程序

試驗2：FastICA算法

一．算法原理：

獨(dú)立分量分析（ICA）的過程如下圖：在信源s(t)中各分量相互獨(dú)立的假設(shè)下，由觀測x(t)通過解混系統(tǒng)B把他們分開開來，使輸出y(t)迫近s(t)！

圖1-ICA的一般過程

ICA算法的研究可分為基于信息論準(zhǔn)則的迭代估計方法和基于統(tǒng)計學(xué)的代數(shù)方法兩大類，從原理上來說，它們都是利用了源信號的獨(dú)立性和非高斯性?；谛畔⒄摰姆椒ㄑ芯恐?，各國學(xué)者從最大熵、最小互信息、最大似然和負(fù)熵最大化等角度提出了一系列估計算法。如FastICA算法,Infomax算法，最大似然估計算法等?；诮y(tǒng)計學(xué)的方法主要有二階累積量、四階累積量等高階累積量方法。本試驗探討FastICA算法。

1.數(shù)據(jù)的預(yù)處理

一般狀況下，所獲得的數(shù)據(jù)都具有相關(guān)性，所以尋常都要求對數(shù)據(jù)進(jìn)行初步的白化或球化處理，由于白化處理可去除各觀測信號之間的相關(guān)性，從而簡化了后續(xù)獨(dú)立分量的提取過程，而且，尋常狀況下，數(shù)據(jù)進(jìn)行白化處理與不對數(shù)據(jù)進(jìn)行白化處理相比，算法的收斂性較好。

若一零均值的隨機(jī)向量Z??Z1,?,ZM?滿足EZZT?I，其中：I為單位矩陣，我

T??們稱這個向量為白化向量。白化的本質(zhì)在于去相關(guān)，這同主分量分析的目標(biāo)是一樣的。在ICA

中，對于為零均值的獨(dú)立源信號S?t???S1?t?,...,SN?t??T，有：

且協(xié)方差矩陣是單位陣cov?S??I，因此，源信號S?t?E?SiSj??E?Si?E?Sj??0,當(dāng)i?j，

是白色的。對觀測信號X?t?，我們應(yīng)當(dāng)尋覓一個線性變換，使X?t?投影到新的子空間后變成白化向量，即：

Z?t??W0X?t?（2.1）

其中，W0為白化矩陣，Z為白化向量。

利用主分量分析，我們通過計算樣本向量得到一個變換

W0???1/2UT

其中U和?分別代表協(xié)方差矩陣CX的特征向量矩陣和特征值矩陣?？梢宰C明，線性變換W0滿足白化變換的要求。通過正交變換，可以保證UTU?UUT?I。因此，協(xié)方差矩陣：EZZT?E??1/2UTXXTU??1/2???1/2UTEXXTU??1/2???1/2???1/2?I（2.2）再將X?t??AS?t?式代入Z?t??W0X?t?，且令W0A?A,有

Z?t??W0AS?t??AS?t?（2.3）

由于線性變換A連接的是兩個白色隨機(jī)矢量Z?t?和S?t?，可以得出A一定是一個正交

??????~~~~變換。假使把上式中的Z?t?看作新的觀測信號，那么可以說，白化使原來的混合矩陣A簡化成一個新的正交矩陣A。證明也是簡單的：EZZT?EASSTAT?AESSTAT?AAT?I（2.4）其實(shí)正交變換相當(dāng)于對多維矢量所在的坐標(biāo)系進(jìn)行一個旋轉(zhuǎn)。

在多維狀況下，混合矩陣A是N?N的，白化后新的混合矩陣A由于是正交矩陣，其

~???~~?~??~~~~自由度降為N??N?1?/2，所以說白化使得ICA問題的工作量幾乎減少了一半。

白化這種常規(guī)的方法作為ICA的預(yù)處理可以有效地降低問題的繁雜度，而且算法簡單，

用傳統(tǒng)的PCA就可完成。用PCA對觀測信號進(jìn)行白化的預(yù)處理使得原來所求的解混合矩陣退化成一個正交陣，減少了ICA的工作量。此外，PCA本身具有降維功能，當(dāng)觀測信號的個數(shù)大于源信號個數(shù)時，經(jīng)過白化可以自動將觀測信號數(shù)目降到與源信號維數(shù)一致。

2.FastICA算法

FastICA算法，又稱固定點(diǎn)(Fixed-Point)算法，是由芬蘭赫爾辛基大學(xué)Hyv?rinen等人提出來的。是一種快速尋優(yōu)迭代算法，與普通的神經(jīng)網(wǎng)絡(luò)算法不同的是這種算法采用了批處理的方式，即在每一步迭代中有大量的樣本數(shù)據(jù)參與運(yùn)算。但是從分布式并行處理的觀點(diǎn)看該算法仍可稱之為是一種神經(jīng)網(wǎng)絡(luò)算法。FastICA算法有基于峭度、基于似然最大、基于負(fù)熵最大等形式，這里，我們介紹基于負(fù)熵最大的FastICA算法。它以負(fù)熵最大作為一個搜尋

方向，可以實(shí)現(xiàn)順序地提取獨(dú)立源，充分表達(dá)了投影追蹤（ProjectionPursuit）這種傳統(tǒng)線性變換的思想。此外，該算法采用了定點(diǎn)迭代的優(yōu)化算法，使得收斂更加快速、穩(wěn)健。

由于FastICA算法以負(fù)熵最大作為一個搜尋方向，因此先探討一下負(fù)熵判決準(zhǔn)則。由信息論理論可知：在所有等方差的隨機(jī)變量中，高斯變量的熵最大，因而我們可以利用熵來度量非高斯性，常用熵的修正形式，即負(fù)熵。根據(jù)中心極限定理，若一隨機(jī)變量X由大量相互獨(dú)立的隨機(jī)變量Si?i?1,2,3,...N?之和組成，只要Si具有有限的均值和方差，則不管其為何種分布，隨機(jī)變量X較Si更接近高斯分布。換言之，Si較X的非高斯性更強(qiáng)。因此，在分開過程中，可通過對分開結(jié)果的非高斯性度量來表示分開結(jié)果間的相互獨(dú)立性，當(dāng)非高斯性度量達(dá)到最大時，則說明已完成對各獨(dú)立分量的分開。

負(fù)熵的定義：

Ng?Y??H?YGaus?s?H?Y?（2.5）式中，YGauss是一與Y具有一致方差的高斯隨機(jī)變量，H???為隨機(jī)變量的微分熵

H?Y???pY???lgpY???d?（2.6）

?根據(jù)信息理論，在具有一致方差的隨機(jī)變量中，高斯分布的隨機(jī)變量具有最大的微分熵。當(dāng)Y具有高斯分布時，Ng?Y??0；Y的非高斯性越強(qiáng)，其微分熵越小，Ng?Y?值越大，所以Ng?Y?可以作為隨機(jī)變量Y非高斯性的測度。由于根據(jù)式（3.6）計算微分熵需要知道Y的概率密度分布函數(shù)，這顯然不切實(shí)際，于是采用如下近似公式：Ng?Y???E?g?Y???E?g?YGauss???2（2.7）

其中，E???為均值運(yùn)算；g???為非線性函數(shù)，可取g1?y??tanh(a1y)，或這里，尋常我們?nèi)1?1。1?a1?2，g2?y??yexp?y2/2或g3?y??y3等非線性函數(shù)，

快速ICA學(xué)習(xí)規(guī)則是找一個方向以便WXY?WX具有最大的非高斯性。這里，

T非高斯性用式（3.7）給出的負(fù)熵Ng(WX)的近似值來度量,WX的方差約束為1，對于

TT??T?T?白化數(shù)據(jù)而言，這等于約束W的范數(shù)為1。FastICA算法的推導(dǎo)如下。首先，WX的負(fù)熵的最大近似值能通過對EGWX??T??進(jìn)行優(yōu)化來獲得。根據(jù)????Kuhn-Tucker條件，在

EWTX

?????W22?1的約束下，EGWTX的最優(yōu)值能在滿足下式的點(diǎn)上獲得。

EXgWTX??W?0（2.8）這里，?是一個恒定值，

??????EW0TXgW0TX,W0是優(yōu)化后的W值。下面我們利用牛

????頓迭代法解方程（3.8）。用F表示式（3.8）左邊的函數(shù)，可得F的雅可比矩陣JF?W?如下：JF?W??EXXTg'WTX??I（2.9）為了簡化矩陣的求逆，可以近似為（3.9）式的第一項。由于數(shù)據(jù)被球化，EXXT?I,所以，EXXTg'WTX???????????E?XX??E?g'?WX???E?g'?WX??I。因而雅可比矩陣變成了

TTT對角陣，并且能比較簡單地求逆。因而可以得到下面的近似牛頓迭代公式：

?W??W?EXgWTX??W/Eg'WTX??W?W/W??????????????（2.10）

這里，W是W的新值，??EWTXgWTX，規(guī)格化能提高解的穩(wěn)定性。簡化后就可以得到FastICA算法的迭代公式：

????W??EXgWTX?Eg'WTXWW?W/W??????????(2.11)

實(shí)踐中，F(xiàn)astICA算法中用的期望必需用它們的估計值代替。當(dāng)然最好的估計是相應(yīng)的樣本平均。理想狀況下，所有的有效數(shù)據(jù)都應(yīng)當(dāng)參與計算，但這會降低計算速度。所以尋常用一部分樣本的平均來估計，樣本數(shù)目的多少對最終估計的確切度有很大影響。迭代中的樣本點(diǎn)應(yīng)當(dāng)分別選取，假使收斂不理想的話，可以增加樣本的數(shù)量。

3.FastICA算法的基本步驟：

1.對觀測數(shù)據(jù)X進(jìn)行中心化，使它的均值為0；2.對數(shù)據(jù)進(jìn)行白化，X?Z。

3.選擇需要估計的分量的個數(shù)m，設(shè)迭代次數(shù)p?14.選擇一個初始權(quán)矢量（隨機(jī)的）Wp。

5.令Wp?EZgWpZ?Eg'WpZW,非線性函數(shù)g的選取見前文。6.Wp?Wp???p?1j?1T?????T????WTpWjWj。

7.令Wp?Wp/Wp。

8.假使Wp不收斂的話，返回第5步。9．令p?p?1，假使p?m，返回第4步。

二．MATLAB源程序及說明：

%下程序為ICA的調(diào)用函數(shù)，輸入為觀測的信號，輸出為解混后的信號functionZ=ICA(X)

%去均值

[M,T]=size(X);%獲取輸入矩陣的行/列數(shù)，行數(shù)為觀測數(shù)據(jù)的數(shù)目，列數(shù)為采樣點(diǎn)數(shù)

average=mean(X')';%均值fori=1:M

X(i,:)=X(i,:)-average(i)*ones(1,T);end

%白化/球化

Cx=cov(X',1);%計算協(xié)方差矩陣Cx

[eigvector,eigvalue]=eig(Cx);%計算Cx的特征值和特征向量W=eigvalue^(-1/2)*eigvector';%白化矩陣Z=W*X;%正交矩陣

%迭代

Maxcount=10000;%最大迭代次數(shù)Critical=0.00001;%判斷是否收斂

m=M;%需要估計的分量的個數(shù)W=rand(m);forn=1:m

WP=W(:,n);%初始權(quán)矢量（任意）%Y=WP'*Z;

%G=Y.^3;%G為非線性函數(shù)，可取y^3等%GG=3*Y.^2;%G的導(dǎo)數(shù)count=0;

LastWP=zeros(m,1);

W(:,n)=W(:,n)/norm(W(