FastICA的原理及实现

为什么是ICA而不是PCA

ICA分离出来的是非高斯分布的信号，而PCA是假设高斯分布，非高斯分布在均值为0方差一样的情况下，信息熵比高斯分布小，所以在0附近有比高斯分布更高的峰值。所以更适合学习稀疏特征。ICA是使分量最大独立化，PCA是使重构误差最小。
PCA与ICA的对比参看https://blog.csdn.net/vendetta_gg/article/details/106521295。

FastICA相比ICA的进步

ICA需要假定si的先验分布函数为sigmoid函数，而FastICA不需要假定先验分布。
ICA使用梯度下降法更新，收敛速度是一次的。FastICA使用牛顿法更新，收敛速度至少是二次的，这就是FastICA名字的来源。
ICA只能一次把所有独立成分全求出来，FastICA可以根据需要只求出1～n个独立成分。

论文的解读

此部分参照独立成分分析FastICA算法原理-知乎，并且加入了我的一些理解。
对于d维随机变量xxx，假设是由相互独立的源sss通过AAA矩阵线性组合产生：
x=Asx=Asx=As
如果sss服从高斯分布，则不能还原出唯一的sss，如果sss非高斯，则可以通过找到WWW使s=Wxs=Wxs=Wx，使sss相互独立从而得到WWW和sss。

为什么ICA可以恢复原始的源？

Darmois - Skitovitch theorem 假设Si是相互独立的源噪声，对于两个随机变量
X=a1S1+…+anSn,Y=b1S1+…+bnSnX=a_1S_1+…+a_nS_n,\\ Y=b_1S_1+…+b_nS_nX=a1S1+…+anSn,Y=b1S1+…+bnSn
若XYXYXY相互独立，则对任意aibi≠0a_ib_i \neq 0aibi=0，必有SiS_iSi为高斯分布。即两个独立的分布，除了高斯分布，两个非零的线性组合一定不相互独立。
根据定理，假设n=2，SSS非高斯，于是
x1=a1s1+a2s2,x2=b1s1+b2s2x_1=a_1s_1+a_2s_2, \\ x_2=b_1s_1+b_2s_2x1=a1s1+a2s2,x2=b1s1+b2s2
且z1=w1TX,z2=w2TXz_1=w_1^T X, z_2=w_2^TXz1=w1TX,z2=w2TX，且相互独立，则
z1=c1s1+c2s2,z2=d1s1+d2s2z_1=c_1s_1+c_2s_2,\\ z_2=d_1s_1+d_2s_2z1=c1s1+c2s2,z2=d1s1+d2s2
根据定理，因为s1,s2s_1,s_2s1,s2非高斯，所以c1d1=0，c2d2=0c_1d_1=0，c_2d_2=0c1d1=0，c2d2=0。但是z1,z2z_1,z_2z1,z2不是零，所以z1,z2z_1,z_2z1,z2分别为c1s1，b2s2c_1s_1，b_2s_2c1s1，b2s2（或者反过来），总之可以分离成两个相互独立的z1,z2z_1,z_2z1,z2。即把x1,x2x_1,x_2x1,x2重新线性组合为s1,s2s_1,s_2s1,s2使其相互独立后，s1,s2s_1,s_2s1,s2就是独立源噪声（可能多一个倍数关系）。

目标函数

我们使用互信息（独立的情况下SSS的信息熵为分量sis_isi信息熵的和，用分量信息熵和减去SSS信息熵，越小说明越接近独立）作为目标函数。
I(s1,…,sn)=∑(H(si))−H(S)I(s_1,…,s_n)=\sum(H(s_i))-H(S)I(s1,…,sn)=∑(H(si))−H(S)
首先计算H(S)H(S)H(S)。因为S=WXS=WXS=WX，
所以logP(S)=logP(X)−log∣detW∣logP(S)=logP(X)-log|detW|logP(S)=logP(X)−log∣detW∣
如果detW=1，log∣detW∣=0detW=1，log|detW|=0detW=1，log∣detW∣=0。此时H(S)H(S)H(S)与WWW无关，只需要最小化每一个H(si)H(s_i)H(si)。

如何保证detW=1detW=1detW=1呢？
首先，根据SSS是相互独立的源，可以假设sisisi的方差是一致的。那么SSS的协方差为III的整数倍。不妨设为III。这时XXX可以看成SSS每个源先经过放大器再线性组合得到，并不影响模型。在XXX接收器，线性组合后，变成X^\hat{X}X^，使其协方差为III。这些改动没有改变XXX是SSS线性组合的本质。但是经过这些假设，从SSS到XXX的中间线性变换WWW的行列式绝对值变为1，因为
SST=WXXTWT,SST=nI,XXT=nI.SS^T=WXX^TW^T,\\SS^T=nI, \\XX^T=nI.SST=WXXTWT,SST=nI,XXT=nI.取行列式发现∣detW∣=1|detW|=1∣detW∣=1。
XXX正交化利用协方差特征分解。
XXT=ASSTAT=AAT=EDET,XX^T=ASS^TA^T=AA^T=EDE^T,XXT=ASSTAT=AAT=EDET,则
ED−1/2ETXXTED−1/2ET=ED−1/2ETEDETED−1/2ET.ED^{-1/2}E^TXX^TED^{-1/2}E^T=ED^{-1/2}E^TEDE^TED^{-1/2}E^T.ED−1/2ETXXTED−1/2ET=ED−1/2ETEDETED−1/2ET.注意EET=IEE^T=IEET=I，所以原式=III。
令X^=ED−1/2ETX，X^X^T=I,\hat{X}=ED^{-1/2}E^TX，\hat{X}\hat{X}^T=I,X^=ED−1/2ETX，X^X^T=I,即X^\hat{X}X^是正交矩阵。

完成以后才发现其实原文定义的III并不是这样。原文令J(y)=H(ygauss )−H(y)J({y})=H\left({y}_{\text {gauss }}\right)-H({y})J(y)=H(ygauss )−H(y)，其中ygaussy_{gauss}ygauss是与yyy同方差的正态分布。I(y1,y2,⋯,yn)=J(y)−∑iJ(yi).I\left(y_{1}, y_{2}, \cdots, y_{n}\right)=J(\mathbf{y})-\sum_{i} J\left(y_{i}\right).I(y1,y2,⋯,yn)=J(y)−i∑J(yi).实际上对SSS在假定独立且协方差为III时，H(Sgauss)=∑iH(sigauss)H(S_{gauss})=\sum_i H(s_{i_{gauss}})H(Sgauss)=∑iH(sigauss)。因此这个互信息和上面的定义等价。因为H(Sgauss)H(S_{gauss})H(Sgauss)也是一个常数，再根据上面已经证明的H(S)H(S)H(S)与WWW无关，J(S)=H(Sgauss)−H(S)J(S)=H(S_{gauss})-H(S)J(S)=H(Sgauss)−H(S)与WWW无关。所以最小化I(s1,…,sn)I(s_1,…,s_n)I(s1,…,sn)，只需要最大化∑iJ(yi)\sum_{i} J\left(y_{i}\right)∑iJ(yi)。根据知乎大佬的解释，最小化H(si)H(s_i)H(si)就是最大化非高斯性，因为高斯分布在同均值方差下有最大熵。J(si)J(s_i)J(si)就是衡量非高斯性大小的。

原文用JJJ定义后用了另一个式子来近似J(yi)≈c[E{G(yi)}−E{G(ν)}]2,J\left(y_{i}\right) \approx c\left[E\left\{G\left(y_{i}\right)\right\}-E\{G(\nu)\}\right]^{2},J(yi)≈c[E{G(yi)}−E{G(ν)}]2,其中ccc是正常数，ν∼N(0,1)\nu\sim N(0,1)ν∼N(0,1)，GGG是非二次函数。
JG(wi)=[E{G(wiTx)}−E{G(ν)}]2J_G(w_i)=[E\left\{G(w_i^Tx)\right\}-E\left\{G(\nu)\right\} ]^2JG(wi)=[E{G(wiTx)}−E{G(ν)}]2
si=wiTxs_i=w_i^Txsi=wiTx.
至此已经找到目标函数了：
最大化 ∑i=1nJG(wi)wrt. wi,i=1,⋯,n\sum_{i=1}^{n} J_{G}\left(\mathbf{w}_{i}\right) \text { wrt. } \mathbf{w}_{i}, i=1, \cdots, n∑i=1nJG(wi) wrt. wi,i=1,⋯,n
其中约束条件为si,sjs_i,s_jsi,sj无关，即
E{(wkTx)(wjTx)}=δjkE\left\{\left(\mathbf{w}_{k}^{T} \mathbf{x}\right)\left(\mathbf{w}_{j}^{T} \mathbf{x}\right)\right\}=\delta_{j k}E{(wkTx)(wjTx)}=δjk，其中δjk\delta_{jk}δjk是示性函数。

GGG的选择

然后考虑非二次函数G(u)G(u)G(u)的选择。论文给了几个定理，分别从www估计的一致性，渐进方差和鲁棒性出发给了几个定理和选择标准。具体可以查看原论文。中间有一个地方比较重要，需要GGG满足：
条件1
E{sig(si)−g′(si)}[E{G(si)}−E{G(ν)}]>0E\left\{s_{i} g\left(s_{i}\right)-g^{\prime}\left(s_{i}\right)\right\}\left[E\left\{G\left(s_{i}\right)\right\}-E\{G(\nu)\}\right]>0E{sig(si)−g′(si)}[E{G(si)}−E{G(ν)}]>0
注意左边第二项是JG(wi)J_{G}(w_i)JG(wi)的绝对值。该论文选择的GGG为logP(s)logP(s)logP(s)，并且在sss是指数族分布下进行一些近似。这个在后面优化会用到。
G1(u)=1a1log⁡cosh⁡(a1u)g1(u)=tanh⁡(a1u)G2(u)=−1a2exp⁡(−a2u2/2)g2(u)=uexp⁡(−a2u2/2)G3(u)=14u4g3(u)=u3\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)=a11logcosh(a1u)=tanh(a1u)=−a21exp(−a2u2/2)=uexp(−a2u2/2)=41u4=u3
其中g(∗)g(*)g(∗)是G(∗)G(*)G(∗)的导数。实际上算法还需要g′(∗)g'(*)g′(∗)但是没有给出。实现的时候还需要再算一下。

牛顿法解优化问题

接下来使用拉格朗日方程去解带约束的优化问题。
当E{G(si)}−E{G(ν)}>0E\left\{G\left(s_{i}\right)\right\}-E\{G(\nu)\}>0E{G(si)}−E{G(ν)}>0时，最大化 ∑i=1nJG(wi)wrt. wi,i=1,⋯,n\sum_{i=1}^{n} J_{G}\left(\mathbf{w}_{i}\right) \text { wrt. } \mathbf{w}_{i}, i=1, \cdots, n∑i=1nJG(wi) wrt. wi,i=1,⋯,n等价于最大化∑i=1nE{G(wiTx)}\sum_{i=1}^nE\left\{G(w_i^Tx)\right\}∑i=1nE{G(wiTx)}。逐个考虑wi\mathbf{w}_iwi，即最大化E{G(wiTx)},∥wi∥2=1E\left\{G(\mathbf{w}_i^T\mathbf{x})\right\}, \left\|\mathbf{w}_ i\right\|^2=1E{G(wiTx)},∥wi∥2=1使用拉格朗日方程：
E{xg(wTx)}−βw=0E\left\{\mathbf{x} g\left(\mathbf{w}^{T} \mathbf{x}\right)\right\}-\beta \mathbf{w}=0E{xg(wTx)}−βw=0
可以得到β=E{w0Txg(w0Tx)}\beta=E\left\{\mathbf{w}_{0}^{T} \mathbf{x} g\left(\mathbf{w}_{0}^{T} \mathbf{x}\right)\right\}β=E{w0Txg(w0Tx)}，其中w0\mathbf{w}_0w0是wi\mathbf{w}_iwi的最优点。
记FFF为左边的式子，则JF(w)J_F(\mathbf{w})JF(w)为雅各比矩阵
JF(w)=E{xxTg′(wTx)}−βI.J F(\mathbf{w})=E\left\{\mathbf{x x}^{T} g^{\prime}\left(\mathbf{w}^{T} \mathbf{x}\right)\right\}-\beta \mathbf{I} .JF(w)=E{xxTg′(wTx)}−βI.
E{xxTg′(wTx)}E\left\{\mathbf{x x}^{T} g^{\prime}\left(\mathbf{w}^{T} \mathbf{x}\right)\right\}E{xxTg′(wTx)}合理近似为E{xxT}E{g′(wTx)}=IE{g′(wTx)}E\left\{\mathbf{x x}^{T} \right\}E\left\{g^{\prime}\left(\mathbf{w}^{T} \mathbf{x}\right) \right\}=IE\left\{g^{\prime}\left(\mathbf{w}^{T} \mathbf{x}\right) \right\}E{xxT}E{g′(wTx)}=IE{g′(wTx)}。
接下来用牛顿法求解拉格朗日方程。
w+=w−[E{xg(wTx)}−βw]/[E{g′(wTx)}−β]w∗=w+/∥w+∥\begin{array}{l} \mathbf{w}^{+}=\mathbf{w}-\left[E\left\{\mathbf{x} g\left(\mathbf{w}^{T} \mathbf{x}\right)\right\}-\beta \mathbf{w}\right] /\left[E\left\{g^{\prime}\left(\mathbf{w}^{T} \mathbf{x}\right)\right\}-\beta\right] \\ \mathbf{w}^{*}=\mathbf{w}^{+} /\left\|\mathbf{w}^{+}\right\| \end{array}w+=w−[E{xg(wTx)}−βw]/[E{g′(wTx)}−β]w∗=w+/∥w+∥
根据条件1左边第一项为正，即分母E{g′(wTx)}−βE\left\{g^{\prime}\left(\mathbf{w}^{T} \mathbf{x}\right)\right\}-\betaE{g′(wTx)}−β为负，所以同乘分母的相反数不改变w+\mathbf{w}^{+}w+的方向，只改变大小，单位化以后不变。即：
w+=E{xg(wTx)}−E{g′(wTx)}ww∗=w+/∥w+∥\begin{aligned} \mathbf{w}^{+} &=E\left\{\mathbf{x} g\left(\mathbf{w}^{T} \mathbf{x}\right)\right\}-E\left\{g^{\prime}\left(\mathbf{w}^{T} \mathbf{x}\right)\right\} \mathbf{w} \\ \mathbf{w}^{*} &=\mathbf{w}^{+} /\left\|\mathbf{w}^{+}\right\| \end{aligned}w+w∗=E{xg(wTx)}−E{g′(wTx)}w=w+/∥∥w+∥∥

有了一个w\mathbf{w}w的估计，就可以进行n次从而得到n个独立成分。为了防止收敛到同一个w\mathbf{w}w，可以使用施密特正交化把wi\mathbf{w}_iwi单位正交化。只需要在第i+1个wi+1\mathbf{w}_{i+1}wi+1去掉其在前i个方向上的分量，然后单位化即可：
wi+1=wi+1−∑j=1iwi+1Twjwjwi+1=wi+1/wp+1Twi+1\mathbf{w}_{i+1}=\mathbf{w}_{i+1}-\sum_{j=1}^{i} \mathbf{w}_{i+1}^{T} \mathbf{w}_{j} \mathbf{w}_{j}\\ \mathbf{w}_{i+1}=\mathbf{w}_{i+1} / \sqrt{\mathbf{w}_{p+1}^{T} \mathbf{w}_{i+1}}wi+1=wi+1−j=1∑iwi+1Twjwjwi+1=wi+1/wp+1Twi+1

还有一种对称去相关的方法，与白化基本一致，不赘述直接贴公式
W=(WWT)−1/2W\mathbf{W}=\left(\mathbf{W W}^{T}\right)^{-1 / 2} \mathbf{W}W=(WWT)−1/2W
对称矩阵的-1/2次幂意味着特征分解后的逆矩阵
(WW)−1/2=ED−1/2ET(\mathbf{W W})^{-1 / 2}=\mathbf{E D}^{-1 / 2} \mathbf{E}^{T}(WW)−1/2=ED−1/2ET
在具体得到n个wi\mathbf{w}_{i}wi的实现，我使用随机选取X的样本点和随机从G1,G2,G3G_1, G_2, G_3G1,G2,G3中选择GGG的方式，防止wi\mathbf{w}_{i}wi相关性太强。

python实现及sklearn包调用

numpy实现

首先使用numpy编写自己的算法，然后与sklearn包比较。

import numpy as np

# 1.白化
def whiten(x):# x是一列一个记录，一行一个特征mean=np.mean(x,axis=1)sd=np.std(x,axis=1)std_x=(x-mean)/sdLambda, Vec=np.linalg.eig(std_x.dot(std_x.T))# C = V*Lambda*V'X_white=(np.dot(Vec.dot(np.diag(1/np.sqrt(Lambda))),Vec.T))X_white=X_white.dot(std_x)# print(1/np.sqrt(Lambda))return X_white# 验证白化效果
X=np.mat(np.random.randn(5,100))
A=np.mat(np.array([[1,0,0,0,0],
[0,2,0,0,1],
[0,0,1,1,0],
[1,1,0,1,0],
[0,1,0,0,1]
]))
X=A.dot(X)
print(A.dot(A.T))X=whiten(X)
print(np.dot(X,X.T))

定义G
def G_1(x,a=1.5):# G_1(u)=1/a1 logcosh(a1u)# G_1'=tanh(a1u)# G_1''=a1(1-tan(a1u)^2)# 返回一阶导数和二阶导数g=np.tanh(a*x)g_dif=a*(1-g**2)return (g,g_dif)def G_2(x,a=1):# G_2(u)=-1/a2*exp(-a2*u^2/2)# G_2'=u*exp(-a2*u^2/2)# G_2''=exp(-a2*u^2/2)-a2*u^2*exp(-a2*u^2/2)g_2=x*np.exp(-(a*x**2/2))g_2_dif=g_2/x - g_2*a*xreturn (g_2,g_2_dif)def G_3(x):# G_3(u)=1/4 * u^4# G_3'=u^3# G_3''=3*u^2g_3=x**3g_3_dif=3*x**2return (g_3,g_3_dif)

牛顿法求单个www

def newton_method(X,G,loops,eps=0.0001):"""参数X: 观测矩阵，每列为一个观测，每行为一个特征G: 函数Gloops: 最大循环次数eps: 判断收敛的参数 当循环后w变化的模小于eps停止循环返回值列向量w"""m,n=X.shapew0=np.random.random([m,1])w0=w0/np.dot(w0.T,w0)w=w0for i in range(loops):E1=np.zeros((m,1))E2=np.zeros((m,1))for k in range(n):x=X[:,k]y=np.dot(w.T,x)g,g_dif=G(y)E1+=g.item()*xE2+=g_dif.item()E1=E1/nE2=E2/nw1=E1-E2*ww1=w1/np.dot(w1.T,w1)dif_w=w1-w0dif_w_size=np.dot(dif_w.T,dif_w)w=w1if dif_w_size<eps:break# print('i',w)return w

def extract_w(X,bs,loops,n=None):"""提取n个w作为矩阵输出参数X: 观测矩阵，每列为一个观测，每行为一个特征bs: 一次提取多少个观测求wloops: 一次牛顿法用至多多少次循环n: 求多少个独立成分。默认X的特征数。返回值w作为列向量构成矩阵W"""if n is None:n=X.shape[0]m=X.shape[1]W=np.mat(np.zeros([X.shape[0],n]))for i in range(n):idx=np.random.choice(range(m),bs,False)X_sample=X[:,idx]if np.random.sample(1)>0.5:G=G_1elif np.random.sample(1)>0.5:G=G_2else:G=G_3w=newton_method(X_sample,G,loops)if i>0:for j in range(i-1):# print(w.T)# print(W[:,j])# print(np.dot(w.T,W[:,j]).item())w=w-np.dot(w.T,W[:,j]).item() * W[:,j]w=w/np.dot(w.T,w)W[:,i]=wreturn W

import matplotlib.pyplot as plt
# 生成波函数1
z0=np.hstack((np.ones(10),0,-np.ones(10)))
z1=np.tile(z0,(10))
z1=z1+np.random.random(len(z1))*0.01
plt.plot(z1)# 生成波函数2
z0=[0.05*i for i in range(21)]
z0=z0-np.mean(z0)
z2=np.tile(z0,(10))
z2=z2+np.random.random(len(z2))*0.01
plt.plot(z2)# 生成波函数3
z0=np.sin(np.linspace(0,2*np.pi,num=21,endpoint=False))
z3=np.tile(z0,(10))*0.5
plt.plot(z3)

# 生成X，A是混合矩阵
A=np.mat([[1,2,0],[0.5,2,0.5],[0,1,1.5]])
z1=z1[:180]
z2=z2[:180]
z3=z3[:180]
Z=np.vstack((z1,z2,z3))
X=A*Z
# 绘制X
for i in range(X.shape[0]):plt.plot(range(X.shape[1]),X[i,:].T)

用FastICA处理信号

X=whiten(X)
W=extract_w_2(X,150,500)
# 重构原始信号S
S=np.dot(W.T,X)
for i in range(3):plt.figure()plt.plot(range(S.shape[1]),S[i,:].T)

sklearn包

用sklearn包直接调用FastICA

from sklearn.decomposition import FastICA
fast_ica=FastICA(n_components=3)
Sr=fast_ica.fit_transform(X)
S=np.dot(Sr.T,X)
for i in range(3):plt.figure()plt.plot(range(S.shape[1]),S[i,:].T)

可以看出我实现的算法效果还是不错的，跟调包差不多。没有进一步增加数据集研究表现。毕竟只是学习笔记。

参照

HYVARINEN A. Fast and robust fixed-point algorithms for independent component analysis[J]. IEEE transactions on Neural Networks, IEEE, 1999, 10(3): 626–634.
小杰-独立成分分析FastICA算法原理
「vendetta_gg」PCA和ICA的对比