独立成分分析FastICA算法原理

独立成分分析FastICA算法原理

首先对于d维的随机变量 x ∈ R d × 1 \displaystyle \mathbf{x} \in R^{d\times 1} x∈Rd×1，我们假设他的产生过程是由相互独立的源 s ∈ R d × 1 \displaystyle \mathbf{s} \in R^{d\times 1} s∈Rd×1，通过 A ∈ R d × d \displaystyle A\in R^{d\times d} A∈Rd×d线性组合产生的

x = A s \mathbf{x} =\mathbf{As} x=As

如果s的服从高斯分布的，那么故事结束，我们不能恢复出唯一的s，因为不管哪个方向都是等价的。而如果s是非高斯的，那么我们希望找到w从而 s = w T x \displaystyle s=\mathbf{w}^{T}\mathbf{x} s=wTx，使得 s \displaystyle s s之间的相互独立就可以恢复出s了，我将在后面指出，这等价于最大化每个 s \displaystyle s s的非高斯性。在这之前先说一下为什么ICA可以恢复出原始的源？

这里涉及到一个重要的定理：

Darmois - Skitovitch theorem 假设 S 1 , S 2 , . . . , S n \displaystyle S_{1} ,S_{2} ,...,S_{n} S1​,S2​,...,Sn​是相互独立的源噪声,对于两个随机变量

X = α 1 S 1 + . . . + α n S n Y = β 1 S 1 + . . . + β n S n X=\alpha _{1} S_{1} +...+\alpha _{n} S_{n}\\ Y=\beta _{1} S_{1} +...+\beta _{n} S_{n} X=α1​S1​+...+αn​Sn​Y=β1​S1​+...+βn​Sn​

若X与Y独立，则一定对于任意的 α i β i ≠ 0 \displaystyle \alpha _{i} \beta _{i} \neq 0 αi​βi​​=0,其 S i \displaystyle S_{i} Si​一定是高斯分布。

这个定理揭示了两个事实：

1. 如果我们知道 S i \displaystyle S_{i} Si​不是高斯分布，那么如果X,Y要独立，只能是 α i β i = 0 \displaystyle \alpha _{i} \beta _{i} =0 αi​βi​=0.

2. 另一个稍微反直觉的事情是，即使X,Y都由相同的高斯分布组成，那么X，Y也是有可能独立的。举个简单的例子，假设 S \displaystyle S S都是标准高斯分布，并设， X = S 1 + S 2 , Y = S 1 + a S 2 \displaystyle X=S_{1} +S_{2} ,Y=S_{1} +aS_{2} X=S1​+S2​,Y=S1​+aS2​,他们的协方差为 cov ( X , Y ) = var ( S 1 ) + a var ( S 2 ) \displaystyle \text{cov}( X,Y) =\text{var}( S_{1}) +a\text{var}( S_{2}) cov(X,Y)=var(S1​)+avar(S2​)，显然当 a = − 1 \displaystyle a=-1 a=−1时，他们的协方差等于0，又因为X,Y是高斯分布，因此X,Y是独立的。

根据这个定理，可以推出ICA的有效性。假设n=2，S是非高斯的，于是，在ICA中

X 1 = α 1 S 1 + α 2 S 2 X 2 = β 1 S 1 + β 2 S 2 X_{1} =\alpha _{1} S_{1} +\alpha _{2} S_{2}\\ X_{2} =\beta _{1} S_{1} +\beta _{2} S_{2} X1​=α1​S1​+α2​S2​X2​=β1​S1​+β2​S2​

接下来将演示，信号源是怎么被“分离”出来的。我们希望找到一个旋转，使得X经过变换后独立，即 Z 1 = w 1 T X 1 ,   Z 2 = w 2 T X 2 \displaystyle Z_{1} =\mathbf{w}^{T}_{1} X_{1} ,\ Z_{2} =\mathbf{w}^{T}_{2} X_{2} Z1​=w1T​X1​, Z2​=w2T​X2​，且 Z 1 ⊥ Z 2 \displaystyle Z_{1} \bot Z_{2} Z1​⊥Z2​，于是变换后的结果为

Z 1 = α 1 ′ S 1 + α 2 ′ S 2 Z 2 = β 1 ′ S 1 + β 2 ′ S 2 Z_{1} =\alpha '_{1} S_{1} +\alpha '_{2} S_{2}\\ Z_{2} =\beta '_{1} S_{1} +\beta '_{2} S_{2} Z1​=α1′​S1​+α2′​S2​Z2​=β1′​S1​+β2′​S2​

根据定理，一定有 α 1 ′ β 1 ′ = α 2 ′ β 2 ′ = 0 \displaystyle \alpha '_{1} \beta '_{1} =\alpha '_{2} \beta '_{2} =0 α1′​β1′​=α2′​β2′​=0，由于S的顺序无法确定，不妨假设 α 1 ′ = 0 \displaystyle \alpha '_{1} =0 α1′​=0，于是 β 1 ′ ≠ 0 \displaystyle \beta '_{1} \neq 0 β1′​​=0，而且 α 2 ′ ≠ 0 \displaystyle \alpha '_{2} \neq 0 α2′​​=0（否则Z1就变成0了）,于是 β 2 ′ = 0 \displaystyle \beta '_{2} =0 β2′​=0，事实上只有这两种可能性：

Z 1 = 0 S 1 + α 2 ′ S 2 Z 2 = β 1 ′ S 1 + 0 S 2 Z_{1} =0S_{1} +\alpha '_{2} S_{2}\\ Z_{2} =\beta '_{1} S_{1} +0S_{2} Z1​=0S1​+α2′​S2​Z2​=β1′​S1​+0S2​

和

Z 1 = α 1 ′ S 1 + 0 S 2 Z 2 = 0 S 1 + β 2 ′ S 2 Z_{1} =\alpha '_{1} S_{1} +0S_{2}\\ Z_{2} =0S_{1} +\beta '_{2} S_{2} Z1​=α1′​S1​+0S2​Z2​=0S1​+β2′​S2​

从而 S 1 \displaystyle S_{1} S1​和 S 2 \displaystyle S_{2} S2​就被奇迹般分离了！唯一缺点就是不知道分离之后的顺序而已！

事实上，如果出现source大于X的维度的情况，我们也能分离，这也称为overcomplete-ICA，只要对于X=AS，A满足列满秩就可以。比如有4个S，2个X，于是我们可以通过对X旋转构造出4个 Z i = α i 1 ′ S 1 + α i 2 ′ S 2 + α i 3 ′ S 3 + α i 4 ′ S 4 Z_i=\alpha'_{i1}S_1+\alpha'_{i2}S_2+\alpha'_{i3}S_3+\alpha'_{i4}S_4 Zi​=αi1′​S1​+αi2′​S2​+αi3′​S3​+αi4′​S4​，于是为了使得这4个Z两两独立，根据DS定理，每个Z都只会出现一个S，从而分离出来。

然而我们实际做的时候其实是去最大化高斯性的，他跟独立性有什么关系？如果要s相互独立，那么首先考虑的是他们的互信息：

I ( s 1 , . . . , s n ) = ∑ i = 1 n H ( s i ) − H ( S ) I( s_{1} ,...,s_{n}) =\sum ^{n}_{i=1} H( s_{i}) -H( S) I(s1​,...,sn​)=i=1∑n​H(si​)−H(S)

因为 s i = w T x \displaystyle s_{i} =\mathbf{w}^{T}\mathbf{x} si​=wTx，因此 log ⁡ P S = log ⁡ P X − log ⁡ ∣ det ⁡ W ∣ \displaystyle \log P_{S} =\log P_{X} -\log |\det W| logPS​=logPX​−log∣detW∣。现在，如果W是一个旋转矩阵（其特征值全是1），那么 det ⁡ W = 1 \displaystyle \det W=1 detW=1，于是 log ⁡ ∣ det ⁡ W ∣ = 0 \displaystyle \log |\det W|=0 log∣detW∣=0，因此， H ( S ) H(S) H(S)与W无关，可以看做是常数，因此，最小化互信息将等价于最小化每个方向下的熵 H ( s i ) \displaystyle H( s_{i}) H(si​)，而这其实本质上就是非高斯性，因为高斯的熵是最大的。

如果接来下有两个问题，1. 如何保证W是旋转矩阵。W是旋转矩阵这个性质很重要，他决定了非高斯的优化方法的正确性。但是如何保证这个W是旋转矩阵？其实我们可以通过的X进行白化，使得这个W变成一个旋转矩阵。2.如何具体算非高斯性？接下来先说下白化操作

白化 (Whitening)

求解ICA的第一步是对x白化得到互不相关的列 x ~ = E D − 1 / 2 E T x \displaystyle \tilde{\mathbf{x}} =\mathbf{ED}^{-1/2}\mathbf{E}^{T}\mathbf{x} x~=ED−1/2ETx，这一步一般是用PCA来完成的

x ~ = E D − 1 / 2 E T A s = A ~ s \tilde{\mathbf{x}} =\mathbf{ED}^{-1/2}\mathbf{E}^{T}\mathbf{As} =\tilde{\mathbf{A}}\mathbf{s} x~=ED−1/2ETAs=A~s

其中EDE其实的协方差x的特征值分解，即 x x T = A A T = E D E T \displaystyle \mathbf{xx}^{T} =\mathbf{AA}^{T} =\mathbf{EDE}^{T} xxT=AAT=EDET。白化的一个重要作用就是这个 A ~ \displaystyle \tilde{\mathbf{A}} A~是一个正交矩阵。利用这个正交性，可以证明 E ( x ~ x ~ T ) = I \displaystyle E\left(\mathbf{\tilde{\mathbf{x}}}\tilde{\mathbf{x}}^{T}\right) =I E(x~x~T)=I，也就是每个维度都是线性无关的：

x ~ x ~ T = E D − 1 / 2 E T A s s T A T E D − 1 / 2 E T = E D − 1 / 2 E T E D E T E D − 1 / 2 E T = I \tilde{\mathbf{x}}\tilde{\mathbf{x}}^{T} =\mathbf{ED}^{-1/2}\mathbf{E}^{T}\mathbf{Ass}^{T}\mathbf{A}^{T}\mathbf{ED^{-1/2} E}^{T} =\mathbf{ED}^{-1/2}\mathbf{E}^{T}\mathbf{EDE^{T} ED^{-1/2} E}^{T} =I x~x~T=ED−1/2ETAssTATED−1/2ET=ED−1/2ETEDETED−1/2ET=I

这个白化使得我们需要学习的参数减少了一半，为什么呢，因为正交矩阵的自由度是 n ( n − 1 ) 2 \displaystyle \frac{n( n-1)}{2} 2n(n−1)​，而原始的矩阵的自由度是 n 2 \displaystyle n^{2} n2.现假设后面的x都是经过白化处理的数据。

FastICA的迭代算法

经过白化后，那么ICA要做的就是就是要找到一个最优方向 w \displaystyle \mathbf{w} w，使得改方向的非高斯性最大

max ⁡ J   ( w T x ) \max J\ \left(\mathbf{w}^{T}\mathbf{x}\right) maxJ (wTx)

因此，这只是一个优化问题。那么怎么衡量非高斯性呢，我们可以用负熵来衡量：

J ( y ) = H ( y g a u s s ) − H ( y ) J(\mathbf{y} )=H(\mathbf{y}_{gauss}) -H(\mathbf{y} ) J(y)=H(ygauss​)−H(y)

其中 y g a u s s \displaystyle \mathbf{y}_{gauss} ygauss​是高斯分布，其协方差矩阵跟 y \displaystyle y y是一样的，显然当y是高斯分布的时候等于0.所以当J越大非高斯性就越大。然而这个东西不好算，所以有很多不同的近似方法,一个经典的方法是

J ( y ) ≈ 1 12 E { y 3 } 2 + 1 48 kurt ⁡ ( y ) 2 J(y)\approx \frac{1}{12} E\left\{y^{3}\right\}^{2} +\frac{1}{48}\operatorname{kurt} (y)^{2} J(y)≈121​E{y3}2+481​kurt(y)2

而在fastICA中，使用的是：

J G ( w ) = [ E { G ( w T x ) } − E { G ( ν ) } ] 2 J_{G} (\mathbf{w} )=\left[ E\left\{G\left(\mathbf{w}^{T}\mathbf{x}\right)\right\} -E\{G(\nu )\}\right]^{2} JG​(w)=[E{G(wTx)}−E{G(ν)}]2

其中 E ( ( w T x ) 2 ) = E ( w T x x T w ) = E ( w T w ) = ∥ w ∥ 2 = 1 \displaystyle E\left(\left(\mathbf{w}^{T}\mathbf{x}\right)^{2}\right) =E\left(\mathbf{w}^{T}\mathbf{xx}^{T}\mathbf{w}\right) =E\left(\mathbf{w}^{T}\mathbf{w}\right) =\| \mathbf{w} \| ^{2} =1 E((wTx)2)=E(wTxxTw)=E(wTw)=∥w∥2=1,这是为了限制目标函数不要到无穷大。那么这个w怎么求呢？我们看看他在约束下的导数：

∂ J G ( w ) − β ( ∥ w ∥ 2 − 1 ) ∂ w = 2 E { x G ′ ( w T x ) } − 2 β w = 0 \frac{\partial J_{G} (\mathbf{w} )-\beta \left( \| \mathbf{w} \| ^{2} -1\right)}{\partial \mathbf{w}} =2E\left\{\mathbf{x} G'\left(\mathbf{w}^{T}\mathbf{x}\right)\right\} -2\beta \mathbf{w} =0 ∂w∂JG​(w)−β(∥w∥2−1)​=2E{xG′(wTx)}−2βw=0

当这个方程等于0的时候就得到最优的w了。此外 β \displaystyle \beta β的值可以用最优的w反推出来，假设 w 0 \displaystyle w_{0} w0​是最优的，一定有 E { x G ′ ( w 0 T x ) } − β w o = 0 ⟹ β = E { w 0 T x G ′ ( w 0 T x ) } \displaystyle E\left\{\mathbf{x} G'\left(\mathbf{w}^{T}_{0}\mathbf{x}\right)\right\} -\beta \mathbf{w}_{o} =0\Longrightarrow \beta =E\left\{\mathbf{w^{T}_{0} x} G'\left(\mathbf{w}^{T}_{0}\mathbf{x}\right)\right\} E{xG′(w0T​x)}−βwo​=0⟹β=E{w0T​xG′(w0T​x)}。于是这其实是一个方程求根的问题，我们可以用牛顿法解决：

x n + 1 = x n − f ( x n ) f ′ ( x n ) x_{n+1} =x_{n} -\frac{f( x_{n})}{f^{\prime }( x_{n})} xn+1​=xn​−f′(xn​)f(xn​)​

简单地介绍一下，我们知道方程的根就是他跟x轴的交点，那如果我们用切线来近似曲线，那么可以“认为“切线与x轴的交点就是根，如果不是的话，我们可以一直重复这个流程直到收敛。

∂ E { x G ′ ( w T x ) } − β w ∂ w = E { x x T G ′ ′ ( w T x ) } − β I \frac{\partial E\left\{\mathbf{x} G'\left(\mathbf{w}^{T}\mathbf{x}\right)\right\} -\beta \mathbf{w}}{\partial w} =E\left\{\mathbf{xx}^{T} G''\left(\mathbf{w}^{T}\mathbf{x}\right)\right\} -\beta \mathbf{I} ∂w∂E{xG′(wTx)}−βw​=E{xxTG′′(wTx)}−βI

其中用一个近似的方法 E { x x T G ′ ′ ( w T x ) } ≈ E { x x T } E { G ′ ′ ( w T x ) } = E { G ′ ′ ( w T x ) } I \displaystyle E\left\{\mathbf{xx}^{T} G''\left(\mathbf{w}^{T}\mathbf{x}\right)\right\} \approx E\left\{\mathbf{xx}^{T}\right\} E\left\{G''\left(\mathbf{w}^{T}\mathbf{x}\right)\right\} =E\left\{G''\left(\mathbf{w}^{T}\mathbf{x}\right)\right\}\mathbf{I} E{xxTG′′(wTx)}≈E{xxT}E{G′′(wTx)}=E{G′′(wTx)}I，最终用牛顿法的公式

w n + 1 = w n − E { x G ′ ( w n T x ) } − β w n E { G ′ ′ ( w n T x ) } − β I \mathbf{w}_{n+1} =\mathbf{w}_{n} -\frac{E\left\{\mathbf{x} G'\left(\mathbf{w}_{n}^{T}\mathbf{x}\right)\right\} -\beta \mathbf{w}_{n}}{E\left\{G''\left(\mathbf{w}^{T}_{n}\mathbf{x}\right)\right\} -\beta \mathbf{I}} wn+1​=wn​−E{G′′(wnT​x)}−βIE{xG′(wnT​x)}−βwn​​

我们用 w n \displaystyle \mathbf{w}_{n} wn​而不是 w 0 \displaystyle \mathbf{w}_{0} w0​来近似 β \displaystyle \beta β就得到了递推公式：

w n + 1 = E { G ′ ′ ( w n T x ) } w n − E { x G ′ ( w n T x ) } E { G ′ ′ ( w n T x ) } − β \mathbf{w}_{n+1} =\frac{\mathbf{E\left\{G''\left( w^{T}_{n} x\right)\right\} w_{n}} -E\left\{\mathbf{x} G'\left(\mathbf{w}_{n}^{T}\mathbf{x}\right)\right\}}{E\left\{G''\left(\mathbf{w}^{T}_{n}\mathbf{x}\right)\right\} -\beta } wn+1​=E{G′′(wnT​x)}−βE{G′′(wnT​x)}wn​−E{xG′(wnT​x)}​

最后两边同时乘以 β − E { G ′ ′ ( w n T x ) \displaystyle \beta -E\{G''\left(\mathbf{w}^{T}_{n}\mathbf{x}\right) β−E{G′′(wnT​x)可以进一步简化？论文里面这一步没看懂，希望看懂的告诉我。。。于是迭代公式就变成

w n + 1 = E { x G ′ ( w n T x ) } − E { G ′ ′ ( w n T x ) } w n w n + 1 = w n / ∥ w n ∥ \mathbf{w}_{n+1} =E\left\{\mathbf{x} G'\left(\mathbf{w}_{n}^{T}\mathbf{x}\right)\right\} -\mathbf{E\left\{G''\left( w^{T}_{n} x\right)\right\} w_{n}}\\ \mathbf{w}_{n+1} =\mathbf{w}_{n} /\| \mathbf{w}_{n} \| wn+1​=E{xG′(wnT​x)}−E{G′′(wnT​x)}wn​wn+1​=wn​/∥wn​∥

这个就是fastICA所使用的迭代公式了。

最后，这个G到底是什么，常用的G:

G 1 ( u ) = 1 a 1 log ⁡ cosh ⁡ ( a 1 u ) G 1 ′ ( u ) = tanh ⁡ ( a 1 u ) G 2 ( u ) = − 1 a 2 exp ⁡ ( − a 2 u 2 / 2 ) G 2 ′ ( u ) = u exp ⁡ ( − a 2 u 2 / 2 ) G 3 ( u ) = 1 4 u 4 G 3 ′ ( u ) = u 3 \begin{aligned} G _ { 1 } ( u ) & = \frac { 1 } { a _ { 1 } } \log \cosh \left( a _ { 1 } u \right) \\ G' _ { 1 } ( u ) & = \tanh \left( a _ { 1 } u \right) \\ G _ { 2 } ( u ) & = - \frac { 1 } { a _ { 2 } } \exp \left( - a _ { 2 } u ^ { 2 } / 2 \right) \\ G' _ { 2 } ( u ) & = u \exp \left( - a _ { 2 } u ^ { 2 } / 2 \right) \\ G _ { 3 } ( u ) & = \frac { 1 } { 4 } u ^ { 4 } \\ G' _ { 3 } ( u ) & = u ^ { 3 } \end{aligned} G1​(u)G1′​(u)G2​(u)G2′​(u)G3​(u)G3′​(u)​=a1​1​logcosh(a1​u)=tanh(a1​u)=−a2​1​exp(−a2​u2/2)=uexp(−a2​u2/2)=41​u4=u3​

其中 1 ≥ a 1 ≤ 2 , a 2 ≈ 1 1\ge a_1 \le2,a_2\approx1 1≥a1​≤2,a2​≈1，一般来说

• G 1 G_1 G1​是一般选择
• 如果独立成分s是非常超高斯的，那么用 G 2 G_2 G2​会好一点
• 如果独立成分是sub-gaussian ，且没有异常值，则用 G 3 G_3 G3​
• 如果要减少计算量，那么或许可以用 G 1 G_1 G1​和 G 2 G_2 G2​的线性组合

参考文献

[1] HYVÄRINEN A, OJA E. Independent component analysis: algorithms and applications[J]. Neural networks, Elsevier, 2000, 13(4–5): 411–430.

[2] HYVARINEN A. Fast and robust fixed-point algorithms for independent component analysis[J]. IEEE transactions on Neural Networks, IEEE, 1999, 10(3): 626–634.