第1页 / 共12页
第2页 / 共12页
第3页 / 共12页
第4页 / 共12页
第5页 / 共12页
第6页 / 共12页
第7页 / 共12页
第8页 / 共12页
史上最直白的ICA教程.pdf
ICA: Independent Component Analysis
December 11, 2015
1 前前前言言言
独立成分分析ICA是一个在多领域被应用的基础算法。ICA是一个不定问
题，没有确定解，所以存在各种不同先验假定下的求解算法。相比其他技
术，ICA的开源代码不是很多，且存在黑魔法–有些步骤并没有在论文里提
到，但没有这些步骤是无法得到正确结果的。
本文给出一个ICA最大似然解法的推导，以及FastICA的python实现，
限于时间和实际需求，没有对黑魔法部分完全解读，只保证FastICA实现能
得到正确结果。
有兴趣的童鞋可以在未来补上相关内容。
ICA问问问题题题表表表述述述
2
设X是随机向量，且X ∈ Rn×1，这也就说，X里有n个成员，每个成员是一
个随机变量：
X =
(1)
x1
x2
...
xi
...
xn
= A
x1
x2
...
xi
...
xn
s1
s2
...
si
...
sn
1
其中，xi是一个随机变量。
随机变量有诸多特性，殆由概率论和数理统计教科书详述备尽，在此不
一一叙述。
X里的n个随机变量是相互非独立的，在一定的假设下，可以用n个相互
独立的随机变量线性组合重新表达X，也就是说：
(2)
s1
s2
...
si
...
sn
S =
X = AS
S = A−1X
W = A−1
S = W X
(3)
(4)
(5)
(6)
(7)
(8)
则：
又有:
令：
则:
其中，W ∈ Rn×n。
其中，si是一个随机变量，且两两相互独立，A是满秩矩阵，且A ∈ Rn×n。
令:
记录随机向量X的值m次，则形成数据集：
d1,1
d1,2
d2,2
d2,1
...
...
dn,1 dn,2
... d1,m
... d2,m
...
...
... dn,m
D =
其中，D ∈ Rn×m
ICA的目标，就是在只知道D的情况下，估算A，W ，S的值。
实例：在一个大厅里，有n个人在随机聊天。在大厅的不同角落，布
置n个麦克风记录大厅的声音，每秒一个记录，一共记录m秒。麦克风记录
的混合声音，多个麦克风记录不同位置的混合声音。ICA的目标，就是从混
声录音中将每个人的声音分离出来。
3 理理理论论论推推推导导导
由式(7)可知:
si = wi,1 wi,2
... wi,j
... wi,n
x1
x2
...
xi
...
xn
(9)
2
wi = wi,1 wi,2
... wi,j
... wi,n
(10)
令:
则:
(11)
设随机变量si概率密度函数是psi(si)，其中p的右下角si表示随机变量标
si = wiX
示，括号中的si表示自变量。
由于S的n个成员si是相互独立的，所以S的概率密度函数为:
pS(s) =
psi(si)
(12)
设X的概率密度函数是pX (x)，如何根据si的概率密度函数求pX (x)呢？
这是可以做到的。
i=1
设随机向量X的概率分布函数是FX (x)，根据概率分布函数和概率密度
函数的关系可知：
(13)
(14)
(15)
(16)
(17)
n
pX (x) = F
X (x)
同理，设随机向量S的概率分布函数是FS(s)，则：
pS(s) = F
S(s)
根据概率分布函数的定义，有:
FX (x) = P (X < x)
FS(s) = P (S < s)
那么：
pX (x) = F
X (x)
= (P (X < x))
= (P (U < u))
= (P (U < s))
= (P (U < W x))
= (P (S < W x))
= (FS(W x))
= F
S(W x)
= pS(W x)(W x)
= pS(W x)W
= pS(W x)W
= W n
psi(wix)
i=1
3
其中，上式的第2个等号是概率密度函数的定义，第3个等号是做变量等
价代换，以免直接从X变换到S导致思维混乱，第4个等号到第6个等号是
逐步将X代换到S，第7个等号是回到S的概率分布函数定义，第8个等号到
第10个等号是求导。
从第5个等号开始，对整个等式取行列式运算，因为pX (x)一定是标量，
对标量做行列式运算是它自身。那么，到了第10个等号，又因为pS(W X)一
定是标量，所以可以从行列式运算拿到外面。这里避免的问题的是，如果
不对整个等式取行列式，得到的结果是矩阵W 而不是W，这是没有道理
的。
注意，在上式中，x是一个向量，且x ∈ Rn×1，wi ∈ R1×n，psi(si)是一
个单自变量的函数，pX (x)是一个多自变量函数，它的自变量是x里的多个
变量，这样等式左右的每一步就清晰了。
下一步是根据数据集计算W 的值，从概率的角度来说，如果数据集已经
记录，那么让这个数据集出现概率最大的W 就是最优值。
L =
数据集也就是式(8)出现的概率是:
m
(W n
其中，表示连乘，di是D的第i列，也就是：
di,1
di =
i=1
j=1
psj (wjdi))
(18)
(19)
di,2
...
di,n
n
di的物理意义，也就是第i次记录随机向量X得到的n个值，这n个值分别
对应n个xi随机变量。注意，不要把di和xi混淆，前者表示D的一列数据，
后者是粗体表示一个随机变量。
式(18)有最大值，当它取最大值时候的W 就是最优解。如果以梯度下降
法求解，需要计算它对W 的偏导，直接求偏导比较复杂，故对它两端取自
然对数，则:
m
m
i=1
n
lnL =
=
(lnW +
(lnpsj (wjdi)))
j=1
lnpsj (wjdi) + mlnW
(20)
i=1
j=1
当上式取最大值的时候，L也同时取最大值，所以求L的最大值等价于求上
式的最大值。
用梯度下降法求解上式，需要计算 ∂lnL
∂W
。这是一个复杂的过程，先从计
4
算 ∂L
∂wu,v
开始，它表示W 的第u行第v列的一个成员：
n
n
j=1
i=1
m
m
m
i=1
∂lnL
∂wu,v
=
=
=
1
∂psj (wjdi)
psj (wjdi)
∂wu,v
1
∂psj (wjdi)
∂W
∂wu,v
m
W
W (−1)u+vMuv
m
+
+
(21)
psj (wjdi)
∂wu,v
j=1
1
∂psu(wudi)
psu(wudi)
∂wu,v
i=1
+
m
W (−1)u+vMuv
其中，(−1)u+vMuv是wu,v的代数余子式， ∂psu (wuxi)
体形式求解。
∂wu,v
的值要根据psi(si)的具
对于psi(si)，如果在没有任何先验信息的情况下，是无法求解的。如果
要求解上式，需要对它做一定的假设，在合理的假设下，可以达到相当不
错的近似结果。
设随机变量xi的概率分布函数是sigmoid函数，因为它是递增，可微，且
最大值不超过1，也就是说：
Fsi(si) =
1
1 + e−si
那么，概率密度函数就是:
psi(si) = F
si(si) =
esi
(1 + esi)2
所以有：
故:
∂psu(wudi)
∂wu,v
psu(wudi) =
ewudi
(1 + ewudi)2 = ewudi(1 + ewudi)−2
= ewudidi,v(1 + ewudi)−2 − 2ewudi(1 + ewudi)−3ewudidi,v
di,vewudi
(1 + ewudi)2 (1 − 2
1 + ewudi
ewudi
=
)
= di,vpsu(wudi)
1 − ewudi
1 + ewudi
其中，di,v是di的第v行的一个成员。
因此:
∂lnL
∂wu,v
i=1
m
m
m
i=1
i=1
=
=
=
1
∂psu(wudi)
psu(wudi)
∂wu,v
+
1
psu(wudi)
di,vpsu(wudi)
m
W (−1)u+vMuv
1 − ewudi
1 + ewudi
m
+
W (−1)u+vMuv
1 − ewudi
1 + ewudi
di,v
+
m
W (−1)u+vMuv
5
(22)
(23)
(24)
(25)
(26)
现在对上式进行矩阵化， 令：
(27)
其中，K ∈ Rn×m，W ∈ Rn×n，D ∈ Rn×m，那么，ku,i就是K的第u行的
第i列的一个成员， 令：
K = W D
1 − ex
1 + ex
g(x) =
g(k1,1)
g(k2,1)
...
g(k1,2)
g(k2,2)
...
...
g(k1,m)
g(k2,m)
g(kn,1) g(kn,2)
... g(kn,m)
令:
Z = g(K) =
那么，就得到:
(28)
(29)
(30)
(31)
∂lnL
∂wu,v
= zT
u dv +
m
W (−1)u+vMuv
其中，zu是Z的第u行，dv是D的第v列。
于是，对W 而言，则有：
∂lnL
∂W
= ZT D +
m
W (W ∗)T
其中，W ∗是W 的伴随矩阵，(W ∗)T 是W ∗的转置，它的第i行第j列的元素
是wi,j的代数余子式，也就是(−1)i+jMi,j。
根据矩阵和它的伴随阵的性质可知:
W W ∗ = WI
其中，I是单位矩阵。 根据上两式可知:
∂lnL
∂W
= ZT D +
m
W (W ∗)T
W (WW −1)T
m
= ZT D +
= ZT D + m(W −1)T
那么，在梯度下降法求解W 的时候，更新公式是:
W = W + α(ZT D + m(W −1)T )
(32)
(33)
(34)
其中，α是学习速率。
最后的结论简洁且美，Verweile doch, du bist so sch¨on。然并卵，按照
这个结果实现代码，计算结果是不合理的，无法恢复原始信号。于是，在
实现FastICA之后，可以认为本推导缺少一些黑魔法，至于到底缺少什么并
不知道，限于时间关系和实际需求，不再继续研究下去。
6
4 FastICA
FastICA计算性能更好。《Indepdent Componet analysis》一书在第8章给
出了FastICA的算法流程，如下:
Figure 1: FastICA
5 白白白化化化
FastICA需要对数据做白化处理。设x是一个随机变量，存在一个线性变
换V 将它变换成z：
且：
z = V x
E{zzT} = I
(35)
(36)
那么，V 就是白化变换矩阵。
x的 协 方 差 阵 是Cx = E{xxT}，Cx = P DP T ，P 是Cx的 单 位 特 征 向
量，D是Cx的特征值组成的对角阵。那么，V 的值就是：
V = D− 1
2 P T
7
(37)
证明如下： 根据相关性质，有P T = P −1，由于D对角阵，则(D− 1
D− 1
2 ， 那么：
2 )T =
E{V x(V x)T} = E{V xxT V T}
(38)
= E{V P DP T V T}
= E{V P DP T V T}
= E{D− 1
= E{D− 1
= E{I}
= I
2 P T P DP T P (D− 1
2}
2 DD− 1
2 )T}
6 代代代码码码实实实现现现
基于python2.7，matplotlib，numpy实现ICA，主要参考sklean的FastICA实
现。
# !/ usr / bin / env python
# FastICA from ICA book , table 8 . 4
import math
import random
import m a t p l o t l i b . pyplot as plt
from numpy import *
n _ c o m p o n e n t s = 2
def f1 (x , period = 4 ) :
return 0 . 5 * ( x - math . floor ( x / period ) * period )
def c r e a t e _ d a t a () :
# 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 ] ] ,
# mix matrix
A = array ( [ [ 0 .8 , 0 . 2 ] , [ - 0 .3 , - 0 . 7 ] ] , float32 )
return T , S , dot (A , S )
float32 )
8
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
