ESL-ICA

Independent Component Analysis

假设x∈ℝp×1是观测直，s∈ℝp×1是源，则有:

x=As     (1)

这里提到一个概念，自由（Independent）与不相关的关系(uncorrelated)：

p(x,y) = p(x)p(y) => E[x,y] = E[x]E[y]

前者是自由的充要条件，后者是uncorrelated的条件。
ICA与PCA的区别在于，对于ICA，s是independent而对于PCA则仅要求s是uncorrelated. 那么ICA要求s是non-Gaussian Distribution，因为虽然Gaussian Dist是uncorrelated的但不满足independent的条件。

ICA & Maximum Likelihood

欲求s，我们需要求解A或者W=A−1，假设s的概率分布是ps(s)，则x的分布是：

px(x)=ps(Ws)|W|     (2)x=[x1,x2,...,xp]Ts=[s1,s2,...,sp]T

记:

W=⎡⎣⎢⎢⎢⎢−wT1−−wT2−......−wTp−⎤⎦⎥⎥⎥⎥   s=Wx

则有：

si=wTi×xipx(x)=ps(wTixi)i=1,2,...p

Sigmoid function

我们并不确定ps(s)的表达形式，可以提出先验的pdf，当然不能用Gaussian pdf, 而是由sigmoid function变形而来：

∫ps(s)=sig(s)=11+e−sps(s)=sig′(s)sig′(s)=sig(s)(1−sig(s))

记data set, X∈ℝn×p 我们列出有关W的，ML的推导过程如下：

L(w)=∏i=1npx(x)=∏i=1n∏j=1ppx(xj)=∏i=1n∏j=1pps(wTix)|W|=∏i=1n∏j=1psig′(wTix)|W|l(w)=∑i=1n(∑j=1pln(sig′(wTix))+ln|W|)

∂l(w)∂wj=∑i=1n[∑j=1p(sig(1−sig))′sig′+|W|−1|W|(W−1)T]=∂l(w)∂wj=∑i=1n[∑j=1p(1−2sig(wTix))+(W−1)T]

SGD迭代的方式为：

For a training examplex(i)∈ℝp Wnew=Wold−α(⎡⎣⎢⎢⎢⎢⎢1−2sig(wT1x(i))1−2sig(wT2x(i))......1−2sig(wTpx(i))⎤⎦⎥⎥⎥⎥⎥(x(i))T+(WT)−1)α>0

FastICA

算法流程如下：

Whiten

白化的定义：
设x是一个随机变量，存在一个线性变换V将它变换成z： z=Vx且E[zzT]=I，V就是白化矩阵。
x协方差covariance是

Cx=E[xxT]=PDP−1

P单位化，则

Cx=E[xxT]=PDPT

那么V=D−1/2PT.

Python Implement

#!/usr/bin/env python#FastICA from ICA book, table 8.4 import mathimport randomimport matplotlib.pyplot as pltfrom numpy import *%matpltlibn_components = 2def 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)

返回的是g(x)和E(g’(x))

---

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

这里用到一些技巧：

WWT=VD2VT

D2是WWT的特征值。

[WWT]−1=VD−2VT[WWT]−12=VD−12VT

---

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))        W = dot(gwtx, X.T) / p - g_wtx[:, newaxis] * W

wi=E[xig(wTixi)]−E[g′(wTIxi)]wi

        W1 = do_decorrelation(W)

W1=(WWT)−1/2W

        lim = max( abs(abs(diag(dot(W1, W.T))) - 1) )        W = W1        if lim < 0.0001:            break    return Wdef 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()