📂統計的分析

リザーバーコンピューティング

定義 ¹

時系列 予測の文脈で、 $t$ 時点で $M$ 次元の入力ベクトル $\mathbf{u} (t) \in \mathbb{R}^{M}$ が与えられたとき $P$ 次元の出力ベクトル $\mathbf{s} (t) \in \mathbb{R}^{P}$ を予測しようとする。これに対して リザバー動力学^{reservoir dynamics} $\mathbf{r} (t) \in \mathbb{R}^{N}$ を導入し、その線形変換として未来の出力ベクトルを予測する方法を リザバーコンピューティング^{reservoir computing}と呼ぶ。

$0 < t \leq T$ までのデータが一定の時間間隔 $\Delta t$ で分割されて $K$ 個のデータポイント $\left\{ \left( \mathbf{u} \left( k \Delta t \right) , \mathbf{s} \left( k \Delta t \right) \right) \right\}_{k=1}^{K}$ が与えられているとしよう。$T$ 以降の予測出力ベクトル $\hat{ \mathbf{s} }$ は次のように計算される。 $$ \hat{ \mathbf{s} } (t + \Delta t) = W_{\text{out}} \left[ \left( 1 - \alpha \right) \mathbf{r} (t) + \alpha \tanh \left( A \mathbf{r} (t) + W_{\text{in}} \mathbf{u} (t) \right) \right] $$ ここで $W_{\text{out}}$ を除いてはハイパーパラメータとして次の通りである。

• $\alpha \in [0, 1]$: リーク率^{leakage rate}
• $A \in \mathbb{R}^{N \times N}$: リザバーレイヤの隣接行列
• $W_{\text{in}} \in \mathbb{R}^{N \times M}$: 入力の重み行列
• $\tanh : \mathbb{R}^{N} \to \mathbb{R}^{N}$: ベクトル化されたハイパボリックタンジェント

説明

リザバーコンピューティングは力学系の研究などで非常によく用いられる機械学習手法で、過去のデータは与えられているが今後の観測が制限される状況で有用だ。

$\hat{ \mathbf{s} }$ の式の形によれば $t$ 時点の入力ベクトル $\mathbf{u} (t)$ が与えられれば、既に学習によって得られた $W_{\text{out}}$ を用いて $t + \delta t$ 時点の出力ベクトル $\hat{ \mathbf{s} } (t + \Delta t)$ を予測できる。リザバーコンピューティングの主要目的は次のように $\mathbf{s}$ を $\mathbf{u}$ として説明する関数の近似 $f$ を見つけることだと要約して差し支えない。 $$ \mathbf{s} \left( t + \Delta t \right) = f \left( \mathbf{u} (t) \right) $$

ロスラーアトラクタ: $$ \begin{align*} {{dx} \over {dt}} =& - y - z \\ {{dy} \over {dt}} =& x + ay \\ {{dz} \over {dt}} =& b + (x-c) z \end{align*} $$

上の図はロスラーアトラクタを学習したリザバーコンピュータが、ただ $x$ 成分のみを観測しても $y, z$ 成分を正しく予測する様子を示している。この場合 $\mathbf{u} = (x)$ は $M=1$ 次元のベクトル、$\mathbf{s} = (y,z)$ は $P=2$ 次元のベクトルになる。

実装

リザバー動力学系

まずリザバーレイヤのノード数 $N \in \mathbb{N}$ が与えられる必要がある。$\mathbf{r} \in \mathbb{R}^{N}$ は次のように写像の形をした力学系を成し、一種の潜在空間^{latent space}の役割を果たすと考えられる。 $$ \mathbf{r} \left( t + \Delta t \right) = \left( 1 - \alpha \right) \mathbf{r} (t) + \alpha \tanh \left( A \mathbf{r} (t) + W_{\text{in}} \mathbf{u} (t) + \xi \mathbf{1} \right) $$ 初期状態 $\mathbf{r} (0)$ はランダムに決めればよく、$\alpha \in [0, 1]$ は入力ベクトル $\mathbf{u}$ の影響をどれだけ反映するかに応じて調整すればよい。$\tanh$ は最小二乗法を使う前に式に非線形性^nonlinearityを与えるための役割を果たし、厳密に同じではないが機能的には活性化関数に類似すると見なせ、必ずしもハイパボリックタンジェントである必要はない。$\xi$ はバイアス^biasの大きさ、$\mathbf{1}$ は一ベクトルだが行列代数的にこのバイアス項は大きな意味を持たないため、冒頭に記した式では省略した。

リザバーレイヤ

$N$ 次元の $\mathbf{r}$ が次のタイムステップへ移る際に再帰的に発生する相互作用は リザバーネットワーク^{reservoir network} によって制御される。$A$ は エルデシュ–レーニィ ネットワーク の隣接行列であり、その ランダムネットワーク の 平均次数 が $D$ になるように サンプリング される。$A$ のうち $0$ でない成分は 一様分布 $U[-1, 1]$ から iid に値を抽出して置換した後、$A$ の 固有値 のうち大きさが最大の値、すなわち スペクトル半径 が $\rho > 0$ になるように $A$ にある スカラー を掛けて補正する。

このようにランダムネットワークでリザバーコンピュータの構造を決める方式を エコステートネットワーク^{echo state network} と呼ぶ²。ご覧の通り $\mathbf{r}$ はネットワークがそのまま $\mathbf{r}$ 自身に影響を与え合い、それ自体として 再帰型ニューラルネットワーク^RNN のアイデアを反映する。

入力重み行列

$W_{\text{in}} \in \mathbb{R}^{N \times M}$ は入力ベクトル $\mathbf{u} (t) \in \mathbb{R}^{M}$ が持つ $M$ 個の成分のうち1つだけを受け渡す役割を持つ。直感的に説明すると、$W_{\text{in}}$ の各行は1つの成分のみが $0$ であってはならず、逆に $N$ 個のリザバーネットワークの各ノードはそれぞれちょうど1つの入力成分だけを受け取る。どの入力成分を受け取るかはランダムに決め、その値も一様分布 $U[-\sigma, \sigma]$ からランダムに抽出する。入力重みスケール^{scale of input weights} $\sigma > 0$ はハイパーパラメータであり、見ての通り $W_{\text{in}}$ は合理的な理由があればかなり恣意的に決めてもよい。

出力重み行列

ここまでで行ったのは $\mathbf{r}$ を $\mathbf{u}$ に対する関数として表現することだ。結果として望むのは $\mathbf{s}$ が $\mathbf{r}$ に対する関数として表現され、$\mathbf{u}$ が入力されると $\mathbf{s}$ が出力されることだ。$\mathbf{s}$ の予測ベクトル $\hat{ \mathbf{s} }$ は次のように表される。 $$ \hat{\mathbf{s}} (t) = W_{\text{out}} \mathbf{r} (t) + \mathbf{c} $$ ここで $\mathbf{c}$ は先に見た $\xi \mathbf{1}$ と同様にバイアス項だが、これも大きな意味を持たないため要約した式では省略した。$W_{\text{out}}$ と $\mathbf{c}$ はハイパーパラメータではなく、実際に計算してすなわち学習を通じて得るべきパラメータである。

まず固定されたハイパーパラメータに対して $\mathbf{u} (t)$ により $\mathbf{r} (t)$ をすべて計算した後、次のような行列 $R, S$ を定義する。 $$ \begin{align*} \overline{\mathbf{r}} =& {\frac{ 1 }{ K }} \sum_{k=1}^{K} \mathbf{r} (k \Delta t) \\ \overline{\mathbf{s}} =& {\frac{ 1 }{ K }} \sum_{k=1}^{K} \mathbf{s} (k \Delta t) \\ \left( R \right)_{:, k} =& \mathbf{r} (k \Delta t) - \overline{\mathbf{r}} \\ \left( S \right)_{:, k} =& \mathbf{s} (k \Delta t) - \overline{\mathbf{s}} \end{align*} $$ ここで $\left( X \right)_{:, j}$ は行列の $j$ 番目の列を表す表記である。言い換えれば $R$ と $S$ はデザイン行列だ。これで最小二乗法を使うための準備が整った。

リッジ回帰の最適解: $$ L \left( \beta \right) = \left\| Y - X \beta \right\|_{2}^{2} + \lambda \left\| \beta \right\|_{2}^{2} $$ $\lambda$ が定数として与えられているとき、リッジ回帰の目的関数 $L$ を上のように表す。リッジ回帰の最適解 $\hat{\beta} = \argmin_{\beta} L \left( \beta \right)$ は次のとおりである。 $$ \hat{\beta} = \left( X^{T} X + \lambda I \right)^{-1} X^{T} Y $$ ここで $A^{T}$ は $A$ の転置行列、$I$ は単位行列、$A^{-1}$ は $A$ の逆行列である。

過学習を防ぐための正則化としてリッジ回帰を用いるとすると、スパース性パラメータ $\lambda \geq 0$ に対して $W_{\text{out}}^{\ast}$ は次のように計算される。 $$ W_{\text{out}}^{\ast} = S R^{T} \left( R R^{T} + \lambda I \right)^{-1} $$ 引用した式と乗算の順序が逆になっているのは、デザイン行列が学習データの量に応じて横に長くなるよう定義されており、すでに転置が取られているからだと考えられる。$\mathbf{c}^{\ast}$ は次のように平均的な意味でのドリフトのような感覚で求められる。 $$ \mathbf{c}^{\ast} = \overline{\mathbf{s}} - W_{\text{out}}^{\ast} \overline{\mathbf{r}} $$

結果的に見るとただの多変量非線形回帰分析ではないかと思うかもしれないが、実際その通りだ。省略せずに実装どおりの式は次のとおりである。 $$ \begin{align*} \hat{ \mathbf{s} } (t + \Delta t) =& W_{\text{out}}^{\ast} \mathbf{r} (t + \Delta t) + \mathbf{c}^{\ast} \\ =& W_{\text{out}}^{\ast} \left[ \left( 1 - \alpha \right) \mathbf{r} (t) + \alpha \tanh \left( A \mathbf{r} (t) + W_{\text{in}} \mathbf{u} (t) + \xi \mathbf{1} \right) \right] + \mathbf{c}^{\ast} \\ =& f \left( \mathbf{u} (t) \right) \end{align*} $$

最終的に学習されたモデルを観察すると、勾配降下法を用いるディープラーニングと比べて圧倒的に軽量であることが分かる。

標準化の問題

論文で著者らが見出した現象として、データを標準化した場合はうまく動作するが、そうでない場合は適切に動作しない例が観察された。

ローレンツアトラクタで学習したリザバーコンピュータのパフォーマンスである。上は標準化していないケース、下は標準化したケースの予測結果だ。