ICA快速算法原理與程序

1、實(shí)驗(yàn)2：FastICA算法一算法原理：獨(dú)立分量分析（ICA）的過程如下圖所示：在信源中各分量相互獨(dú)立的假設(shè)下，由觀察通過結(jié)婚系統(tǒng)把他們分離開來，使輸出逼近。圖1-ICA的一般過程ICA算法的研究可分為基于信息論準(zhǔn)則的迭代估計(jì)方法和基于統(tǒng)計(jì)學(xué)的代數(shù)方法兩大類，從原理上來說，它們都是利用了源信號的獨(dú)立性和非高斯性?；谛畔⒄摰姆椒ㄑ芯恐校鲊鴮W(xué)者從最大熵、最小互信息、最大似然和負(fù)熵最大化等角度提出了一系列估計(jì)算法。如FastICA算法, Infomax算法，最大似然估計(jì)算法等?；诮y(tǒng)計(jì)學(xué)的方法主要有二階累積量、四階累積量等高階累積量方法。本實(shí)驗(yàn)主要討論FastICA算法。1. 數(shù)據(jù)的預(yù)處理 一般情

2、況下，所獲得的數(shù)據(jù)都具有相關(guān)性，所以通常都要求對數(shù)據(jù)進(jìn)行初步的白化或球化處理，因?yàn)榘谆幚砜扇コ饔^測信號之間的相關(guān)性，從而簡化了后續(xù)獨(dú)立分量的提取過程，而且，通常情況下，數(shù)據(jù)進(jìn)行白化處理與不對數(shù)據(jù)進(jìn)行白化處理相比，算法的收斂性較好。若一零均值的隨機(jī)向量滿足，其中：為單位矩陣，我們稱這個(gè)向量為白化向量。白化的本質(zhì)在于去相關(guān)，這同主分量分析的目標(biāo)是一樣的。在ICA中，對于為零均值的獨(dú)立源信號，有：，且協(xié)方差矩陣是單位陣，因此，源信號是白色的。對觀測信號，我們應(yīng)該尋找一個(gè)線性變換，使投影到新的子空間后變成白化向量，即： （2.1）其中，為白化矩陣，為白化向量。利用主分量分析，我們通過計(jì)算樣本向量得

3、到一個(gè)變換其中和分別代表協(xié)方差矩陣的特征向量矩陣和特征值矩陣?？梢宰C明，線性變換滿足白化變換的要求。通過正交變換，可以保證。因此，協(xié)方差矩陣： （2.2）再將式代入，且令,有 （2.3）由于線性變換連接的是兩個(gè)白色隨機(jī)矢量和，可以得出一定是一個(gè)正交變換。如果把上式中的看作新的觀測信號，那么可以說，白化使原來的混合矩陣簡化成一個(gè)新的正交矩陣。證明也是簡單的： （2.4）其實(shí)正交變換相當(dāng)于對多維矢量所在的坐標(biāo)系進(jìn)行一個(gè)旋轉(zhuǎn)。在多維情況下，混合矩陣是的，白化后新的混合矩陣由于是正交矩陣，其自由度降為，所以說白化使得ICA問題的工作量幾乎減少了一半。白化這種常規(guī)的方法作為ICA的預(yù)處理可以有效地降低問

4、題的復(fù)雜度，而且算法簡單，用傳統(tǒng)的PCA就可完成。用PCA對觀測信號進(jìn)行白化的預(yù)處理使得原來所求的解混合矩陣退化成一個(gè)正交陣，減少了ICA的工作量。此外，PCA本身具有降維功能，當(dāng)觀測信號的個(gè)數(shù)大于源信號個(gè)數(shù)時(shí)，經(jīng)過白化可以自動(dò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ò)算法。

5、FastICA算法有基于峭度、基于似然最大、基于負(fù)熵最大等形式，這里，我們介紹基于負(fù)熵最大的FastICA算法。它以負(fù)熵最大作為一個(gè)搜尋方向，可以實(shí)現(xiàn)順序地提取獨(dú)立源，充分體現(xiàn)了投影追蹤（Projection Pursuit）這種傳統(tǒng)線性變換的思想。此外，該算法采用了定點(diǎn)迭代的優(yōu)化算法，使得收斂更加快速、穩(wěn)健。因?yàn)镕astICA算法以負(fù)熵最大作為一個(gè)搜尋方向，因此先討論一下負(fù)熵判決準(zhǔn)則。由信息論理論可知：在所有等方差的隨機(jī)變量中，高斯變量的熵最大，因而我們可以利用熵來度量非高斯性，常用熵的修正形式，即負(fù)熵。根據(jù)中心極限定理，若一隨機(jī)變量由許多相互獨(dú)立的隨機(jī)變量之和組成，只要具有有限的均值和方差

6、，則不論其為何種分布，隨機(jī)變量較更接近高斯分布。換言之，較的非高斯性更強(qiáng)。因此，在分離過程中，可通過對分離結(jié)果的非高斯性度量來表示分離結(jié)果間的相互獨(dú)立性，當(dāng)非高斯性度量達(dá)到最大時(shí)，則表明已完成對各獨(dú)立分量的分離。負(fù)熵的定義： （2.5）式中，是一與具有相同方差的高斯隨機(jī)變量，為隨機(jī)變量的微分熵 （2.6） 根據(jù)信息理論，在具有相同方差的隨機(jī)變量中，高斯分布的隨機(jī)變量具有最大的微分熵。當(dāng)具有高斯分布時(shí)，；的非高斯性越強(qiáng)，其微分熵越小，值越大，所以可以作為隨機(jī)變量非高斯性的測度。由于根據(jù)式（3.6）計(jì)算微分熵需要知道的概率密度分布函數(shù)，這顯然不切實(shí)際，于是采用如下近似公式： （2.7）其中，為均值

7、運(yùn)算；為非線性函數(shù)，可取，或或等非線性函數(shù)，這里，通常我們?nèi)?。快速ICA學(xué)習(xí)規(guī)則是找一個(gè)方向以便具有最大的非高斯性。這里，非高斯性用式（3.7）給出的負(fù)熵的近似值來度量,的方差約束為1，對于白化數(shù)據(jù)而言，這等于約束的范數(shù)為1。FastICA算法的推導(dǎo)如下。首先，的負(fù)熵的最大近似值能通過對進(jìn)行優(yōu)化來獲得。根據(jù)Kuhn-Tucker條件，在的約束下，的最優(yōu)值能在滿足下式的點(diǎn)上獲得。 （2.8）這里，是一個(gè)恒定值， ,是優(yōu)化后的值。下面我們利用牛頓迭代法解方程（3.8）。用表示式（3.8）左邊的函數(shù)，可得的雅可比矩陣如下： （2.9）為了簡化矩陣的求逆，可以近似為（3.9）式的第一項(xiàng)。由于數(shù)據(jù)被球化

8、，,所以，。因而雅可比矩陣變成了對角陣，并且能比較容易地求逆。因而可以得到下面的近似牛頓迭代公式： （2.10）這里，是的新值，規(guī)格化能提高解的穩(wěn)定性。簡化后就可以得到FastICA算法的迭代公式： (2.11)實(shí)踐中，F(xiàn)astICA算法中用的期望必須用它們的估計(jì)值代替。當(dāng)然最好的估計(jì)是相應(yīng)的樣本平均。理想情況下，所有的有效數(shù)據(jù)都應(yīng)該參與計(jì)算，但這會(huì)降低計(jì)算速度。所以通常用一部分樣本的平均來估計(jì)，樣本數(shù)目的多少對最后估計(jì)的精確度有很大影響。迭代中的樣本點(diǎn)應(yīng)該分別選取，假如收斂不理想的話，可以增加樣本的數(shù)量。3. FastICA算法的基本步驟：1. 對觀測數(shù)據(jù)進(jìn)行中心化，使它的均值為0；2. 對

9、數(shù)據(jù)進(jìn)行白化，。3. 選擇需要估計(jì)的分量的個(gè)數(shù)，設(shè)迭代次數(shù)4. 選擇一個(gè)初始權(quán)矢量（隨機(jī)的）。5. 令,非線性函數(shù)的選取見前文。6. 。7. 令。8. 假如不收斂的話，返回第5步。9令，如果，返回第4步。二MATLAB源程序及說明：%下程序?yàn)镮CA的調(diào)用函數(shù)，輸入為觀察的信號，輸出為解混后的信號function Z=ICA(X)%-去均值-M,T = size(X); %獲取輸入矩陣的行/列數(shù)，行數(shù)為觀測數(shù)據(jù)的數(shù)目，列數(shù)為采樣點(diǎn)數(shù) average= mean(X')' %均值for i=1:M X(i,:)=X(i,:)-average(i)*ones(1,T); end%-白

10、化/球化-Cx = cov(X',1); %計(jì)算協(xié)方差矩陣Cxeigvector,eigvalue = eig(Cx); %計(jì)算Cx的特征值和特征向量W=eigvalue(-1/2)*eigvector' %白化矩陣Z=W*X; %正交矩陣 %-迭代-Maxcount=10000; %最大迭代次數(shù)Critical=0.00001; %判斷是否收斂m=M; %需要估計(jì)的分量的個(gè)數(shù)W=rand(m);for n=1:m WP=W(:,n); %初始權(quán)矢量（任意）% Y=WP'*Z;% G=Y.3;%G為非線性函數(shù)，可取y3等% GG=3*Y.2; %G的導(dǎo)數(shù) count=0

11、; LastWP=zeros(m,1); W(:,n)=W(:,n)/norm(W(:,n); while abs(WP-LastWP)&abs(WP+LastWP)>Critical count=count+1; %迭代次數(shù) LastWP=WP; %上次迭代的值 % WP=1/T*Z*(LastWP'*Z).3)'-3*LastWP; for i=1:m WP(i)=mean(Z(i,:).*(tanh(LastWP)'*Z)-(mean(1-(tanh(LastWP)'*Z).2).*LastWP(i); end WPP=zeros(m,1);

12、 for j=1:n-1 WPP=WPP+(WP'*W(:,j)*W(:,j); end WP=WP-WPP; WP=WP/(norm(WP); if count=Maxcount fprintf('未找到相應(yīng)的信號); return; end end W(:,n)=WP;endZ=W'*Z;%以下為主程序，主要為原始信號的產(chǎn)生，觀察信號和解混信號的作圖clear all;clc;N=200;n=1:N;%N為采樣點(diǎn)數(shù)s1=2*sin(0.02*pi*n);%正弦信號t=1:N;s2=2*square(100*t,50);%方波信號a=linspace(1,-1,25)

13、;s3=2*a,a,a,a,a,a,a,a;%鋸齒信號s4=rand(1,N);%隨機(jī)噪聲S=s1;s2;s3;s4;%信號組成4*NA=rand(4,4);X=A*S;%觀察信號 %源信號波形圖figure(1);subplot(4,1,1);plot(s1);axis(0 N -5,5);title('源信號');subplot(4,1,2);plot(s2);axis(0 N -5,5);subplot(4,1,3);plot(s3);axis(0 N -5,5);subplot(4,1,4);plot(s4);xlabel('Time/ms');%觀察信

14、號（混合信號）波形圖figure(2);subplot(4,1,1);plot(X(1,:);title('觀察信號（混合信號）');subplot(4,1,2);plot(X(2,:);subplot(4,1,3);plot(X(3,:);subplot(4,1,4);plot(X(4,:); Z=ICA(X); figure(3);subplot(4,1,1);plot(Z(1,:);title('解混后的信號');subplot(4,1,2);plot(Z(2,:);subplot(4,1,3);plot(Z(3,:);subplot(4,1,4);plo

15、t(Z(4,:);xlabel('Time/ms'); 三實(shí)驗(yàn)結(jié)果： 實(shí)驗(yàn)結(jié)果如下所示：其中圖2為源信號的波形圖，圖3為觀察信號（混合信號）波形圖，圖4為解混后的信號波形圖。從圖4可以看出，執(zhí)行ICA后，可以將圖2中的4種信號分離出來，且誤差較小（個(gè)別信號發(fā)生反相）。與實(shí)驗(yàn)一相比，分離效果要提高很多。圖2-源信號圖3-觀察信號（混合信號）圖4-解混后的信號四源信號的分布特性對分離效果的影響：下列3個(gè)圖為源信號中含有不同個(gè)數(shù)的高斯白噪聲分解后得到的相應(yīng)解混信號：其中圖5源信號只含一個(gè)隨機(jī)噪聲，圖6源信號含兩個(gè)隨機(jī)噪聲，圖7源信號含三個(gè)隨機(jī)噪聲。從圖5,6,7可以看出，源信號所含的

16、高斯白噪聲越多，分離后得到的信號與源信號相比誤差越大，效果越差；所含高斯白噪聲越少，分離效果越好。圖5-源信號只含一個(gè)隨機(jī)噪聲分離后得到的波形圖圖6-源信號含兩個(gè)隨機(jī)噪聲分離后得到的波形圖圖7-源信號含三個(gè)隨機(jī)噪聲分離后得到的波形圖下圖8,9,10分別為方波信號，正弦信號和鋸齒波3種信號與同一個(gè)噪聲混合后，經(jīng)ICA解混后得到的結(jié)果。為便于比較，源信號中3種信號所含的能量相等。信號的非高斯性可用四階累積量（峰度）進(jìn)行描述，>0稱為超高斯型，<0稱為亞高斯型，可用的大小作為信號距離高斯型程度的度量。越大，表明信號距高斯型越遠(yuǎn)，即信號的非高斯性越強(qiáng)。下圖中3種信號在源信號中的四階累積量分別為：方波為-8.0，正弦波為-6.0，鋸齒波為-4.8154，表明方波