ICA的著名应用是盲源分离,于是这里就以盲源分离为例子进行说明。
题目
假设n个人面前有n个话筒,然后这n个人说话时这n个话筒进行录音,这n个人说了m句话,最后从这n个话筒中收集一些录音,目标:从这些录音中分离出每个人的声音。
如下图所示:
下面开始解题。
题目整理
首先将信息整理如下:
源信号,如下图所示:
PS:源信号是未知的(知道还求啥),但为了分析还是需要把源信号表示出来。
第一个人是s1,他第一个时刻发送的信号是s11,第二个时刻发送的信号是s12,以此类推第m个时刻发送的时刻是s1m。第二个人,第三个人,….,第n个人同理。
因此列向量中的每个元素就代表每个人,行向量的每个元素代表每个时刻。
观测信号,如下图所示:
解读方式和源信号一样。
这里我们需要的是源信号和观测信号行向量,因为我们只知道每个时刻大伙说了什么话。于是sm和xm都是n维的向量。
然后,我们假设源信号经过线性加权得到了观测信号, 即:sm乘上一个权值A后得到了xm,即:xm = Asm。接下来为了表示所有的sm和xm,就用s代表sm组成的矩阵,x同理,于是式子更新为x = A·s。
到此,模型出来了,即:
x =A·s
后面要做的是如何求解该模型,但因为除了观测信号x之外,A未知(不知道怎么混合的),s也未知(源信号不知道)。
别忘了我们的目标:分离源信号。所以要把上面的式子转换成源信号=什么什么,即,转换成:
s = A-1·x,
为了方便书写,用W表示A的逆矩阵(或广义逆矩阵),于是上式变成:
s = W·x
之后,我们的核心目标就是如何求W(x已知,s是目标)。
到此题目整理完毕。
PS:这里可以看出ICA的一个不足:顺序和振幅不稳定,即:对于ICA模型x = A·s,甲求得A=A1,S=S1;乙求得A=-A1/2,S=-2S1,那甲乙的结果都能得出源信号,但甲和乙求得的结果的顺序就是反的(有个负号),振幅甲是乙的1/2。
进行假设
首先,先做两个假设:
假设1:源信号彼此间独立。
假设2:源信号是非高斯分布。
上面的假设1不用解释,假设2是的根据是中心极限定理,即:一组有确定方差的独立随机变量的和趋近于高斯分布。
不过假设2可以多解释些,嗯….如果打个比方的话:假设源信号是一个个颜色时,高斯分布就是颜色混合后的那一坨脏色,任何颜色(源信号)一旦和其他颜色混合(成为高斯分布),那该颜色就再也无法被识别。换句话说:给定随机变量A和B时,A+B比A或B更接近高斯分布,再换句话说:如果能够找到一组独立信号,或者说找到了一组“最不像高斯分布”的信号,则它们极有可能是源信号。
开始推导ICA
推导的工具就是最大似然估计。
假定第i个源信号的概率密度函数为pi(s),第j时刻的n个源信号记做向量sj,则在j时刻向量sj 的联合密度为:
根据xj = A·sj和http://blog.csdn.net/xueyingxue001/article/details/51915390中“复合函数的概率密度”的知识,得:
从而得到似然函数:
进一步得到对数似然:
从而,目标函数为:
但是,对于目标函数,除了我们相求的W外,我们还不知道每个人说话的概率密度,即:pi(WiT·xj),这点比较郁闷,不过我们可以做一个有很强根据的假定:
假定源信号的累积概率密度函数是Logistics/Sigmoid函数。
为什么说这个假定有很强的依据?
如下图所示:
这是Sigmoid的函数图像,从图像中可以很清楚的看到,在负无穷时Sigmoid函数接近0,正无穷时接近1,所以,我们就可以将其看成是某个随机变量的累积概率。所以在没有任何先验信息时,我们不妨认为源信号的累积概率密度大体符合Sigmoid函数。(当然也可以假设成其他函数,最后在对比效果。)
接下来给出Sigmoid函数的公式和求导(即:概率密度):
公式:F(x) = 1/(1+e-x)
概率密度/求导:f(x) = F’(x) =F(x)(1-F(x))
为了之后的步骤方便,再求次导数:
概率密度在求导:f’(x) = F’’(x) = f(x)(1-2F(x))
那按照极大似然估计的知识,接下来我们就直接对J(W)对wij求偏导:
PS:上面第二行利用了http://blog.csdn.net/xueyingxue001/article/details/51915390中“标量对矩阵求导”的知识
既然对wij求偏导是上面的结果,那对整个W求偏导就是把上面的结果写成向量的形式,即:
于是梯度下降的公式就有了,即:
PS1:如果W不是方阵,则只要把W-1替换成W+就可以了。
PS2:α是学习率。
使用时:
虽然教科书中说在使用ica时,对于一个数据要经过如下步骤:
原始数据->中心化去均值->白化->PCA降维->ICA
但实时上:
若噪声不太强,PCA可以忽视
对于白化,有些数据不做此步骤会更好。
最后,从矩阵的角度总结ICA:
对于ICA的模型:x = A·s
假设x是个m*n的矩阵,A是个m*k的矩阵,那s就需要是个k*n的矩阵,于是ICA就是把x分解成2个矩阵的乘积,这两个矩阵就是求出来的源信号。
代码:
#-*-coding:utf-8-*- import sys import time import math import datetime import threading import numpy asnp importmatplotlib.pyplot as plt def show_data(x,y): num_list_x = np.arange(len(x)) plt.figure(figsize=(20, 10)) plt.xlim(0, len(num_list_x)+50) plt.plot(num_list_x, x, color='b',linestyle='-', label='x') plt.legend() plt.show() if y != None: plt.figure(figsize=(20, 10)) num_list_y = np.arange(len(y)) plt.xlim(0, len(num_list_y)+50) plt.plot(num_list_y, y, color='r',linestyle='-', label='y') plt.legend() plt.show() def logistic(t): return 1.0/(1+np.exp(-t)) defn_multiply(t, x): return t * x deftrans_inverse(x): if not isinstance(x, np.ndarray): x = np.array(x) return np.linalg.inv(x.transpose()) def ica(x): m = len(x) # 样本数目,这里是 2 个 n = len(x[0]) # mic 数目,这里是 1000 个 # 建立个 1000*1000 的列表 w = [[0.0]*n for t in range(n)] # 建立个 1000*1000 的列表 iw = [[0.0]*n for t in range(n)] # 建立个 1000*1000 的列表 w1 = [[0.0]*n for t in range(n)] # 将对角线初始化为1 for i in range(n): w[i][i] = 1 # 初始化学习率为 alpha = 0.001 # shuffle(x) # 迭代次数最多不超过 200 次 for time in range(200): # 分离出"样本数目"个样本 for i in range(m): for j in range(n): t = np.dot(w[j], x[i]) t = 1 - 2*logistic(t) w1[j] = n_multiply(t,x[i]) # w1[j] = t*x[i] iw = trans_inverse(w) # iw = w^T^(-1) iw = w1 + iw w1 = np.add(w1, iw) #w1 += iw w1 = np.dot(w1, alpha) w = np.add(w, w1) print time, ":\t", w return w def mix(A, x,y): mix = np.dot(A, np.array([x, y])) for i in xrange(len(mix[0])): mix[0][i] = mix[0][i] + mix[1][i] #show_data(mix[0], None) return mix def decode(w,x): ps = np.dot(x, w) return ps[0], ps[1] if __name__ =="__main__": s1 =[math.sin(float(x)/20) for x in range(0, 1000, 1)] s2 = [float(x)/50 for x in range(0, 50, 1)]* 20 #s1 = [math.sin(float(x)/20) for x inrange(0, 100, 1)] #s2 = [float(x)/50 for x in range(0, 50,1)] * 2 #show_data(s1, s2) # 假定有N个人在说话,则在任何一个时刻,原始的源信号是N维的,混合信号也是N维的。从而,混合矩阵A和解混矩阵W都是N*N维的。——这个维度,和观测时间M无关。 # 即:混合矩阵A必须是N*N维,才能使得n维的原始信号乘上A后可以得到个n维的混合信号,解混矩阵同理。 A = [[0.6, 0.4], [0.45, 0.55]] # 混合矩阵 x = mix(A, s1, s2) #s1/s2线性加权得到输入数据x # 去均值,这个老师的效果不错,但我这里反而变差了.... # 实际运用中还是要根据实际场合看看是否选取 x_mean = np.mean(x) for i in xrange(x.shape[0]): for j in xrange(x.shape[1]): x[i][j] = x[i][j] - x_mean # ica w = ica(x) # 解混 [ps1, ps2] = decode(w, x) show_data(ps1, ps2)