1、 实验 2： FastICA算法 一算法原理 ： 独立分量分析 （ ICA） 的过程如下图所示：在信源 ()st 中各分量相互独立的假设下，由观察 ()xt 通过结婚系统 B 把他们分离开来，使输出 ()yt 逼近 ()st 。 图 1-ICA的一般过程 ICA 算法的研究可分为基于信息论准则的迭代估计方法 和基于统计学的代数方法两大类，从原理上来说，它们都是利用了源信号的独立性和非高斯性。 基于信息论的方法研究中，各国学者从最大熵、最小互信息、最大似然和负熵最大化等角度提出了一系列估计算法。如FastICA算法 , Infomax 算法 ，最大似然估计算法等 。基于统计学的方法主要有二阶累积
2、量、四阶累积量等高阶累积量方法。 本实验主要讨论 FastICA算法 。 1. 数据的预处理 一般情况下，所获得的数据都具有相关性，所以通常都要求对数据进行初步的白化或球化处理，因为白化处理可去除各观测信号之间的相关性，从而简化 了后续独立分量的提取过程，而且，通常情况下，数据进行白化处理与不对数据进行白化处理相比，算法的收敛性较好。 若一零均值的随机向量 TMZZZ ,1 满足 IZZE T ，其中： I 为单位矩阵，我们称这个向量为白化向量。白化的本质在于去相关，这同主分量分析的目标是一样的。在ICA 中 ， 对 于 为 零 均 值 的 独 立 源 信 号 TN tStStS ,.,1 ，
3、有： jiSESESSE jiji 当,0，且协方 差矩阵是单位阵 IS cov ，因此，源信号 tS是白色的。对观测信号 tX ，我们应该寻找一个线性变换，使 tX 投影到新的子空间后变成白化向量，即： tXWtZ 0 （ 2.1） 其中， 0W 为白化矩阵， Z 为白化向量。 利用主分量分析，我们通过计算样本向量得到一个变换 TUW 2/10 其中 U 和 分别代表协方差矩阵 XC 的特征向量矩阵和特征值矩阵。可以证明，线性变换0W 满足白化变换的要求。通过正交变换，可以保证 IUUUU TT 。因此，协方差矩阵： IUXXEUUXXUEZZE TTTTT 2/12/12/12/12/12
4、/1 （ 2.2） 再将 tAStX 式代入 tXWtZ 0 ，且令 AAW 0 ,有 tSAtASWtZ 0 （ 2.3） 由于线性变换 A 连接的是两个白色随机矢量 tZ 和 tS ，可以得出 A 一定是一个正交变换。如果把上式中的 tZ 看作新的观测信号，那么可以说，白化使原来的混合矩阵 A 简化成一个新的正交矩阵 A 。证明也是简单的： IAAASSEAASSAEZZE TTTTTT （ 2.4） 其实正交变换相当于对多维矢量所在的坐标系进行一个旋转。 在多维情况下，混合矩阵 A 是 NN 的，白化后新的混合矩阵 A 由于是正交矩阵，其自由度降为 2/1 NN ，所以说白化使得 ICA
5、问题的工作量几乎减少了一半。 白化这种常规的方法作为 ICA的预处理可以有效地降低问题的复杂度，而且算法简单，用传统的 PCA就可完成。用 PCA对观测信号进行白 化的预处理使得原来所求的解混合矩阵退化成一个正交阵，减少了 ICA的工作量。此外， PCA本身具有降维功能，当观测信号的个数大于源信号个数时，经过白化可以自动将观测信号数目降到与源信号维数相同。 2. FastICA算法 FastICA 算法，又称固定点 (Fixed-Point)算法， 是由芬兰赫尔辛基大学 Hyvrinen 等人提出来的。 是一种快速寻优迭代算法，与普通的神经网络算法不同的是这种算法采用了批处理的方式，即在每一步
6、迭代中有大量的样本数据参与运算。但是从分布式并行处理的观点看该算法仍可称之为是 一种神经网络算法。 FastICA算法有基于峭度、基于似然最大、基于负熵最大等形式，这里，我们介绍基于负熵最大的 FastICA算法。它以负熵最大作为一个搜寻方向，可以实现顺序地提取独立源，充分体现了投影追踪（ Projection Pursuit）这种传统线性变换的思想。此外，该算法采用了定点迭代的优化算法，使得收敛更加快速、稳健。 因为 FastICA 算法以负熵最大作为一个搜寻方向，因此先讨论一下负熵判决准则。由信息论理论可知：在所有等方差的随机变量中，高斯变量的熵最大，因而我们可以利用熵来度量非高斯性，常用
7、熵 的修正形式，即负熵。根据中心极限定理，若一随机变量 X 由许多相互独立的随机变量 NiSi ,.3,2,1 之和组成，只要 iS 具有有限的均值和方差，则不论其为何种分布，随机变量 X 较 iS 更接近高斯分布。换言之， iS 较 X 的非高斯性更强。因此，在 分离过程中，可通过对分离结果的非高斯性度量来表示分离结果间的相互独立性，当非高斯性度量达到最大时，则表明已完成对各独立分量的分离。 负熵的定义： YHYHYN G a u s sg （ 2.5） 式中， GaussY 是一与 Y 具有相同方差的高斯随机变量， H 为随机变量的微分熵 dppYHYY lg（ 2.6） 根据信息理论，在
8、具有相同方差的随机变量中，高斯分布的随机变量具有最大的微分熵。当 Y 具有高斯分布时， 0YNg ； Y 的非高斯性越强，其微分熵越小， YNg 值越大，所以 YNg 可以作为随机变量 Y 非高斯性的测度。由于根据式（ 3.6）计算微分熵需要知道 Y 的概率密度分布函数，这显然不切实际，于是采用如下近似公式 ： 2G a u s sg YgEYgEYN （ 2.7） 其中， E 为均值运算； g 为非线性函数，可取 )tanh( 11 yayg ，或 2/e x p 22 yyyg 或 33 yyg 等非线性函数，这里， 21 1a ，通常我们取 11a 。 快速 ICA学习规则是找一个方向以
9、便 XWYXW TT 具有最大的非高斯性。这里，非高斯性用式（ 3.7）给出的负熵 )( XWN Tg 的近似值来度量 , XWT 的方差 约束为 1，对于白化数据而言，这等于约束 W 的范数为 1。 FastICA算法的推导如下。首先， XWT 的负熵的最大近似值能通过对 XWGE T 进行优化来获得。根据 Kuhn-Tucker 条件，在 122 WXWE T 的约束下， XWGE T 的最优值能在满足下式的点上获得。 0 WXWXgE T （ 2.8） 这里， 是一个恒定值， XWXgWE TT 00 , 0W 是优化后的 W 值。下面我们利用牛顿迭代法解方程（ 3.8）。用 F 表示式
10、（ 3.8）左边的函数，可得 F 的雅可比矩阵 WJF 如下： IXWgXXEWJF TT （ 2.9） 为了简化矩阵的求逆，可以近似为（ 3.9）式的第一项。由于数据被球化， IXXE T ,所以， IXWgEXWgEXXEXWgXXE TTTTT 。因而雅可比矩阵变成了对角阵，并且能比较容易地求逆。因而可以得到下面的近似牛顿迭代公式： WWWXWgEWXWXgEWW TT/ （ 2.10） 这里， W 是 W 的新值， XWXgWE TT ，规格化能提高解的稳定性。简化后就可以得到 FastICA算法的迭代公式： WWWWXWgEXWXgEW TT/ (2.11) 实践中， FastICA
11、算法中用的期望必须用它们的估计值代替。当然最好的估计是相应的样本平均。理想情况下，所有的有效 数据都应该参与计算，但这会降低计算速度。所以通常用一部分样本的平均来估计，样本数目的多少对最后估计的精确度有很大影响。迭代中的样本点应该分别选取，假如收敛不理想的话，可以增加样本的数量。 3. FastICA算法的基本步骤： 1. 对观测数据 X 进行中心化，使它的均值为 0； 2. 对数据进行白化， ZX 。 3. 选择需要估计的分量的个数 m ，设迭代次数 1p 4. 选择一个初始权矢量（随机的） pW 。 5. 令 WZWgEZWZgEW TpTpp ,非线性函数 g 的选取见前文。 6. jp
12、j jTppp WWWWW 11。 7. 令ppp WWW /。 8. 假如 pW 不收敛的话，返回第 5 步。 9令 1pp ，如果 mp ，返回第 4 步。 二 MATLAB 源程序及说明： %下程序为 ICA的调用函数，输入为观察的信号，输出为解混后的信号 function Z=ICA(X) %-去均值 - M,T = size(X); %获取输入矩阵的行 /列数，行数为观测数据的数目，列数为采样点数 average= mean(X); %均值 for i=1:M X(i,:)=X(i,:)-average(i)*ones(1,T); end %-白化 /球化 - Cx = cov(X,
13、1); %计算协方差矩阵 Cx eigvector,eigvalue = eig(Cx); %计算 Cx的特征值和特征向量 W=eigvalue(-1/2)*eigvector; %白化矩阵 Z=W*X; %正交矩阵 %-迭代 - Maxcount=10000; %最大迭代次数 Critical=0.00001; %判断是否收敛 m=M; %需要估计的分量的个数 W=rand(m); for n=1:m WP=W(:,n); %初始权矢量（任意） % Y=WP*Z; % G=Y.3;%G为非线性函数，可取 y3等 % GG=3*Y.2; %G的导数 count=0; LastWP=zeros(
14、m,1); W(:,n)=W(:,n)/norm(W(:,n); while abs(WP-LastWP) %迭代次数 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); for j=1:n-1 WPP=WPP+(WP*W(:,j)*W(:,j); end WP=WP-WPP; WP=WP/(norm(WP); if count=
15、Maxcount fprintf(未找到相应的信号 ); return; end end W(:,n)=WP; end Z=W*Z; %以下为主程序，主要为原始信号的产生，观察信号和解混 信号的作图 clear all;clc; N=200;n=1:N;%N为采样点数 s1=2*sin(0.02*pi*n);%正弦信号 t=1:N;s2=2*square(100*t,50);%方波信号 a=linspace(1,-1,25);s3=2*a,a,a,a,a,a,a,a;%锯齿信号 s4=rand(1,N);%随机噪声 S=s1;s2;s3;s4;%信号组成 4*N A=rand(4,4); X=
16、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); %观察信号（混合信号）波形图 figure(2);subplot(4,1,1);plot(X(1,:);title(观察信号（混合信号） ); subplot(4,1,2);plot(X(2,:);
17、 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);plot(Z(4,:);xlabel(Time/ms); 三 实验结果： 实验结果如下 所示：其中图 2为源信号的波形图， 图 3为观察信号（混合信号）波形图，图 4为 解混后的信号波形图。 从图 4可以看出， 执行 I
18、CA后，可以将图 2中的 4种信号 分离出来，且误差较小 （个别信号发生反相） 。 与实验一相比，分离效果要提高很多。 图 2-源信号 图 3-观察信号（混合信号） 图 4-解混后的信号 四源信号的分布特性对分离效果的影响 ： 下列 3个图为源信号中含有不同个数的高斯白噪声分解后得到的相应解混信号 ：其中图5源信号只含一个随机噪声 ，图 6源信号 含两 个随机噪声 ，图 7源信号 含三 个随机噪声 。 从图 5,6,7可以看出， 源信号所含的 高斯白噪声越多， 分离 后得到的信号与源信号相比误差越大， 效果越差；所含高斯白噪声越少，分离效果越好。 图 5-源信号只含一个随机噪声 分离 后得到的
19、波形图 图 6-源信号 含两 个随机噪声 分离 后得到的波形图 图 7-源信号 含三 个随机噪声 分离 后得到的波形图 下图 8,9,10分别为方波信号，正弦信号和锯齿 波 3种 信号与 同一个噪声混合后，经 ICA解混后得到的 结果。为便于比较，源信号中 3种信号所含的能量相等。信号的非高斯性可用四阶累积量（峰度） 进行描述， 0称为超高斯型， 0称为亚高 斯型，可用 的大小作为信号距离高斯型程度的度量。 越大，表明信号距高斯型越远，即信号的非高斯性越强。 下图中 3种信号在源信号中的四阶累积量 分别为：方波为 -8.0，正弦波为 -6.0，锯齿波为 -4.8154，表明方波的非高斯性最强，锯齿波德非高斯性最弱。 经 ICA分离后 得到的 方波、正弦波和锯齿波 与源信号中对应的三种信号的相关系数依次为： -0.9994， -0.9987， -0.9925。可以看出，分解后得到的三种信号中方波与源信号最为接近， 锯齿波与源信号差距最大。 结合四阶累积量，可以得出：在同一个 ICA系统中，信 号的非高斯性越强，分离出来的信号越接近源信号，分离效果越好；反之，分离效果越差。 图 8-方波信号分离 后得到的波形图 图 9-正弦信号分离 后得到的波形图 图 10-锯齿波信号分离 后得到的波形图
