独立成分分析（Independent Component Analysis）

简介

综合转载于：

• 独立成分分析（Independent Component Analysis）
• ICA扩展描述
• 详解独立成分分析
• 独立成分（ICA）分析
• 独立成分分析Independent component analysis(ICA)

问题

上节提到的PCA是一种数据降维的方法，但是只对符合高斯分布的样本点比较有效，那么对于其他分布的样本，有没有主元分解的方法呢？

经典的鸡尾酒宴会问题（cocktail party problem）。假设在party中有 $n$ 个人，他们可以同时说话，我们也在房间中一些角落里共放置了 $n$ 个声音接收器（Microphone）用来记录声音。宴会过后，我们从 $n$ 个麦克风中得到了一组数据 $\{x^{(i)}(x_1^{(i)},x_2^{(i)},\dots,x_n^{(i)}),i=1,2,\dots,m\}$ ，$i$ 表示采样的时间顺序，也就是说共得到了 $m$ 组采样，每一组采样都是 $n$ 维的。我们的目标是单单从这 $m$ 组采样数据中分辨出每个人说话的信号。

将第二个问题细化一下，有 $n$ 个信号源 $\text{s}(s_1,s_2,\dots,s_n)^T$，$s\in \mathbb R^n$，每一维都是一个人的声音信号，每个人发出的声音信号独立。$\rm A$ 是一个未知的混合矩阵（mixing matrix），用来组合叠加信号 $\rm s$，那么

$$\rm x=As$$

$\rm x$ 的意义在上文解释过，这里的 $\rm x$ 不是一个向量，是一个矩阵。其中每个列向量是 $x^{(i)}$，$x^{(i)}=\text{A}s^{(i)}$ 表示成图就是

$x^{(i)}$ 的每个分量都由 $s^{(i)}$ 的分量线性表示。$\rm A$ 和 $\rm s$ 都是未知的，$\rm x$ 是已知的，我们要想办法根据 $\rm x$ 来推出 $\rm s$。这个过程也称作为盲信号分离。

令 $\rm W=A^{-1}$，那么 $s^{(i)}=\text A^{-1}x^{(i)}=\rm Wx^{(i)}$，将 $\rm W$ 表示成

其中 $w_i\in\mathbb R^n$，其实就是将 $w_i$ 写成行向量形式。那么得到：

$$s_j^{(i)}=w_j^Tx^{(i)}$$

ICA的不确定性

由于 $\rm W$ 和 $\rm s$ 都不确定，那么在没有先验知识的情况下，无法同时确定这两个相关参数。比如上面的公式 $\rm s=Wx$。当 $\rm W$ 扩大两倍时，$\rm s$ 只需要同时扩大两倍即可，等式仍然满足，因此无法得到唯一的 $\rm s$。同时如果将人的编号打乱，变成另外一个顺序，如上图的蓝色节点的编号变为3,2,1，那么只需要调换 $\rm A$ 的列向量顺序即可，因此也无法单独确定 $\rm s$。这两种情况称为原信号不确定。

还有一种ICA不适用的情况，那就是信号不能是高斯分布的。假设只有两个人发出的声音信号符合多值正态分布，$\text{s}\sim N(0,1)$，I是2*2的单位矩阵，$\rm s$ 的概率密度函数就不用说了吧，以均值0为中心，投影面是椭圆的山峰状（参见多值高斯分布）。因为$ \rm x = As$，因此，$\rm x$ 也是高斯分布的，均值为0，协方差为 $\rm E[xx^T]=E[Ass^TA^T]=AA^T$。

令 $\rm R$ 是正交阵（$\rm RR^T=R^TR=I$），$\rm A’=AR$。如果将 $\rm A$ 替换成 $\rm A’$ 。那么 $\rm x’ =A’s$。$\rm s$ 分布没变，因此 $\rm x’$ 仍然是均值为0，协方差。因此$\rm E[x’x’^T]=E[A’ss^TA’^T]=E[ARss(AR)^T]=ARR^TA^T=AA^T$，不管混合矩阵是 $\rm A$ 还是 $\rm A’$，$\rm x$ 的分布情况是一样的，那么就无法确定混合矩阵，也就无法确定原信号。

密度函数和线性变换

在讨论ICA具体算法之前，我们先来回顾一下概率和线性代数里的知识。

假设我们的随机变量 $\rm s$ 有概率密度函数 $p_s(s)$（连续值是概率密度函数，离散值是概率）。为了简单，我们再假设$\rm s$ 是实数，还有一个随机变量 $\rm x=As$，$\rm A$ 和 $\rm x$ 都是实数。令 $p_x$ 是 $\rm x$的概率密度，那么怎么求 $p_x$ ？

令 $\rm W=A^{-1}$，首先将式子变换成 $\rm s = Wx$，然后得到 $p_x(x)=p_s(\text Ws)$，求解完毕。可惜这种方法是错误的。比如 $\rm s$ 符合均匀分布的话（$\rm s\sim Uniform[0,1]$），那么 $\rm s$ 的概率密度是 $p_s(s)=1\{0\le s\le 1\}$，现在令 $\rm A=2$，即 $\rm x=2s$，也就是说 $\rm x$ 在[0,2]上均匀分布，可知 $p_x(x)=0.5$。然而，前面的推导会得到 $p_x(x)=p_s(0.5s)=1$。正确的公式应该是

$$p_x(x)=p_s(\text Wx)|\text W|$$

推导方法：

$$\text F_{\text X}(\text a)=\text P(\text X\le\text a)=\text P(\text A\text S\le \text x)=\text P(\text S\le \text W\text x) = \text F_{\text s}(\text W\text x)\\ p_x(x)=\text F'_{\text X}(\text W\text x)=p_s(\text W\text x)|\text W|$$

更一般地，如果 $\rm s$ 是向量，$\rm A$ 可逆的方阵，那么上式子仍然成立。

ICA算法

ICA算法归功于Bell和Sejnowski，这里使用最大似然估计来解释算法，原始的论文中使用的是一个复杂的方法Infomax principal。

我们假定每个 $s_i$ 有概率密度 $p_s$，那么给定时刻原信号的联合分布就是

$$p(s) = \prod_{i=1}^n p_s(s_i)$$

这个公式代表一个假设前提：每个人发出的声音信号各自独立。有了 $p(s)$，我们可以求得 $p(x)$：

$$p(x)=p_s(\text Wx)|\text W| = |\text W|\prod_{i=1}^np_s(w_i^Tx_i)$$

左边是每个采样信号 $x$（$n$ 维向量）的概率，右边是每个原信号概率的乘积的 $|\text W|$ 倍。

前面提到过，如果没有先验知识，我们无法求得 $\rm W$ 和 $\rm s$。因此我们需要知道 $p_s(s_i)$，我们打算选取一个概率密度函数赋给s，但是我们不能选取高斯分布的密度函数。在概率论里我们知道密度函数 $p(x)$ 由累计分布函数（cdf）$F(x)$ 求导得到。$F(x)$ 要满足两个性质是：单调递增和在[0,1]。我们发现sigmoid函数很适合，定义域负无穷到正无穷，值域0到1，缓慢递增。我们假定 $\rm s$ 的累积分布函数符合sigmoid函数。

$$g(\text s) = \frac{1}{1+e^{-s}}$$

求导后：

$$p_s(s) = g'(s)=\frac{e^s}{(1+e^s)^2}$$

这就是 $\rm s$ 的密度函数。这里 $\rm s$ 是实数。

如果我们预先知道 $\rm s$ 的分布函数，那就不用假设了，但是在缺失的情况下，sigmoid函数能够在大多数问题上取得不错的效果。由于上式中 $p_s(s)$ 是个对称函数，因此 $\text E[\rm s]=0$ （$\rm s$ 的均值为0），那么 $\rm E[x]=E[As]=0$，$\rm x$ 的均值也是0。

知道了 $p_s(s)$，就剩下 $\rm W$ 了。给定采样后的训练样本 $\{x^{(i)}(x_1^{(i)},x_2^{(i)},\dots,x_n^{(i)}),i=1,2,\dots,m\}$，样本对数似然估计如下：

使用前面得到的 $\text x$ 的概率密度函数，得

$$l(\text W)=\sum_{i=1}^m\Big(\sum_{j=1}^n\log g'(w_j^Tx^{(i)})+\log|\text W|\Big)$$

大括号里面是 $p(x^{(i)})$。

接下来就是对 $\rm W$ 求导了，这里牵涉一个问题是对行列式 $|\rm W| $ 进行求导的方法，属于矩阵微积分。这里先给出结果，在文章最后再给出推导公式。

$$\rm \nabla_W|W|=|W|(W^{-1})^T$$

最终得到的求导后公式如下，$\log g’(s)$ 的导数为 $1-2g(s)$：

其中 $\alpha$ 是梯度上升速率，人为指定。

当迭代求出 $\rm W$ 后，便可得到 $s^{(i)} = \text Wx^{(i)}$ 来还原出原始信号。

注意：我们计算最大似然估计时，假设了 $x^{(i)}$ 与 $x^{(j)}$ 之间是独立的，然而对于语音信号或者其他具有时间连续依赖特性（比如温度）上，这个假设不能成立。但是在数据足够多时，假设独立对效果影响不大，同时如果事先打乱样例，并运行随机梯度上升算法，那么能够加快收敛速度。

回顾一下鸡尾酒宴会问题，$\rm s$ 是人发出的信号，是连续值，不同时间点的s不同，每个人发出的信号之间独立（$s_i$ 和 $s_j$ 之间独立）。$\rm s$ 的累计概率分布函数是 sigmoid 函数，但是所有人发出声音信号都符合这个分布。 $\rm A$ 代表了 $\rm s$ 相对于 $\rm x$ 的位置变化，$\rm x$ 是 $\rm s$ 和 $\rm A$ 变化后的结果。

实例

$\rm s =2$ 时的原始信号

观察到的 $\rm x$ 信号

使用ICA还原后的 $\rm s$ 信号

行列式的梯度

对行列式求导，设矩阵 $\rm A$ 是 $n×n$ 的，我们知道行列式与代数余子式有关，

$$|\text A| = \sum_{i=1}^n(-1)^{i+j}\text A_{ij}|A_{\backslash i,\backslash,j}|\space (\text{for any }j \in 1,\dots,n)$$

$\text A_{\backslash i,\backslash j}$是去掉第 $i$ 行，第 $j$ 列后的余子式，那么对 $a_{k,l}$ 求导得：

$$\frac{\partial}{\partial a_{k,l}}|\text A|=\frac{\partial}{\partial a_{k,l}}\sum_{i=1}^n(-1)^{i+j}a_{ij}|A_{\backslash i, \backslash j}| = (-1)^{k+l}|\text A_{\backslash k, \backslash j}| = (A^\ast)_{lk}$$

$\rm A^\ast$ 是伴随矩阵，因此

$$\rm \nabla_A|A| = (A^\ast)^T = |A|(A^{-1})^T$$

ICA算法扩展描述

上面介绍的内容基本上是讲义上的，与我看的另一篇《Independent Component Analysis: Algorithms and Applications》（Aapo Hyvärinen and Erkki Oja）有点出入。下面总结一下这篇文章里提到的一些内容（有些我也没看明白）。

首先里面提到了一个与“独立”相似的概念“不相关（uncorrelated）”。Uncorrelated属于部分独立，而不是完全独立，怎么刻画呢？

如果随机变量 $y_1$ 和 $y_2$ 是独立的，当且仅当 $p(y_1,y_2)=p(y_1)p(y_2)$。

如果随机变量 $y_1$ 和 $y_2$ 是不相关的，当且仅当 $\text E(y_1, y_2) = E(y_1)E(y_2)$

第二个不相关的条件要比第一个独立的条件“松”一些。因为独立能推出不相关，不相关推不出独立。

证明如下：

$$p_1(y_1) = \int p(y_1, y_2)dy_2,\\ p(y_1, y_2) = p_1(y_1)p_2(y_2).$$ $$\begin{aligned} &E\left\{h_{1}\left(y_{1}\right) h_{2}\left(y_{2}\right)\right\}\\ =& \iint h_{1}\left(y_{1}\right) h_{2}\left(y_{2}\right) p\left(y_{1}, y_{2}\right) d y_{1} d y_{2}\\ =&\iint h_{1}\left(y_{1}\right)p_{1}\left(y_{1}\right) h_{2}\left(y_{2}\right) p_{2}\left(y_{2}\right) d y_{1} d y_{2}\\ =&\int h_{1}\left(y_{1}\right) p_{1}\left(y_{1}\right) d y_{1} \int h_{2}\left(y_{2}\right) p_{2}\left(y_{2}\right) d y_{2}\\ =&E\left\{h_{1}\left(y_{1}\right)\right\} E\left\{h_{2}\left(y_{2}\right)\right\} \end{aligned}$$

反过来不能推出。

比如，$y_1$ 和 $y_2$ 的联合分布如下(0,1)，(0, -1)，(1,0)，(-1,0)。

$$E(y_1,y_2) = E(y_1)E(y_2)=0$$

因此 $y_1$ 和 $y_2$ 不相关，但是

$$E(y_1^2y_2^2)=0\neq \frac{1}{4}=E(y_1^2)E(y_2^2)$$

因此 $y_1$ 和 $y_2$ 不满足上面的积分公式，$y_1$ 和 $y_2$ 不是独立的。

上面提到过，如果 $s^{(i)}$ 是高斯分布的，$\rm A$ 是正交的，那么 $x^{(i)}$ 也是高斯分布的，且 $x^{(i)}$ 与 $x^{(j)}$ 之间是独立的。那么无法确定 $\rm A$，因为任何正交变换都可以让 $x^{(i)}$ 达到同分布的效果。但是如果 $s^{(i)}$ 中只有一个分量是高斯分布的，仍然可以使用ICA。

那么ICA要解决的问题变为：如何从 $\rm x$ 中推出 $\rm s$ ，使得 $\rm s$ 最不可能满足高斯分布？中心极限定理告诉我们：大量独立同分布随机变量之和满足高斯分布。

我们一直假设的是 $x^{(i)}$ 是由独立同分布的主元 $s^{(i)}$ 经过混合矩阵 $\rm A$ 生成。那么为了求 $s^{(i)}$，我们需要计算 $s^{(i)}$ 的每个分量 $y_j^{(i)}=w_j^Tx^{(i)}$。定义 $z_j=A^Tw_j$，那么 $y_j^{(i)}=w_j^Tx^{(i)}=w_j^T\text A s^{(i)}=z_j^Ts^{(i)}$，之所以这么麻烦再定义 $\rm z$ 是想说明一个关系，我们想通过整出一个 $w_j$ 来对 $x^{(i)}$ 进行线性组合，得出 $\rm y$ 。而我们不知道得出的 $\rm y$ 是否是真正的 $\rm s$ 的分量，但我们知道 $\rm y$ 是 $\rm s$ 的真正分量的线性组合。由于我们不能使 $\rm s$ 的分量成为高斯分布，因此我们的目标求是让 $\rm y$（也就是 $w_j^Tx^{(i)}$）最不可能是高斯分布时的 $\rm w$。那么问题递归到如何度量 $\rm y$ 是否是高斯分布的了。

一种度量方法是kurtosis方法，公式如下：

$$\text{kurt}(y)=E\{y^4\}-3(E\{y^2\})^2$$

如果 $y$ 是高斯分布，那么该函数值为0，否则绝大多数情况下值不为0。但这种度量方法不怎么好，有很多问题。看下一种方法：

负熵（Negentropy）度量方法。我们在信息论里面知道对于离散的随机变量 $Y$，其熵是

$$H(Y)=-\sum_i P(Y=a_i)\log P(Y=a_i)$$

连续值是：

$$H(y)=-\int f(y)\log f(y)dy$$

在信息论里有一个强有力的结论是：高斯分布的随机变量是同方差分布中熵最大的。也就是说对于一个随机变量来说，满足高斯分布时，最随机。定义负熵的计算公式如下：

$$J(y)=H(y_{\text{gauss}})-H(y)$$

也就是随机变量 $\rm y$ 相对于高斯分布时的熵差，这个公式的问题就是直接计算时较为复杂，一般采用逼近策略。

$$J(y) \approx \frac{1}{12}E\{y^3\}^2+\frac{1}{48}\text{kurt}(y)^2$$

这种逼近策略不够好，作者提出了基于最大熵的更优的公式：

$$J(y)\approx \sum_{i=1}^p k_i[E\{G_i(y)\}-E\{G_i(v)\}]^2$$

之后的FastICA就基于这个公式。另外一种度量方法是最小互信息方法：

$$I(y_1,y_2,\dots,y_m)=\sum_{i=1}^mH(y_i)-H(\rm y)$$

这个公式可以这样解释，前一个 $H$ 是 $y_i$ 的编码长度（以信息编码的方式理解），第二个 $H$ 是 $\rm y$ 成为随机变量时的平均编码长度。之后的内容包括FastICA就不再介绍了，我也没看懂。

ICA的投影追踪解释

投影追踪在统计学中的意思是去寻找多维数据的“interesting”投影。这些投影可用在数据可视化、密度估计和回归中。比如在一维的投影追踪中，我们寻找一条直线，使得所有的数据点投影到直线上后，能够反映出数据的分布。然而我们最不想要的是高斯分布，最不像高斯分布的数据点最interesting。这个与我们的ICA思想是一致的，寻找独立的最不可能是高斯分布的 $\rm s$。

在下图中，主元是纵轴，拥有最大的方差，但最interesting的是横轴，因为它可以将两个类分开（信号分离）。

ICA算法的前处理步骤

1. 中心化：也就是求 $\rm x$ 均值，然后让所有 $\rm x$ 减去均值，这一步与PCA一致。

2. 白化：目的是将 $\rm x$ 乘以一个矩阵变成 $\tilde{\text x}$，使得 $\tilde {\text x}$ 的协方差矩阵是 $I$。原始的向量是 $\rm x$ 。转换后的是 $\tilde {\text x}$。$\tilde {\text x}$ 的协方差矩阵是 $I$，即

   $$E\{\tilde {\text x}\tilde {\text x}^T\}=I$$

   我们只需用下面的变换，就可以从 $\rm x$ 得到想要的 $\tilde {\text x}$。

   $$\tilde {\text x}=\text E\text D^{-1/2}\text E^T\text x$$

   其中使用特征值分解来得到 $\rm E$（特征向量矩阵）和$ \rm D$（特征值对角矩阵），计算公式为

   $$E\{\text x\text x^T\}=\text E\text D\text E^T$$

   下面用个图来直观描述一下：

   假设信号源 $s_1$ 和 $s_2$ 是独立的，比如下图横轴是 $s_1$，纵轴是 $s_2$，根据 $s_1$ 得不到 $s_2$。

   

   我们只知道他们合成后的信号 $\rm x$，如下

   

   此时 $x_1$ 和 $x_2$ 不是独立的（比如看最上面的尖角，知道了 $x_1$ 就知道了 $x_2$ ）。那么直接代入我们之前的极大似然概率估计会有问题，因为我们假定 $\rm x$ 是独立的。

   因此，漂白这一步为了让 $\rm x$ 独立。漂白结果如下：

   

   可以看到数据变成了方阵，在 $\tilde {\text x}$ 的维度上已经达到了独立。

   然而这时 $\text x$ 分布很好的情况下能够这样转换，当有噪音时怎么办呢？可以先使用前面提到的PCA方法来对数据进行降维，滤去噪声信号，得到k维的正交向量，然后再使用ICA。

3. 进一步的预处理：给定数据集的ICA成功与否可能会跟特定应用的预处理步骤有关。比如数据中包含时间信号，那么带通滤波也将会是很有用的。如果我们对观测信号 $x_i(t)$ 进行线性滤波得到新的信号 $x_i^\ast(t)$，具有相同混合矩阵的ICA模型仍然适用于 $x_i^\ast(t)$。观测 $\text x(1),\dots,\text x(T)$ 是矩阵 $\rm X$ 的列，对于 $\rm S$ 也是如此，那么ICA模型可被表示为：

   $$\rm X=AS$$

   现在，$\rm X$ 的时间滤波相当于 $\rm X$ 从右边乘以一个矩阵 $\rm M$，如下：

   $$\rm X^\ast = XM = ASM = AS^\ast$$

   这表明ICA模型依然有效。

FastICA算法

在前面的小节中，介绍了非高斯性的不同度量方式，也就是ICA估计的目标函数。在实际中，还需要一种最大化对比函数的算法。这一小节介绍一种非常有效的最大化方法，这里假设数据都是经过预处理的。

单元FastICA

首先来看单元FastICA，通过“单元”，我们指的是计算单元，最终是人工神经元，具有神经元能够通过学习规则更新的权值向量 $\rm w$。FastICA学习规则寻找方向，也就是使投影 $\rm w^Tx$ 最大化非高斯性的单位向量 $\rm w$，非高斯性通过负熵 $J(\rm w^Tx)$ 的近似度量。回想一下 $\rm w^Tx$ 的方差必须是统一的，对于白化的数据，这相当于将 $\rm w$ 的范数统一。FastICA是基于定点迭代的方案的寻找 $\rm w^Tx$ 非高斯性的最大值。这也可以用近似牛顿迭代法推导出来。

FastICA算法的基本形式如下：

1. 选择初始的权值向量 $\rm w$ （随机选择）
2. 令 $\mathbf{w}^{+}=E\left\{\mathbf{x} g\left(\mathbf{w}^{T} \mathbf{x}\right)\right\}-E\left\{g^{\prime}\left(\mathbf{w}^{T} \mathbf{x}\right)\right\} \mathbf{w}$
3. 令 $\mathbf{w}=\mathbf{w}^{+} /\left|\mathbf{w}^{+}\right|$
4. 如果不收敛，就返回步骤2

收敛是指 $\rm w$ 的新旧值指向同一个方向，也就是，它们的点积等于1，向量不一定收敛到同一个点，因为 $\rm w$ 和 $\rm -w$ 定义的方向是相同的，这也是因为独立分量只能定义为乘法符号。

FastICA推导如下，首先 $\rm w^Tx$ 的负熵的近似的最大值是从 $E\left\{G\left(\mathbf{w}^{T} \mathbf{x}\right)\right\}$ 获得的，在约束 $E\left\{\left(\mathbf{w}^{T} \mathbf{x}\right)^{2}\right\}=|\mathbf{w}|^{2}=1$ 下的 $E\left\{G\left(\mathbf{w}^{T} \mathbf{x}\right)\right\}$的最优值是从下式中获得的：

$$E\left\{\mathbf{x} g\left(\mathbf{w}^{\mathrm{T}} \mathbf{x}\right)\right\}-\beta \mathbf{w}=0$$

用牛顿法来解上面的等式，令等式左边为 $F$，然后可以获得雅可比矩阵 $J F(\mathbf{w})$：

$$J F(\mathbf{w})=E\left\{\mathbf{x} \mathbf{x}^{\mathrm{T}} g^{\prime}\left(\mathbf{w}^{\mathrm{T}} \mathbf{x}\right)\right\}-\beta \mathbf{I}$$

为了简化矩阵，我们对上式中的第一部分取近似，由于数据是球形的，一个合理的近似是 $E\left\{\mathbf{x} \mathbf{x}^{\mathrm{T}} g^{\prime}\left(\mathbf{w}^{\mathrm{T}} \mathbf{x}\right)\right\} \approx E\left\{\mathbf{x} \mathbf{x}^{\mathrm{T}}\right\} E\left\{g^{\prime}\left(\mathbf{w}^{\mathrm{T}} \mathbf{x}\right)\right\}=E\left\{g^{\prime}\left(\mathbf{w}^{\mathrm{T}} \mathbf{x}\right)\right\} \mathbf{I}$ 因此雅可比矩阵变成了对角矩阵，并且很容易求其逆，从而得到下面的近似牛顿迭代：

$$\mathbf{w}^{+}=\mathbf{w}-\left[E\left\{\mathbf{x} g\left(\mathbf{w}^{\mathrm{T}} \mathbf{x}\right)\right\}-\beta \mathbf{w}\right] /\left[E\left\{g^{\prime}\left(\mathbf{w}^{\mathrm{T}} \mathbf{x}\right)\right\}-\beta\right]$$

通过对上式两边同时乘以 $\beta-E\left\{g^{\prime}\left(\mathbf{w}^{\mathrm{T}} \mathbf{x}\right)\right\}$ 可以进一步简化。实际上，FastICA的期望必须由他们的估计来代替。 自然估计当然是相应的样本均值。 理想情况下，应该使用所有可用数据，但这通常不是一个好主意，因为计算可能变得过于苛刻。 然后可以使用较小的样本估计平均值，其大小可能对最终估计的准确性具有相当大的影响。 应在每次迭代时单独选择样本点。 如果收敛不令人满意，则可以增加样本量。

多元FastICA

前面讲的单单元算法只是估计一个独立成分或者一个投影追踪的方向，为了估计几个独立成分，我们需要使用具有权重向量 $\mathbf{w}$ 的若干单元（例如，神经元）来运行单单元FastICA算法。为了防止不同的向量收敛到相同的最大值，需要在每次迭代后对输出 $\mathbf{w}_{1}^{T} \mathbf{x}, \ldots, \mathbf{w}_{n}^{T} \mathbf{x}$ 去相关，有3种方法可以做到这一点，对于白化的 $\mathbf x$ 这种去相关相当于正交化。

实现去相关的简单方法是基于类似Gram-Schmidt的去相关的放缩方案。这意味着需要一个个的估计独立成分，当我们已经估计了 $p$ 个独立成分或者 $p$ 个向量 $\mathbf{w}_{1}, \dots, \mathbf{w}_{p}$，对 $\mathbf{w}_{p+1}$ 运行单元定点算法，在每一次迭代之后，从 $\mathbf{w}_{p+1}$ 中减去投影 $\mathbf{w}_{p+1}^{T} \mathbf{w}_{j} \mathbf{w}_{j}$ ，其中 $j=1, \dots, p$，最后重新标准化 $\mathbf{w}_{p+1}$：

1. $\mathbf{w}_{p+1}=\mathbf{w}_{p+1}-\sum_{j=1}^{r} \mathbf{w}_{p+1}^{\mathrm{T}} \mathbf{w}_{j} \mathbf{w}_{j}$
2. $\mathbf{w}_{p+1}=\mathbf{w}_{p+1} / \sqrt{\mathbf{w}_{p+1}^{\mathrm{T}} \mathbf{w}_{p+1}}$

然而，在某些应用程序中，可能需要使用对称的去相关，其中没有向量比其他向量具有“特权”，这可以通过矩阵平方根法做到，令：

$$\mathbf{W}=\left(\mathbf{W} \mathbf{W}^{T}\right)^{-\frac{1}{2}} \mathbf{W}$$

其中 $\mathbf{W}$ 是向量 $\left(\mathbf{w}_{1}, \dots, \mathbf{w}_{n}\right)^{T}$ 组成的矩阵，平方根的倒数 $\left(\mathbf{W} \mathbf{W}^{T}\right)^{-\frac{1}{2}}$ 可以从 $\mathbf{W} \mathbf{W}^{T}=\mathbf{F} \Lambda \mathbf{F}^{T}$ 特征值分解得到，$\left(\mathbf{W} \mathbf{W}^{T}\right)^{-\frac{1}{2}}=\mathbf{F} \Lambda^{\frac{1}{2}} \mathbf{F}^{T}$，一个更简单的迭代方法如下：

1. $\mathbf{W}=\mathbf{W} / \sqrt{\left|\mathbf{W} \mathbf{W}^{\mathrm{T}}\right|}$
2. 重复：$\mathbf{W}=\frac{3}{2} \mathbf{W}-\frac{1}{2} \mathbf{W} \mathbf{W}^{\mathrm{T}} \mathbf{W}$，直至收敛

步骤1中的范数可以使用矩阵的任意范数，比如2-范数。

FastICA 和最大似然

最后，给出FastICA和最大似然估计的联系。以矩阵形式写，看到FastICA采用以下形式：

$$\mathbf{W}^{+}=\mathbf{W}+\operatorname{diag}\left(\alpha_{i}\right)\left[\operatorname{diag}\left(\beta_{i}\right)+E\left\{g(\mathbf{y}) \mathbf{y}^{\mathrm{T}}\right\}\right] \mathbf{W}$$

其中 $\mathbf{y}=\mathbf{W} \mathbf{x}$，$\beta_{i}=-E\left\{y_{i} g\left(y_{i}\right)\right\}$，$\alpha_{i}=-\frac{1}{\beta_{i}+E\left\{g^{\prime}\left(y_{i}\right)\right\}}$，矩阵 $\mathbf{W}$ 在每一步后都需要正交化。上述版本的FastICA可以与随机梯度法的最大似然比较：

$$\mathbf{W}^{+}=\mathbf{W}+\mu\left[\mathbf{I}+g(\mathbf{y}) \mathbf{y}^{\mathrm{T}}\right] \mathbf{W}$$

其中 $\mu$ 是学习速率。$g$是独立成分的概率密度函数的函数：$g=\frac{f_{i}^{\prime}}{f_{i}}$，其中 $f_{i}$ 是独立成分的概率密度函数。FastICA可以作为ICA数据模型最大似然估计的定点算法。在FastICA中，通过选择矩阵 $\operatorname{diag}\left(\alpha_{i}\right)$ 和 $\operatorname{diag}\left(\beta_{i}\right)$ 来优化收敛速度。 FastICA的另一个优点是它可以估计子和超高斯独立分量，这与普通ML算法形成对比，普通ML算法仅适用于给定的分布类别。

FastICA的属性

与现有的ICA方法相比，FastICA算法和底层对比度函数具有许多所需的特性。

1. 在ICA数据模型的假设下，收敛是立方的（或至少是二次的）。 这与基于（随机）梯度下降方法的普通ICA算法形成对比，其中收敛仅是线性的。 这意味着非常快速的收敛，正如通过对真实数据的模拟和实验所证实的那样。
2. 与基于梯度的算法相反，没有选择步长参数。这意味着该算法易于使用。
3. 该算法使用任何非线性 $g$ 直接找到（实际上）任何非高斯分布的独立分量。 这与许多算法形成对比，其中必须首先获得概率分布函数的一些估计，并且必须相应地选择非线性。
4. 可以通过选择合适的非线性 $g$ 来优化该方法的性能。 特别地，可以获得稳健和/或最小方差的算法。 实际上，两个非线性具有一些最佳性质。
5. 可以逐个估计独立分量，这大致相当于进行投影追踪。 这在探索性数据分析中很有用，并且在仅需要估计某些独立组件的情况下减少了该方法的计算负荷。
6. FastICA具有神经算法的大部分优点：它是并行的，分布式的，计算简单的，并且需要很少的存储空间。 只有在不断变化的环境中需要快速适应性时，随机梯度法似乎才是优选的。

原论文中还涉及到了ICA算法的应用：在脑磁图去除数据伪迹中的应用；在金融数据中发现隐藏因子；为自然图像去噪。有兴趣的可以去阅读原文。

MATLAB实践

%FastICA的matlab代码
function Z = ICA( X )
 
%去均值
[M,T]=size(X);   %获取输入矩阵的行列数，行数为观测数据的数目，列数为采样点数
average=mean(X')';  %均值
for i=1:M
    X(i,:)=X(i,:)-average(i)*ones(1,T);
end
 
%白化/球化
Cx=cov(X',1);  %计算协方差矩阵Cx
[eigvector,eigvalue]=eig(Cx);   %计算Cx的特征值和特征向量
W=eigvalue^(-1/2)*eigvector';   %白化矩阵
Z=W*X;    %正交矩阵
 
%迭代
Maxcount=10000;  %最大迭代次数
Critical=0.00001;  %判断是否收敛
m=M;
W=rand(m);
for n=1:m
    WP=W(:,n);  %初始权矢量（任意）
    %Y=WP'*Z;
    %G=Y.^3;%G为非线性函数，可取y^3等
    %GG=3*Y.^2；   %G的导数
    count=0;
    LastWP=zeros(m,1);
    W(:,n)=W(:,n)/norm(W(:,n));
    while abs(WP-LastWP)&abs(WP+LastWP)>Critical
        count=count+1;  %迭代次数
        LastWP=WP;      %上次迭代的值
        %WP=1/T*Z*((LastWP'*Z).^3)'-3*LastWP;
        for i=1:m
            WP(i)=mean(Z(i,:).*(tanh((LastWP)'*Z)))-(mean(1-(tanh((LastWP))'*Z).^2)).*LastWP(i);
        end
        WPP=zeros(m,1);
        for j=1:n-1
            WPP=WPP+(WP'*W(:,j))*W(:,j);
        end
        WP=WP-WPP;
        WP=WP/(norm(WP));
        if count==Maxcount
            fprintf('未找到相应的信号');
            return;
        end
    end
    W(:,n)=WP;
end
Z=W'*Z;
end

Python实践

#coding:utf-8
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from sklearn.decomposition import FastICA, PCA
# 生成观测模拟数据
np.random.seed(0)
n_samples = 2000
time = np.linspace(0, 8, n_samples)
s1 = np.sin(2 * time)  # 信号源 1 : 正弦信号
s2 = np.sign(np.sin(3 * time))  # 信号源 2 : 方形信号
s3 = signal.sawtooth(2 * np.pi * time)  # 信号源 3: 锯齿波信号
S = np.c_[s1, s2, s3]
S += 0.2 * np.random.normal(size=S.shape)  # 增加噪音数据
S /= S.std(axis=0)  # 标准化
# 混合数据
A = np.array([[1, 1, 1], [0.5, 2, 1.0], [1.5, 1.0, 2.0]])  # 混合矩阵
X = np.dot(S, A.T)  # 生成观测信号源
# ICA模型
ica = FastICA(n_components=3)
S_ = ica.fit_transform(X)  # 重构信号
A_ = ica.mixing_  # 获得估计混合后的矩阵
# PCA模型
pca = PCA(n_components=3)
H = pca.fit_transform(X)  # 基于PCA的成分正交重构信号源
# 图形展示
plt.figure()
models = [X, S, S_, H]
names = ['Observations (mixed signal)',
         'True Sources',
         'ICA recovered signals',
         'PCA recovered signals']
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.subplots_adjust(0.09, 0.04, 0.94, 0.94, 0.26, 0.46)
plt.show()

小结

ICA的盲信号分析领域的一个强有力方法，也是求非高斯分布数据隐含因子的方法。从之前我们熟悉的样本-特征角度看，我们使用ICA的前提条件是，认为样本数据由独立非高斯分布的隐含因子产生，隐含因子个数等于特征数，我们要求的是隐含因子。

而PCA认为特征是由k个正交的特征（也可看作是隐含因子）生成的，我们要求的是数据在新特征上的投影。同是因子分析，一个用来更适合用来还原信号（因为信号比较有规律，经常不是高斯分布的），一个更适合用来降维（用那么多特征干嘛，k个正交的即可）。有时候也需要组合两者一起使用。这段是我的个人理解，仅供参考。