ICA独立成分分析

ICA的数学推导

ICA算法的思路比较简单，但是推导过程比较复杂，本文只是梳理了推理路线。

假设我们有n个混合信号源$X\subset{R^{n}}$和n个独立信号$S\subset{R^{n}}$，且每个混合信号可以由n个独立信号的线性组合产生，即：$X=\left[ \begin{matrix} x_1&\\ x_2&\\ …&\\ x_n& \end{matrix} \right]S=\left[ \begin{matrix} s_1&\\ s_2&\\ …&\\ s_n& \end{matrix} \right]X=AS => S=WX,W=A^{-1}$

假设我们现在对于每个混合信号，可以取得m个样本，则有如下n*m的样本矩阵：$D=\left[ \begin{matrix} d_{11}&d_{12}&…&d_{1m}&\\ …&\\ d_{n1}&d_{n2}&…&d_{nm}\\ \end{matrix} \right]$

由于S中的n个独立信号是相互独立的，则它们的联合概率密度为：$$p_S(s)=\Pi_{i=1}^{n}p_{s_i}(s_i)$$

由于$s=Wx$，因此我们可以得出：$p_X(x)=F^{‘}_{X}(x)=|\frac{\partial{s}}{\partial{x}}|*p_S(s(x))=|W|*\Pi_{i=1}^{n}p_{s_i}(w_ix)$

考虑目前有m个样本，则可以得到所有样本的似然函数：$L=\Pi_{i=1}^{m}(|W|*\Pi_{j=1}^{n}p_{s_j}(w_{j\cdot}d_{\cdot{i}}))$

取对数之后，得到：$lnL=\Sigma_{i=1}^{m}\Sigma_{j=1}^{n}lnp_{s_j}(w_{j\cdot}d_{\cdot{i}})+mln|W|$

之后只要通过梯度下降法对lnL求出最大值即可，即求使得该样本出现概率最大的参数W。
此时假设我们上面的各个独立信号的概率分布函数为sigmoid函数，但是不确定这里的g函数和下面fastICA中的g函数是否有关联）：$F_{s_i}(s_i)=\frac{1}{1+e^{-s_i}}$

最终，我们求得：$\frac{\partial{lnL}}{\partial{W}}=Z^TD+\frac{m}{|W|}(W^*)^T$

其中：$$Z=g(K)=\left[ \begin{matrix} g(k_{11})&g(k_{12})&…&g(k_{1m})&\\ …&\\ g(k_{n1})&g(k_{n2})&…&g(k_{nm})\\ \end{matrix} \right]g(x)=\frac{1-e^x}{1+e^x}K=WDD=\left[ \begin{matrix} d_{11}&d_{12}&…&d_{1m}&\\ …&\\ d_{n1}&d_{n2}&…&d_{nm}\\ \end{matrix} \right]$$

由于伴随矩阵具有以下性质：$$WW^*=|W|I$$

因此我们可以求出：$\frac{\partial{lnL}}{\partial{W}}=Z^TD+m(W^{-1})^T$

因此可以得到梯度下降更新公式：$$W=W+\alpha(Z^TD+m(W^{-1})^T)$$

至此，ICA的基本推理就此结束。下面我们来看一下fastICA的算法过程（没有数学推理）。

fastICA的算法步骤

观测信号构成一个混合矩阵，通过数学算法进行对混合矩阵A的逆进行近似求解分为三个步骤：

去均值。去均值也就是中心化，实质是使信号X均值是零。
白化。白化就是去相关性。
构建正交系统。

在常用的ICA算法基础上已经有了一些改进，形成了fastICA算法。fastICA实际上是一种寻找$w^Tz(Y=w^Tz)$的非高斯最大的不动点迭代方案。具体步骤如下：

观测数据的中心化(去均值)
数据白化(去相关)，得到z
选择需要顾及的独立源的个数n
随机选择初始权重W（非奇异矩阵）
选择非线性函数g
迭代更新：
- $w_i \leftarrow E\{zg(w_i^Tz)\}-E\{g^{‘}(w_i^Tz)\}w$
- $W \leftarrow (WW^T)^{-1/2}W$
判断收敛，是下一步，否则返回步骤6
返回近似混合矩阵的逆矩阵

代码实现
基于python2.7，matplotlib，numpy实现ICA，主要参考sklean的FastICA实现。

import math
import random
import matplotlib.pyplot as plt
from numpy import *

n_components = 2

def f1(x, period = 4):
return 0.5(x-math.floor(x/period)period)

def create_data():
#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]], float32)
#mix matrix
A = array([[0.8, 0.2], [-0.3, -0.7]], float32)
return T, S, dot(A, S)

def whiten(X):
#zero mean
X_mean = X.mean(axis=-1)
X -= X_mean[:, newaxis]
#whiten
A = dot(X, X.transpose())
D , E = linalg.eig(A)
D2 = linalg.inv(array([[D[0], 0.0], [0.0, D[1]]], float32))
D2[0,0] = sqrt(D2[0,0]); D2[1,1] = sqrt(D2[1,1])
V = dot(D2, E.transpose())
return dot(V, X), V

def _logcosh(x, fun_args=None, alpha = 1):
gx = tanh(alpha * x, x); g_x = gx ** 2; g_x -= 1.; g_x *= -alpha
return gx, g_x.mean(axis=-1)

def do_decorrelation(W):
#black magic
s, u = linalg.eigh(dot(W, W.T))
return dot(dot(u * (1. / sqrt(s)), u.T), W)

def do_fastica(X):
n, m = X.shape; p = float(m); g = _logcosh
#black magic
X *= sqrt(X.shape[1])
#create w
W = ones((n,n), float32)
for i in range(n):
for j in range(i):
W[i,j] = random.random()

#compute W
maxIter = 200
for ii in range(maxIter):
    gwtx, g_wtx = g(dot(W, X))
    W1 = do_decorrelation(dot(gwtx, X.T) / p - g_wtx[:, newaxis] * W)
    lim = max( abs(abs(diag(dot(W1, W.T))) - 1) )
    W = W1
    if lim < 0.0001:
        break
return W

def show_data(T, S):
plt.plot(T, [S[0,i] for i in range(S.shape[1])], marker=”*”)
plt.plot(T, [S[1,i] for i in range(S.shape[1])], marker=”o”)
plt.show()

def main():
T, S, D = create_data()
Dwhiten, K = whiten(D)
W = do_fastica(Dwhiten)
#Sr: reconstructed source
Sr = dot(dot(W, K), D)
show_data(T, D)
show_data(T, S)
show_data(T, Sr)

if name == “main“:
main()

参考：

http://skyhigh233.com/blog/2017/04/01/ica-math/

ICA的数学推导

fastICA的算法步骤

相关文章：

发表评论 取消回复

发表评论取消回复