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

Transcript

1. 並列化時代の乱数生成 ~ 軽量 で 再現可能 な枠組みを求めて ~ @abap34 2024/08/31 1 / 78
2. @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
3. 乱数 ⇩ 今日のおはなし 3 / 78
4. 現代の 擬似 乱数 あれこれ 4 / 78
5. [1] イントロダクション ▶︎ 乱数とは？ ▶︎ 真の乱数 ▶︎ 乱数生成と再現性 [2] 擬似乱数生成と検定 ▶︎ 擬似乱数 ▶︎ 線形合同法 ▶︎ スペクトル検定 [3] 並列実行と再現性 ▶︎ 並列処理入門 ▶︎ 再現性が失われる例 [4] Julia における乱数生成 ▶︎ 並列計算でも再現性を保つ工夫 [6] まとめ ▶︎ 参考文献 5 / 78
6. このスライドは、おもに科学技術計算、機械学習 etc... における擬似乱数生成に ついて 「再現性」 の観点からプログラミング言語がどのような工夫をしているの か？というのを書いたものです！ なるべく、擬似乱数や OS などの前提知識なしでも読めるように書かれています が、網羅的・体系的ではないです 例えば、基本的なのに扱っていない生成器がたくさんあります 他にも、乱数検定などの重要なトピックも触れる程度です これらについては、末尾の参考文献などを参照してください。 (+ 各ページで、重要だが飛ばしたものについて参考文献を書いています) ソースコード: https://github.com/abap34/juliatokyo12/tree/main/src このスライドについて 6 / 78
7. [1] イントロダクション ▶︎ 乱数とは？ ▶︎ 真の乱数 ▶︎ 乱数生成と再現性 ▶︎ 擬似乱数 7 / 78
8. [定義: 乱数列 (ゆるふわ)] 規則性がない、 つぎの値が予測できない列 乱数とは？ 8 / 78
9. 世の中は乱数 (のように見えるもの) で溢れている ギャンブル 気象現象 株価 (?) などなど... ⇨ 計算機上で表現するときに乱数が 必要 乱数を使うとき 9 / 78
10. モンテカルロ法 ... 乱数を使った数値計算 モンテカルロ法で の近似値を求め る 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
11. 焼きなまし法 (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
12. 1. ミニバッチのサンプリング 2. パラメータの初期化 3. Dropout 4. Augmentation などなど... 乱数だらけ！ 例) 深層学習での乱数たち 12 / 78
13. 世の中と計算機は乱数に溢れてる ⇨ 実際どうやるか考えよう 乱数生成のことを考えよう 13 / 78
14. ここからは から一様に値を取ってくる方法を考えます ⇨ 一様分布に従う乱数から 所望の分布 • • • • • に従う乱数に変換可能 why? (の方法はこの資料の本筋ではないので割愛します。 興味がある方は、 "逆関数サンプリング", "棄却サンプリング" などを検索するか、参考文献 [6] などを見てみてください) ここからの目標 14 / 78
15. ユーザ目線... たいていのプログラミング言語の 標準ライブラリで 「乱数を生成する機能」 が提供されている 例1) Python >>> import random >>> random.random() 0.2176257387331907 例2) C++ #include <random> 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
16. なんなら OS や CPU命令のレベルでも 提供されている /dev/random /dev/urandom getrandom システムコール RDRAND 命令 例) getrandom(2) を使った乱数生成 #include <sys/random.h> int main() { int rand; getrandom(&rand, sizeof(rand), 0); printf("random number: %d", rand); } https://onecompiler.com/c/42pn5drwx 乱数の作り方 16 / 78
17. [定義 (再掲)] 乱数列 (ゆるふわ) 規則性がない、 つぎの値が予測できない列 ⇩ 決定的な動作しかしない計算機上では、「規則性が ない • • 」列は生成できない. 乱数を計算機で作るには？ 17 / 78
18. アプローチ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
19. これらの方法は (とってくる情報がちゃんとしていれば) 実際予測できない ⇩ 「真の乱数」 (True Random Number, TRNG) とよばれる 乱数生成のアプローチ1. 19 / 78
20. 真の乱数は「予測できない」という乱数の本来の性質を忠実に持っている！ ⇨ 例えば、パスワードの生成などではとても便利・たぶん安全 ➤ cat /dev/random | head -c 10 | base64 UacNyS/a8zKKnw== 真の乱数の困るところ 20 / 78
21. プログラムの 再現性 < 同じプログラムを実行したら、同じ結果になってほしい。 乱数生成と再現性 21 / 78
22. 例) 乱数を生成する関数を受けとって、モンテカルロ法で の近似値を求める関数 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
23. 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
24. 実行のたび、計算結果は変わる。 ( 乱数 • • だからあたりまえ) 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
25. 「パスワードの生成」が再現できる必要はないが、例えば 今のようなモンテカルロ法での計算 ランダムテスト 機械学習モデルの訓練 で、「たまたまうまくいって二度と同じ結果を出せない」のはとても困る (検証がとても大変) 乱数生成と再現性 25 / 78
26. 最近では非常に大量の乱数が必 要な計算が増えている！ (例: モンテカルロ積分は で確率収束するので がたくさんあると気持ちが安心) 乱数生成がボトルネックになる ような計算もたくさん (言語比較のベンチマーク失敗あ るある) 乱数生成のパフォーマンス 26 / 78
27. システムコールをたくさん呼ぶの で、真の乱数生成はパフォーマンス が悪いがち 例) で円周率計算 生成 実行時間 [ms] rand() 2.289 rand(RandomDevice()) 536.367 乱数生成のパフォーマンス 27 / 78
28. rand() ← これ何？ ん？ 28 / 78
29. [2] 擬似乱数生成と検定 ▶︎ 擬似乱数生成 ▶︎ 線形合同法 ▶︎ スペクトル検定 29 / 78
30. アプローチ2: 擬似乱数生成器 (Pseudo Random Number Generator, PRNG) 決定的な動作のみで、 乱数っぽいもの • • • • • • • を作る たいていのプログラミングの標準ライブラリで 「乱数生成器」 として提供される そこで、擬似乱数 暗号論の文脈では、もう少し厳密に擬似乱数の定義について議論ができますが、 (筆者があまり詳しくないので) この資料ではあまり触れません。 興味がある方は https://www.ieice-hbkb.org/files/01/01gun_03hen_11.pdf などを参考にすると良さそうです。 30 / 78
31. 決定的なアルゴリズムのみによって、「乱数」 (のように見えるもの) を作る。 擬似乱数生成 31 / 78
32. この例で伝えたいこと: 擬似乱数生成の基本的な枠組み 乱数の質、検定の考え方 線形合同法の網羅的な説明 一つめの例: 線形合同法 32 / 78
33. どうも乱数っぽい配列 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
34. 実は、 を として計算したときの (を で割って に収めたもの) が先ほどの配列 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
35. [線形合同法] 乱数シード: とパラメータ: のもとで として得られる によって乱数を生成するアルゴリズム 線形合同法 35 / 78
36. 1. 単に決定的なアルゴリズムを実 行しているだけなので、 どの計算機でも 「シード ( ) 」を固定すれば 常に同じ結果が得られる！ 2. とても高速！ (真の乱数比) 擬似乱数生成によるメリット 36 / 78
37. とても 高速 • • (加算、乗算、剰余演算を 1 回ずつ) とても 省メモリ • • • • (生成器の状態としては だけあればよい) とても 実装が楽！ • • • • • (rand(n, x, a, c, M) = [x = (a * x + c) % M for _ in 1:n] だけ) 線形合同法のいいところ 37 / 78
38. 分布が偏ることがある 多次元結晶構造 ができる 線形合同法のやばいところ 38 / 78
39. (見やすさのために) パラメータとして を使って , の中に点を 打ってみると... sample2d!(rng::LCGPRNG) = sample!(rng, 2) sample3d!(rng::LCGPRNG) = sample!(rng, 3) 多次元結晶構造 39 / 78
40. 満遍なく値をとらず、規則的な並びになってしまった！ ⇩　なぜ？ 多次元結晶構造 40 / 78
41. は のみによって一意に定まる. ここから次のことが分かる. [線形合同法の基本的な性質] (1) 同一周期で同じ数は二度と出現しない (2) 得られる列は長さが高々 の周期を持つ 多次元結晶構造の原因 41 / 78
42. さらに、ここから (1) 座標が等しい相異なる 点が得られることはない (2) 座標が等しい相異なる 点が得られることはない も分かる！ ( をよく見ると確かにそうなって いる) 多次元結晶構造の原因 42 / 78
43. ⇩ < 理想的には全ての点から平等に選ばれてほしい！　　が、 -クイーン問題 (の斜めなしver) の制約を満たす点のみから選ばれることになり 一様な分布が得られない可能性がある. 多次元結晶構造の原因 より進んだ結果として、線形合同法を使って に点を打つと、全ての点は 個の並行な超平面に乗ることが知られています。 (Marsagliaの定理) 43 / 78
44. 一方で、 のプロットのように 各超平面間の間隔が小さければ実用 上は問題なさそう？ ⇩ 超平面の間隔を指標として使うこと で良いパラメータかを評価できそう ⇩ スペクトル検定 スペクトル検定 最初に書いたようにこの資料の本筋ではないので、具体的な定義やアルゴリズムは 割愛します。興味があるかたは資料末尾の参考文献 [1] などが詳しいです。 44 / 78
45. (とても時間があるなら、) 自分の使いたい擬似乱数に求め られている性質を考える その性質を乱数検定でチェック (検定にパス いい乱数 は偽だが 検定に落ちる だめな乱数 は真) 乱数検定による乱数の評価 右の画像は、 "RANDU" と呼ばれる、かつて19７０年代ごろに広く使われていた線形 合同法のパラメータを使って三次元空間に点を打ったものです。かなり偏っているの がわかります。 これをスペクトル検定で調べると実際かなり大きい値が得られます。 45 / 78
46. 線形合同法自体の性質を まとめると... 実装はお手軽で 空間計算量も軽い が、質の悪い乱数を生みがち (紹介したもの意外でも致命的な弱点がたくさん... 参考文献の [3] などをどうぞ） ⇨ 本当にメモリの制約が厳しい場合をのぞき、現代で使う場面はあまりない 線形合同法 46 / 78
47. 決定的な動作によって、「乱数っぽいもの」を生成できる それがちゃんと「欲しい乱数」であるかは乱数検定をすることで 検討 できる この例で言いたかったこと 47 / 78
48. 二つめの例: 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
49. 色々なところでデフォルトの乱数生成器として使われている / いた. 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
50. 色々なところでデフォルトの乱数生成器として使われている / いた 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
51. 疑問: なぜ Julia は Mersenne Twister をやめたのか？ ⇨ ヒントは、 並列処理 • • • • Mersenne Twister ? 51 / 78
52. [3] 並列実行と再現性 ▶︎ 並列計算入門 ▶︎ 再現性が失われる例 52 / 78
53. 重要な事実: 現代の科学技術計算は 並列化 抜きには考えられない 並列計算 53 / 78
54. [並列処理 (Parallel Processing)] 同時点 で複数の計算を行うこと とてもとても単純に考えれば、 2つのこ とを同時にやれば2倍速くなるので、処理 を早く終わらすことを目指せる 並列計算入門 逆に、特定の期間で複数の計算を (実際の実行形態はともかく) 行うことを並行処理 (Concurrent Processing) と呼び ます。 54 / 78
55. プロセス (Process) 実行中のプログラムのインスタンス OS は各プロセスにメモリ空間を割 り当てたりして、各プロセスが独立 して動作するようにしている プロセスとスレッド 55 / 78
56. スレッド (Thread) 「軽量プロセス」とも 同じプロセスから作られた他の スレッドとメモリ空間などを共 有している プロセスとスレッド 56 / 78
57. 1. マルチプロセス (Multiprocessing) プロセスを複数立ち上げて、それぞれのプロセスで計算を行う ので、メモリ空間は独立で、競合しづらい 逆に、データのやりとりのコストが大きい プロセスの作成自体も割と高コスト 2. マルチスレッド (Multithreading) プロセス内で複数のスレッドを立ち上げて、それぞれのスレッドで計算を行う メモリ空間は共有なので、データをメモリに置くだけで「通信」できる 逆に、競合が起きて大変なことになる可能性がある プロセスよりは作成が軽い 並列計算入門 57 / 78
58. 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
59. で計算してみると、 ➤ 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
60. 約3倍の高速化！ マルチプロセス シングルプロセス 実際の実行時間 12.71 secs 37.31 secs ユーザー時間 78.52 secs 37.32 secs 並列計算の実行 1. マルチプロセス 60 / 78
61. 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
62. ➤ 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
63. 並列計算をするときの再現性を考える. 並列計算と再現性 63 / 78
64. プログラムの 再現性 < 同じプログラムを実行したら、同じ結果になってほしい。 復習: 再現性 64 / 78
65. 擬似乱数を使えば、 シード を固定することで、再現性を保った上で乱数が使えた. 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
66. 状況: マルチスレッドでモンテカルロ法による円周率の高速化を考えてみる. # 内部的に何かしらの状態をもつ、擬似乱数生成器 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
67. 復習: 各スレッドではメモリ空間は共通だった。 なので、特に何もしなくても「グローバルに一つの乱数生成器を参照する」 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
68. 復習: 各スレッドではメモリ空間は共通だった。 なので、特に何もしないと 「グローバルに一つの乱数生成器を参照する」 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
69. それぞれのスレッドで乱数を要求すると、 C D [1, 2] [3, 4, 5] [1, 4] [2, 3, 5] ... ⇨ シードを固定することによって 「乱数生成器が出すもの • • • • • • • • • • 」 は変わらないが、各スレッドで得られるものは変わ ってしまう！ 再現性が失われる例 69 / 78
70. 乱数を使った計算を並列にやるのは頻出のパターン 機械学習などまさにそう！ 安全に、パフォーマンスよくやる需要がとてもある ⇨ 実は 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
71. [4] Julia における乱数生成 ▶︎ 並列計算でも再現性を保つ工夫 71 / 78
72. 共通の乱数生成器を使うと、乱数の要求タ イミング次第で得られる乱数が異なる. ⇩ 新しいスレッドで計算を行うとき、 各スレッド用の乱数生成器を決定論的 なアルゴリズムで作成して (task.c)、 そこから取り出す ことで、 タイミン グに依らず同じ結果になる！ 並列計算でも再現性を保つ工夫 　 この「作成」の方法については、参考文献 [2] などを参照してください. 72 / 78
73. 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
74. 内部的には、 xoshiro256++ というアルゴリズムを使っている. (Xorshift という著名なアルゴリズムの派生) julia> copy(Random.default_rng()) Xoshiro(0xf41c78cbc898a2cf, 0xbe2efdab9d7feb15, 0x84443883cae5402c, 0xdd254c2b8a0684b1, 0xd003a1925dc045c1) ⇨ Mersenner Twister と比較して、 圧倒的に内部状態が小さい。 why? Task 生成レベルまで介入するので、パフォーマンスに非常に注意する必要があ る！ xoshiro256++ それはそれとして、乱数の質なども含めて Mersenner Twister はどちらかというと古くなりつつあるようです。詳しくは参考文献の [4] や [5] を見てください. 74 / 78
75. s1, s2, s3, s4 だけで OK. xoshiro256++ 75 / 78
76. 言語レベルで対応しているおかげで、特に意識しなくても　 並列計算でも再現性を保つことができる！ Julia の並列の乱数生成まとめ 76 / 78
77. 擬似乱数生成器 (PRNG) を使うことで高速に再現性を持った「乱数」を生成で き、さまざまな計算に活用できる。 さまざまな擬似乱数生成器が存在して、乱数検定などでその質を考察できる。 並列計算をするときは再現性に注意が必要。 Julia を使うことでそこまで意識せ ずとも再現性を持った乱数生成ができる! (おまけ...この置き換えまわりの Issue, PR の議論はとても参考になるし面白いのでおすすめです: https://github.com/JuliaLang/julia/pull/34852 からだいたい辿れるはず) 全体のまとめ 77 / 78
78. 参考文献 [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