原 线性独立成分分析(ICA)与鸡尾酒会问

对于鸡尾酒会问题，一种简单的情况如下：有n个人在同时说话，同时又m个声音接收器捕捉到了信号之间的线性组合，于是我们可以得到m组声音数据。那么，如何利用这m组接收到的声音信号恢复成原来的n组独立信号呢？

鸡尾酒会问题示意图

在上个世纪末，鸡尾酒会问题催生了各种盲源分离问题，也有不少机器学习算法被应用于此问题。其中，独立成分分析是简便且有效的一种。

##问题分析 设原始信号为S1、S2、...、Sn，一定存在一个mn的权重矩阵A，使得每个声音接收器对应A中一行。如m=n=2时，令A=[[0.6, 0.4], [0.45, 0.55]]，则第一个麦克风接收到的信号为x1=0.6S1+0.4S2，第二个麦克风接收到的信号为x2=0.45S1+0.55S2。即$$AS=x$$ 将该公式左右同乘以A的逆矩阵$$W=A^{-1}$$，得$$A^{-1}AS=A^{-1}x$$，即$$S=Wx$$，这样问题就解决了（当m≠n时，A不是方阵，此时可以乘A的伪逆）。于是问题来了：如何求A的逆矩阵W？

##独立成分分析(Independent Component Analysis) ICA算法始于Bell和Sejnowski，下面用极大似然估计简单推导一下ICA的求解过程：

###引理：已知x的概率分布$$f_{x}(x)$$，若$$Y=Ax$$，则对其累积分布函数有$$F_{Y}(y)=P{Y\leq y}=P{Ax\leq y}=P{x\leq y/A}=F_{X}(y/A)$$。求导得$$f_{Y}(y)=frac{f_{X}(x/A)}{|A|}$$。

每个信号源$$S_{i}$$都有其概率密度$$p_{i}(s)$$，在第j时刻，向量s的联合概率密度为$p(s_{j})=\prod_{i=1}^{n} p_{i}(s_{i,j})$。通过引理，我们可以将该式转化为$$L(x_{j})=|W|p_{i}(Wx_{j})=|W|\prod_{i=1}^{n} p_{i}(w_{i}^{T}x_{j})$$。 于是可以得到似然函数$$L(x)=\prod_{j=1}^{m} p(x_{j})=\prod_{j=1}^{m} [|W|\prod_{i=1}^{n} p_{i}(w_{i}^{T}x_{j})]$$ 现在似然函数还是个连乘的形式，不便计算。我们将其取对数，得到对数似然函数$$l(x)=\sum_{j=1}^{m} [\sum_{i=1}^{n} ln p_{i}(w_{i}^{T}x_{j}) + ln|W|]$$ 这种形式的似然函数就可以当做我们的成本函数cost function了，即$$J(W)=\sum_{j=1}^{m} [\sum_{i=1}^{n} ln p_{i}(w_{i}^{T}x_{j}) + ln|W|]$$ 但是，式中的$$p_{i}$$，即先验概率$$p(s)$$是未知的，但是我们可以假定它服从某种已知的分布。一般的，我们可以假设$$p(s)$$的累积分布函数$$F(s)$$是sigmoid函数。当然了，假设它们服从正态分布(erf函数)，甚至可以使用tanh、arctan当做s的累积分布函数，这些先验假设都是可行的。 以sigmoid函数为例，令$$g(s)=frac{1}{1+e^{-s}}$$，易知$$g'(s)=g(s)(1-g(s))$$、$$g''(s)=g'(s)(1-2g(s))$$。带入对数似然函数得$$J(W)=\sum_{i=1}^{m}[\sum_{j=1}^{n} ln g'(wx) + ln|W|]$$。 将对数似然函数对$$W$$求偏导，得$$\frac{\partial J(W)}{\partial W}=[1-2g(Wx)]x + (W^{-1})^{T}$$。 到了这一步，相信大家已经知道下面怎么做了：梯度上升到最大值，求出此时的$$W$$即可。也就是$$W += \alpha * (1-2g(W*x) + (W^{-1})^{T})$$。这时候一般梯度上升、随机梯度上升和mini-batch算法就可以派上用场了。 用python实现出来的代码如下：

def sigmoid(x):
    return 1.0 / (1 + np.exp(-x))

def ica(x):
    # x = np.array(x)
    n = len(x[0])
    m = len(x)
    w = np.eye(n)
    iw = np.zeros([n, n])
    w1 = np.zeros([n, n])
    alpha = 0.001

    for time in range(200):
        for i in range(m):
            for j in range(n):
                # 使用sigmoid函数作为CDF
                t = 1 - 2 * sigmoid(np.dot(w[j], x[i]))
                w1[j] = t * x[i]
                # 使用arctan函数作为CDF
                # t = -2 * np.dot(w[j], x[i]) / (1 + np.dot(w[j] ** 2, x[i] ** 2))                
                # w1[j] = t * x[i]
                # 使用tanh函数作为CDF
                # t = -2 * x[i] * np.tanh(np.dot(w[j], x[i]))
                # w1[j] = t
                # 使用erf函数(正态分布)的CDF
                # t = -2 * np.dot(w[j], x[i])
                # w1[j] = t * x[i]
            iw = np.linalg.inv(w.transpose())
            w1 += iw * alpha
            w1 = w1.transpose() * alpha
            w += w1
        # print time, ":\t"
        # print w
    return w

##仿真结果 首先，使用python模拟，生成一个扫描波、一个正弦波如下；使用权重矩阵$$A=\begin{pmatrix} 0.6 & 0.4 \ 0.45 & 0.55 \end{pmatrix}\$$对原始信号进行加权，得到的混合信号如下：

原始信号：

原始信号

混合信号：

混合信号

import numpy as np
from matplotlib import pyplot as plt

def show(s1, s2, label1, label2):
    x = [i for i in range(1000)]
    plt.plot(x, s1, 'r', label = label1)
    plt.hold(True)
    plt.plot(x, s2, 'g', label = label2)
    plt.hold()
    plt.legend()
    plt.show()

if __name__ == "__main__":
    s1 = np.array([np.sin(float(x) / 20) for x in range(1000)])
    s2 = np.array([float(x) / 50 for x in range(50)] * 20)
    show(s1, s2, 'Signal 1', 'Signal 2')

    A = np.array([[0.6, 0.4], [0.45, 0.55]])
    x = np.dot(A, np.array([s1, s2])).transpose()
    show(x[:,0], x[:,1], 'Mic 1', 'Mic 2')
    # x -= x.mean()
    w = ica(x)
    [ps1, ps2] = np.dot(x, w).transpose()
    show(ps1, ps2, 'Component 1', 'Component 2')

将混合信号x输入给ICA算法，并使用得到的W进行还原，得到的结果如下：

ICA_result

虽然有些局部失真，但得到这种结果已经很不错了。下面是使用其他CDF函数得到的结果：

arctan as CDF

erf as CDF

tanh as CDF

陈政/arc001 原创作品转载请注明出处

我的博客即将搬运同步至腾讯云+社区，邀请大家一同入驻：https://cloud.tencent.com/developer/support-plan?invite_code=6klbojckkbuw