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