1.概述

        ICA (IndependentComponent Analysis) 又名独立成分分析，是20世纪90年代发展起来的一种新的信号处理技术，它是从多维统计数据中找出隐因子或独立成分的方法。从线性变换和线性空间的角度来看，源信号为相互独立的非高斯信号，可以看作为线性空间的基信号，而观测信号则为源信号的线性组合,ICA就是在源信号和线性变换均不可知的情况下，从观测的混合信号中估计出数据空间的基本结构或者源信号。

2.问题引入

        考虑下面的情况，也就是经典的鸡尾酒宴会问题（cocktail party problem）。假设你在一个非常拥挤的房间，房间里很多人都在讲话。你的耳朵实际上就像两个麦克风，它们在听房间里不同语音信号的线性组合。你的目标是将混合信号分解为它们对应的组成部分，也就是从采样数据中分辨出每个人说话的信号。这就是鸡尾酒会上的问题，也是盲信号分离(BSS)或盲源分离的一个例子，盲意味着我们对信号的来源一无所知。除了在语音信号处理方面的明显应用外，在分析脑电图EEG和脑磁图信号MEG、金融数据以及其他任何潜在源或因子以线性方式混合在一起的数据集(不一定是时序数据集)时，也会出现这个问题。

下面对上述问题形式化的定义如下。令是在t时刻的观测信号，代表的是源信号。我们假设观测信号与源信号的关系符合如下表达式：

                                        ​​​​​​​        ​​​​​​​        

其中是一个的矩阵，且。在这一节中，假定每一个时间点的数据都是相互独立的，因而我们无需建立时间相关模型（所以我们可以用  替换  下标，但是我们坚持用 ，为了和其它ICA文献保持一致），也就是说我们的目的就是推理源信号，下图就是一个具体的例子：

在这里， 被称作混合矩阵，如果，意味着，传感器的数目与源的数目相等，矩阵就是方阵，为了简单，我们假设噪声的量级是0，即  ，也就是噪声是零均值的白化噪声.

目前看来，这个模型和FA模型好像没什么区别（PCA就是对应于没有噪声，除了我们没有要求是标准正交的）。但是我们在这里会用一个不一样的先验，在PCA里面，我们假设每一个源或者成分是独立的，其先验都是高斯的，即

        ​​​​​​​        ​​​​​​​        ​​​​​​​                ​​​​​​​        

但是这里假设的源分布不在约束为高斯分布，可以放松为是任意的非高斯分布。

不失一般性，我们可以将源分布的方差约束为1，因为任何其他方差都可以通过 的行实现适当的缩放而得到方差为1的结果，由此得到的模型称为独立成分分析(independent component analysis, ICA)。

高斯分布在ICA中不作为先验源的原因是它不能实现源的唯一恢复，具体我们参照上图的图(c)可以很明显的看到。这是因为PCA似然对于源  和混合矩阵  的任何正交变换具有不变性，PCA可以恢复信号所在的最佳线性子空间，但不能唯一地恢复信号本身。

因此在PCA下，我们缺乏一个旋转信息，所以我们能估计出子空间，但是不能唯一确定模型的参数，但是ICA可以。但是ICA需要 是一个方阵，这样它才可逆，如果不是方阵，那么也就不能唯一的恢复出最终的结果，但我们可以去计算出后验。当然不管对于哪一种情况，我们都需要去估计 以及源的分布​​​​​​​。

3.最大似然估计

模型对应的混合矩阵（方阵） 以及白化噪声的ICA模型。和之前一样，可以假设观测值是经过中心化处理的零均值结果，因此可以假设也是0均值的。另外，我们假设观测值已被白化处理，关于白化的解释可参见

所以数据经过中心化处理后被白化的数据，必然有 ,同时我们根据白化噪声处理，

        ​​​​​​​        ​​​​​​​                

所以 必须是标准正交的。这样需要估计的参数就从减少到了，这样就简化了算法的复杂度。

令,称之为recognition weights，而相应的称作generative weights。

令，根据jacobian Matrix的计算可得

        ​​​​​​​        ​​​​​​​        

那么假定有  个独立同分布的数据样本，根据上式得对应的log-likelihood函数如下：

        ​​​​​​​        ​​​​​​​       

其中是  的第  行。因为我们知道是标准正交的，所以式子最左边第一项是可以扔掉的（是个常数），那么我们的NLL(negative log-likelihood）就可以写成如下的形式（针对所有可能的数据，写成期望的形式）：

其中,并且 。由此得到在的行是标准正交的约束条件下，最小化NLL。关于这个的求解可以用梯度法，牛顿法，或者是EM算法。不同的算法各有各的优势，相比能用EM，肯定是EM更好。

4. FastICA 算法

以下将详述fast ICA算法，在这里我们只讲利用近似的牛顿法解ICA问题，主要是基于（Hyvarinen and Oja 2000）的研究结果。

为了简单起见，假设上述模型的隐变量只有一个，并且我们假设源的分布是已知的并且是一样的（这里只有一个隐变量不就意味着我们的源只有一个？）。因此，令，同时再令,那么目标函数以及其梯度和Hessian矩阵如下：

        ​​​​​​​        ​​​​​​​        

其中 是拉格朗日乘子，我们做如下的近似处理，

        ​​​​​​​        ​​​​​​​        

做了这样的近似，Hessian矩阵就是一个常数乘以单位阵，所以我们的牛顿法更新step如下：

也可以写成如下形式

        ​​​​​​​        ​​​​​​​        ​​​​​​​        

在实践中，期望可以被训练集的蒙特卡洛估计所取代，这提供了一个有效的在线学习算法）。执行此更新之后，应该使用以下公式将其投射回约束表面，即

通过上述迭代表达式不停的迭代，迭代算法直至收敛，但这里也可能不收敛（因为的符号的任意性），但是我们向量定义的方向一定会收敛，因此可以假定收敛条件就是会收敛于1。我们可以顺序化的学习，并把当前学习的在之前的子空间中投影的部分给去掉，我们也可以并行的去学习。后面这个可能会更常用，因为不像PCA，这里的特征并没有什么顺序，所以用相同的处理方式处理它们可能会更好。
 

5.代码实现

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

from sklearn.decomposition import FastICA, PCA

""" First generate simulated data """

np.random.seed(0)  # set seed for reproducible results
n_samples = 2000
time = np.linspace(0, 8, n_samples)

s1 = np.sin(2 * time)  # Signal 1 : sinusoidal signal
s2 = np.sign(np.sin(3 * time))  # Signal 2 : square signal
s3 = signal.sawtooth(2 * np.pi * time)  # Signal 3: sawtooth signal

S = np.c_[s1, s2, s3]
S += 0.2 * np.random.normal(size=S.shape)  # Add noise

S /= S.std(axis=0)  # Standardize data
# Mix data
A = np.array([[1, 1, 1], [0.5, 2, 1.0], [1.5, 1.0, 2.0]])  # Mixing matrix
X = np.dot(S, A.T)  # Generate observations

"""Now try to recover the sources"""
# compute ICA
ica = FastICA(n_components=3)
S_ = ica.fit_transform(X)  # Get the estimated sources
A_ = ica.mixing_  # Get estimated mixing matrix

# compute PCA
pca = PCA(n_components=3)
H = pca.fit_transform(X)  # estimate PCA sources

plt.figure(figsize=(9, 6))

models = [X, S, S_, H]
names = ['Observations (mixed signal)',
         'True Sources',
         'ICA estimated sources',
         'PCA estimated sources']
colors = ['red', 'steelblue', 'orange']

for ii, (model, name) in enumerate(zip(models, names), 1):
    plt.subplot(4, 1, ii)
    plt.title(name)
    for sig, color in zip(model.T, colors):
        plt.plot(sig, color=color)

plt.tight_layout()
plt.show()