「並列化時代の乱数生成」

Slide 1

Slide 1 text

並列化時代の乱数生成 ~ 軽量で再現可能な枠組みを求めて ~ @abap34 2024/08/31 1 / 78

Slide 2

Slide 2 text

@abap34 @abap34 https://abap34.com 所属東京工業大学情報理工学院情報工学系 B3 趣味個人開発機械学習最近は自炊再ブーム Most Used Languages Python 835,745 Bytes Julia 726,183 Bytes TeX 372,733 Bytes C++ 210,563 Bytes C 92,641 Bytes JavaScript 62,902 Bytes Racket 41,235 Bytes Java 35,466 Bytes Go 15,696 Bytes TypeScript 14,456 Bytes 自己紹介 2 / 78

Slide 3

Slide 3 text

乱数 ⇩ 今日のおはなし 3 / 78

Slide 4

Slide 4 text

現代の擬似乱数あれこれ 4 / 78

Slide 5

Slide 5 text

[1] イントロダクション ▶︎ 乱数とは？ ▶︎ 真の乱数 ▶︎ 乱数生成と再現性 [2] 擬似乱数生成と検定 ▶︎ 擬似乱数 ▶︎ 線形合同法 ▶︎ スペクトル検定 [3] 並列実行と再現性 ▶︎ 並列処理入門 ▶︎ 再現性が失われる例 [4] Julia における乱数生成 ▶︎ 並列計算でも再現性を保つ工夫 [6] まとめ ▶︎ 参考文献 5 / 78

Slide 6

Slide 6 text

このスライドは、おもに科学技術計算、機械学習 etc... における擬似乱数生成について「再現性」の観点からプログラミング言語がどのような工夫をしているのか？というのを書いたものです！なるべく、擬似乱数や OS などの前提知識なしでも読めるように書かれていますが、網羅的・体系的ではないです例えば、基本的なのに扱っていない生成器がたくさんあります他にも、乱数検定などの重要なトピックも触れる程度ですこれらについては、末尾の参考文献などを参照してください。 (+ 各ページで、重要だが飛ばしたものについて参考文献を書いています) ソースコード: https://github.com/abap34/juliatokyo12/tree/main/src このスライドについて 6 / 78

Slide 7

Slide 7 text

[1] イントロダクション ▶︎ 乱数とは？ ▶︎ 真の乱数 ▶︎ 乱数生成と再現性 ▶︎ 擬似乱数 7 / 78

Slide 8

Slide 8 text

[定義: 乱数列 (ゆるふわ)] 規則性がない、つぎの値が予測できない列乱数とは？ 8 / 78

Slide 9

Slide 9 text

世の中は乱数 (のように見えるもの) で溢れているギャンブル気象現象株価 (?) などなど... ⇨ 計算機上で表現するときに乱数が必要乱数を使うとき 9 / 78

Slide 10

Slide 10 text

モンテカルロ法 ... 乱数を使った数値計算モンテカルロ法での近似値を求める julia> monte_carlo_pi(n) = 4 * mean(rand()^2 + rand()^2 < 1 for _ in 1:n) monte_carlo_pi (generic function with 1 method) julia> monte_carlo_pi(10^9) 3.141590948 乱数を使ったさまざまな計算例: モンテカルロ法 10 / 78

Slide 11

Slide 11 text

焼きなまし法 (Simulated Annealing) による最適化 function anealing(f::Function, init::Float64) xᵢ₋₁ = init fxᵢ₋₁ = f(xᵢ₋₁) while T > ε T *= α xᵢ = neighbor(xᵢ₋₁, T) fxᵢ = f(xᵢ) if fxᵢ < fxᵢ₋₁ || random_accept(fxᵢ, fxᵢ₋₁, T) xᵢ₋₁, fxᵢ₋₁ = xᵢ, fxᵢ end end return xᵢ₋₁, fxᵢ₋₁ end 例) 焼きなまし法 https://ja.wikipedia.org/wiki/焼きなまし法 11 / 78

Slide 12

Slide 12 text

1. ミニバッチのサンプリング 2. パラメータの初期化 3. Dropout 4. Augmentation などなど... 乱数だらけ！例) 深層学習での乱数たち 12 / 78

Slide 13

Slide 13 text

世の中と計算機は乱数に溢れてる ⇨ 実際どうやるか考えよう乱数生成のことを考えよう 13 / 78

Slide 14

Slide 14 text

ここからはから一様に値を取ってくる方法を考えます ⇨ 一様分布に従う乱数から所望の分布 ● ● ● ● ● に従う乱数に変換可能 why? (の方法はこの資料の本筋ではないので割愛します。興味がある方は、 "逆関数サンプリング", "棄却サンプリング" などを検索するか、参考文献 [6] などを見てみてください) ここからの目標 14 / 78

Slide 15

Slide 15 text

ユーザ目線... たいていのプログラミング言語の標準ライブラリで「乱数を生成する機能」が提供されている例1) Python >>> import random >>> random.random() 0.2176257387331907 例2) C++ #include int main() { std::mt19937 engine(34); std::uniform_real_distribution<> dist(0, 1); std::cout << dist(engine) << std::endl; } 例3) Julia julia> rand() 0.028646780286208817 乱数の作り方 15 / 78

Slide 16

Slide 16 text

なんなら OS や CPU命令のレベルでも提供されている /dev/random /dev/urandom getrandom システムコール RDRAND 命令例) getrandom(2) を使った乱数生成 #include int main() { int rand; getrandom(&rand, sizeof(rand), 0); printf("random number: %d", rand); } https://onecompiler.com/c/42pn5drwx 乱数の作り方 16 / 78

Slide 17

Slide 17 text

[定義 (再掲)] 乱数列 (ゆるふわ) 規則性がない、つぎの値が予測できない列 ⇩ 決定的な動作しかしない計算機上では、「規則性がない ● ● 」列は生成できない. 乱数を計算機で作るには？ 17 / 78

Slide 18

Slide 18 text

アプローチ1: 外から持ってくる. 例1) /dev/random : キーボード入力のような環境からの入力をもとにした予測困難な情報を貯めておいて、 (エントロピープール) そこから計算をして乱数を生成する例2) BRNG (Banana Random Number Generator): バナナに含まれる放射性カリウムの崩壊を観測して乱数を生成する乱数生成のアプローチ一部の OS では /dev/random でも擬似乱数 (後述) を使っていたようです。(例えば FreeBSD は動作環境が特定されやすいので環境ノイズが予測されやすく、よく設計された擬似乱数生成器の方が安全という判断によるものらしいです。)　さらに Linux でも /dev/random 周りではいろいろと変化があり、多くの情報が古くなっていそうです。例えば Linux Kernel v5.6 以降では /dev/random はエントロピープールが枯渇してもブロックされなくなりました。情報源として、 Linux ではネットワークの情報を使ったり使わなかったりとかいろいろと変化や議論があったみたいです。例えば、 2008年の記事ですが以下のようなものがあります: https://atmarkit.itmedia.co.jp/flinux/rensai/watch2008/watch01b.html 18 / 78

Slide 19

Slide 19 text

これらの方法は (とってくる情報がちゃんとしていれば) 実際予測できない ⇩ 「真の乱数」 (True Random Number, TRNG) とよばれる乱数生成のアプローチ1. 19 / 78

Slide 20

Slide 20 text

真の乱数は「予測できない」という乱数の本来の性質を忠実に持っている！ ⇨ 例えば、パスワードの生成などではとても便利・たぶん安全 ➤ cat /dev/random | head -c 10 | base64 UacNyS/a8zKKnw== 真の乱数の困るところ 20 / 78

Slide 21

Slide 21 text

プログラムの再現性 < 同じプログラムを実行したら、同じ結果になってほしい。乱数生成と再現性 21 / 78

Slide 22

Slide 22 text

例) 乱数を生成する関数を受けとって、モンテカルロ法での近似値を求める関数 using Random function monte_carlo_pi(sampler, N) inside = 0 for i in 1:N x = sampler() y = sampler() if x^2 + y^2 < 1 inside += 1 end end return (4 * inside) / N end 乱数生成と再現性 C言語版: https://onecompiler.com/c/42prbqnjz 22 / 78

Slide 23

Slide 23 text

rng = RandomDevice() として /dev/random から乱数を生成することができる julia> rng = RandomDevice() RandomDevice() julia> rand(rng) 0.3809207951383946 julia> rand(rng) # 引数で RNG を指定すること以外は見た目はまったく rand() と同じ 0.5802374206262335 乱数生成と再現性内部的には Libc.getrandom! が呼ばれます。 Libc.getrandom! は　https://docs.libuv.org/en/v1.x/misc.html#c.uv_random　を呼んでいて、これが getrandom(2) を使うか、直接 /dev/random を読みます。 23 / 78

Slide 24

Slide 24 text

実行のたび、計算結果は変わる。 ( 乱数 ● ● だからあたりまえ) julia> N = 10^6 1000000 julia> device_rng = RandomDevice() RandomDevice() julia> monte_carlo_pi( () -> rand(device_rng), N ) 3.13876 julia> monte_carlo_pi( () -> rand(device_rng), N ) 3.1438 乱数生成と再現性 24 / 78

Slide 25

Slide 25 text

「パスワードの生成」が再現できる必要はないが、例えば今のようなモンテカルロ法での計算ランダムテスト機械学習モデルの訓練で、「たまたまうまくいって二度と同じ結果を出せない」のはとても困る (検証がとても大変) 乱数生成と再現性 25 / 78

Slide 26

Slide 26 text

最近では非常に大量の乱数が必要な計算が増えている！ (例: モンテカルロ積分はで確率収束するのでがたくさんあると気持ちが安心) 乱数生成がボトルネックになるような計算もたくさん (言語比較のベンチマーク失敗あるある) 乱数生成のパフォーマンス 26 / 78

Slide 27

Slide 27 text

システムコールをたくさん呼ぶので、真の乱数生成はパフォーマンスが悪いがち例) で円周率計算生成実行時間 [ms] rand() 2.289 rand(RandomDevice()) 536.367 乱数生成のパフォーマンス 27 / 78

Slide 28

Slide 28 text

rand() ← これ何？ん？ 28 / 78

Slide 29

Slide 29 text

[2] 擬似乱数生成と検定 ▶︎ 擬似乱数生成 ▶︎ 線形合同法 ▶︎ スペクトル検定 29 / 78

Slide 30

Slide 30 text

アプローチ2: 擬似乱数生成器 (Pseudo Random Number Generator, PRNG) 決定的な動作のみで、乱数っぽいもの ● ● ● ● ● ● ● を作るたいていのプログラミングの標準ライブラリで「乱数生成器」として提供されるそこで、擬似乱数暗号論の文脈では、もう少し厳密に擬似乱数の定義について議論ができますが、 (筆者があまり詳しくないので) この資料ではあまり触れません。興味がある方は https://www.ieice-hbkb.org/files/01/01gun_03hen_11.pdf などを参考にすると良さそうです。 30 / 78

Slide 31

Slide 31 text

決定的なアルゴリズムのみによって、「乱数」 (のように見えるもの) を作る。擬似乱数生成 31 / 78

Slide 32

Slide 32 text

この例で伝えたいこと: 擬似乱数生成の基本的な枠組み乱数の質、検定の考え方線形合同法の網羅的な説明一つめの例: 線形合同法 32 / 78

Slide 33

Slide 33 text

どうも乱数っぽい配列 100-element Vector{Float64}: 0.24924475536681712 0.36249492410570383 0.0996150195132941 0.9264233387075365 0.043930134968832135 ⋮ 0.643334295367822 0.2540650968439877 0.9413922114763409 0.10687562916427851 0.392702643526718 0.603784283157438 擬似乱数生成 33 / 78

Slide 34

Slide 34 text

実は、をとして計算したときの (をで割ってに収めたもの) が先ほどの配列 100-element Vector{Float64}: 0.24924475536681712 0.36249492410570383 0.0996150195132941 0.9264233387075365 0.043930134968832135 ⋮ 0.643334295367822 0.2540650968439877 0.9413922114763409 0.10687562916427851 0.392702643526718 0.603784283157438 擬似乱数生成パラメータは C++ の std::minstd_rand と同じものを使っています。 34 / 78

Slide 35

Slide 35 text

[線形合同法] 乱数シード: とパラメータ: のもとでとして得られるによって乱数を生成するアルゴリズム線形合同法 35 / 78

Slide 36

Slide 36 text

1. 単に決定的なアルゴリズムを実行しているだけなので、どの計算機でも「シード ( ) 」を固定すれば常に同じ結果が得られる！ 2. とても高速！ (真の乱数比) 擬似乱数生成によるメリット 36 / 78

Slide 37

Slide 37 text

とても高速 ● ● (加算、乗算、剰余演算を 1 回ずつ) とても省メモリ ● ● ● ● (生成器の状態としてはだけあればよい) とても実装が楽！ ● ● ● ● ● (rand(n, x, a, c, M) = [x = (a * x + c) % M for _ in 1:n] だけ) 線形合同法のいいところ 37 / 78

Slide 38

Slide 38 text

分布が偏ることがある多次元結晶構造ができる線形合同法のやばいところ 38 / 78

Slide 39

Slide 39 text

(見やすさのために) パラメータとしてを使って , の中に点を打ってみると... sample2d!(rng::LCGPRNG) = sample!(rng, 2) sample3d!(rng::LCGPRNG) = sample!(rng, 3) 多次元結晶構造 39 / 78

Slide 40

Slide 40 text

満遍なく値をとらず、規則的な並びになってしまった！ ⇩　なぜ？多次元結晶構造 40 / 78

Slide 41

Slide 41 text

はのみによって一意に定まる. ここから次のことが分かる. [線形合同法の基本的な性質] (1) 同一周期で同じ数は二度と出現しない (2) 得られる列は長さが高々の周期を持つ多次元結晶構造の原因 41 / 78

Slide 42

Slide 42 text

さらに、ここから (1) 座標が等しい相異なる点が得られることはない (2) 座標が等しい相異なる点が得られることはないも分かる！ ( をよく見ると確かにそうなっている) 多次元結晶構造の原因 42 / 78

Slide 43

Slide 43 text

⇩ < 理想的には全ての点から平等に選ばれてほしい！　　が、 -クイーン問題 (の斜めなしver) の制約を満たす点のみから選ばれることになり一様な分布が得られない可能性がある. 多次元結晶構造の原因より進んだ結果として、線形合同法を使ってに点を打つと、全ての点は個の並行な超平面に乗ることが知られています。 (Marsagliaの定理) 43 / 78

Slide 44

Slide 44 text

一方で、のプロットのように各超平面間の間隔が小さければ実用上は問題なさそう？ ⇩ 超平面の間隔を指標として使うことで良いパラメータかを評価できそう ⇩ スペクトル検定スペクトル検定最初に書いたようにこの資料の本筋ではないので、具体的な定義やアルゴリズムは割愛します。興味があるかたは資料末尾の参考文献 [1] などが詳しいです。 44 / 78

Slide 45

Slide 45 text

(とても時間があるなら、) 自分の使いたい擬似乱数に求められている性質を考えるその性質を乱数検定でチェック (検定にパスいい乱数は偽だが検定に落ちるだめな乱数は真) 乱数検定による乱数の評価右の画像は、 "RANDU" と呼ばれる、かつて19７０年代ごろに広く使われていた線形合同法のパラメータを使って三次元空間に点を打ったものです。かなり偏っているのがわかります。これをスペクトル検定で調べると実際かなり大きい値が得られます。 45 / 78

Slide 46

Slide 46 text

線形合同法自体の性質をまとめると... 実装はお手軽で空間計算量も軽いが、質の悪い乱数を生みがち (紹介したもの意外でも致命的な弱点がたくさん... 参考文献の [3] などをどうぞ） ⇨ 本当にメモリの制約が厳しい場合をのぞき、現代で使う場面はあまりない線形合同法 46 / 78

Slide 47

Slide 47 text

決定的な動作によって、「乱数っぽいもの」を生成できるそれがちゃんと「欲しい乱数」であるかは乱数検定をすることで検討できるこの例で言いたかったこと 47 / 78

Slide 48

Slide 48 text

二つめの例: Mersenne Twister Julia 1.6 まではデフォルトに採用されていた、とても実用的・主要な乱数生成器 https://docs.julialang.org/en/v1.6/stdlib/Random/ ↑ Julia の Doc に広島大のページが貼られていて凄い多次元でも均等に分布することが保証されている！周期がとても長い！ ( ) 線形合同法では、高々 SIMD で高速化可能！ (http://www.math.sci.hiroshima-u.ac.jp/m- mat/MT/SFMT/index.html) Mersenne Twister 48 / 78

Slide 49

Slide 49 text

色々なところでデフォルトの乱数生成器として使われている / いた. julia> versioninfo() Julia Version 1.6.7 ... LIBM: libopenlibm LLVM: libLLVM-11.0.1 (ORCJIT, westmere) julia> Random.default_rng() MersenneTwister(0xb64166ab5b12df4dc1df3351babcb816, (0, 1002, 0, 1)) Nobody ever got fired for choosing Mersenne Twister. ─ Chris Wellons (https://nullprogram.com/blog/2017/09/21/) Mersenne Twister 49 / 78

Slide 50

Slide 50 text

色々なところでデフォルトの乱数生成器として使われている / いた julia> versioninfo() Julia Version 1.11.0-rc2 Commit 34c3a63147b (2024-07-29 06:24 UTC) ... LLVM: libLLVM-16.0.6 (ORCJIT, apple-m2) Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores) julia> copy(Random.default_rng()) Xoshiro(0xd61a5db6fa36f012, 0xfe333f52e0297386, 0x303b3ad67aa60728, 0xb67bb83a380206b7, 0xde1ca16d107d0c59) Mersenne Twister ? 50 / 78

Slide 51

Slide 51 text

疑問: なぜ Julia は Mersenne Twister をやめたのか？ ⇨ ヒントは、並列処理 ● ● ● ● Mersenne Twister ? 51 / 78

Slide 52

Slide 52 text

[3] 並列実行と再現性 ▶︎ 並列計算入門 ▶︎ 再現性が失われる例 52 / 78

Slide 53

Slide 53 text

重要な事実: 現代の科学技術計算は並列化抜きには考えられない並列計算 53 / 78

Slide 54

Slide 54 text

[並列処理 (Parallel Processing)] 同時点で複数の計算を行うこととてもとても単純に考えれば、 2つのことを同時にやれば2倍速くなるので、処理を早く終わらすことを目指せる並列計算入門逆に、特定の期間で複数の計算を (実際の実行形態はともかく) 行うことを並行処理 (Concurrent Processing) と呼びます。 54 / 78

Slide 55

Slide 55 text

プロセス (Process) 実行中のプログラムのインスタンス OS は各プロセスにメモリ空間を割り当てたりして、各プロセスが独立して動作するようにしているプロセスとスレッド 55 / 78

Slide 56

Slide 56 text

スレッド (Thread) 「軽量プロセス」とも同じプロセスから作られた他のスレッドとメモリ空間などを共有しているプロセスとスレッド 56 / 78

Slide 57

Slide 57 text

1. マルチプロセス (Multiprocessing) プロセスを複数立ち上げて、それぞれのプロセスで計算を行うので、メモリ空間は独立で、競合しづらい逆に、データのやりとりのコストが大きいプロセスの作成自体も割と高コスト 2. マルチスレッド (Multithreading) プロセス内で複数のスレッドを立ち上げて、それぞれのスレッドで計算を行うメモリ空間は共有なので、データをメモリに置くだけで「通信」できる逆に、競合が起きて大変なことになる可能性があるプロセスよりは作成が軽い並列計算入門 57 / 78

Slide 58

Slide 58 text

Distributed パッケージを使ってマルチプロセスによる並列計算ができる using Distributed using Random N_WORKERS = 8 @everywhere using Statistics @everywhere monte_carlo_pi(n) = 4 * mean(rand()^2 + rand()^2 < 1 for _ in 1:n) function monte_carlo_distibuted(; N, N_WORKERS) tasks = [remotecall(monte_carlo_pi, i, N ÷ N_WORKERS) for i in workers()] return mean(fetch.(tasks)) end println("Estimated: ", monte_carlo_distibuted(N=10^10, N_WORKERS=N_WORKERS)) 並列計算の実行 1. マルチプロセス 58 / 78

Slide 59

Slide 59 text

で計算してみると、 ➤ time julia -p 8 --project="." src/3_parallel/mp_montecarlo.jl Estimated: 3.1415537572 ________________________________________________________ Executed in 12.71 secs fish external usr time 78.52 secs 0.09 millis 78.52 secs sys time 0.95 secs 1.25 millis 0.95 secs ➤ time julia --project="." src/3_parallel/montecarlo.jl Estimated: 3.1415815196 ________________________________________________________ Executed in 37.31 secs fish external usr time 37.32 secs 19.64 millis 37.30 secs sys time 0.33 secs 4.81 millis 0.32 secs 並列計算の実行 1. マルチプロセス 59 / 78

Slide 60

Slide 60 text

約3倍の高速化！マルチプロセスシングルプロセス実際の実行時間 12.71 secs 37.31 secs ユーザー時間 78.52 secs 37.32 secs 並列計算の実行 1. マルチプロセス 60 / 78

Slide 61

Slide 61 text

Base.Threads を使ってマルチスレッドによる並列計算ができる using Base.Threads using Statistics N_THREADS = 8 N = 10^10 monte_carlo_pi(N) = 4 * mean(rand()^2 + rand()^2 < 1 for _ in 1:N) function monte_carlo_threaded(; N, N_THREADS) tasks = [Threads.@spawn monte_carlo_pi(N ÷ N_THREADS) for _ in 1:N_THREADS] return mean(fetch.(tasks)) end println("Estimated: ", monte_carlo_threaded(N=N, N_THREADS=N_THREADS)) 並列計算の実行 2. マルチスレッド 61 / 78

Slide 62

Slide 62 text

➤ time julia --threads 8 --project="." src/3_parallel/mt_montecarlo.jl Estimated: 3.1415765036 ________________________________________________________ Executed in 7.77 secs fish external usr time 50.13 secs 0.10 millis 50.13 secs sys time 0.45 secs 1.34 millis 0.45 secs マルチスレッドマルチプロセスシングルプロセス・スレッド実際の実行時間 7.77 secs 12.71 secs 37.31 secs ユーザー時間 50.13 secs 78.52 secs 37.32 secs 並列計算の実行 2. マルチスレッド 62 / 78

Slide 63

Slide 63 text

並列計算をするときの再現性を考える. 並列計算と再現性 63 / 78

Slide 64

Slide 64 text

プログラムの再現性 < 同じプログラムを実行したら、同じ結果になってほしい。復習: 再現性 64 / 78

Slide 65

Slide 65 text

擬似乱数を使えば、シードを固定することで、再現性を保った上で乱数が使えた. using Random Random.seed!(34) rand(3) # 3-element Vector{Float64}: # 0.5142608579591283 # 0.7836918551121808 # 0.5253260048512097 Random.seed!(34) rand(3) # 3-element Vector{Float64}: # 0.5142608579591283 # 0.7836918551121808 # 0.5253260048512097 復習: 再現性 65 / 78

Slide 66

Slide 66 text

状況: マルチスレッドでモンテカルロ法による円周率の高速化を考えてみる. # 内部的に何かしらの状態をもつ、擬似乱数生成器 struct PRNG state end monte_carlo_pi(rng::PRNG, n)　= 4 * mean(rand(rng)^2 + rand(rng)^2 < 1 for _ in 1:n) monte_carlo_pi(PRNG(), 10^6) # -> 3.1415... ⇩　マルチスレッドで計算並列計算と再現性 66 / 78

Slide 67

Slide 67 text

復習: 各スレッドではメモリ空間は共通だった。なので、特に何もしなくても「グローバルに一つの乱数生成器を参照する」 function monte_carlo_threaded(rng::PRNG; N, N_THREADS) tasks = [Threads.@spawn monte_carlo_pi(rng, N ÷ N_THREADS) for _ in 1:N_THREADS] return mean(fetch.(tasks)) end 並列計算と再現性 67 / 78

Slide 68

Slide 68 text

復習: 各スレッドではメモリ空間は共通だった。なので、特に何もしないと「グローバルに一つの乱数生成器を参照する」 function monte_carlo_threaded(rng::PRNG; N, N_THREADS) tasks = [Threads.@spawn monte_carlo_pi(rng, N ÷ N_THREADS) for _ in 1:N_THREADS] return mean(fetch.(tasks)) end 並列計算と再現性 68 / 78

Slide 69

Slide 69 text

それぞれのスレッドで乱数を要求すると、 C D [1, 2] [3, 4, 5] [1, 4] [2, 3, 5] ... ⇨ シードを固定することによって「乱数生成器が出すもの ● ● ● ● ● ● ● ● ● ● 」は変わらないが、各スレッドで得られるものは変わってしまう！再現性が失われる例 69 / 78

Slide 70

Slide 70 text

乱数を使った計算を並列にやるのは頻出のパターン機械学習などまさにそう！安全に、パフォーマンスよくやる需要がとてもある ⇨ 実は Julia はかなり低いレベルでこれに対応しているそこで < It feels strange to me to make RNGs this low-level, but it might be necessary. 対策の需要 https://github.com/JuliaLang/julia/pull/34852#issuecomment-590670749 70 / 78

Slide 71

Slide 71 text

[4] Julia における乱数生成 ▶︎ 並列計算でも再現性を保つ工夫 71 / 78

Slide 72

Slide 72 text

共通の乱数生成器を使うと、乱数の要求タイミング次第で得られる乱数が異なる. ⇩ 新しいスレッドで計算を行うとき、各スレッド用の乱数生成器を決定論的なアルゴリズムで作成して (task.c)、そこから取り出すことで、タイミングに依らず同じ結果になる！並列計算でも再現性を保つ工夫　この「作成」の方法については、参考文献 [2] などを参照してください. 72 / 78

Slide 73

Slide 73 text

Task のレベルで乱数生成に介入している！！ (乱数生成が今どきどれくらい重要かを示唆していて、めちゃくちゃ面白い) void jl_rng_split(uint64_t dst[JL_RNG_SIZE], uint64_t src[JL_RNG_SIZE]) JL_NOTSAFEPOINT { // load and advance the internal LCG state uint64_t x = src[4]; src[4] = dst[4] = x * 0xd1342543de82ef95 + 1; // high spectrum multiplier from https://arxiv.org/abs/2001.05304 // random xor constants static const uint64_t a[4] = { ... c *= m[i]; c ^= c >> 43; dst[i] = c; } } (https://github.com/JuliaLang/julia/blob/4c2f728a9976a5651acfe2f7eba703e6d0b64562/src/task.c#L1094) 並列計算でも再現性を保つ工夫 73 / 78

Slide 74

Slide 74 text

内部的には、 xoshiro256++ というアルゴリズムを使っている. (Xorshift という著名なアルゴリズムの派生) julia> copy(Random.default_rng()) Xoshiro(0xf41c78cbc898a2cf, 0xbe2efdab9d7feb15, 0x84443883cae5402c, 0xdd254c2b8a0684b1, 0xd003a1925dc045c1) ⇨ Mersenner Twister と比較して、圧倒的に内部状態が小さい。 why? Task 生成レベルまで介入するので、パフォーマンスに非常に注意する必要がある！ xoshiro256++ それはそれとして、乱数の質なども含めて Mersenner Twister はどちらかというと古くなりつつあるようです。詳しくは参考文献の [4] や [5] を見てください. 74 / 78

Slide 75

Slide 75 text

s1, s2, s3, s4 だけで OK. xoshiro256++ 75 / 78

Slide 76

Slide 76 text

言語レベルで対応しているおかげで、特に意識しなくても　並列計算でも再現性を保つことができる！ Julia の並列の乱数生成まとめ 76 / 78

Slide 77

Slide 77 text

擬似乱数生成器 (PRNG) を使うことで高速に再現性を持った「乱数」を生成でき、さまざまな計算に活用できる。さまざまな擬似乱数生成器が存在して、乱数検定などでその質を考察できる。並列計算をするときは再現性に注意が必要。 Julia を使うことでそこまで意識せずとも再現性を持った乱数生成ができる! (おまけ...この置き換えまわりの Issue, PR の議論はとても参考になるし面白いのでおすすめです: https://github.com/JuliaLang/julia/pull/34852 からだいたい辿れるはず) 全体のまとめ 77 / 78

Slide 78

Slide 78 text

参考文献 [1] 伏見正則. 乱数. 筑摩書房, 2023. [2] L’Ecuyer, Pierre, et al. "Random numbers for parallel computers: Requirements and methods, with emphasis on GPUs." Mathematics and Computers in Simulation 135 (2017): 3-17. [3] 松本眞. あなたの使っている乱数、大丈夫？ -危ない標準乱数と、メルセンヌ・ツイスター開発秘話-. http://www.math.sci.hiroshima- u.ac.jp/m-mat/TEACH/ichimura-sho-koen.pdf. [4] Vigna, Sebastiano. "It is high time we let go of the Mersenne Twister." arXiv preprint arXiv:1910.06437 (2019). [5] Chris Wellons. "Finding the Best 64-bit Simulation PRNG". https://nullprogram.com/blog/2017/09/21/ [6] 杉山将. 統計的機械学習: 生成モデルに基づくパターン認識. オーム社, 2009. ISBN 9784274502484. Google Books, https://books.google.co.jp/books?id=yK-IQgAACAAJ. ソースコード https://github.com/abap34/juliatokyo12 Appendix 78 / 78