独立成分分析の実装

はじめに

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

独立成分分析

確率変数 \textbf{s} = {s_1, s_2, ..., s_N} が線形に"mix"された観測 \bf{x} が得られたとき，元の観測をどうにかして復元したい． 独立成分分析では，元の確率変数 \textbf{s}の要素が独立であるという仮定から，復元を試みる．

"mix"行列を Aとすると，

\textbf{x} = A\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}を決める．

中心極限定理と独立

ここで \textbf{y} = [y_1, y_2]とすると，

y_1 = \textbf{w}^\top A \textbf{s} ここで， \textbf{w}は Wの行ベクトル

ここで， \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)がある． \textbf{y} の分散共分散行列が単位行列のとき， \mathrm{kurt} ( \textbf{y} ) = \mathbb{ E } [ y^4 ]-3

尖度はガウス分布に近いときに0となる． よって，尖度の絶対値を最大化するように 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では不動点反復法(?)で最適化するらしい． 制約付きの最適化問題なので，リーマン多様体上の最適化としても解ける気がする．

参考