スパコンで脳を再現する
山﨑 匡 (@neuralgorithm)

ようこそ！

電気通信大学 大学院 情報理工学研究科
NumericalBrain.Org

BLSCスプリングスクール 2018年3月28,29日 電気通信大学

スケジュール
13:00- 開会式 & 研究室へ移動
13:30- ごあいさつ & 資料配付 & 自己紹介
13:45- クラスタマシンへのログイン

満員御礼

14:00- 一時間目: 高校生のための計算神経科学入門
14:50- 休憩
15:00- 二時間目: 神経回路シミュレーション事始め
15:50- 休憩
16:00- 三時間目: ランダムネットワークの並列計算
16:50- まとめ ＆ 計算機見学 ＆ 閉会式へ移動
17:15- 閉会式
17:45 終了

資料配布

・授業テキスト (当日配布版)・スライドの縮小コピーを配ります。

今日やらないこと

・人工知能
・深層学習

ノートを取る必要はありません。メモはテキストやスライドに書き込んで下さい。話を聞く方が重要。

もしこういうのを期待して来てたらごめんなさい

自己紹介

クラスタマシンへのログインクラスタマシンへのログイン
1. 各自のユーザ名を配ります (guest##, ##は数字)
2. 端末エミュレータでssh -Y <hostname> -l guest##と入力してください
3. ログインできるとguest##@node00:~$ □一時間目

高校生のための計算神経科学入門と表示されると思います
4. プロンプトに続けてslと入力して何が起きるかを確認してください

ニューロンの数式dV= −gleak (V (t) − Eleak ) + Isyn (t) + Iext (t)dt!
"Isyn (t) = −gsyn (t) V (t) − Esyn
Cτsyndgsyn
dt= −gsyn (t) +!

入力スパイクwj δj (t)j

シナプス結合重みV (t) > θ ⇒ δ(t) = 1, V (t) = Vreset

ニューロンの計算

神経回路シミュレーションとは人工の脳を作ることに相当

実験可能

脳のシステム的理解
・神経回路の挙動の再現
・仮説の検証・実験結果の解釈と予言実験ではできない操作が可能

行動
実験可能？

実験可能

生理

実験可能

解剖

EUのThe Human Brain Project
ヒト全脳シミュレーション
10年で10億ユーロ

JUQUEEN (IBM, 5.9PFs)

Markram et al. Cell (2015)

何をどう計算するか？

1. 何を計算するか？ → 先ほどの式を計算する
2. どう計算するか？ → 数値計算法が必要

高校物理の超簡単なおさらい

高校物理の超簡単なおさらい

問1

問2

時刻0秒で原点(x=0)に静止している車が、速度v=1メートル/秒

右に一定速度vメートル/秒で走っている車があり、時刻1秒でで右に移動を開始した。

原点から位置x(1)メートル、1秒後の時刻2秒で位置x(2)メート
ルにいたとする。vをxで表せ。

問1a: 移動を開始してから1秒後の位置 x(1) を答えよ。
問1b: 2秒後、3秒後、... a秒後の位置をそれぞれ答えよ。
答え

答えx(1) = 1 メートル、x(a) = a メートルv = ( x(2) - x(1) ) / 1 メートル/秒

高校物理の超簡単なおさらい
問3

高校物理の超簡単なおさらい問4

右に一定速度vメートル/秒で走っている車があり、時刻t秒で原点から位置x(t)メートル、Δt秒後の時刻t+Δt秒で位置

問3の設定の元で、x(t+Δt)をx(t), v, Δtで表せ。

x(t+Δt)メートルにいたとする。vをxで表せ。

v = ( x(t+Δt) - x(t) ) / Δt

Δ: 「デルタ」と読む

答え

答え

v = ( x(t+Δt) - x(t) ) / Δt メートル/秒

x(t+Δt) = x(t) + v × Δt

(×は普通のかけ算)

この式の意味を考えよう

もう一歩進めて考えよう

x(t+Δt) = x(t) + v × Δtx(t+Δt) = x(t) + v × Δt・ x(0)とvからx(Δt)がわかる

・時刻tでの位置 x(t)・そのときの速度 v
がわかれば、Δt 後の位置 x(t+Δt) がわかる

・x(Δt)とvからx(Δt+Δt) = x(2Δt)がわかる
・x(2Δt)とvからx(2Δt+Δt) = x(3Δt)がわかる
・ :つまり、最初の位置x(0)と速度vがわかれば、未来の位置 x(t) は(Δt毎に)全て順番にわかるということを言ってるんだけど、わかりますか？

数値計算法を導入するということも言えるんだけど、わかりますか？

問1の答えを数値的に求めよう

x(t+Δt) = x(t) + v × Δtx(0)とvからx(Δt), x(2Δt), x(3Δt), …を順番に求める

以下はクラスタマシンで作業する→ 前方オイラー法 (大学2年後期の「数値計算」で学ぶ)

・lsと入力するとファイルの一覧が表示される
・car.c というファイルがあるはず

ここからクラスタマシンで実際にプログラムを実行します

ちなみに v は一般に v(t) でかまわない・nano car.c
として、ファイルを開く(nanoはエディタ)

car.cの説明

プログラムを実行する#include<stdio.h>

int main ( void ){double t = 0; // 時刻: 最初は0秒から
double x = 0; // 位置: 最初は原点から0メートルから
double v = 1.0; // 速度: 1メートル/秒
double dt = 1.0; // 時間の刻み幅: 1秒ずつ進める
while ( t < 10.0 ) { // 10 秒間繰り返し

以下はクラスタマシンで作業する
・nanoで^X (Controlを押しながらXを押す)→ nano が終了する
・gcc -std=c99 -O3 -o car car.c(TODO: おまじないの説明をする)

printf ( "%f %f\n", t, x ); // 今の時刻と位置を表示
x = x + v * dt; // 次の時刻の位置を計算
t = t + dt; // 時間をdt秒進める・./car
計算結果の数字列が表示されるはず}}

return 0;

計算結果を表示する

こんなのが表示されるはず・./car > car.dat
計算結果をファイルに保存
・gnuplot

位置 (メートル)

以下はクラスタマシンで作業する

gnuplotはグラフを表示するアプリ
・plot ’car.dat’

時間 (秒)

・終わったら Ctrl + D でgnuplotを終了する
・時間があったら、vをtに書き換えて等加速度運動にする
書き直したら^X → Y → リターン

一時間目まとめ

一時間目まとめ・脳はニューロンがシナプスで繋がったネットワークである

・脳はニューロンがシナプスで繋がったネットワークである

・1個のニューロンとシナプスの振る舞いは数式で書ける

・1個のニューロンとシナプスの振る舞いは数式で書ける

・その数式は数値的に答えを求めることができる

・その数式は数値的に答えを求めることができる

・ので、全ニューロン・全シナプスの数式を全てプログラム

・ので、全ニューロン・全シナプスの数式を全てプログラムして数値的に解けば、(原理的には)脳と同じ動作をする

して数値的に解けば、(原理的には)脳と同じ動作をする

プログラム (いわばコンピュータ上の脳) を作成できるプログラム (いわばコンピュータ上の脳) を作成できる
二時間目は実際にニューロンとシナプスの式を解きます

ちなみに、( x(t+Δt) - x(t) ) / Δtは、Δtが十分小さいときdx/dt (もしくは xʼ)と書く。これをxの (時間に関する)微分という。もちろん v(t) = dx/dt である。

10分休憩！

二時間目

神経回路シミュレーション事始め

内容

ニューロン1個のシミュレーション

外部電流 Iext (t)
膜抵抗 R・ニューロン1個のシミュレーション・ニューロン2個 (繋がない) のシミュレーション・ニューロン2個 (繋ぐ) のシミュレーションをやります。

膜電位 v (t)積分発火型ニューロンモデル

ニューロン1個のシミュレーション
0’neuron1.dat’

膜電位 v (ミリボルト)-10dv= − (v (t) − Vleak ) + RIext (t),dtv (t) > θ ⇒ Spike (t) = 1, v (t) ← Vreset ,-20

τ-30-40v (0) = Vinit .

-50-60-70

何だこれ？と思っても心配無用！説明します！

0 100 200 300 400 500 600 700 800 900 1000
時間 (ミリ秒)

積分発火型ニューロンモデル

dv= − (v (t) − Vleak ) + RIext (t),dtv (t) > θ ⇒ Spike (t) = 1, v (t) ← Vreset ,τ

積分発火型ニューロンモデル

τdv= − (v (t) − Vleak ) + RIext (t),dt()
Spike (t) 1 ( )v (0) = Vinit .

積分発火型ニューロンモデル
τdv= − (v (t) − Vleak ) + RIext (t),dt()
Spike (t) τ(タウ),1 ( )Vleak, R は定数
Δt も定数式変形をすると、積分発火型ニューロンモデル

dv= − (v (t) − Vleak ) + RIext (t),dtv (t) > θ ⇒ Spike (t) = 1, v (t) ← Vreset ,τv (0) = Vinit .
v (0) と Iext (t) を決めればあとは順番に計算できる何だこれ？と思っても心配無用！説明します！

積分発火型ニューロンモデルdv= − (v (t) − Vleak ) + RIext (t),dtv (t) > θ ⇒ Spike (t) = 1, v (t) ← Vreset ,τv (0) = Vinit .

θ (シータ): 定数(閾値)Vreset: 定数

何だこれ？と思っても心配無用！説明します！

積分発火型ニューロンモデル

dv= − (v (t) − Vleak ) + RIext (t),dtv (t) > θ ⇒ Spike (t) = 1, v (t) ← Vreset ,τv (0) = Vinit .

何だこれ？と思っても心配無用！説明します！

コードを試すneuron1.c

1. #include <stdio.h>2.
3. #define TAU 20.0
4. #define V_LEAK -65.0
5. #define V_INIT (V_LEAK)
6. #define V_RESET (V_LEAK)
7. #define THETA -55.0
8. #define R 1.0
9. #define DT 1.0
10. #define T 1000.0
11. #define I_EXT 12.0
12.
13. void loop ( void )
14. {
15.
double t = 0;
16.
double v = V_INIT;
17.
int spike = 0;18.
while ( t < T ) {19.
printf ( "%f %f\n", t, (spike ? 0.0 : v ) );20.
double dv = ( DT / TAU ) * ( - ( v - V_LEAK ) + R * I_EXT );21.
spike = ( v > THETA ) ? 1 : 0;22.
v = ( v > THETA ) ? V_RESET : v + dv;
23.
t = t + DT;
24.
}
25. }
26.
27. int main ( void )
28. {
29.
loop ( );30.
return 0;31. }初期値の設定・gcc -std=c99 -O3 -o neuron1 neuron1.c・./neuron1
数字が表示されるのを確認
・./neuron1 > neuron1.dat
・gnuplot・plot ’neuron1.dat’ with line
・[課題] 外部電流の強さを10.0 (ナノアンペア)や15.0(ナノアンペア)に変更して、発射されるスパイク数の変化を調べてください

要説明

ニューロン2個のシミュレーション
neuron2.c(抜粋)1. void loop ( void )2. {3.
double t = 0;
4.
double v [ 2 ] = { V_INIT, V_INIT - 10.0 };5.
int spike [ 2 ] = { 0, 0 };6.
while ( t < T ) {7.
printf ( "%f %f %f\n", t, ( spike [ 0 ] ? 0.0 : v [ 0 ] ),8.
( spike [ 1 ] ? 0.0 : v [ 1 ] ) );9.
double dv [ 2 ] = { 0.0, 0.0 };10.
for ( int i = 0; i < 2; i++ ) {11.
dv [ i ] = ( DT / TAU ) * ( - ( v [ i ] - V_LEAK ) + R * I_EXT );
12.
}
13.
for ( int i = 0; i < 2; i++ ) {14.
spike [ i ] = ( v [ i ] > THETA ) ? 1 : 0;15.
v [ i ] = ( v [ i ] > THETA ) ? V_RESET : v [ i ] + dv [ i ];
16.
}
17.
t = t + DT;
18.
}
19. }

要説明

コードを試す・gcc -std=c99 -O3 -o neuron2 neuron2.c・./neuron2
数字が表示されるのを確認
・./neuron2 > neuron2.dat
・gnuplot・plot ’neuron2.dat’ using 1:2 with line,’neuron2.dat’ using 1:3 with line

ニューロン2個のネットワーク

指数関数型シナプス
1. void loop ( void )2. {double t = 0;
3.
double g [ 2 ] = { 0.0, 0.0 };4.
double v [ 2 ] = { V_INIT, V_INIT - 10.0 };5.
int spike [ 2 ] = { 0, 0 };6.
7.
while ( t < T ) {8.
printf ( "%f %f %f\n", t, ( spike [ 0 ] ? 0.0 : v [ 0 ] ),9.
( spike [ 1 ] ? 0.0 : v [ 1 ] ) );10.
double dg [ 2 ] = { 0.0, 0.0 };11.
double dv [ 2 ] = { 0.0, 0.0 };12.
for ( int i = 0; i < 2; i++ ) {13.
dg [ i ] = ( DT / TAU_SYN ) * ( - g [ i ] + W * spike [ ( i + 1 ) % 2 ] );14.
dv [ i ] = ( DT / TAU ) * ( - ( v [ i ] - V_LEAK ) + g [ i ] + R * I_EXT );
15.
}
16.
for ( int i = 0; i < 2; i++ ) {17.
spike [ i ] = ( v [ i ] > THETA ) ? 1 : 0;18.
g [ i ] = g [ i ] + dg [ i ];19.
v [ i ] = ( v [ i ] > THETA ) ? V_RESET : v [ i ] + dv [ i ];
20.
}
21.
t = t + DT;
22.
}
23.
24. }

network2.c(抜粋)スパイク {0,1}

シナプス後膜電位10.90.8

重み w

0.7
0.6
0.5
0.4
0.3
0.2
0.10(時間)

τsyn

0100

200300

400

500

600

700

800

900

1000

(時間)

dg
= −g(t) + w · Spike(t)dt

要説明

二時間目まとめ課題・gcc -std=c99 -O3 -o network2 network2.c
・./network2
数字が表示されるのを確認
・./network2 > network2.dat

・ニューロン1個のシミュレーションから始めて、ニューロン2個のネットワークまで作成した・シナプスの結合重みを変えるだけでネットワークの挙動は・gnuplot・plot ’network2.dat’ using 1:2 with line,

変化した・脳はネットワークの繋ぎ方がとても重要である’network2.dat’ using 1:3 with line
・[課題] Wの符号を負にして再度試してください

10分休憩！

三時間目

ランダムネットワークの並列計算

ニューロン4000個のネットワーク
1. void loop ( void )random.c
2. {(抜粋)
double t = 0.;
3.
timer_start ();4.
ニューロン数 (4000)
while ( t < T ) {5.
for ( int i = 0; i < N; i++ ) {6.
calculateSynapse ( i );7.
updateMembranePotential ( i );8.
}9.
outputSpike ( t );10.
t = t + DT;
11.
12. }
13. double elapsedTime = timer_elapsed ();14. printf ( "Elapsed time = %f sec.\n", elapsedTime);
15. }
16.
要説明プログラムを実行してみる

計算が遅すぎる問題とその解決策

・make random
自動的にコンパイルされる・計算時間が結構かかりますよね (16秒くらい)・./random

・ここまで1個のCPU (コア) しか使っていません

・spike.datができる

・今使っているクラスタマシンは全部で192コアです

・gnuplot・ので、192コア全部使わないともったいない→ 並列計算・plot ’spike.dat’ with dots
計算時間も測定する

例: 8ニューロンを独立の4CPUで計算
123

4

5678t

8ニューロンを4CPUで計算

CPU 1

CPU 2

CPU 3

CPU 4

1357246

スパイク交換 = MPIによる通信
ネットワークを介して、異なるCPU間で情報を転送する
MPI (Message Passing Interface) というライブラリがある

8計算

12345678

12345678

今日は一番簡単なMPI̲Allgatherを使います
スパイク交換t+Δt

spike [ ]計算1234567CPU 28
スパイク交換t+2Δt1234567CPU 18CPU 3

1

spike [ ]234MPI̲Allgather
56CPU 478

12345678

12345678

12345678

12345678(時間)random_mpi.c
1. void loop ( const int mpi_size, const int mpi_rank )
(抜粋)
2. {
3.
const int n_each = ( N + mpi_size - 1 ) / mpi_size;4.
int spike_local [ n_each ];4000 / コア数
5.
double t = 0.;
timer_start ();6.
(例えば100コアあれば40)
while ( t < T ) {7.
for ( int n = 0; n < n_each; n++ ) {8.
calculateSynapse ( n, n_each * mpi_rank );9.
updateMembranePotential ( n, n_each * mpi_rank,10.
spike_local );
11.
}
12.
MPI_Allgather ( spike_local, n_each, MPI_INT, spike, n_each,MPI_INT, MPI_COMM_WORLD );13.
if ( mpi_rank == 0 ) { outputSpike ( t ); }t = t + DT;
14.
15. }
16. double elapsedTime = timer_elapsed ();17. if ( mpi_rank == 0 ) { printf ( "Elapsed time = %f sec.\n",elapsedTime); }
18. }
要説明

計算時間はどう変化するか？

コードを試す・make random_mpi
・mpirun -hostfile hostfile -np 16 ./random_mpi
・計算時間を測定する
・-np の後ろの数字を変えて測定する
10人いるから、1, 2, 4, 8, 16, 32, 64, 128, 160, 192をそれぞれ1人ずつ順番に試そう

三時間目まとめもし、CPU数を2倍にすると計算時間が1/2になったら、理想的な並列化がなされている → 強スケーリング
・4000ニューロンからなるランダムネットワークのシミュプログラムの全ての箇所が並列化できるわけではない
並列化できない箇所に計算時間がかかっていたら、できる箇所を頑張っても意味が無い → アムダールの法則

なので、並列化は必ずしもあらゆる問題に対する解決策(いわゆる銀の弾丸) ではない

レーションをおこなった
・MPIを使った並列計算で、計算を高速化したまとめ・脳は作れる！

全体まとめ

・ただし、計算機を上手に使う必要あり・I類ならそのような使い方をばっちり学べます

質問？

閉会式に移動

おしまい！