P&Dアドベントカレンダー2018の23日目担当ぽこつんです。
22日に作業はじめてプログラムがうまくいかないまま徹夜で書いているので間違いがあるかもしれません。¹遠慮なくマサカリ飛ばしてください。
23日だからまだ大丈夫って余裕こきすぎた...

1. はじめに

Reservoir Computingというものをご存知でしょうか？
Qiitaにおいても過去に紹介されていたりしますが、あまりメジャーじゃないようでQiitaの記事以外は一口解説程度の情報しか落ちていません。
以前から興味があったので論文を読んでみたところ、もしかしたら結構突っ込んだところまで紹介できるかもと思い今回記事にすることにしました。
人工知能学会などもリザーバーと読んでますが正確にはレザバーらしいです。
またどうせやるならと名著ゼロから作るDeep Learningにあやかって下手なライブラリを使わずnumpyだけでの実装に挑戦してみました。

2. 環境

python 3.7
- matplotlib 3.0.2
- numpy 1.15.4
- scipy 1.2.0

3. Reservoir Computingとは

RNN(Recurrent Neural Network: 過去の出力を現在の入力として扱えるようにすることで時系列などの順番を扱えるようにしたNN)の一種として認識されています。
Echo State Network(以下ESN)とLiquid State MACHINE(以下LSM)というほぼ同時期に考えられた2つのモデルが始まりと言われています。
現状ESNのほうが知られているかもしれません。少なくともGoogleの検索に引っかかりやすい気がします。
ここで述べる内容も主にEcho State Networkについて考えていきたいと思います。

Reservoir Computingは入力層、中間(Reservoir)層、出力層の3層からなっています。
Reservoir特有の特徴としては入力層とReservoir層において重みの更新が行われない点です。
最初にランダムな値で決め打ちしてしまったら、あとのことはReservoir層のノードの内部状態と出力層の重みに仕事を丸投げして出力層の重みを調整することで学習していきます。
Reservoir層でのみ同層間の結合を許しReservoir層のノードが内部状態を持つことでモデル内に短期記憶を用意し、時系列データの持つ因果関係を扱えるようにしているみたいです。
Reservoir(貯水池)の名前自体、この中間層が入力データを貯めているようにみえることから来ているみたいなので、だいたいこんな認識でいいんじゃないでしょうか。

これを使って何が嬉しいか

RCの過去紹介記事でも触れられている通りRNNの場合、過去の情報についても遡って逆伝搬で学習するため回帰なしのNeural Networkに比べて学習コストがとんでもないことになるという問題がありました。
Reservoir Computingの場合、学習の部分を出力層のみに任せることでこの大量の学習コストの問題を解消できるようです。
ただまだまだこれからの分野ということもあり実験の結果とりあえずそうなるんや!!と言った部分も論文の所々で見え隠れしているため今後に期待と言った状況らしいです。

また学習方法として出力層の重みは最小二乗誤差をベースに学習するため従来のRNNだと時間を遡りすぎたことにより起こる伝搬の消失や局所解に嵌ったりといったことが起こりにくくなるようです。隠れ層などから伝搬せずに出力層のみのfittingで一気にカタをつけようとしているからですかね...?

4. モデルの説明

ここから本格的になかでどんなことが起きているのかを説明していこうと思います。

まず各種変数の定義です。

$u(t)$: 時刻tにおける入力
$W_{in}$: 入力層の重み
$K$: 入力層のノード数
$x(t)$: 時刻tにおけるreservoir層の内部状態
$X$: $x(t)$の行列記法
$W_{res}$ : Reservoir層の重み
$N$: Reservoir層のノード数
$W_{out}$: 出力層の重み
$z(t)$: 出力
$K$: 出力層のノード数

とりあえずずっと出てくる奴らはこのくらいです。
他の奴は出てきしだい紹介していきます。といっても核となる数式は3つくらいなのでぱぱっと理解していきましょう。

1つめはreservoir層の内部状態$x(t)$の式です。

x(t) = f((1 - \delta)x(t-1) + \delta(W_{in}u(t-1) + W_{res}x(t-1)))

この式では大きくモデルの性能に関係してくる変数が2つあります。1つめは$\delta$です。$\delta$は漏れ率を表します。$\delta$が小さければ小さいほど第一項の寄与率が大きくなり第二項の寄与率が小さくなります。つまり時刻 $t-1$ までの関数の状態を保存している$x(t-1)$の寄与率が大きくなり、直前の入力を表す$u(t-1)$の寄与率が小さくなるため、$\delta$が小さいほど過去の状態を長く覚えていられるようになるが直前の入力の影響が小さくなります。逆に$\delta$が大きいほど直前の入力はよく覚えているが、過去の状態については極短期間しか覚えられない鳥頭になります。

2つめの性能に大きく関係してくる変数は入力$u(t)$です。式の全体に活性化関数$f$がかかっていますがESNでは$\tanh$を使用します。

ここで$tanh$の数式を考えてみるとわかりますが入力が0に近い場合はほぼ線形をとり、0から離れるほど非線形になります。ものすごく入力が大きかったら出力は1か-1に偏ることになります。よって複雑な式の場合、入力の大きさは一度正規化して1未満にしておいたほうが学習がうまくいきやすくなることがわかります。逆に単純な形を学習させたいのなら入力 $u(t)$ に下駄をはかせて$x(t)$を1か-1に偏らせたほうがすぐに学習が終わると思います(未検証)。

2つめの式は$W_{out}$の学習に関係する式です。

W_{out} = (X^TX + \lambda I)^{-1}X^{T}Y_{target}\\

と言ってもリッジ回帰で正則化した行列$X$と出力としてほしい理想値(教師データ)$Y_{target}$の最小二乗誤差が0となる値を$W_{out}$に代入しているだけです。
行列の最小二乗誤差についてはこちらの方がすごく丁寧にまとめてくださっています。
これでReservoirの学習が完了するというのが正直未だに信じられないのですができるみたいです。

最後に出力$z(t)$に関係する式です。

z(t) = f^{out}(W_{out}x(t))\\

$f^{out}$も$tanh$を使用します。reservoir層の出力と出力層の重みの積和を活性化関数にかけるだけです。まんまですね。

他に初期条件として、
$W_{in}$: -0.1 or 0.1の一様分布
$W_{res}$: 平均0, 分散1の正規分布
$x(0)$: ゼロ行列

となっています。これらは実験の結果これが良さそうだとなったハイパーパラメータの類です。

以上でモデルの計算に関係してくる式は終わりです。

手順のおさらい

手順としては以下のようになります。

1. 各種重みを初期化

2. 訓練データを流し込む

3. 流し込まれたデータに対し全時刻での内部状態を記憶

4. 記憶した内部状態を基にここで初めて出力層の重みを更新

5. 更新データと内部状態からデータを予測

5. その他諸々の知識

学習モデルの精度に関わってくるパラメータの部分をもう少し書いていきます。少し込み入った話になるので興味ない人は飛ばしましょう。

Echo State Propertyとは

Echo State Property(以下ESP)とは与えられたタスクをReservoir Computingが達成するために満足すべき条件のことです。ざっくり言ったらfading Memoryという現在の入力は現在の内部状態に過去の入力や状態よりも大きな影響を及ぼすという機能を満たす必要があります。
つまりESPを持つということはESN は最終的に安定状態になった時に初期状態の影響を全て洗い流してしまえるということみたいです。
ここでESPを満たす場合、$W_{res}$の固有値の絶対値が最大のもの(spectral radius)を1ではないが1に限りなく近い値として持たなければいけないようです。
このあたりは正直自分も全然理解できていない場所で又別の論文を読む必要があるようです。
プログラム的にはspectral radiusで$W_{res}$を割ったあとに0.99をかけておけばなんとかなります。

ここでspectral radiusと$x(t)$の式との関係を考えた人は鋭いです。$W_{res}$が変化するということは$u(t)$の大きさがモデルに与える影響とほとんど同じ影響をspectral radiusはモデルに与えます。
違う点としては、

$u(t)$の大きさは$x(t+1)$における$u(t)$項に影響を与える -> 直近の入力$u(t)$の$x(t+1)$に対する寄与率に影響
spectral radiusは$x(t+1)$における$x(t)$に影響を与える -> 過去のReservoir層の状態すべての$x(t+1)$に対する寄与率に影響
という点です。

6. 実装

実装です。正直コレをやってみたいがために書いていたようなもんです。~~numpyのみを武器に戦っていきます~~ ~~負けました間に合いませんでした~~
ある程度の予測ができるようになりました。

大きくわけて3つのファイルに分けています。²

ファイル名	用途
reservoir_computing.py	ESNの実装をしています train()で学習, predict()で予測
input_generator.py	入力のジェネレータクラスとりあえずsin関数をnumpy配列で作ってくれる
main.py	main関数ここがいちばん汚い

reservoir_computing.py

import numpy as np
from scipy import linalg

class ReservoirNetWork:

    def __init__(self, inputs, num_input_nodes, num_reservoir_nodes, num_output_nodes, leak_rate=0.1, activator=np.tanh):
        self.inputs = inputs # 学習で使う入力
        self.log_reservoir_nodes = np.array([np.zeros(num_reservoir_nodes)]) # reservoir層のノードの状態を記録

        # init weights
        self.weights_input = self._generate_variational_weights(num_input_nodes, num_reservoir_nodes)
        self.weights_reservoir = self._generate_reservoir_weights(num_reservoir_nodes)
        self.weights_output = np.zeros([num_reservoir_nodes, num_output_nodes])

        # それぞれの層のノードの数
        self.num_input_nodes = num_input_nodes
        self.num_reservoir_nodes = num_reservoir_nodes
        self.num_output_nodes = num_output_nodes

        self.leak_rate = leak_rate # 漏れ率
        self.activator = activator # 活性化関数

    # reservoir層のノードの次の状態を取得
    def _get_next_reservoir_nodes(self, input, current_state):
        next_state = (1 - self.leak_rate) * current_state
        next_state += self.leak_rate * (np.array([input]) @ self.weights_input
            + current_state @ self.weights_reservoir)
        return self.activator(next_state)

    # 出力層の重みを更新
    def _update_weights_output(self, lambda0):
        # Ridge Regression
        E_lambda0 = np.identity(self.num_reservoir_nodes) * lambda0 # lambda0
        inv_x = np.linalg.inv(self.log_reservoir_nodes.T @ self.log_reservoir_nodes + E_lambda0)
        # update weights of output layer
        self.weights_output = (inv_x @ self.log_reservoir_nodes.T) @ self.inputs

    # 学習する
    def train(self, lambda0=0.1):
        for input in self.inputs:
            current_state = np.array(self.log_reservoir_nodes[-1])
            self.log_reservoir_nodes = np.append(self.log_reservoir_nodes,
                [self._get_next_reservoir_nodes(input, current_state)], axis=0)
        self.log_reservoir_nodes = self.log_reservoir_nodes[1:]
        self._update_weights_output(lambda0)

    # 学習で得られた重みを基に訓練データを学習できているかを出力
    def get_train_result(self):
        outputs = []
        reservoir_nodes = np.zeros(self.num_reservoir_nodes)
        for input in self.inputs:
            reservoir_nodes = self._get_next_reservoir_nodes(input, reservoir_nodes)
            outputs.append(self.get_output(reservoir_nodes))
        return outputs

    # 予測する
    def predict(self, length_of_sequence, lambda0=0.01):
        predicted_outputs = [self.inputs[-1]] # 最初にひとつ目だけ学習データの最後のデータを使う
        reservoir_nodes = self.log_reservoir_nodes[-1] # 訓練の結果得た最後の内部状態を取得
        for _ in range(length_of_sequence):
            reservoir_nodes = self._get_next_reservoir_nodes(predicted_outputs[-1], reservoir_nodes)
            predicted_outputs.append(self.get_output(reservoir_nodes))
        return predicted_outputs[1:] # 最初に使用した学習データの最後のデータを外して返す

    # get output of current state
    def get_output(self, reservoir_nodes):
        # return self.activator(reservoir_nodes @ self.weights_output) 修正前
        return reservoir_nodes @ self.weights_output # 修正後

    #############################
    ##### private method ########
    #############################

    # 重みを0.1か-0.1で初期化したものを返す
    def _generate_variational_weights(self, num_pre_nodes, num_post_nodes):
        return (np.random.randint(0, 2, num_pre_nodes * num_post_nodes).reshape([num_pre_nodes, num_post_nodes]) * 2 - 1) * 0.1

    # Reservoir層の重みを初期化
    def _generate_reservoir_weights(self, num_nodes):
        weights = np.random.normal(0, 1, num_nodes * num_nodes).reshape([num_nodes, num_nodes])
        spectral_radius = max(abs(linalg.eigvals(weights)))
        return weights / spectral_radius

input_generator.py

import numpy as np
import signalz

class InputGenerator:

    def __init__(self, start_time, end_time, num_time_steps):
        self.start_time = start_time
        self.end_time = end_time
        self.num_time_steps = num_time_steps 

    def generate_sin(self, amplitude=1.0):
        return np.sin( np.linspace(self.start_time, self.end_time, self.num_time_steps) ) * amplitude

    def generate_mackey_glass(self, a=0.2, b=1, c=0.9, d=17, e=10, initial=0.1):
        return (signalz.mackey_glass(self.num_time_steps+200, a=a, b=b, c=c, d=d, e=e, initial=initial) - 0.8)[200:]

main.py

import numpy as np
from input_generator import InputGenerator
from reservoir_network import ReservoirNetWork
import matplotlib.pyplot as plt

T = np.pi * 16 
RATIO_TRAIN = 0.6
dt = np.pi * 0.01
AMPLITUDE = 0.9
LEAK_RATE=0.02
NUM_INPUT_NODES = 1
NUM_RESERVOIR_NODES = 150
NUM_OUTPUT_NODES = 1
NUM_TIME_STEPS = int(T/dt)

# example of activator
def ReLU(x):
    return np.maximum(0, x)

def main():
    i_gen = InputGenerator(0, T, NUM_TIME_STEPS)
    data = i_gen.generate_sin(amplitude=AMPLITUDE)
    num_train = int(len(data) * RATIO_TRAIN)
    train_data = data[:num_train]

    model = ReservoirNetWork(inputs=train_data,
        num_input_nodes=NUM_INPUT_NODES, 
        num_reservoir_nodes=NUM_RESERVOIR_NODES, 
        num_output_nodes=NUM_OUTPUT_NODES, 
        leak_rate=LEAK_RATE)

    model.train() # 訓練
    train_result = model.get_train_result() # 訓練の結果を取得

    num_predict = int(len(data[num_train:]))
    predict_result = model.predict(num_predict)

    t = np.linspace(0, T, NUM_TIME_STEPS)
    ## plot
    plt.plot(t, data, label="inputs")
    plt.plot(t[:num_train], train_result, label="trained")
    plt.plot(t[num_train:], predict_result, label="predicted")
    plt.axvline(x=int(T * RATIO_TRAIN), label="end of train", color="green") # border of train and prediction
    plt.legend()
    plt.title("Echo State Network Sin Prediction")
    plt.xlabel("time[ms]")
    plt.ylabel("y(t)")
    plt.show()

if __name__=="__main__":
    main()

結果

~~予測がうまくいってません~~ 修正しました。

2019/7/31追記

更に修正しました。出力層の活性化関数を消しました。

Sin関数(振幅0.9)の学習結果

Mackey Glass方程式の学習結果(LEAK_RATE=0.05, 最初の200ステップはイレギュラーな値のため学習から外している)

青色の線が学習データです。
縦線の左側までが学習に使ったものです。
縦線の右側緑色の線が学習済みモデルに予測させた出力です。
多少のズレが生じていますが,だいたいの波形は学習できていると思います。
mackey glassの方は...うーんって感じですね.
~~予想していたよりは微妙な結果³になってしまいましたがまだまだ改善の余地はあると思います。~~

2019/9/29追記

比較しやすいように教師データをモデルに予測させている時刻においても描画するようにしました.
またもう少し複雑なモデルとしてMackey Glass方程式も学習させてみました.

最後に

徹夜と相まってかなりきつかったです。
ただこんなことがないとなかなかここまで根詰めて調べないと思うのでいい機会だったと思います。
あと遅れてすみません。
P&Dアドベントカレンダーも残るは2日です。@AK-10, @arabian9ts 素晴らしい記事を楽しみにしています!

参考文献

12/25修正しました ↩
AM 3:00くらいまでは綺麗にかけていました。 ↩
予定ではやっべ、これkerasにpull request送れるわ!!ってレベルでした ↩

ゼロから作るReservoir Computing