其他分享
首页 > 其他分享> > 独立成分(ICA)分析 原理及公式推导 示例

独立成分(ICA)分析 原理及公式推导 示例

作者:互联网

独立成分分析(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的算法流程,如下:

白化

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()    
   

在这个实现中,创建了两个数据源,一个是正弦函数,一个是线性周期函数,它们的图形如下: 
 
将这两个数据源混合成两个新数据源,也就是“可观测”的数据,它们的图像如下: 

经过FastICA处理后,重建数据源。注意,此时的数据源在图形形状上跟初始数据源具有相似性,但幅度是不一样的,且可能会发生翻转,这是因为ICA是一个不定问题,有多个解符合假设,不是唯一解。 

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

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

       



来自为知笔记(Wiz)



标签:...,xi,推导,示例,ICA,wu,si,ewudi,dot
来源: https://www.cnblogs.com/jins-note/p/10835491.html