ICA使用牛顿迭代法对FastICA算法经行改进

ICA用牛顿迭代法改进的FastICA算法

ICA算法原理：

独立分量分析（ICA）的过程如下图所示：在信源()st中各分量相互独立的假设下，由观察xt通过结婚系统B把他们分离开来，使输出yt逼近st。

图1-ICA的一般过程

ICA算法的研究可分为基于信息论准则的迭代估计方法和基于统计学的代数方法两大类，从原理上来说，它们都是利用了源信号的独立性和非高斯性。基于信息论的方法研究中，各国学者从最大熵、最小互信息、最大似然和负熵最大化等角度提出了一系列估计算法。如FastICA算法, Infomax算法，最大似然估计算法等。基于统计学的方法主要有二阶累积量、四阶累积量等高阶累积量方法。本实验主要讨论FastICA算法。

1. 数据的预处理

一般情况下，所获得的数据都具有相关性，所以通常都要求对数据进行初步的白化或球化处理，因为白化处理可去除各观测信号之间的相关性，从而简化了后续独立分量的提取过程，而且，通常情况下，数据进行白化处理与不对数据进行白化处理相比，算法的收敛性较好。

若一零均值的随机向量

满足

，

其中：I为单位矩阵，我们称这个向量为白化向量。白化的本质在于去相关，这同主分量分析的目标是一样的。在ICA中，对于为零均值的独立源信号

，

有：

，

且协方差矩阵是单位阵cov( S ) = I，因此，源信号

S( t )是白色的。对观测信号X( t )，我们应该寻找一个线性变换，使X( t )投影到新的子空间后变成白化向量，即：

其中，W0为白化矩阵，Z为白化向量。

利用主分量分析，我们通过计算样本向量得到一个变换

其中U和 分别代表协方差矩阵XC的特征向量矩阵和特征值矩阵。可以证明，线性变换W0满足白化变换的要求。通过正交变换，可以保证

因此，协方差矩阵：

再将

代入

且令

有

由于线性变换A~连接的是两个白色随机矢量Z( t )和S( t )，可以得出A~ 一定是一个正交变换。如果把上式中的Z( t )看作新的观测信号，那么可以说，白化使原来的混合矩阵A简化成一个新的正交矩阵A~。证明也是简单的：

其实正交变换相当于对多维矢量所在的坐标系进行一个旋转。

在多维情况下，混合矩阵A是N*N 的，白化后新的混合矩阵A~

由于是正交矩阵，其自由度降为N*(N-1)/2，所以说白化使得ICA问题的工作量几乎减少了一半。

白化这种常规的方法作为ICA的预处理可以有效地降低问题的复杂度，而且算法简单，用传统的PCA就可完成。用PCA对观测信号进行白化的预处理使得原来所求的解混合矩阵退化成一个正交阵，减少了ICA的工作量。此外，PCA本身具有降维功能，当观测信号的个数大于源信号个数时，经过白化可以自动将观测信号数目降到与源信号维数相同。

2. FastICA算法

FastICA算法，又称固定点(Fixed-Point)算法，是由芬兰赫尔辛基大学Hyv?rinen等人提出来的。是一种快速寻优迭代算法，与普通的神经网络算法不同的是这种算法采用了批处理的方式，即在每一步迭代中有大量的样本数据参与运算。但是从分布式并行处理的观点看该算法仍可称之为是一种神经网络算法。FastICA算法有基于峭度、基于似然最大、基于负熵最大等形式，这里，我们介绍基于负熵最大的FastICA算法。它以负熵最大作为一个搜寻方向，可以实现顺序地提取独立源，充分体现了投影追踪（Projection Pursuit）这种传统线性变换的思想。此外，该算法采用了定点迭代的优化算法，使得收敛更加快速、稳健。

因为FastICA算法以负熵最大作为一个搜寻方向，因此先讨论一下负熵判决准则。由信息论理论可知：在所有等方差的随机变量中，高斯变量的熵最大，因而我们可以利用熵来度量非高斯性，常用熵的修正形式，即负熵。根据中心极限定理，若一随机变量X由许多相互独立的随机变量 NiSi,...3,2,1 之和组成，只要iS具有有限的均值和方差，则不论其为何种分布，随机变量X较iS更接近高斯分布。换言之，iS较X的非高斯性更强。因此，在分离过程中，可通过对分离结果的非高斯性度量来表示分离结果间的相互独立性，当非高斯性度量达到最大时，则表明已完成对各独立分量的分离。

负熵的定义：

式中，GaussY是一与Y具有相同方差的高斯随机变量，H为随机变量的微分熵

根据信息理论，在具有相同方差的随机变量中，高斯分布的随机变量具有最大的微分熵。当Y具有高斯分布时，Ng(Y) = 0 ；Y的非高斯性越强，其微分熵越小，Ng(Y)值越大，所以Ng(Y)可以作为随机变量Y非高斯性的测度。由于计算微分熵需要知道Y的概率密度分布函数，这显然不切实际，于是采用如下近似公式：

其中，E为均值运算；g为非线性函数，可取 ，或

等非线性函数，这里，通常

我们取。快速ICA学习规则是找一个方向以便具有最大的非高斯性。这里，

非高斯性用式（3.7）给出的负熵的近似值来度量,的方差约束

为1，对于

白化数据而言，这等于约束W的范数为1。FastICA算法的推导如下。首先，XWT 的负熵的最大近似值能通过对

进行优化来获得。根据Kuhn-Tucker条件，在

的约束下

的最优值能在满足下式的点上获得。

以下是牛顿迭代法在生物医学工程专业中的实际应用，其中选择了3篇有关于牛顿迭代法在ICA算法中优化的应用？。

ICA 去除EEG 中眼动伪差和工频干扰方法研究

眼动伪差和工频干扰是临床脑电图(EEG)中常见噪声 , 严重影响其有用信息提取 . 本文尝试采用独立分量分析 (Independent Component Analysis , ICA) 方法分离 EEG 中此类噪声通过对早老性痴呆症( Alzheimer disease , AD)患者临床 EEG 信号( 含眼动伪差和混入工频干扰 , 信噪比仅 0dB)作 ICA 分析 , 比较了最大熵(Infomax)和扩展最大熵(Extended Infomax) ICA 算法的分离效果 , 证实虽然最大熵算法可以分离出眼动慢波 , 但难以消除工频干扰 , 为此需采用扩展的最大熵算法 ; 并知 ICA 方法在极低信噪比时也有较好的抗干扰性 , 且在处理非平稳信号时有好的鲁棒性 ; 文中还结合近似熵 (approximate entropy , ApEn )分析说明利用 ICA 去除干扰后有助于恢复和保持原始 EEG 信号的非线性特征 .

ICA 算法原理

假设存在 N 道相互独立的未知源信号 X =[ x1 , x2 , … ,xN] , 通过某个未知的待检测分析系统后线性叠加为在 M 个传感器上接收到的观测信号 U =[ u1 , u2 , … , uM] . 已经证明在 M ≥N 条件下 , 如果源信号 X 不含一个以上的高斯过程则存在混合矩阵 A , 使得 U =AX . ICA 算法的基本思路在于求解一矩阵 W , 使其作用于观测信号 U 所得估计信号X ′ =WU , 在统计独立的意义下最逼近于未知源信号 X , 即要找到矩阵 W 使得信号 X ′的各个分量尽可能地相互独立 . 为此需建立一个合适的代价函数 , 再采用某种优化算法分离源信号 . 对于 EEG 信号 , 代价函数可取各导联信号间的高阶统计量、互信息熵或最大似然估计 , 优化算法则可采用牛顿迭代法、基于神经网络的自适应算法、随机梯度法、自然梯度法等本文所采用的最大熵( Infomax) 算法是根据常规的随机梯度法求解矩阵 W , 算法的目标判据为独立分量通过非线性环节后信息熵极大 , 其调节公式为ΔW =μ [ ( WT) -1 + (1 +2y)xT] ( 1)由于算法中涉及矩阵求逆 , 故又提出了用自然梯度代替常规梯度的改进算法 , 以避免求逆过程 , 减少计算量 . 如此 , 矩阵调节改为下式ΔW =μ [ (WT ) -1 + ( 1 +2y)xT] WTW =μ [ I + ( 1 -2y)uT] W其计算量和收敛速度均得到较大改进 .但上述算法仅适用于超高斯信号 , 而眼动肌电和工频干扰之类噪声多属亚高斯信号 . 为了使算法能同时适用于亚高斯信号 , 又推导出所谓扩展的最大熵(Extended

Infomax)算法 , 其调节公式为ΔW =μ [ I -Ktanh (u )uT -uuT ] W ,kii =1 (超高斯信号 ) , kii =-1 (亚高斯信号) 式中K 为对角阵 . 算法通过其主对角元素 Kii取不同值来区分超高斯和亚高斯信号 , 以便将脑电数据中 EEG 有用成分和肌电、工频等干扰信号分解开来 . 而简单的最大熵算法则难以做到这点 .

改进的Fast ICA 算法在事件相关电位提取中的应用研究

事件相关电位的特征提取分析在大脑认知的神经生理基础和临床应用研究中起着非常重要的作用。独立分量分析( Independent Co mpo nent A naly sis, ICA)作为目前比较流行的盲源信号处理算法, 其特性非常适合应用于事件相关电位的提取。文章主要讨论了独立分量分析的基本原理、判决条件和算法, 针对快速定点算法( F ast Independent Component Analy sis alg orithm , Fa st ICA) 迭代次数较多和对初始权值敏感的缺点, 我们引入利用梯度法改进的修正因子, 在此基础上对Fast ICA 进行优化, 得到改进算法, 改进算法降低了对初始权值的敏感性, 减少了迭代次数, 从而提高了算法的收敛速度。最后将其应用于ERP 信号的提取当中, 实验表明, 在分离效果相当的前提下, 收敛速度得到了较大的提高。

ICA 的基本理论

ICA 的线性模型可表示为:

x i( t)=Asi( t) ( 1)

ICA 的目的就是寻求一优化的分离矩阵W , 通过它能由观测信号xi( t)近似地恢复相互独立的源信号si( t):

u i( t)= W x i( t)= W As i( t) ( 2)

基于负熵判踞的Fast ICA 算法

Fast ICA 是一种基于固定点( fixed-point) 迭代理论来寻求非高斯性最大值的方法[ 4] 。由中心极限定理可知,非高斯性可以作为随机信号相互依赖的度量, 所以当非高斯性达到最大时, 表明已完成对各独立分量的分离。由信息论可知负熵可以度量信号的非高斯程度, 因此采用负熵作为独立性判据,可以从观测信号中分离出独立分量。

负熵判据

对于一概率密度函数为p( y) 的随机变量y , 负熵定义为:

J( y)=H( ygauss)-H( y)

H( y)=-∫ p( y) lgp( y) dy

为方便计算,一种较好的负熵近似是:

J( y)∞[ E{ G( y) }-E{ G( y ga us s) } ] 2 ( 5)

Fast ICA 算法实现

预处理

在运用ICA 方法之前, 适当的对原始观测信号进行一些预处理是非常必要的, 它可使ICA 的工作量大大减小,从而有利于提高ICA 算法的效率, 也能使问题更符合前述约束条件。信号预处理包括中心化和白化处理。中心化( 去均值) 是为了使实际的盲源分离问题能够符合ICA 数学模型,而白化则是对去均值后信号向量X 施加一个线性变换,使得新向量X 的各个分量互不相关, 同时X 的协方差矩阵为单位阵。

独立分量提取

对X 进一步处理, 即依据负熵判据来寻找矩阵以实现独立分量的分离。依据牛顿迭代

定理[ 5] , Fast ICA 算法的调整公式为:

w i →E{ zg( w iTz) }-E{ g′ ( w i Tz) }

式中: w i 为某一次牛顿迭代的结果; z 为预处理后的信号; w T i z 为z 在w i 上的投影; w 为调整权值。式只估计出了一个独立分量, 若要估计n 个独立分量,在每次提取一个分量之后, 要从观测信号中减去该独立分量,如此重复直到所有分量都被提取出来为止。去掉已提取的独立分量的方法如下( 假设已

估计了p 个分量):

w p +1 =w p +1 -∑w p T+1w j w j

w p +1 =w p +1/ w T p +1 w p +1

改进的Fast ICA 算法

快速分离算法Fast ICA 的迭代依赖于初始权值。算法的初始权值一般是随机选取的, 因此就会因为权值的值不同而导致每次迭代的效率不同, 得到的独立分量也会稍有不同, 特别是对源信号独立性不太好的情况,对权值的选取比较敏感, 不同的初始权值可能导致收敛性能的不同, 甚至会出现收敛和不收敛两种极端的情况。虽然可以选取白化过程中获得的主成分作为初始权值的值, 但算法很容易收敛到该白化初始值,尽管有良好的收敛性能,但分离效果不佳,与主成分分析效果相当。为了解决该问题, 就应该放宽算法对初始权值的要求, 也即实现大范围收敛。为了改善ICA 算法对初始权值的要求, 决定在Fast ICA 迭代公式中引入松弛因子[ 6] ,迭代公式如下:

w p +1 =w p -α p

[ E{ xg( wp T x) }-βw p]

[ E{ g′ ( w p T x) }-β]

通过引入松弛因子,保证从w k 某个开始进入牛顿迭代法的收敛区域,使得算法在任何情况下均可以达到收敛的效果。当α p =1 时, 该算法就退化到Fast ICA 算法了。带松弛因子的Fast ICA 算法对初始权值不敏感, 但松弛因子的选择却直接影响算法的收敛速度, 由多元微分学可知, 负梯度的方向就是函数值下降最快的方向,因此我们可以选择将梯度值作为松弛因子。计算E{ xg( w P Tx) }在w 处的梯度值:

其中

,调整迭代公

式,得出改进后算法如下:

w′ i →E{ zg( w i Tz) }+α E{ g′ ( w i Tz) } w

w″ i →E{ zg( w i Tz) }+E{ g′ ( w′ i Tz) } w ( 11)

w i →w″ i/ ‖w″ i ‖

改进算法实现的具体步骤为:

( 1) 对观测信号x 进行去均值处理;

( 2) 对去均值后的x 进行白化处理;

( 3)初始化随机权矢量w 0 ,设置收敛误差0 <ε<1 ;

( 4) 按照公式( 10) 计算出松弛因子α;

( 5) 依照公式( 11) 调整w p +1 , 计算出w p +1后,进行去相关和归一化;

( 6) 如果w p +1 -w p <ε ,则算法收敛,估计出独

立分量,

改进的I C A算法及其在 f MR I 信号上的应用

针对目前广泛使用的两种独立成分分析( I C A) 算法(f i x e d -p o i n t 算法和i n f o ma x算法) 在处理功能磁共振成像( f MRI ) 数据时速度较慢的特点, 给出了独立成分分析的一个优化模型,在此基础上,提出了一种快速的牛顿型迭代算法. 该算法采用修正后的牛顿迭代形式,使收敛速度达到三阶. 将文中算法与其它两种算法应用于实际 f MRI 数据,实验结果表明,文中算法能够很好地分离出任务成分,同时大大减少了运算量,提高了运算速度,在处理大数据量的f MRI 信号方面有明显的优势.

I C A的优化模型

1. 1 问题描述

标准独立成分分析的数学模型为

x=A s ( 1)式中:x=( x 1 , x 2 , …, x N)T为N维随机混合信号, 即观测信号; A为N×N常数混合矩阵; s=( s 1 , s 2, …,s N) T为分量彼此独立的源信号. I C A的任务是在只知道观测信号x 的情况下, 寻找N×N阶解混矩阵

W=( w 1 , w 2, …, w N) T,

使得恢复信号y=Wx的各个分量尽可能统计独立,则可以认为

y=( y 1, y 2, …,y N) T就是要找的源信号.

1. 2 优化模型

ICA的实现可以采用不同的估计准则, 包括信息极大化、非高斯性极大化、极大似然估计方法和互信息极小化等,但它们在某些条件下彼此是等价的,可以相互转化. 文中采用非高斯性极大化估计方法,利用负熵来度量随机变量的非高斯性.密度函数为p ( y ) 的随机变量y 的熵定义为

H( y )=-∫ p ( y ) l g p ( y ) d y ( 2)

而负熵J 定义为

J ( y )=H( y g a u s s)-H( y ) ( 3)

其中y g a u s s 表示服从高斯分布的随机变量, 且与随机变量y 具有相同的方差. y 的非高斯性越强, J ( y ) 值越大, 即最大化非高斯性就是最大化负熵J ( y ). 但负熵的计算较为复杂, 难以直接应用, Hy v ?r i n e n等[ 7] 给出了负熵的一个近似表达式( 假设随机变量为零均值且具有单位方差) :

J ( y ) ∝[ E( G (y ) )-E( G( υ )) ] 2 ( 4)

式中: υ 为具有零均值且与y 有相同单位方差的高斯随机变量; G( ·) 为非二次函数, 实际计算中可以选择增长速度不太快的函数,以获得稳定的估计. 由式( 4) 可知, J (y ) 最大即E (G ( y ) ) 最大.在I C A算法中, 需要对观测信号进行预处理,包括中心化和白化处理.

中心化的目的是使观测信号的均值为零,而白化则是将观测信号x 线性变换为具有单位方差且各分量互不相关的z ,即E( z z T)=I . 数据经过白化处理后, 只需要在正交矩阵空间中寻找解混矩阵W, 从而减少了变量维数, 简化了问题的求解. 这样, 得到如下约束优化问题: ma xψ ( w i)=E( G (y ))=E( G( w T iz ) ) ( 5)

s . t . w i 2 =1, i =1, 2, …, N.

I C A的一种牛顿迭代算法

2. 1 算法的导出

从式( 5) 可知,问题( 5) 的约束是等式且较简单,因此可采用在约束集上投影的方法,即先用牛顿法来求解问题( 5) 相应的无约束优化问题, 然后在每一步迭代后,将解向量w i投影到单位球面上以满足约束.由式( 5) 可知, 对ψ(w i)=E( G(w T iz ) ) 求极值, 即

E(G ′ ( w T iz ) )=E(z g( w T iz ))=0, g ( ·)=G′ ( ·) 是G( ·) 的一阶导数. 由牛顿迭代定理可以得到:

这里g ′ ( ·) 是g ( ·) 的一阶导数. 在算法( 6) 中需要计算矩阵的逆,使得实际计算中算法的收敛速度很慢. 而对于白化后的信号z , 可以将其作简化处理,一个合理的估计式如下: E(z z Tg ′ ( w T iz ) ) ≈E( z z T)E ( g ′ (w T iz ) )=E( g ′ ( w T iz ) ) I ( 7)这样矩阵E(z z Tg ′ ( w T iz )) 就变成了对角矩阵,可以很方便地求其逆矩阵. 在牛顿迭代中,由于初始值选择的不同,经常会造成算法不能正确收敛. 为克服此问题,在算法( 6) 中引入步长参数μ 得到: w i→w i-μ E( z g (w T iz ) ) E (g ′ ( w T i z ))

w i→w i w i( 8)

μ 在迭代过程中可以变化, 当μ取足够小( 如为0.10或0. 01) 时, 可大大增强算法的稳定性.

2. 2 算法步骤

牛顿迭代法

x n+1 =x n -f ( x n)/f ′ ( x n) , n=0, 1, 2, … ( 9)

是求解非线性方程f ( x )=0 的重要方法. 算法( 6)就是根据式( 9) 得出的. 可以证明, 当,式( 9) 为二阶收敛. 为加速收敛, 对牛

顿迭代进行了如下修正:

可以证明修正后的牛顿迭代法为三阶收敛[ 8] . 根据

修正的形式,文中对算法( 8) 进行改进得到:

( 11)

算法( 11) 可以提取一个独立成分. 若要提取N个独

立成分,将算法( 11) 写成矩阵形式:

( 12)

其中g ( y )=[ g (y 1), g ( y 2 ) , …, g ( y n) ] T, 每次迭代后, 需使矩阵W成为正交矩阵, 这可以通过对称正交化来实现,即W=(WW T) -1 2 W.

综上所述,文中算法的具体步骤如下: ( 1) 中心化数据x , 使其均值为零; ( 2) 白化数据得到z ; ( 3)选择一个正交矩阵作为初始矩阵; ( 4) 按照式( 12)进行迭代; ( 5) 对称正交化矩阵W=(WW T) -1 2 W;( 6) 如果算法未达到终止条件,返回步骤( 4) .