Shoubu の猫小脳シミュレーションは，スーパーコンピュータのアプリケーションの大部分には参考にならない

1. 違う実験結果を比較すること

この記事では，前の記事に対して [Y+17] の著者（以下著者）が比較できないとコメントした計算時間の結果を比較しています．ただし，著者による論文 [Y+17] では [K+14] と（違うニューロンモデルであるにも関わらず）ニューロン数という基準で揃えるだけで比較しています．さらに，著者による別の論文 [L+14] では [YI13] と（顆粒層と小脳全体という異なる脳の部分であるにも関わらず）ニューロン数という基準で揃えるだけで比較しています．
特に [Y+17] と [K+14] の比較では，入力がおよそ 10,000 倍となる [K+14] に対して，数千倍の時間差 (0.783秒対2481秒) がついたことで [K+14] はとても遅いと [Y+17] は評価しています．

また，コメントでは猫の小脳全体のシミュレーションだけに限定することを強調していますが，「小脳，大脳皮質，猫限定，猫以外」，どの脳シミュレーションもスパイキングニューラルネットワークという同一の枠に入っています．どれも微分方程式の計算が根幹であることに変わりなく，コンピュータの計算性能を比較する上で大きな違いはありません．

コメントに沿って今後これらの論文が修正されるかもしれませんが，以前の記事とこの記事で比較に用いた基準は[Y+17], [L+15] と同程度か，それより細かくなっています．

著者の論文が修正され，論文でも「猫の小脳全体のシミュレーションは猫の顆粒層や猫の大脳皮質のシミュレーションと比較できない」と書き換えられたと仮定します．それでも，これらの記事のように比較することは，スパイキングニューラルネットワークという同一モデルのシミュレーションを比較するという点で，スーパーコンピュータの利用者にとって参考になる情報だと考えます．その際，同一モデルだが対象が違うという著者の主張を，利用者が判断できるように，前の記事とこの記事では比較する基準を明確に示しています．以下では，[Y+17] の比較と同様の粗い基準と [Y+17] より細かい基準の2通りで比較します．そして比較する対象を明確にするために，図で表します.

2. 比較基準

次の例を考えます．辺に重みがついたグラフ構造 G を入力に取って，G に含まれる辺の重みの総和を求める 2つのプログラム A, B があり，異なるコンピュータ C, D でそれぞれ実行するとします．

A を C で実行したとき，2,500秒かかり, B を D で実行したとき，0.7秒だった，という報告がありました．(A, C) の組と (B, D) の組はどちらが速いでしょうか?

ただし，報告をよく読むと入力が違いました．(A, C) の組の入力 G1 は頂点数 10^9, 辺の本数 10^13 で, (B, D) の組の入力 G2 は頂点数 10^9, 辺の本数 10^9 であるときの結果でした．(B, D) の組はそれでも速いでしょうか?

この例を単純化した場合と，それに対する論文 [Y+17] の内容と [Y+17] の著者からのコメントを次の図で示します．ただし，図に含まれるグラフの重みの値は意味を持ちません．

見た目の速さと本質的な速さを区別するために，まず入力長（+α）を基準に取ります．

2.1 粗い基準（シナプス数（入力長）と時間ステップ）

脳シミュレーションのプログラムの入力は，グラフ構造です．グラフ構造の頂点がニューロン，辺がシナプスにそれぞれ対応します．グラフ構造を隣接リストで表すと，その大きさ（入力長）はシナプス数 n に比例します．

シナプスによるニューロン間の情報伝達を計算するとき，すべてのシナプスについて計算しなければなりません．10^9 ニューロン，10^13 シナプスの脳を対象とするとき，10^13 に比例する時間をかけなければ，全部のシナプスを処理することはできません．このため，シナプス数 n を基準にするのはニューロン数より妥当です．

そして，脳シミュレーションにおいて，時間ステップ u を1ミリ秒から 0.1ミリ秒に変更すれば10倍の計算時間がかかります．

そこで, (n/u) を基準にして，1秒あたりのシミュレーションにかかる時間 t を (n/u) で割った値 ut/n を比較します.

2.2 細かい基準（微分方程式の演算コスト）

2.1 では [Y+17] と同レベルとなる，入力に関する値だけを基準にしました．ここではそれより細かく，脳シミュレーションの大部分を占める微分方程式の演算コスト L を基準にします．
L = ([ニューロン数 m] / [時間ステップ u]) x ([微分方程式で使う値の参照コスト L1] + [微分方程式の解を求めるコスト L2]),
L1 = [隣接ニューロン数]x[メモリ参照コスト L3] = ([シナプス数]/[ニューロン数])x[メモリ参照コスト L3],
L3 = [プロセッサ内メモリの参照コスト]x[同じプロセッサのメモリに値が入っている確率 p] + [プロセッサ外メモリの参照コスト] x (1-p).
L2 は [I04] に示された，ニューロンモデルに対する演算コストを用います．
L3 において，プロセッサ内メモリの参照コストは 1 とし，プロセッサ外メモリの参照コストは（[渡辺14] からプロセッサ内メモリ参照は最大100サイクルかかることから）100サイクルの2倍の200サイクルを推定コストとして用います.

L2 は入力長に依存しない定数であるのに対して，L1 は入力長であるシナプス数に比例する項です．このため，十分大きな入力長で，シナプス数がニューロン数より十分大きいとき，L = (m/u) x (n/m x L3 + 定数) であり n/u x L3 とほぼ等しいとみなすことができます．さらに L3 を定数とみなすと，シナプス数と時間ステップによる基準で比較することは，十分大きな入力が与えられたときの，微分方程式演算コストを基準としたときの比較の近似になっています．

3. 比較

A. 猫の小脳モデル

最も緻密
[Y+17] は Runge-Kutta の計算を2次で打ち切っていて，[YN12] は4次まで計算しています．[YN12] の方が近似誤差が小さいので，緻密になっています.
最速
猫の顆粒層のシミュレータ [L+15] と猫の小脳全体のシミュレータの違いは，扱える細胞種（関連するパラメータ）であり，[L+15] にパラメータを追加すれば，本質的に違いはありません．このため，[L+15] のシミュレータを拡張した場合を C. 猫の顆粒層を含むシミュレーションで，計算速度を推定します．

B. 10^9 ニューロン

[Y+17] にはシナプス数に関する記述が 4. Discussion における「シナプス数は 10^9 のオーダーである」しかないため，[Y+17] のシナプス数をニューロン数 x 4 と多めに推定します.

粗い基準

細かい基準

プロセッサ内メモリに値が含まれる確率 p

メモリ参照コスト L3 に含まれる p を推定します．1ニューロンあたりの平均シナプス数（隣接するニューロン数），各プロセッサに均等にニューロンを割り振ったときの1プロセッサあたりのニューロン数 m/q を表にまとめます．

1プロセッサに含まれるニューロン数 m/q より，1ニューロンあたりの平均シナプス数が大きければ，プロセッサ外にあるニューロンを参照する必要があります．参照するプロセッサ外のニューロン数が7番目の項目（2行目 3,896.5）になります．1行目，4行目では，隣接するニューロン数は1プロセッサに含まれるニューロン数より小さくなっているので，ニューロンのプロセッサへの理想的な割り当てが存在するなら，それぞれのニューロンに対して，隣接するニューロンはすべて同じプロセッサに収まる可能性があります．

[K+14] はランダム結合なので，別の推定が必要になります．10^9 ニューロンから同一プロセッサ内の 2,600 ニューロンを選ぶ割合は 2600/1.11x10^9 < 0.001 になるため，ここでは p = 0.001 とおきます.

微分方程式解を求める演算コスト L2

[I04] に示された演算コストを使います.
[A+09] が用いた Izhikevich モデルでは L2=13, [K+14] が用いた Integrate-and-Fire モデル (IF) では L2=5です. [Y+17] が用いた conductance-based モデルは Hodgkin-Huxley モデルに対応します．しかし, [I03] では 0.1ミリ秒ごとに微分方程式を解いて1ミリ秒のコストを 1,200 としている一方, [Y+17] は時間ステップが1ミリ秒なので，演算コストは 1/10 になります．さらに，[Y+17] は Runge-Kutta の計算を2次で打ち切っていて，4次まで計算した [YN12] の 50% の計算時間だと [YI13] で報告しているので，演算コストはさらに半分になり，最終的なコストは L2=60 になります．

より細かく見るならば，[Y+17] ではシナプス入力の種類は AMPA しかなく，
シナプスの可塑性も実装されていませんが，[A+09] ではシナプス入力について4種類 AMPA, NMDA, GABA A, GABA B を取っていて， スパイクタイミング依存シナプス可塑性 (STDP) を実装している，などの違いがあります（この点では [A+09] は [Y+17] より緻密で計算コストがかかります）．しかし，この比較ではこれらの違いまでは考慮せず，同一とみなして（[Y+17] に有利な条件で）比較します．

これらの値を用いて演算コストを求め，それを基準にした計算時間を比較した結果が次の表です．

[K+14] の演算コストは [Y+17] の演算コストのおよそ30,000倍です．このことから，[Y+17] の計算時間が [K+14] より短いのは，Shoubu の計算が速いのではなく，演算コストの低い問題を扱っただけであることがわかります．そして，演算コストに対する計算時間において [Y+17] は [A+09], [K+14] に対して10倍多く時間がかかっています．

計算時間よりも深刻なのは，[Y+17] がバンド幅に負荷を掛けない，スーパーコンピュータのアプリケーションの大部分とは異なる性質を持っていることです．

ネットワーク構造とプロセッサ割り当て

[Y+17] のシミュレータに入力されるネットワーク構造は，規則的な構造で頂点あたりの辺本数は平均すれば数本です（プルキンエ細胞のように辺本数の多い細胞は少数あります）．このとき，プロセッサ間通信が少なくなるように，ニューロンをプロセッサに割り当てることが可能で，プロセッサ間通信を含む演算コストは小さく抑えられます．

一方，[K+14] のシミュレータに入力されるネットワーク構造は，ランダム構造で，頂点あたりの辺本数はおよそ 10^4 本です．このとき，プロセッサ間通信が少なくするように，ニューロンをプロセッサに割り当てることは難しく，どのように割り当てても，プロセッサ外への通信がほとんどになります．つまり，プロセッサ間通信を含む演算コストは大きくなります．

この [Y+17] と [K+14] の差は，そのままバンド幅への負荷の違いになります．スーパーコンピュータにおける大部分のアプリケーションは，バンド幅への負荷が掛かるものが多くを占めます．これより，バンド幅への負荷が軽い [Y+17] の結果は，スーパーコンピュータで実行する大多数のアプリケーションの参考になりません．

招待講演の聴衆はスーパーコンピュータを利用する研究者でしたが，講演者は脳の小脳に限定した話をしているので，聴衆である研究者の大多数の参考にならなくても仕方がありません．

C. 猫の顆粒層を含むネットワークシミュレーション

猫の顆粒層から猫の小脳へのシミュレータの拡張
[L+15] は, [YI13], [Y+17] と同じニューロンモデルを用いています．[Y+17] の表1, 2, 3 にあるように，細胞種の違いはパラメータの違いなので，[L+15] の脳シミュレータが猫の顆粒層以外の小脳の細胞を扱うには，基本的にパラメータを追加すれば良いことになります．その他の細かい違いを考慮する必要はありますが，本質的には，[L+15] の脳シミュレータを猫の小脳を扱えるようにするための，技術的困難はありません．

ニューロン数を基準とする比較
[L+15] のニューロン数は，顆粒細胞が 102,400, ゴルジ細胞が 1,024 です．さらに，10倍のニューロン数でも計算時間はほぼ同じという推定が図11に示されています．[Y+17] (細胞数の近いシミュレーション, 3.3. Computer simulation in a realistic scenario) のニューロン数は，顆粒細胞が 1,048,576, ゴルジ細胞が 1,024, その他が 66 です. 

小脳全体の全細胞における顆粒層（顆粒細胞, ゴルジ細胞）の細胞数の割合 99.99% 以上なので，[Y+17], [L+15] の主張（ニューロンモデルが違う [K+14], [YI13] でも比較できる）に従えば，[L+15] の計算時間0.028秒は [Y+17] の計算時間 0.72 秒程度より十分短い，というのは妥当です．

シナプス数を基準とする比較
[Y+17] の主張（ニューロン数を基準に比較する）ではなく，シナプス数を基準に [L+15] と [Y+17] を比較します．

仮定1 [L+15] のシナプス数 n が少ないから計算時間が短いだけで，シナプス数が多ければ [Y+17] の計算時間と変わらない.
[Y+17] のシナプス数を 4x[ニューロン数] = 4m と仮定しているので，顆粒層のシナプス数を n とすると，2つのシナプス数基準の計算時間が等しいということは，

(0.0256 / n) = (0.72 / 4m),
n = (0.0256/0.72)x4m < 0.142m < m

となります．この仮定の下では，シナプス数がニューロン数より大幅に少なくなり, 脳を表すグラフ構造が非連結になります．つまり，シミュレーションする1つの脳として意味を成しません．これは仮定1が誤っていたからです．
よって，[L+15] の計算時間が短いのは，シナプス数が少ないからではなく，[L+15] が [Y+17] よりも同じ計算コストを早く計算するからだ，ということがわかります．

仮定2 [L+15] の顆粒層を表すグラフ構造のシナプス数は，連結になる最低限の数 n = m-1 しかない．

これは [Y+17] と [L+15] の比較において，[Y+17] に有利な仮定です．さらに，小脳全体を表すには [Y+17] に対する仮定と同様に 4m のシナプスが必要だと仮定します．このとき，[L+15] のシミュレータを拡張して小脳全体を扱えるようにしたとき，1秒当たりのシミュレーションに必要な計算時間は 0.028 の4倍である 0.112 秒で行えると予想されます．これは [Y+17] における1秒あたりのシミュレーションに必要な計算時間 0.72 秒前後より，十分短い予想時間です．

D. 小脳全体

実験内容は次の通り. 
[B+05] 8の字型の軌跡の学習
[C+08] 実ロボットアームの制御
[YI13] パブロフの遅延型瞬目反射条件づけ, バッティングロボット
[L+16] フィードバック制御ループ, フィードフォワード制御ループ
[Y+17] 視運動性反応

著者の研究以外では，異なる研究グループでも，論文中にシナプス数（の概数）がわかりやすく示されています．そして，シナプス数はニューロン数の20倍から100倍超になっています．

[L+16] を参照すると，小脳の研究において大きなニューロン数や実時間計算が必須でないことがわかります．特に, 実験において [C+08] や [YI13] のようにコンピュータの外側とやり取りしない限り，実時間計算でなくても研究を実施できます．

10^9 ニューロンの研究 [K+14] で用いた汎用脳シミュレータ NEST 用の小脳シミュレーション [L+16] のコードが公開されています [G17]. よって現在では, [K+14] の計算環境で小脳のシミュレーションを行うことは容易なので，小脳研究に関しても [K+14] と [Y+17] を比較することは参考になるでしょう．NEST は C++ が動作する環境なら実行できるので，手元のコンピュータで小さなモデルを試してから，大規模なモデルをスーパーコンピュータでシミュレーションする，という柔軟な対応が可能です．NEST を用いた（小脳に限らないさまざまな）脳シミュレーションの論文一覧はこちらにあります．

4. 比較結果

猫の小脳
最大 [Y+17]
最速 [Y+17] 0.72秒前後（ただし，[L+15] のシミュレータを拡張すれば 0.112 秒以内と大幅に短縮されると推定される）
最も精緻 [YN12] (Runge-Kutta 4次) （[Y+17] ではない）

10^9 ニューロン
最大 [K+14] シナプス数が最大
最速 [A+09] 演算コストあたりの計算時間が最速
最も精緻 [K+14] （最大の演算コスト）

5. まとめ

この記事では2つの基準で脳シミュレーションの計算時間を比較しました．
1つは，シミュレータの入力長を表すシナプス数とシミュレーションの時間ステップを基準に取りました．
もう1つは，より細かく，シミュレーションの根幹である，微分方程式を解く演算コストを基準に取りました．
そして，1つ目の基準は，シナプス数が十分大きいときに2つ目の細かい基準の近似になっていることを示しました．
これらの基準では [Y+17] が過去の研究より遅いことに加えて，[Y+17] はバンド幅に負荷を掛けないアプリケーションなので，[Y+17] の結果はスーパーコンピュータのアプリケーションの大部分には参考にならないこともわかりました．

PEZY Computing 社長らによる論文では，

特に，猫の小脳をリアルタイムで実現できたことは ZS-1 の有効性を端的に示した成果といえる．

と主張していますが，バンド幅を必要とするようなスーパーコンピュータのアプリケーションの大部分には参考にならず，有効性はごく限られたアプリケーションについてしか示されていません．

参考文献

[A+09] R. Ananthanarayanan et al. (2009) The Cat is Out of the Bag: Cortical Simulations with 10^9 Neurons, 10^{13} Synapses. In Supercomputing 09: Proceedings of the ACM/IEEE SC2009 Conference on High Performance Networking and Computing (Portland, OR). https://doi.org/10.1145/1654059.1654124
[B+05] C. Christian Boucheny et al. (2005) Real-Time Spiking Neural Network: An Adaptive Cerebellar Model. Proceedings of the 8th international conference on Artificial Neural Networks: computational Intelligence and Bioinspired Systems (IWANN'05), 136-144. https://doi.org/10.1007/11494669_18
[C+08] R. R. Carrillo et al. (2008) A real-time spiking cerebellum model for learning robot control. Bio Systems 94(1-2):18-27. http://dx.doi.org/10.1016/j.biosystems.2008.05.008
[G17] J. Garrido (2017) SpikingCerebellum, Neuron-model implementation for cerebellum learning. Code for NEST Simulator https://github.com/jgarridoalcazar/SpikingCerebellum
[I04] E. M. Izhikevich (2004) Which Model to Use for Cortical Spiking Neurons? IEEE Transactions on Neural Networks, 15, 1063-1070.
https://doi.org/10.1109/TNN.2004.832719
[K+14] S. Kunkel et al. (2014) Spiking network simulation code for petascale computers. Frontiers in Neuroinformatics 8: 78. https://doi.org/10.3389/fninf.2014.00078
[L+15] J. Luo et al. (2015) Real-Time Simulation of Passage-of-Time Encoding in Cerebellum Using a Scalable FPGA-Based System. IEEE Transactions on Biomedical Circuits and Systems 10(3): 742-753. October 2015. http://dx.doi.org/10.1109/TBCAS.2015.2460232
[L+16] N. Luque et al. (2016) Distributed Cerebellar Motor Learning: A Spike-Timing-Dependent Plasticity Model. Frontier Computational Neuroscience. 2016; 10:17. Published online 2016 March 2. http://dx.doi.org/10.3389/fncom.2016.00017
[渡辺13] 渡辺宙志 (2013) 計算科学技術特論A 第9回 高速化チューニングとその関連技術２
http://www.aics.riken.jp/jp/course/%E7%AC%AC9%E5%9B%9E+%E9%AB%98%E9%80%9F%E5%8C%96%E3%83%81%E3%83%A5%E3%83%BC%E3%83%8B%E3%83%B3%E3%82%B0%E3%81%A8%E3%81%9D%E3%81%AE%E9%96%A2%E9%80%A3%E6%8A%80%E8%A1%93%EF%BC%92.html
[YN12] T. Yamazaki and S. Nagao (2012) A Computational Mechanism for Unified Gain and Timing Control in the Cerebellum. PLoS ONE 7: e33319. https://doi.org/10.1371/journal.pone.0033319
[YI13] T. Yamazaki and J. Igarashi (2013) Real-time cerebellum: a large-scale spiking network model of the cerebellum that runs in real time using a graphics processing unit. Neural Network 47: 103-111. https://doi.org/10.1016/j.neunet.2013.01.019
[Y+17] T. Yamazaki et al. (2017) Real-time simulation of a cat-scale artificial cerebellum on PEZY-SC processors. The International Journal of High Performance Computing Applications. First Published June 6, 2017. http://journals.sagepub.com/doi/full/10.1177/1094342017710705