独立成分分析の実装

はじめに

統計的因果探索を読んでいる． この中で，独立成分分析が出てくるが，そういえばあんまり理解できていなかったので，実装して理解する． 日本語の本も見たけど，結局 Independent Component Analysis: Algorithms and Applicationsがわかりやすかった．

独立成分分析

独立成分分析では，観測 \textbf{x} = A\textbf{s}が与えられたとき，元の確率変数 \textbf{s} = {s_1, s_2, ..., s_N}を復元することを試みる． ここで，独立成分分析では，確率変数 \textbf{s}の要素が独立であるという仮定をおいて，信号を復元する．

以下では，矩形波と制限波からなる \textbf{s}を例にプログラムを書く．

import pylab as plt
%matplotlib inlinee

import numpy as np
from scipy import signal

from sklearn.decomposition import PCA

import torch
from torch.autograd import Variable

np.random.seed(0)

n_samples = 5000
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

S = np.c_[s1,  s2]

S /= S.std(axis=0)  # Standardize data
A = np.array([[1, 1], 
              [0.5, 2],
             ])  # Mixing matrix

X = np.dot(S, A.T)  # Generate observations

元の信号は以下．

plt.plot(S)
plt.title("S")
plt.show()

f:id:ksknw:20180519231649p:plain

この信号に対して，以下のような観測が得られたとする．

plt.plot(X)
plt.title("X")
plt.show()

f:id:ksknw:20180519231705p:plain

独立成分分析では， \textbf{s}の各要素が独立であるという仮定から，復元行列 Wをいい感じに求めて，復元値 \textbf{y} = W\textbf{x}を決める．

中心極限定理と独立

2次元の系列データ \textbf{y} = [y_1, y_2]を復元するとする． このとき， [tex: y_1 = \textbf{w}1^\top A \textbf{s}] ここで， \textbf{w}1は Wの行ベクトル(1行目)である． ここで， \textbf{z}=A\textbf{w} とすると y_1 = \textbf{z}^\top \textbf{s}であり，復元値 y_1, y_2 は独立な確率変数 \textbf{s}を重み付きで足したものであるといえる．

中心極限定理から独立な確率変数の和は，ガウス分布に近づくことが多い(参考文献では"usually")． このため，復元がうまく行っていない場合には， \textbf{y}は \textbf{s}に比べてガウス分布に近づくことが多い． 独立成分分析では，この性質を利用して \textbf{y}がよりガウス分布から離れるように Wを求める．

ガウス分布度合いを測る指標として尖度(kurtosis)がある．尖度はガウス分布に近いときに0となる． \textbf{y} の分散共分散行列が単位行列のとき，尖度は \mathrm{kurt} ( \textbf{y} ) = \mathbb{ E } [ y^4 ]-3 として求まる． 独立成分分析では，尖度の絶対値を最大化するような Wを勾配法によって求める． ただし， Wの大きさを大きくすると，どこまでの誤差関数を小さくできてしまうので， Wのノルムを1に制約する． 今回は適当に毎ステップでノルムが１になるようにリスケールすることで，制約を満たす解が求まるようにした．

pca = PCA(n_components = X.shape[1], whiten=True)
X = pca.fit_transform(X)

def kurt(Y):
    return torch.mean(Y**4, dim=0) - 3
def kurt_loss(Y):
    # Y の分散共分散が単位行列であることが前提
    return -torch.mean(torch.abs(kurt(Y)))

def ICA(X, loss_func = kurt_loss):
    X = Variable(torch.from_numpy(X), requires_grad =False).float()
    W = Variable((torch.randn((2, 2))), requires_grad =True)
    W.data = W.data / torch.norm(W.data)
    hist_loss = []
    
    for i in range(1000):
        Y = torch.mm(X, W)
        loss = loss_func(Y)
        hist_loss.append(loss.data.numpy())
        loss.backward()

        W.data = W.data - W.grad.data * 0.1
        W.grad.data.zero_()
        W.data = W.data / torch.norm(W.data)
    return Y, W, hist_loss

Y, W, hist_loss = ICA(X, loss_func=kurt_loss)
plt.plot(Y.data.numpy())
plt.show()
plt.plot(hist_loss)

f:id:ksknw:20180519231733p:plain

f:id:ksknw:20180519231743p:plain

以上より，いい感じに復元することができた．

Negentropyを用いた誤差関数

上ではガウス分布度合いを測る指標として，尖度を用いた． しかし，尖度は外れ値に対して敏感であることが知られている． そこで、ガウス分布度合いを測る指標として，Negentropyを用いる方法が提案されている．

J(\textbf{y}) = H(\textbf{y}_\mathrm{gauss}) - H(\textbf{y})

ここで Hはエントロピー． ガウス分布はエントロピーが最も高い確率分布であるため，その差を使ってガウス分布度合いを評価することができるということ． しかし，結局，この計算は確率密度関数が必要なので，データ点だけから求めるのは一般には難しい． そこで，以下のように近似的にNegentropyを求める．

J(\textbf{y}) \approx \sum _{i=1}^p k_i { \mathbb{E} [G_i(y) ] - \mathbb{E} [G_i(\nu) ] }^2 期待値をサンプルから近似することで，この式は計算できる． ただし \nu \sim \mathcal{N}(0,\mathbb{I}), k_i>0 ． また， Gには色々な関数が設定できるが，例えば， G_1(u) = 1/a_1 \log \mathrm{cosh} a_1 uがある．

def negentropy_loss_logcosh(Y, a=1):
    return -torch.mean(torch.log(torch.cosh(a * Y))/a )**2 

Y, W, hist_loss = ICA(X, loss_func=negentropy_loss_logcosh)
plt.plot(Y.data.numpy())
plt.show()
plt.plot(hist_loss)

f:id:ksknw:20180519231759p:plain

f:id:ksknw:20180519231815p:plain

おわりに

独立成分分析の理解のために，自分で実装してみた． 今回は，天下り的に独立とガウス分布らしさの関係について考えたが，論文では同じ最適化問題が相互情報量を最小化するという考えからも導出されている． また，今回は雑に最適化問題を解いたが，FastICAでは不動点反復法(?)で最適化するらしい． 制約付きの最適化問題なので，リーマン多様体上の最適化としても解ける気がする．

参考