概要

動機

実装を通してイジングモデルを理解する
GPU における計算手法を確立する

GPU 計算については、最高のパフォーマンスを求めるというよりも、できる限り楽に、そこそこのパフォーマンスを出せる方法を模索したい。

イジングモデルについて

特に参考にさせてもらったサイト

Pythonで粒子のスピンを操ろう！: 磁石とはなんぞや？
Ising Modelを平易に解説してみる
 Gen Kuroki さんの実装

結果

画像クリックか下のアドレスでGIF動画（12MB）へ
https://imgur.com/8yqpnzt

コード（Julia Notebook）
https://gist.github.com/Lirimy/127d6a557ba10d549d17c6687ee9c5c2

１フレームですべてのピクセルを平均１回評価する。
サイズ 640×360, 450 フレーム (30 fps * 15 秒)
計算だけなら 3 秒、動画を出力しても 11 秒で完了する。

サイズ	640x360	1920x1080
計算のみ	2.8 sec	20 sec
計算＋画像出力	7.0 sec	51 sec
動画化 (gif)	3.8 sec	34 sec
動画化 (mp4)	5.8 sec	52 sec
動画容量 (gif)	12 MB	103 MB
動画容量 (mp4)	33 MB	298 MB
GPU 使用率	30 %	50 %

イジングモデル動画化についての知見

当初は　mp4 形式で YouTube にアップロードしようと目論んで、実際やったのだが、再エンコードされてしまって残念な画質になってしまっている。アップロード前の元動画は 15秒 33 MB でおよそ 18 Mbps を必要とするが、 YouTube の推奨ビットレートが 1 Mbps なので、当然の成り行きである。

そもそもイジングモデルではピクセルがランダムにフリップするので、空間的に相関が少なく、いわゆるエンコ殺しの画像になっており、mp4とは相性が悪いようだ。一方で、色的に見ると２状態に対応する２色しか使われていないので、gifのパレットが有効に機能するようだ。gif動画はmp4と同条件で 12 MB であり、およそ　1/3 の容量に収まっている。動画変換も速い。

実装の準備

GPU ライブラリの選択

可視化に優れた Julia を使いたいので、今回は ArrayFire のラッパー ArrayFire.jl を選択した。

高レベルなので書きやすい
バックエンドを CUDA と OpenCL から選べる

今回は見送ったが、他の候補たち

ArrayFire (C++) を Julia から呼び出し
CUDA ラッパー
CUDA Fortran を Julia から呼び出し

ArrayFire.jl の使い方

ライブラリを使う場合、普通は公式のドキュメントを読めばなんとかなるのだが、残念ながらまだドキュメントが整っていないようなのでどうすべきか、という話。

ArrayFire.jl の Supported Functions で欲しい関数が実装されているか探す
wrap.jl で検索して、関数の定義を確認する
ArrayFire 本家で渡すべき値を確認する

今回使った convolve2 を例にすると、

convolve2(signal::AFArray,filter::AFArray,mode::af_conv_mode,domain::af_conv_domain)

が関数の定義で、４つの引数（デフォルトも設定されていない）が必要だとわかる。
convolve2 を見ると、signal, filter のほかに、 convMode と convDomain を設定すればよいとわかる。
以上をまとめて、次のように呼び出せば GPU での２次元畳込みができる。

convolve2(s, kernel, AF_CONV_EXPAND, AF_CONV_SPATIAL)

可視化

状態の経時変化はなんとしても見てみたかった。
既存のライブラリだと Plots.jl がアニメーション機能を有するが、プロットするのに heatmap を使用することになる。しかし、可能ならば１セルを１ピクセルに対応させたいという希望があったので、自分でモジュールを書いてみた。

コードは　Gist においてある。
中身は animation.jl をもとにいじっただけだが、速度的にずっと改善してある（３倍くらい？）。
Plots.jl では、中間ファイルに PNG を使用しており、実は PNG のエンコーディングに時間がかかっているようなので、中間ファイル形式を BMP に変更した。

Julia におけるアニメーションのしくみ

ArrayFire.jl を用いた実装

GPU に合わせたコーディング

GPU でやると極端に遅い処理

条件分岐・個別の処理
インデックス操作（読むのはOK、書くのがダメ）

これらをアルゴリズム上うまく回避して書かねばならない。

他の注意としては、特に GeForce などのコンシューマ向け GPU は Float32 に特化しており、うっかり Float64 にキャストされると遅くなってしまう。

s = rand(AFArray{Float32}, 100, 100)

s / 10   # OK
0.1 * s  # NG

Float32(0.1) * s
0.1f0 * s

隣接するセルの和

カーネルを畳み込むことで実現した。

kernel = AFArray(Float32.([0 1 0; 1 0 1; 0 1 0]))

周期境界条件

行列積で実現した。

畳み込み時のパディングを工夫してみる

比較

比較演算の結果 (Bool) を 1/0 と見て、素直にキャストなしに演算できたのでそれでも良いが、今回は signbit を主に利用してみた。

signbit は sign() のラップで、負のときのみ１、他の場合は０を出力する。特筆すべきは、出力が Float32 なので、一貫して同じ型で計算できることである。が、速度にどの程度影響するかは確かめられていない。

a > b         # Bool
signbit(b-a)  # Float32

逐次処理

イジングモデルの計算の本来の形は

どこか１つのセルに注目して評価する
セルの値を更新して、次のセルへ

であるが、GPU だと素直に実装できない。

そこで、操作を少々変更して、全体のうちランダムな 1/10 のセルの評価・更新を繰り返すようにした。結果を見るとうまくいっていそうな気はするが、問題があれば教えてほしい。

格子模様を描き、色ごとに計算すればよい。
https://qiita.com/Lirimy/items/e25efa016c54d07368c8

メモリ管理

メモリ管理をせずに計算していると、唐突にこんなエラーが出現する。

GPU is out of memory, to avoid this in the future you can:
  @afgc function f(input)   # free all temporary variables inside the function scope
  @afgc a = b + c           # free all temporary arrays inside the assignment scope
  @afgc a .= b + c          # replace with a new value, free the old and all temps
  afgc(threshold)           # garbage collect after GPU memory usage reaches threshold
  finalize(array)           # manually free GPU memory

しかし、何も考えず頻繁にガベージコレクションを呼び出しているとパフォーマンスが大きく低下する。
下手なことをやっていたときは　gc (CPU側) が処理時間の 99% を占めていた。
あまりにもひどいので gc_enable(false) などと禁断の呪文を呼び出していたが、無事解決策が見つかったので、これからやる人は参考にしてほしい。

メインループ内で afgc()

手軽さを求めるならこれが良さそう。

const MB = UInt64(1024*1024);

for i in 1:100 # メインループ

    # 処理

    afgc(3200MB)
end

nvidia-smi -l 1 でメモリ使用量を見て、できるだけ大きいしきい値を設定しよう。
afgc() は CPU/GPU 両方のガベージコレクションを呼び出す。

関数定義に　`@afgc` をつける

この方式だと、 @afgc が　256 回呼ばれるごとに afgc() が呼ばれる。
上の方式でうまくいったので採用しなかったが、やり方を残しておく。

@noinline @afgc function update!(s::AFArray{Float32, 2})

    # s を更新する

    return nothing
end

戻り値はなし
念のため、明示的にスコープをわける

CPU 側の gc

基本的には afgc() 内で呼ばれるので問題はないはずなのだが、画像出力せずにベンチマークしていると、なぜか画像出力したときの数倍遅くなってしまった。よくわからないが、とりあえず gc(false) を入れておけば解決はした。false はフルで gc しないという意味らしい。

まとめ

成果

イジングモデルを大まかに理解できた
GPU の計算手法を確立できた
- GPU の様々な制約への対応
- メモリ管理
２次元配列の可視化手法を確立できた

課題

Julia の gc() の挙動を理解する
そもそもゴミを出さない、または減らす手法を確立する
パッケージの取り扱いを理解する（よく名前が衝突するため）

JuliaGPU / ArrayFire でイジングモデルを計算する

概要

動機