《ica快速算法原理和matlab算法程序》由会员分享，可在线阅读，更多相关《ica快速算法原理和matlab算法程序（10页珍藏版）》请在金锄头文库上搜索。
1、 实验实验 2：FastICA 算法算法 一算法原理一算法原理： 独立分量分析（ICA）的过程如下图所示：在信源( )s t中各分量相互独立的假设下，由 观察( )x t通过结婚系统B把他们分离开来，使输出( )y t逼近( )s t。 图 1-ICA 的一般过程 ICA 算法的研究可分为基于信息论准则的迭代估计方法和基于统计学的代数方法两大 类， 从原理上来说， 它们都是利用了源信号的独立性和非高斯性。 基于信息论的方法研究中， 各国学者从最大熵、最小互信息、最大似然和负熵最大化等角度提出了一系列估计算法。如 FastICA 算法, Infomax 算法，最大似然估计算法等。基于统计学的方法
2、主要有二阶累积量、 四阶累积量等高阶累积量方法。本实验主要讨论 FastICA 算法。 1. 数据的预处理数据的预处理 一般情况下， 所获得的数据都具有相关性， 所以通常都要求对数据进行初步的白化或球 化处理， 因为白化处理可去除各观测信号之间的相关性， 从而简化了后续独立分量的提取过 程，而且，通常情况下，数据进行白化处理与不对数据进行白化处理相比，算法的收敛性较 好。 若一零均值的随机向量T M ZZZ, 1 满足IZZE T ，其中：I为单位矩阵，我 们称这个向量为白化向量。白化的本质在于去相关，这同主分量分析的目标是一样的。在 ICA中 ， 对 于 为 零 均 值 的 独 立 源 信
3、号 T N tStStS,., 1 ， 有 ： jiSESESSE jiji 当, 0， 且协方差矩阵是单位阵 IS cov， 因此， 源信号 tS 是白色的。对观测信号 tX，我们应该寻找一个线性变换，使 tX投影到新的子空间后变 成白化向量，即： tXWtZ 0 （2.1） 其中， 0 W为白化矩阵，Z为白化向量。 利用主分量分析，我们通过计算样本向量得到一个变换 T UW 2/1 0 其中U和分别代表协方差矩阵 X C的特征向量矩阵和特征值矩阵。可以证明，线性变换 0 W满足白化变换的要求。通过正交变换，可以保证IUUUU TT 。因此，协方差矩阵： IUXXEUUXXUEZZE TTT
4、TT 2/12/12/12/12/12/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，所以说白化使得
5、ICA 问题的工作量几乎减少了一半。 白化这种常规的方法作为 ICA 的预处理可以有效地降低问题的复杂度，而且算法简单， 用传统的 PCA 就可完成。 用 PCA 对观测信号进行白化的预处理使得原来所求的解混合矩阵 退化成一个正交阵，减少了 ICA 的工作量。此外，PCA 本身具有降维功能，当观测信号的 个数大于源信号个数时，经过白化可以自动将观测信号数目降到与源信号维数相同。 2. FastICA 算法算法 FastICA 算法，又称固定点(Fixed-Point)算法，是由芬兰赫尔辛基大学 Hyv rinen 等人 提出来的。 是一种快速寻优迭代算法， 与普通的神经网络算法不同的是这种算法
6、采用了批处 理的方式， 即在每一步迭代中有大量的样本数据参与运算。 但是从分布式并行处理的观点看 该算法仍可称之为是一种神经网络算法。FastICA 算法有基于峭度、基于似然最大、基于负 熵最大等形式，这里，我们介绍基于负熵最大的 FastICA 算法。它以负熵最大作为一个搜寻 方向，可以实现顺序地提取独立源，充分体现了投影追踪（Projection Pursuit）这种传统线 性变换的思想。此外，该算法采用了定点迭代的优化算法，使得收敛更加快速、稳健。 因为 FastICA 算法以负熵最大作为一个搜寻方向，因此先讨论一下负熵判决准则。由 信息论理论可知：在所有等方差的随机变量中，高斯变量的熵
7、最大，因而我们可以利用熵 来度量非高斯性，常用熵的修正形式，即负熵。根据中心极限定理，若一随机变量X由许 多相互独立的随机变量NiSi,.3 , 2 , 1之和组成， 只要 i S具有有限的均值和方差， 则不论 其为何种分布，随机变量X较 i S更接近高斯分布。换言之， i S较X的非高斯性更强。因 此，在分离过程中，可通过对分离结果的非高斯性度量来表示分离结果间的相互独立性， 当非高斯性度量达到最大时，则表明已完成对各独立分量的分离。 负熵的定义： YHYHYN G a u s sg （2.5） 式中， Gauss Y是一与Y具有相同方差的高斯随机变量，H为随机变量的微分熵 dppYH YY
8、 lg （2.6） 根据信息理论，在具有相同方差的随机变量中，高斯分布的随机变量具有最大的微分 熵。当Y具有高斯分布时， 0YNg；Y的非高斯性越强，其微分熵越小， YNg值越 大，所以 YNg可以作为随机变量Y非高斯性的测度。由于根据式（3.6）计算微分熵需要 知道Y的概率密度分布函数，这显然不切实际，于是采用如下近似公式： 2 G a u s sg YgEYgEYN （2.7） 其 中 ，E为 均 值 运 算 ；g为 非 线 性 函 数 ， 可 取 )tanh( 11 yayg， 或 2/e x p 2 2 yyyg或 3 3 yyg等非线性函数， 这里，21 1 a， 通常我们取1 1
9、a。 快速 ICA 学习规则是找一个方向以便XWYXW TT 具有最大的非高斯性。这里， 非高斯性用式（3.7）给出的负熵)(XWN T g 的近似值来度量,XW T 的方差约束为 1，对于 白化数据而言，这等于约束W的范数为 1。FastICA 算法的推导如下。首先，XW T 的负熵 的最大近似值能通过对XWGE T 进行优化来获得。根据 Kuhn-Tucker 条件，在 1 22 WXWE T 的约束下，XWGE T 的最优值能在满足下式的点上获得。 0 WXWXgE T （2.8） 这里，是一个恒定值， XWXgWE TT 00 , 0 W是优化后的W值。下面我们利用牛 顿迭代法解方程（
10、3.8） 。用F表示式（3.8）左边的函数，可得F的雅可比矩阵 WJF如 下： IXWgXXEWJF TT （2.9） 为了简化矩阵的求逆，可以近似为（3.9）式的第一项。由于数据被球化，IXXE T ,所 以，IXWgEXWgEXXEXWgXXE TTTTT 。 因而雅可比矩阵变成了 对角阵，并且能比较容易地求逆。因而可以得到下面的近似牛顿迭代公式： WWW XWgEWXWXgEWW TT / / （2.10） 这里， W是W的新值，XWXgWE TT ，规格化能提高解的稳定性。简化后就可 以得到 FastICA 算法的迭代公式： WWW WXWgEXWXgEW TT / (2.11) 实践
11、中，FastICA 算法中用的期望必须用它们的估计值代替。当然最好的估计是相应的 样本平均。理想情况下，所有的有效数据都应该参与计算，但这会降低计算速度。所以通 常用一部分样本的平均来估计，样本数目的多少对最后估计的精确度有很大影响。迭代中 的样本点应该分别选取，假如收敛不理想的话，可以增加样本的数量。 3. FastICA 算法的基本步骤：算法的基本步骤： 1. 对观测数据X进行中心化，使它的均值为 0； 2. 对数据进行白化，ZX 。 3. 选择需要估计的分量的个数m，设迭代次数1p 4. 选择一个初始权矢量（随机的） p W。 5. 令WZWgEZWZgEW T p T pp ,非线性函
12、数g的选取见前文。 6. j p j j T ppp WWWWW 1 1 。 7. 令 ppp WWW/。 8. 假如 p W不收敛的话，返回第 5 步。 9令1 pp，如果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); en
13、d %-白化/球化- Cx = cov(X,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的导数 coun
14、t=0; LastWP=zeros(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=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*
