2017-11-25

これからスパイキングニューラルネットワークを始める人へ2 : ニューロンモデルの実装

スパイキングニューラルネットワークについて日本語の資料が少ないと思ったので、あんまり詳しくない人に向けて関連情報をまとめてみようかと思います。(説明の記事はいっぱいあるけど、実践的な内容が少ない)
多分、ソースコードとか調べても見つけにくいと思うので、サンプルとして最後にC++のコードもつけようと思います。
なるべく簡単にまとめようと思いますが難しかったらコメントくだせぇ。

※注意：この記事は公開しながら編集していますので、ところどころ変かもしれないですが、よろしくです。

今回は実装です。基本はC++で書きますが、他のソースコードへのリンクも載せるので多分参考になると思います。(リンクはたぶんMatlabのコードが多い)
この記事だけでも、内容がわかるようにはしてますが、一応前回からの続きです。↓
これからスパイキングニューラルネットワークを始める人へ1 : ニューロンのモデルについて - オープンメモ置き場

ケッコー見てくれてる人がいるようで更新が遅くてすみません。。。頑張ります！
追記2018/8/10
修論終わったら頑張って書きます！！

数値解析

前回の記事でニューロンのモデルの式を解説して、その特性を常微分方程式で表しました。その常微分方程式をコンピューターで解くために、数値解析の知識が必要です。
なんで必要かというと、微分の定義を思い出して欲しいんですが

$\displaystyle{ \frac{df(t)}{dt}=\lim_{\Delta h\to0}{\frac{f(t+\Delta h)-f(t)}{\Delta h}} }$

微分 - Wikipedia
こんな感じだったと思います。
げげっ、無限に小さい $\displaystyle\Delta h$ がいます。こいつはコンピューターでは扱えません。コンピューターでは全ての変数(時間も)が一定の幅を持って進めていくしかないですから。この $\displaystyle\Delta h$ を無限に小さくせずに、1とか0.1とかの一定の数値で近似して計算するしかないのです。こいつを数値解析でなんとかします。

f:id:Nonatonic:20171201002519j:plain:w500 — 常微分方程式の可視化の一例

以下では、代表的な方法であるオイラー法とルンゲ=クッタ法を紹介します。
各方法の説明の前に、これらの方法が常微分方程式から何にたどり着ければゴールなのかを説明します。

スタート

まぁスタートは常微分方程式です。
状態変数をxとして、こんな式をコンピュータ内で扱えるようにします。

$\displaystyle{ \frac{dx}{dt}=f(t,x) }$

ゴール

状態変数の初期値を決めます。そしてそこから時間を進めていったら状態がどう変化するかがわかればいいわけです。
ニューロンなら今何mVで、次に何mVになるのかを知りたいわけです。
これはつまり、時間 $\displaystyle t_{n}$ の状態を $\displaystyle x_{n}$ 、時間を $\displaystyle h$ 進めた時間 $\displaystyle t_{n+1}$ 時の状態を $\displaystyle x_{n+1}$ として、以下の式における $\displaystyle g$ を定義することができたらゴールです。

$\displaystyle{ \begin{align} t_{n+1}&=t_{n}+h\\ x_{n+1}&=x_{n}+g(t_{n},x_{n}, h,f) \end{align} }$

オイラー法

一番単純な方法がオイラー法で、こんなの

$\displaystyle{ x_{n+1}=x_{n}+hf(t_{n},x_{n}) }$

こいつは $\displaystyle t$ から $\displaystyle t+h$ の時間の間 $\displaystyle x$ の変化率が一定であることを仮定しています。
方法としては単純で、とりあえずシミュレーションをやるにはこれを使えば何の問題もないです。
変化率が一定であると仮定してるので、一定じゃない時に誤差が大きくなります。
オイラー法 - Wikipedia

ルンゲ=クッタ法

正確には4次ルンゲ=クッタ法です。1次も2次もあるんですけど、4次が1番計算コストと正確性のバランスが良いらしく、4次しか使われないので省略されがちです。こんな感じ

$\displaystyle{ \begin{align} x_{n+1}&=x_{n}+\frac{h}{6}\left( k_{1}+2k_{2}+2k_{3}+k_{4}\right)\\ k_{1}&=f\left(t_{n},x_{n}\right)\\ k_{2}&=f\left(t_{n}+\frac{h}{2},x_{n}+\frac{h}{2}k_{1}\right)\\ k_{3}&=f\left(t_{n}+\frac{h}{2},x_{n}+\frac{h}{2}k_{2}\right)\\ k_{4}&=f\left(t_{n}+h,x_{n}+hk_{3}\right) \end{align} }$

$\displaystyle k_{1}$ から $\displaystyle k_{4}$ は全てグラフの傾きです。複雑に見えますが、実際は精度のために初期値を変えて4回同じ計算をしてるだけです。
証明は難しいため省略しますが、イメージはこんな感じです。

初期値を変化させて4回計算した結果を加重平均します。
式を見てもなんじゃこりゃだと思うので、プログラミングしてみて慣れてください。数値解析でかなりよく使われる常套手段です。
みんな最初は名称に疑問を持って「昨日ルンゲ、食い過ぎたわ」みたいなボケを言うという噂があります。
誤差の比較の図です。

ルンゲ＝クッタ法 - Wikipedia

数値解析の実装

なるべくはそのままコンパイルすれば動くプログラムを見せたいので、ちょっと長くなるかもしれないですけどすみません。
出力はcsvにします。
csvはExcelとかで読み込んで可視化してみて下さい。
できる人はRかPythonで可視化して下さい。
あんまり式の数が多いと分かりづらいのでフィッツフュー・南雲モデルを実装します。

$\displaystyle{ \frac{dv}{dt}=v-\frac{v^{3}}{3}-w+I\\ \frac{dw}{dt}=0.08\left(v+0.7-0.8w\right) }$

フィッツフュー-南雲モデル - Wikipedia
FitzHugh-Nagumo model - Scholarpedia

そのまま書く

上で紹介した方法をそのまま書いていきます。

オイラー法

#include<iostream>
#include<fstream>

int main(){

    //initialize
    const double h = 1.; // 1ステップの時間[ms]
    const double simulationEnd = 100.; // シミュレーション時間[ms]
    double v = 0.; // 状態変数
    double w = 0.;
    double I = 4.;

    //output
    ofstream output("Euler_Result.csv"); // 出力ファイル

    //simulation
    for(double t=0; t<simulationEnd; t=t+h){
        output << t << "," << v <<std::endl;
        v += h * (v - v*v*v/3 - w + I);
        w += h * (0.08 * (v + 0.7 - 0.8*w));
    }

    output.close();
    return 0;
}

ルンゲ=クッタ法

#include<iostream>
#include<fstream>

int main(){

    //initialize
    const double h = 1.; // 1ステップの時間[ms]
    const double simulationEnd = 100.; // シミュレーション時間[ms]
    double v = 0.; // 状態変数
    double w = 0.;
    double I = 4.;

    //output
    ofstream output("RK4_Result.csv"); // 出力ファイル

    //simulation
    for(double t=0; t<simulationEnd; t=t+h){
        output << t << "," << v <<std::endl;
        v_calc = v; // 初期値を変えて計算を行うので
        w_calc = w; // 計算用の変数を用意します
        k_1 = v_calc - v_calc*v_calc*v_calc/3 - w_calc + I;
        l_1 = 0.08 * (v_calc + 0.7 - 0.8*w_calc);

        v_calc = v_calc + h * k1 / 2;
        w_calc = w_calc + h * l1 / 2;
        k_2 = v_calc - v_calc*v_calc*v_calc/3 - w_calc + I;
        l_2 = 0.08 * (v_calc + 0.7 - 0.8*w_calc);

        v_calc = v_calc + h * k2 / 2;
        w_calc = w_calc + h * l2 / 2;
        k_3 = v_calc - v_calc*v_calc*v_calc/3 - w_calc + I;
        l_3 = 0.08 * (v_calc + 0.7 - 0.8*w_calc);

        v_calc = v_calc + h * k3;
        w_calc = w_calc + h * l3;
        k_4 = v_calc - v_calc*v_calc*v_calc/3 - w_calc + I;
        l_4 = 0.08 * (v_calc + 0.7 - 0.8*w_calc);

        v += h * (k_1 + 2 * k_2 + 2 * k_3 + k_4) / 6;
        w += h * (l_1 + 2 * l_2 + 2 * l_3 + l_4) / 6;
    }

    output.close();
    return 0;
}

行列で書く

boost/odeint

次回へ続く！

これからスパイキングニューラルネットワークを始める人へ3 : シナプスについて - オープンメモ置き場

リンク集

Izhikevich Model

Hodgkin-Huxley Model

FitzHugh-Nagumi Model

Hindmarsh-Rose Model

Wilson

2017-11-25

これからスパイキングニューラルネットワークを始める人へ1 : ニューロンのモデルについて

スパイキングニューラルネットワークニューロン

※注意：この記事は公開しながら編集していますので、ところどころ変かもしれないですが、よろしくです。

目次
- テレビでみるAI
- もっと凄そうなことがしてみたいぜ！
スパイキングニューラルネットワークとは
- てかニューロンって何？
ニューロンの活動電位を数式で再現すること、それがニューロンのモデルの使命です。
次回に続く！

テレビでみるAI

今巷で大人気のあの人工知能さんは、大体のものがいわゆるディープラーニングという技術を応用した何かです。ディープラーニングっていうのはニューラルネットワークっていうものを、ネットワークが深くなってもうまくいくようにした感じのやつなのですが、まぁこいつが脳的システムにおいて何か本質的なことを教えてくれるかというと正直よくわかりません。
あいつはニューロンをめっちゃ簡略化してるし、時間的なダイナミクスも持たないし、生命とか意識とかの問題を本気で考えたりする自分のような人間にはちょっと疑わしいものなのです。

もっと凄そうなことがしてみたいぜ！

そう思ったあなたには、スパイキングニューラルネットワークをおすすめします。僕と一緒に深遠なる脳の世界を探検しましょう！

スパイキングニューラルネットワークとは

「スパイキングニューラルネットワークってなんだよ(怒)」って人はこの記事にたどり着かないかもしれませんけど一応説明。
スパイキングニューラルネットワーク(Spiking Neural Network: SNN)は、「脳ってニューロンとシナプスでできてるんだから、ニューロンを数式にして、シナプスを数式にして、プログラムにしてくっつけたら脳になるんじゃね(適当)」って感じの分野です。
説明は適当ですけど出発点は大体こんな感じです。この状態から脱却していくために、とにかく頑張って勉強しましょうねってことで、この記事ではニューロンを数式にすることからコツコツ始めましょう。

↓Wikipediaにもちょっとだけ書いてあります
ニューラルネットワーク - Wikipedia

てかニューロンって何？

ニューロン(Neuron、神経細胞)は細胞の一種です。
生物って細胞がいっぱい集まってできてます。細胞たちは生物がうまく生きていくためにいろんな機能を分担してやりくりしてます。

f:id:Nonatonic:20171210031015p:plain:w400

ニューロンは生物の中でも脳や神経をうまく機能させるために頑張っているのです。

f:id:Nonatonic:20171210031027p:plain:w300

そんなニューロンさんは脳とか神経とかを運営するのに電気を使っています。それが活動電位って呼ばれるものです。ニューロンさんの活動電位は一瞬バンって電圧が上がってすぐ下がるのでスパイク(Spike)とも呼ばれています。
波形はこんな感じです。↓

f:id:Nonatonic:20180129162810p:plain:w200

ニューロンさんのスパイクをちゃんと再現してネットワークにしてあげるので、スパイキングニューラルネットワークって呼ばれてます。

↓これはちゃんとWikipediaに乗ってる！
神経細胞 - Wikipedia
活動電位 - Wikipedia
http://www.brain.riken.jp/jp/aware/neurons.html

※注意：ニューロンのことを調べているあなたは、もしかしたらどこかで「ニューロンは情報を伝達している」みたいな文を見たかもしれません。僕はこれは半分本当だけど半分嘘だと思ってます。研究分野における「物理」と「情報」の間は、あなたが思うより若くて、確定的でなく、わかりにくい分野です。ニューロンが活動電位を持つことは証明されていますが、ニューロンが情報を扱っていることを、みんなが納得する方法で立証できた人はいないように思います。まぁ話半分でいいですけど、みんなが言うから正しいみたいな思い込みには注意で。

ニューロンの活動電位を数式で再現すること、それがニューロンのモデルの使命です。

ということでニューロンは活動電位が特徴で、その電位でなんかの機能を果たしているようです。
ニューロンの特徴をシミュレーションするために、この活動電位を数式にするのです！

とその前にモデルの分類から

活動電位を数式にすることが目的なのはわかってもらえたと思うのですが、個々のモデルの説明の前にとりあえず代表的なモデルたちの名前を紹介しちゃいまっくす。

ホジキン-ハクスレーモデル(Hodgkin-Huxley Model)
フィッツフュー-南雲モデル(FitHugh-Nagumo Model)
インテグレートアンドファイアモデル(Integrate and Fire Model: IF Model)
イジケヴィッチモデル(Izhikevich Model)
ハインドマーシュ-ローズモデル(Hindmarsh-Rose Model)
ウィルソンモデル(Wilson Model)

↓こんなただの羅列でも日本語だと情報がないのです。また今度詳しくまとめますよ。たぶん。
Biological neuron model - Wikipedia
Which Model to Use for Cortical Spiking Neurons?

※参考：英語でもいい人は最初はIzhikevichの論文を読むのが一番だと思う。Izhikevichはロシアのタフガイで、プログラムとかいろんなじゃあ実際どうなのよを教えてくれる。ただ上の論文はIzhikevichが自分のモデルの良さをアピールするためのやつなので参考程度で。

とりあえずホジキン-ハクスレーモデル

一番古くて全部入りなのが、ホジキンさんとハクスレーさんがイカのニューロンを測定して作ったこのモデル。イカのニューロンは巨大で測定しやすい上、非常時には食料にもなるため、イカちゃんは研究対象として優秀！
モデルとしては、生理学的に妥当で、ニューロンの様々な特徴を再現でき、計算コストが多分一番高い。
式はこんな感じ

$\displaystyle{ \begin{align} C_{m}\frac{dV}{dt}&=-\overline{g}_{Na}m^{3}h(V-V_{Na})-\overline{g}_{K}n^{4}(V-V_{K})-\overline{g}_{l}(V-V_{l})+I \\ \frac{dm}{dt}&=\alpha_{m}(V)(1-m)-\beta_{m}(V)m \\ \frac{dh}{dt}&=\alpha_{h}(V)(1-h)-\beta_{h}(V)h \\ \frac{dn}{dt}&=\alpha_{n}(V)(1-n)-\beta_{n}(V)n \end{align} }$

Hodgkin–Huxley model - Wikipedia

結構難しく見えますが、そんなことはないです。
モデルの基本はまず内部状態を表している変数がなんなのかを知ることです。
ニューロンの場合、一番大事な状態を表している変数は膜電位 $\displaystyle V$ です。
膜電位は細胞を包んでいる細胞膜の中と外の電位差で、この膜電位がスパイクを起こします。そして、この膜電位は様々な要因により変動します。
まず、細胞膜はコンデンサーになっています。このコンデンサーのキャパシタンスが膜容量 $\displaystyle C_{m}$ です。

$\displaystyle{ I=C_{m} \frac{dV}{dt} \\ }$

f:id:Nonatonic:20180506233447p:plain:w400

ニューロンには内側の電流 $\displaystyle I_{i}$ と外側の電流 $\displaystyle I$ がありまして、式はこうなります。

$\displaystyle{ I=C_{m} \frac{dV}{dt}+I_{i} \\ }$

そして、ホジキン-ハクスレーモデルにおいては電流が4種類あります。(種類数はモデル次第)
イオン電流3種類(合計 $\displaystyle I_{i}$ )

ナトリウムイオン電流 $\displaystyle I_{Na}$
カリウムイオン電流 $\displaystyle I_{K}$
その他のイオン電流 $\displaystyle I_{l}$

それと

他のニューロンなど外部からの電流 $\displaystyle I$

よってさっきの式はこうなる

$\displaystyle{ I=C_{m} \frac{dV}{dt}+I_{Na}+I_{K}+I_{l} }$

[画像(仮):plain:w400]

イオン電流は細胞膜にあるイオンチャネルというものを通じて細胞内外を出入りしてるおり、基本的にはこういう感じの式に従ってる。

$\displaystyle{ I_{i}=g_{i}(V-V_{i}) }$

f:id:Nonatonic:20180506233035p:plain:w400

イオンチャネルのイオンの透過性が $\displaystyle g_{i}$ であり、膜電位 $\displaystyle V$ と各イオンの平衡電位 $\displaystyle V_{i}$ との差に従ってイオン電流が流れている。まぁ普通にオームの法則な感じ。

詳細は各イオンによって違っていてホジキン-ハクスレーモデルの奴らはこう

$\displaystyle{ \begin{align} I_{Na}&=\overline{g}_{Na}m^{3}h(V-V_{Na})\\ I_{K}&=\overline{g}_{K}n^{4}(V-V_{K})\\ I_{l}&=\overline{g}_{l}(V-V_{l}) \end{align} }$

[画像(仮):plain:w400]

ここで $\displaystyle m、h、n$ はイオンチャネルの開き具合を表していて、0から1の間をとる無次元量。毛穴の開き具合みたいな。こいつらが残りの状態変数。さっきの式に代入して

$\displaystyle{ I=C_{m} \frac{dV}{dt}+\overline{g}_{Na}m^{3}h(V-V_{Na})+\overline{g}_{K}n^{4}(V-V_{K})+\overline{g}_{l}(V-V_{l}) }$

[画像(仮):plain:w400]

あとは $\displaystyle m、h、n$ の時間変化が存在して、全部こんな感じの式

$\displaystyle{ \frac{dx}{dt}=\alpha_{x}(V)(1-x)-\beta_{x}(V)x }$

こいつの両辺を $\displaystyle{\alpha_{x}(V)+\beta_{x}(V)}$ で割ってやると

$\displaystyle{ \frac{1}{\alpha_{x}(V)+\beta_{x}(V)}\frac{dx}{dt}=\frac{\alpha_{x}(V)}{\alpha_{x}(V)+\beta_{x}(V)}-x }$

なんだこの式はってことで

$\displaystyle{ \tau_{x}(V)=\frac{1}{\alpha_{x}(V)+\beta_{x}(V)}\\ x_{\infty}(V)=\frac{\alpha_{x}(V)}{\alpha_{x}(V)+\beta_{x}(V)} }$

こいつらで置き換えると

$\displaystyle{ \tau_{x}(V)\frac{dx}{dt}=x_{\infty}(V)-x }$

簡単のために $\displaystyle V$ を固定で考えてみると、この式では $\displaystyle x$ の時間変化量は $\displaystyle x_{\infty}$ と $\displaystyle x$ の差に依存していていることがわかる。これって $\displaystyle x$ は時間を進めていくと $\displaystyle x_{\infty}$ に近づいていくことを表している。そして $\displaystyle \tau_{x}$ はその変化量がどのぐらいの強さで効いてくるかを決めている。いわゆる時定数ってやつ。
$\displaystyle V$ は時間で変化するので、 $\displaystyle m、h、n$ は、 $\displaystyle V$ を使って計算した目標値 $\displaystyle x_{\infty}$ に向かって、 $\displaystyle V$ を使って計算した時定数 $\displaystyle \tau_{x}$ の分だけ遅れて追従するって感じ。めっちゃ振り回されてる。都合のいい女みたいな。
で、 $\displaystyle m、h、n$ の振り回され方を決めている $\displaystyle \alpha_{x}(V)$ と $\displaystyle \beta_{x}(V)$ が実際はどうなっているかというとこんな感じ

$\displaystyle{ \begin{align} \alpha_{m}(V)&=\frac{0.1(25-V)}{\exp\left(\frac{25-V}{10}\right)-1} &\beta_{m}(V)&=4\exp\left(\frac{-V}{20}\right)\\ \alpha_{h}(V)&=0.07\exp\left(\frac{-V}{80}\right) &\beta_{h}(V)&=\frac{1}{\exp\left(\frac{30-V}{10}\right)-1}\\ \alpha_{n}(V)&=\frac{0.01(10-V)}{\exp\left(\frac{10-V}{10}\right)-1} &\beta_{n}(V)&=0.125\exp\left(\frac{-V}{80}\right) \end{align} }$

というわけで最初の式に戻って来る。

[画像(仮):plain:w400]

これで一応説明終わり。とりあえず役者は揃った。状態変数は $\displaystyle V、m、h、n$ の4つ！
そんなに難しくないって言ったけど、微分がわからない人にはハードルが高くなってしまった。すまぬ。
多分ホジキン-ハクスレーモデルがニューロンのモデルの中で一番式が多いので、これ以上めんどいのはないと思う。たぶん。

常微分方程式

解説だけで結構大変だったけど、これから上の式たちをプログラムで解いていかなきゃならない。こういう形で表された式には名前があって、それが常微分方程式(Ordinary Differential Equation: ODE)っていうのだけれど、コンピューターの中でこいつらを解くためには数値解析を用いないといけないのである。まぁそれは次回やるんだけど、自分で調べたい人のために名前だけ紹介。次回はオイラー法(Euler method)とルンゲ＝クッタ法(Runge-Kutta method: RK4)から始めます。

常微分方程式 - Wikipedia
数値解析 - Wikipedia
オイラー法 - Wikipedia
ルンゲ＝クッタ法 - Wikipedia

次回に続く！

次回は実際にプログラムを書いてとにかく1ニューロンのシミュレーションができるようになるとこまで行きたいですな。
これからスパイキングニューラルネットワークを始める人へ2 : ニューロンモデルの実装 - オープンメモ置き場

目次

数値解析

スタート

ゴール

オイラー法

ルンゲ=クッタ法

数値解析の実装

そのまま書く

オイラー法

ルンゲ=クッタ法

行列で書く

boost/odeint

次回へ続く！

リンク集

Izhikevich Model

Hodgkin-Huxley Model

FitzHugh-Nagumi Model

Hindmarsh-Rose Model

Wilson

目次

テレビでみるAI

もっと凄そうなことがしてみたいぜ！

スパイキングニューラルネットワークとは

てかニューロンって何？

ニューロンの活動電位を数式で再現すること、それがニューロンのモデルの使命です。

とその前にモデルの分類から

とりあえずホジキン-ハクスレーモデル

常微分方程式

次回に続く！