独立成分（ICA）分析 原理及公式推导 示例

2024-02-29 19:19:22

独立成分分析（Independent component analysis）  

前言

独立成分分析ICA是一个在多领域被应用的基础算法。ICA是一个不定问题，没有确定解，所以存在各种不同先验假定下的求解算法。相比其他技术，ICA的开源代码不是很多，且存在黑魔法–有些步骤并没有在论文里提到，但没有这些步骤是无法得到正确结果的。

本文给出一个ICA最大似然解法的推导，以及FastICA的python实现，限于时间和实际需求，没有对黑魔法部分完全解读，只保证FastICA实现能得到正确结果。

有兴趣的童鞋可以在未来补上相关内容。

ICA问题表述

设XX是随机向量，且X∈Rn×1X∈Rn×1，这也就说，XX里有nn个成员，每个成员是一个随机变量： 

X=⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜x1x2...xi...xn⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟X=(x1x2...xi...xn)


其中， xi xi是一个随机变量。

 

随机变量有诸多特性，殆由概率论和数理统计教科书详述备尽，在此不一一叙述。

XX里的nn个随机变量是相互非独立的，在一定的假设下，可以用nn个相互独立的随机变量线性组合重新表达XX，也就是说： 

⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜x1x2...xi...xn⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟=A⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜s1s2...si...sn⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟(x1x2...xi...xn)=A(s1s2...si...sn)


其中，sisi是一个随机变量，且两两相互独立，AA是满秩矩阵，且A∈Rn×nA∈Rn×n。 
令: 

S=⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜s1s2...si...sn⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟S=(s1s2...si...sn)


则： 

X=ASX=AS


又有: 

S=A−1XS=A−1X


令： 

W=A−1W=A−1


则: 

S=WXS=WX


其中，W∈Rn×nW∈Rn×n。

 

记录随机向量XX的值mm次，则形成数据集： 

D=⎛⎝⎜⎜⎜⎜d1,1d2,1...dn,1d1,2d2,2...dn,2............d1,md2,m...dn,m⎞⎠⎟⎟⎟⎟D=(d1,1d1,2...d1,md2,1d2,2...d2,m............dn,1dn,2...dn,m)


其中，D∈Rn×mD∈Rn×m

 

ICA的目标，就是在只知道DD的情况下，估算AA，WW，SS的值。

实例：在一个大厅里，有nn个人在随机聊天。在大厅的不同角落，布置nn个麦克风记录大厅的声音，每秒一个记录，一共记录m秒。麦克风记录的混合声音，多个麦克风记录不同位置的混合声音。ICA的目标，就是从混声录音中将每个人的声音分离出来。

理论推导

由前可知: 

 si=(wi,1wi,2...wi,j...wi,n)⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜x1x2...xi...xn⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟ si=(wi,1wi,2...wi,j...wi,n)(x1x2...xi...xn)


令: 

wi=(wi,1wi,2...wi,j...wi,n)wi=(wi,1wi,2...wi,j...wi,n)


则: 

si=wiXsi=wiX

 

设随机变量sisi概率密度函数是psi(si)psi(si)，其中pp的右下角sisi表示随机变量标示，括号中的sisi表示自变量。

由于SS的nn个成员sisi是相互独立的，所以SS的概率密度函数为: 

pS(s)=∏i=1npsi(si)pS(s)=∏i=1npsi(si)

 

设XX的概率密度函数是pX(x)pX(x)，如何根据sisi的概率密度函数求pX(x)pX(x)呢？这是可以做到的。

设随机向量XX的概率分布函数是FX(x)FX(x)，根据概率分布函数和概率密度函数的关系可知： 

pX(x)=F′X(x)pX(x)=FX′(x)

 

同理，设随机向量SS的概率分布函数是FS(s)FS(s)，则： 

pS(s)=F′S(s)pS(s)=FS′(s)

 

根据概率分布函数的定义，有: 

FX(x)=P(X<x)FX(x)=P(X<x)

 

FS(s)=P(S<s)FS(s)=P(S<s)


那么： 

pX(x)=F′X(x)=(P(X<x))′=(P(U<u))′=(P(U<s))′=(∥P(U<Wx)∥)′=(∥P(S<Wx)∥)′=(∥FS(Wx)∥)′=∥F′S(Wx)∥=∥pS(Wx)(Wx)′∥=∥pS(Wx)W∥=pS(Wx)∥W∥=∥W∥∏i=1npsi(wix)pX(x)=FX′(x)=(P(X<x))′=(P(U<u))′=(P(U<s))′=(‖P(U<Wx)‖)′=(‖P(S<Wx)‖)′=(‖FS(Wx)‖)′=‖FS′(Wx)‖=‖pS(Wx)(Wx)′‖=‖pS(Wx)W‖=pS(Wx)‖W‖=‖W‖∏i=1npsi(wix)


其中，上式的第2个等号是概率密度函数的定义，第3个等号是做变量等价代换，以免直接从X变换到S导致思维混乱，第4个等号到第6个等号是逐步将X代换到S，第7个等号是回到SS的概率分布函数定义，第8个等号到第10个等号是求导。

 

从第5个等号开始，对整个等式取行列式运算，因为pX(x)pX(x)一定是标量，对标量做行列式运算是它自身。那么，到了第10个等号，又因为pS(WX)pS(WX)一定是标量，所以可以从行列式运算拿到外面。这里避免的问题的是，如果不对整个等式取行列式，得到的结果是矩阵WW而不是∥W∥‖W‖，这是没有道理的。

注意，在上式中，xx是一个向量，且x∈Rn×1x∈Rn×1，wi∈R1×nwi∈R1×n，psi(si)psi(si)是一个单自变量的函数，pX(x)pX(x)是一个多自变量函数，它的自变量是xx里的多个变量，这样等式左右的每一步就清晰了。

下一步是根据数据集计算WW的值，从概率的角度来说，如果数据集已经记录，那么让这个数据集出现概率最大的WW就是最优值。

前述数据集出现的概率是: 

L=∏i=1m(∥W∥∏j=1npsj(wjdi))L=∏i=1m(‖W‖∏j=1npsj(wjdi))


其中，∏∏表示连乘，didi是DD的第ii列，也就是： 

di=⎛⎝⎜⎜⎜⎜di,1di,2...di,n⎞⎠⎟⎟⎟⎟di=(di,1di,2...di,n)

 

didi的物理意义，也就是第ii次记录随机向量XX得到的nn个值，这nn个值分别对应nn个xixi随机变量。注意，不要把didi和xixi混淆，前者表示DD的一列数据，后者是粗体表示一个随机变量。

上式有最大值，当它取最大值时候的WW就是最优解。如果以梯度下降法求解，需要计算它对WW的偏导，直接求偏导比较复杂，故对它两端取自然对数，则: 

lnL=∑i=1m(ln∥W∥+∑j=1n(lnpsj(wjdi)))=∑i=1m∑j=1nlnpsj(wjdi)+mln∥W∥lnL=∑i=1m(ln‖W‖+∑j=1n(lnpsj(wjdi)))=∑i=1m∑j=1nlnpsj(wjdi)+mln‖W‖


当上式取最大值的时候，LL也同时取最大值，所以求LL的最大值等价于求上式的最大值。

 

用梯度下降法求解上式，需要计算∂lnL∂W∂lnL∂W。这是一个复杂的过程，先从计算∂L∂wu,v∂L∂wu,v开始，它表示WW的第uu行第vv列的一个成员： 

∂lnL∂wu,v=∑i=1m∑j=1n1psj(wjdi)∂psj(wjdi)∂wu,v+m∥W∥∂∥W∥∂wu,v=∑i=1m∑j=1n1psj(wjdi)∂psj(wjdi)∂wu,v+m∥W∥(−1)u+vMuv=∑i=1m1psu(wudi)∂psu(wudi)∂wu,v+m∥W∥(−1)u+vMuv∂lnL∂wu,v=∑i=1m∑j=1n1psj(wjdi)∂psj(wjdi)∂wu,v+m‖W‖∂‖W‖∂wu,v=∑i=1m∑j=1n1psj(wjdi)∂psj(wjdi)∂wu,v+m‖W‖(−1)u+vMuv=∑i=1m1psu(wudi)∂psu(wudi)∂wu,v+m‖W‖(−1)u+vMuv


其中，(−1)u+vMuv(−1)u+vMuv是wu,vwu,v的代数余子式，∂psu(wuxi)∂wu,v∂psu(wuxi)∂wu,v的值要根据psi(si)psi(si)的具体形式求解。

 

对于psi(si)psi(si)，如果在没有任何先验信息的情况下，是无法求解的。如果要求解上式，需要对它做一定的假设，在合理的假设下，可以达到相当不错的近似结果。

设随机变量xixi的概率分布函数是sigmoid函数，因为它是递增，可微，且最大值不超过1，也就是说： 

Fsi(si)=11+e−siFsi(si)=11+e−si


那么，概率密度函数就是: 

psi(si)=F′si(si)=esi(1+esi)2psi(si)=Fsi′(si)=esi(1+esi)2


所以有： 

psu(wudi)=ewudi(1+ewudi)2=ewudi(1+ewudi)−2psu(wudi)=ewudi(1+ewudi)2=ewudi(1+ewudi)−2


故: 

∂psu(wudi)∂wu,v=ewudidi,v(1+ewudi)−2−2ewudi(1+ewudi)−3ewudidi,v=di,vewudi(1+ewudi)2(1−2ewudi1+ewudi)=di,vpsu(wudi)1−ewudi1+ewudi∂psu(wudi)∂wu,v=ewudidi,v(1+ewudi)−2−2ewudi(1+ewudi)−3ewudidi,v=di,vewudi(1+ewudi)2(1−2ewudi1+ewudi)=di,vpsu(wudi)1−ewudi1+ewudi


其中，di,vdi,v是didi的第vv行的一个成员。

 

因此: 

∂lnL∂wu,v=∑i=1m1psu(wudi)∂psu(wudi)∂wu,v+m∥W∥(−1)u+vMuv=∑i=1m1psu(wudi)di,vpsu(wudi)1−ewudi1+ewudi+m∥W∥(−1)u+vMuv=∑i=1mdi,v1−ewudi1+ewudi+m∥W∥(−1)u+vMuv∂lnL∂wu,v=∑i=1m1psu(wudi)∂psu(wudi)∂wu,v+m‖W‖(−1)u+vMuv=∑i=1m1psu(wudi)di,vpsu(wudi)1−ewudi1+ewudi+m‖W‖(−1)u+vMuv=∑i=1mdi,v1−ewudi1+ewudi+m‖W‖(−1)u+vMuv

 

现在对上式进行矩阵化， 
令： 

K=WDK=WD


其中，K∈Rn×mK∈Rn×m，W∈Rn×nW∈Rn×n，D∈Rn×mD∈Rn×m，那么，ku,iku,i就是KK的第uu行的第ii列的一个成员， 
令： 

g(x)=1−ex1+exg(x)=1−ex1+ex


令: 

Z=g(K)=⎛⎝⎜⎜⎜⎜g(k1,1)g(k2,1)...g(kn,1)g(k1,2)g(k2,2)g(kn,2).........g(k1,m)g(k2,m)g(kn,m)⎞⎠⎟⎟⎟⎟Z=g(K)=(g(k1,1)g(k1,2)...g(k1,m)g(k2,1)g(k2,2)...g(k2,m)...g(kn,1)g(kn,2)...g(kn,m))

 

那么，就得到: 

∂lnL∂wu,v=zTudv+m∥W∥(−1)u+vMuv∂lnL∂wu,v=zuTdv+m‖W‖(−1)u+vMuv


其中，zuzu是ZZ的第uu行，dvdv是DD的第vv列。

 

于是，对WW而言，则有： 

∂lnL∂W=ZTD+m∥W∥(W∗)T∂lnL∂W=ZTD+m‖W‖(W∗)T


其中，W∗W∗是WW的伴随矩阵，(W∗)T(W∗)T是W∗W∗的转置，它的第ii行第jj列的元素是wi,jwi,j的代数余子式，也就是(−1)i+jMi,j(−1)i+jMi,j。

 

根据矩阵和它的伴随阵的性质可知: 

WW∗=∥W∥IWW∗=‖W‖I


其中，II是单位矩阵。 
根据上两式可知: 

∂lnL∂W=ZTD+m∥W∥(W∗)T=ZTD+m∥W∥(∥W∥W−1)T=ZTD+m(W−1)T∂lnL∂W=ZTD+m‖W‖(W∗)T=ZTD+m‖W‖(‖W‖W−1)T=ZTD+m(W−1)T

 

那么，在梯度下降法求解WW的时候，更新公式是: 

W=W+α(ZTD+m(W−1)T)W=W+α(ZTD+m(W−1)T)


其中，αα是学习速率。

 

最后的结论简洁且美，Verweile doch, du bist so schön。然并卵，按照这个结果实现代码，计算结果是不合理的，无法恢复原始信号。于是，在实现FastICA之后，可以认为本推导缺少一些黑魔法，至于到底缺少什么并不知道，限于时间关系和实际需求，不再继续研究下去。

FastICA

FastICA计算性能更好。《Indepdent Componet analysis》一书在第8章给出了FastICA的算法流程，如下:

独立成分（ICA）分析 原理及公式推导 示例

白化

FastICA需要对数据做白化处理。设xx是一个随机变量，存在一个线性变换VV将它变换成zz： 

z=Vxz=Vx


且： 

E{zzT}=IE{zzT}=I


那么，VV就是白化变换矩阵。

 

xx的协方差阵是Cx=E{xxT}Cx=E{xxT}，Cx=PDPTCx=PDPT，PP是CxCx的单位特征向量，DD是CxCx的特征值组成的对角阵。那么，VV的值就是： 

V=D−12PTV=D−12PT


证明如下： 
根据相关性质，有PT=P−1PT=P−1，由于DD对角阵，则(D−12)T=D−12(D−12)T=D−12， 
那么： 

E{Vx(Vx)T}=E{VxxTVT}=E{VPDPTVT}=E{VPDPTVT}=E{D−12PTPDPTP(D−12)T}=E{D−12DD−12}=E{I}=IE{Vx(Vx)T}=E{VxxTVT}=E{VPDPTVT}=E{VPDPTVT}=E{D−12PTPDPTP(D−12)T}=E{D−12DD−12}=E{I}=I

 

代码实现

基于python2.7，matplotlib，numpy实现ICA，主要参考sklean的FastICA实现。

#!/usr/bin/env python

#FastICA from ICA book, table 8.4 

import math
import random
import matplotlib.pyplot as plt
from numpy import *

n_components = 2

def f1(x, period = 4):
    return 0.5*(x-math.floor(x/period)*period)

def create_data():
    #data number
    n = 500
    #data time
    T = [0.1*xi for xi in range(0, n)]
    #source
    S = array([[sin(xi)  for xi in T], [f1(xi) for xi in T]], float32)
    #mix matrix
    A = array([[0.8, 0.2], [-0.3, -0.7]], float32)
    return T, S, dot(A, S)

def whiten(X):
    #zero mean
    X_mean = X.mean(axis=-1)
    X -= X_mean[:, newaxis]
    #whiten
    A = dot(X, X.transpose())
    D , E = linalg.eig(A)
    D2 = linalg.inv(array([[D[0], 0.0], [0.0, D[1]]], float32))
    D2[0,0] = sqrt(D2[0,0]); D2[1,1] = sqrt(D2[1,1])
    V = dot(D2, E.transpose())
    return dot(V, X), V

def _logcosh(x, fun_args=None, alpha = 1):
    gx = tanh(alpha * x, x); g_x = gx ** 2; g_x -= 1.; g_x *= -alpha
    return gx, g_x.mean(axis=-1)

def do_decorrelation(W):
    #black magic
    s, u = linalg.eigh(dot(W, W.T))
    return dot(dot(u * (1. / sqrt(s)), u.T), W)

def do_fastica(X):
    n, m = X.shape; p = float(m); g = _logcosh
    #black magic
    X *= sqrt(X.shape[1])
    #create w
    W = ones((n,n), float32)
    for i in range(n): 
        for j in range(i):
            W[i,j] = random.random()

    #compute W
    maxIter = 200
    for ii in range(maxIter):
        gwtx, g_wtx = g(dot(W, X))
        W1 = do_decorrelation(dot(gwtx, X.T) / p - g_wtx[:, newaxis] * W)
        lim = max( abs(abs(diag(dot(W1, W.T))) - 1) )
        W = W1
        if lim < 0.0001:
            break
    return W

def show_data(T, S):
    plt.plot(T, [S[0,i] for i in range(S.shape[1])], marker="*")
    plt.plot(T, [S[1,i] for i in range(S.shape[1])], marker="o")
    plt.show()

def main():
    T, S, D = create_data()
    Dwhiten, K = whiten(D)
    W = do_fastica(Dwhiten)
    #Sr: reconstructed source
    Sr = dot(dot(W, K), D)
    show_data(T, D)
    show_data(T, S)
    show_data(T, Sr)

if __name__ == "__main__":
    main()
          
x
85       1 
#!/usr/bin/env python
2 




3


#FastICA from ICA book, table 8.4 




4






5


import math




6


import random




7


import matplotlib.pyplot as plt




8


from numpy import *




9






10


n_components = 2




11






12


def f1(x, period = 4):




13


    return 0.5*(x-math.floor(x/period)*period)




14






15


def create_data():




16


    #data number




17


    n = 500




18


    #data time




19


    T = [0.1*xi for xi in range(0, n)]




20


    #source




21


    S = array([[sin(xi)  for xi in T], [f1(xi) for xi in T]], float32)




22


    #mix matrix




23


    A = array([[0.8, 0.2], [-0.3, -0.7]], float32)




24


    return T, S, dot(A, S)




25






26


def whiten(X):




27


    #zero mean




28


    X_mean = X.mean(axis=-1)




29


    X -= X_mean[:, newaxis]




30


    #whiten




31


    A = dot(X, X.transpose())




32


    D , E = linalg.eig(A)




33


    D2 = linalg.inv(array([[D[0], 0.0], [0.0, D[1]]], float32))




34


    D2[0,0] = sqrt(D2[0,0]); D2[1,1] = sqrt(D2[1,1])




35


    V = dot(D2, E.transpose())




36


    return dot(V, X), V




37






38


def _logcosh(x, fun_args=None, alpha = 1):




39


    gx = tanh(alpha * x, x); g_x = gx ** 2; g_x -= 1.; g_x *= -alpha




40


    return gx, g_x.mean(axis=-1)




41






42


def do_decorrelation(W):




43


    #black magic




44


    s, u = linalg.eigh(dot(W, W.T))




45


    return dot(dot(u * (1. / sqrt(s)), u.T), W)




46






47


def do_fastica(X):




48


    n, m = X.shape; p = float(m); g = _logcosh




49


    #black magic




50


    X *= sqrt(X.shape[1])




51


    #create w




52


    W = ones((n,n), float32)




53


    for i in range(n): 




54


        for j in range(i):




55


            W[i,j] = random.random()




56






57


    #compute W




58


    maxIter = 200




59


    for ii in range(maxIter):




60


        gwtx, g_wtx = g(dot(W, X))




61


        W1 = do_decorrelation(dot(gwtx, X.T) / p - g_wtx[:, newaxis] * W)




62


        lim = max( abs(abs(diag(dot(W1, W.T))) - 1) )




63


        W = W1




64


        if lim < 0.0001:




65


            break




66


    return W




67






68


def show_data(T, S):




69


    plt.plot(T, [S[0,i] for i in range(S.shape[1])], marker="*")




70


    plt.plot(T, [S[1,i] for i in range(S.shape[1])], marker="o")




71


    plt.show()




72






73


def main():




74


    T, S, D = create_data()




75


    Dwhiten, K = whiten(D)




76


    W = do_fastica(Dwhiten)




77


    #Sr: reconstructed source




78


    Sr = dot(dot(W, K), D)




79


    show_data(T, D)




80


    show_data(T, S)




81


    show_data(T, Sr)




82






83


if __name__ == "__main__":



 
 

84


    main()    







 
 





在这个实现中，创建了两个数据源，一个是正弦函数，一个是线性周期函数，它们的图形如下： 
独立成分（ICA）分析 原理及公式推导 示例 
将这两个数据源混合成两个新数据源，也就是“可观测”的数据，它们的图像如下： 
独立成分（ICA）分析 原理及公式推导 示例


经过FastICA处理后，重建数据源。注意，此时的数据源在图形形状上跟初始数据源具有相似性，但幅度是不一样的，且可能会发生翻转，这是因为ICA是一个不定问题，有多个解符合假设，不是唯一解。 
独立成分（ICA）分析 原理及公式推导 示例

来源： https://blog.csdn.net/lizhe_dashuju/article/details/50263339

<wiz_tmp_tag id="wiz-table-range-border" contenteditable="false" style="display: none;">

 
 
 
 




来自为知笔记(Wiz)





 

上一篇：ANSA二次开发之第三方包导入-通过conda

下一篇：随笔-安装theano的那些故事（亲测有效，附安装包）