跳慮跋考

興味も思考も行先不明

Reading "Topographic Independent Component Analysis"

感覚信号の効率的な(即ちスパースな)符号化に関して「視覚野・聴覚野地図の同一適応アルゴリズムによる解釈」に興味を持ったので、論文で使用されているトポグラフィック独立成分分析(TICA)の論文「Topographic Independent Component Analysis」を読もうとしたものの正直よく解らなかった為、Hyvarinen 先生の著「詳解 独立成分分析―信号解析の新しい世界」で勉強してからリベンジした次第。以下に論文の要諦を纏めてみる。

いい加減 blog というメディアに適合していない気がするが……。

モデル

共通の分散

通常の ICA では信号源たる確率変数 si が互いに独立であることを仮定するが、現実には全く独立である場合と言うのは少なくて ICA を行っても幾許かの従属性が「残存」してしまうので、それを使って独立成分の二次元構造を決定しようというのが TICA の趣旨である。
特にこの論文ではエネルギー、即ち si2 の相関に着目し、その相関が強いほど近い位置にあるとする。これは近い位置の成分ほどより共起しやすい(が値自体一方から一方を予測しがたい)と考えている事となる。
例えば変数 σ により成分 si, sj (i ≠ j) に共通な分散が与えられるとし、zi, zj を白色*1な確率変数として
 \begin{cases} s_i = z_i \sigma \\ s_j = z_j \sigma \end{cases}
とおく。ここで σ や zi, zj は平均 0 で互いに独立とする。
この時、si と sj の値は互いに無相関だがエネルギーはそうではない。実際、独立性と白色性の仮定から
 E\{s_is_j\} = E\{ z_i \sigma \cdot z_j \sigma \} = E\{ z_i z_j \} E\{ \sigma^2 \} = 0
 E\{ s_i^2 s_j^2 \} - E\{ s_i^2 \} E\{ s_j^2 \} = E\{ z_i^2 z_j^2 \} E\{ \sigma^4 \} - E\{z_i^2\} E\{z_j^2\} E\{ \sigma^2 \}^2
  = E\{ \sigma^4 \} - E\{ \sigma^2 \}^2
となり、最後の式は σ2 の分散なので非負である。

成分の空間構造

si の分散 σi2 を確率変数とし、各 si はその分散が与えられた下で互いに独立とする。
σi2 の従属性は近傍函数 h(i,j) により記述する。通常近傍函数はある種の距離について単調減少する函数と定義され、対称性 h(i,j) = h(j,i) や、任意の i について h(i,i) = const. が要求される。
例えば、i 番目の成分の座標が (xi, xj) であるとすると
 h(i,j) = \begin{cases} 1 & (|x_i - x_j| \le 1\ \textrm{or}\ |y_i - y_j| \le 1) \\ 0 & (\textrm{otherwise}) \end{cases}
とすれば自身と上下左右及び斜めの9近傍が取れる。
これにより、非線形函数 φ を*2用いて
 \sigma_i = \phi\left( \sum_{k=1}^{n} h(i,k) u_k \right)
で分散を定義する。(正直この φ が何を表すのかよく分からないのだが……)
以上の結果として TICA のモデルは、
 s_i = \phi\left( \sum_{k=1}^{n} h(i,k) u_k \right) z_i
と書くことができる。

尤度函数

導出

まずトポグラフィックな成分 s と「高次」の独立成分 u の結合密度分布
 p(\textbf{s}, \textbf{u}) = \prod_i p(s_i | u_i) p(u_i)
  = \prod_i \frac{1}{ \phi\left( \sum_{k=1}^{n} h(i,k) u_k \right) } p^{s|u}_i \left( \frac{s_i}{ \phi\left( \sum_{k=1}^{n} h(i,k) u_k \right) } \right) \cdot p^u_i(u_i)
から(ps|ui は pzi なので確率密度の変換をしている)、s の周辺密度分布が
 p_s(\textbf{s}) = \int \prod_i p^{s|u}_i\left( \frac{s_i}{ \phi\left( \sum_{k=1}^{n} h(i,k) u_k \right) } \right) \frac{ p^u_i(u_i) }{ \phi\left( \sum_{k=1}^{n} h(i,k) u_k \right) } d\textbf{u}
で得られる。W = (w1, ..., wn)T として、 ICA の基本モデル s = Wx より
 p_x\left( \textbf{x}(t) \right) = p_s\left( W \textbf{x}(t) \right) |\det W|
となるから、尤度函数
 L(W) = \prod_{t=1}^{T} p_x\left( \textbf{x}(t) \right)
  = \prod_{t=1}^{T} \int \prod_i p^{s|u}_i\left( \frac{ \textbf{w}_i^T \textbf{x}(t) }{ \phi\left( \sum_{k=1}^{n} h(i,k) u_k \right) } \right) \frac{ p^u_i(u_i) }{ \phi\left( \sum_{k=1}^{n} h(i,k) u_k \right)} |\det W| d\textbf{u}
となる。

近似

上の尤度函数をそのまま用いるのは計算量が多く賢明でないので、扱い易い近似を求める。
まず ps|ui や pui は全ての i について等しいとし、ps|u = ps|uiガウス分布とする。これは条件付きでない s の分布が優ガウス*3である事を意味する。また非線形函数 φ を
 \phi\left( \sum_{k=1}^{n} h(i,k) u_k \right) = \left( \sum_{k=1}^{n} h(i,k) u_k \right)^{-1/2}
で定める。*4
これらの近似より s の周辺密度分布は
 p_s(\textbf{s}) = \int \frac{1}{(2\pi)^{n/2}} \exp\left[ -\frac{1}{2} \sum_i s_i^2 \left( \sum_k h(i,k) u_k \right) \right] \prod_i p_u(u_i) \sqrt{ \sum_k h(i,k) u_k } d\textbf{u}
となり、和の順序を交換すると
 p_s(\textbf{s}) = \int \frac{1}{(2\pi)^{n/2}} \exp\left[ -\frac{1}{2} \sum_k u_k \left( \sum_i h(i,k) s_i^2 \right) \right] \prod_i p_u(u_i) \sqrt{ \sum_k h(i,k) u_k } d\textbf{u}
と書ける。
ここで更に
 \sqrt{ \sum_k h(i,k) u_k } \approx \sqrt{ h(i,i) u_i }
なる近似によって(実際にはこれは下界になっている)
 \widetilde{p}_s(\textbf{s}) = \int \frac{1}{(2\pi)^{n/2}} \exp\left[ -\frac{1}{2} \sum_k u_k \left( \sum_i h(i,k) s_i^2 \right) \right] \prod_i p_u(u_i) \sqrt{ h(i,i) u_i } d\textbf{u}
となる。h(i,i) は定数であったので h0 とおき、この積分を各 uk について分解すると
 \widetilde{p}_s(\textbf{s}) = \prod_{k=1}^{n} \int \frac{1}{\sqrt{2\pi}} \exp\left[ -\frac{1}{2} u_k \left( \sum_i h(i,k) s_i^2 \right) \right] p_u(u_k) \sqrt{ h_0 u_k } du_k
だが、uk についての対称性から
 G(y) = \log \int \frac{1}{\sqrt{2\pi}} \exp\left( -\frac{1}{2} uy \right) p_u(u) \sqrt{ h_0 u } \,du
という函数を定義すれば
 \widetilde{p}_s(\textbf{s}) = \prod_{k=1}^{n} \exp\left[ G\left( \sum_i h(i, k) s_i^2 \right) \right]
と簡潔な形になる。

すると、対数尤度函数
 \log \widetilde{L}(W) = \sum_{t=1}^{T} \sum_{j=1}^{n} G\left( \sum_{i=1}^{n} h(i, j) (\textbf{w}_i^T \textbf{x}(t))^2 \right) + T \log|\det W|
により近似できる。*5
ICA に於いて独立成分の確率分布(の対数)にあたる函数 G は、概形さえ合っていれば正しく収束する事が知られており、著者は指数分布に関連した函数として
 G_1(y) = -\alpha_1 \sqrt{y + \epsilon} + \beta_1
を挙げている。ε は数値的な安定性の為に加えられる定数である。係数は si の分散が 1 であるという条件を満たす様に決められる定数だが、以下で W を直交行列とする制約を課す為にこれらは対数尤度を線形変換するだけなので、α1=1, β1=0 などと適当に決めても尤度の最大点には影響しない。

学習則

さて尤度が得られたので、勾配法によりこれを最大化しよう。まず両辺を T で割り
 \frac{1}{T} \log \widetilde{L}(W) = E\left\{ \sum_{j=1}^{n} G\left( \sum_{i=1}^{n} h(i, j) (\textbf{w}_i^T \textbf{x})^2 \right) \right\} + \log|\det W|
と書き直す。(右辺第1項は標本平均なので、厳密にはこう書くべきでない。)すると wk に関する勾配は
 \frac{\partial}{\partial w_{kl}} G\left( \sum_{i=1}^{n} h(i, j) (\textbf{w}_i^T \textbf{x})^2 \right) = G'\left( \sum_{i=1}^{n} h(i, j) (\textbf{w}_i^T \textbf{x})^2 \right) \cdot 2 h(k, j) (\textbf{w}_k^T \textbf{x}) x_l
となるから、ベクトル勾配に直すと
 \frac{\partial}{\partial \textbf{w}_k} G\left( \sum_{i=1}^{n} h(i, j) (\textbf{w}_i^T \textbf{x})^2 \right) = 2 G'\left( \sum_{i=1}^{n} h(i, j) (\textbf{w}_i^T \textbf{x})^2 \right) h(k, j) \left( \textbf{w}_k^T \textbf{x}(t) \right) \textbf{x}
となる。
観測値は z = Vx により白色化されているとし、W は直交行列とする。この時 det W は定数なので、結局
 \frac{\partial}{\partial \textbf{w}_k} \frac{1}{T} \log \widetilde{L}(W) = 2 E\left\{ r_k \left( \textbf{w}_k^T \textbf{z} \right) \textbf{z} \right\}
となる。ここで
 r_k = \sum_{j=1}^{n} h(k, j) g\left( \sum_{i=1}^{n} h(i, j) (\textbf{w}_i^T \textbf{z})^2 \right)
であり、g は G の導関数である。また W は仮定を満たす為に
 W \leftarrow (WW^T)^{-1/2} W
で毎回正規直交化されなければならない。

毎回標本平均を取るのは実用的でない場面が多いので、確率的勾配降下法に基づくオンライン学習則として
 \Delta \textbf{w}_k \propto r_k \left( \textbf{w}_k^T \textbf{z} \right) \textbf{z}
を用いる。これは単に期待値計算を外しただけだが、学習率などに適当な条件を課せばちゃんと収束する事が証明されている。

以上により取り敢えず学習則が導出できたので、一旦筆を擱くとする。やっと原理が解ってきたので愈々実装と行きたいところ。
ただ φ はどうも納得がいかない……。

*1:共分散行列が単位行列に等しい。

*2:不幸にも環境によっては \phi ではなく \varphi が表示されているかもしれないが大目に見て欲しい。

*3:正の尖度を持つ。ガウス分布(正規分布)よりも平均近くの密度が高い。

*4:この函数が単調「減少」である事はどうにも不思議なのだが、私の「u_i が大きいほど活性が大きい」という認識が間違っているのだろうか?

*5:尤度が近傍のエネルギーの和のみの函数でよいというのは Kohonen の不変特徴部分空間の「部分空間の射影のノルムが不変特徴を定める」という話と関連があるらしい。