KEMBAR78
Implementing sobol's quasirandom sequence generator | PDF
ALGORITHM 659

Implementing Sobol’s
Quasirandom Sequence
Generator

PAUL BRATLEY and BENNETT L. FOX

ACM Transactions on Mathematical Software

Vol. 14 (1988) 88;

1
読んできた論文

● Ilya M. Sobolさんという方が考えた Quasi-random Sequence Generatorの実装
方法について書かれた論文
○ “Quasi” とついているように、通常の乱数とは少し違い、Low Discrepancyと
いう性質を持つ。
● Sobolさんの論文ではなく、アルゴリズムの実装にフォーカスして解説してくれてい
る論文 (引用数: 961件)
○ Theoretical underpinningsにはあんまりちゃんと触れない。
● 今回はGo言語で実装してGoptunaにも導入
2
Low discrepancy

3
Low discrepancy

4
ハイパーパラメーター最適化を例に考えると、
似たようなパラメーターは、 (経験的に) 似たような
評価値になることが多いので、散らばってくれると
嬉しい。
アルゴリズムの用途

● Optimization

○ BBO Challengeで自分たちも使用 

○ “We employ Quasi Monte Carlo sampling
based on Sobol sequences instead of Latin
hypercube sampling.” 

● Numerical Integration 

○ 疑似乱数(pseudorandom)を使用する
Monte-Carlo Integrationと区別して
”Quasi-Monte Carlo Integration” と呼ばれ
る

5
アルゴリズム

6
アルゴリズム

大きく2つのフェーズに分かれる

1. Direction numbersの生成

2. Sobol sequenceの生成

7
Direction numbersの生成

8
Direction numbers の形式と性質

9
i 1 2 3 4 5 6
mi 1 3 7 5 7 43
vi 0.1 0.11 0.111 0.0101 0.00111 0.101011
Direction numbers の形式と性質

10
i 1 2 3 4 5 6
mi 1 3 7 5 7 43
vi 0.1 0.11 0.111 0.0101 0.00111 0.101011
こういう形式のdirection numbersを生成。
※ 論文上はこのような浮動小数点数です
が、実装の都合上プログラムでは uint32など
で表現されています。
Direction numbers の形式と性質

11
i 1 2 3 4 5 6
mi 1 3 7 5 7 43
vi 0.1 0.11 0.111 0.0101 0.00111 0.101011
direction numbersは、miという値から計
算されるので、miの作り方を知る必要が
あります。
Direction numbers の形式と性質

12
i 1 2 3 4 5 6
mi 1 3 7 5 7 43
vi 0.1 0.11 0.111 0.0101 0.00111 0.101011
miはこのような値を持ちます。
Direction numbers の形式と性質

13
i 1 2 3 4 5 6
mi 1 3 7 5 7 43
vi 0.1 0.11 0.111 0.0101 0.00111 0.101011
miの生成

14
mi のつくりかた ①

15
次のようなPrimitive Polynomial (原子多項式)をφ(2^d-1)/d個用意。

各係数は0 or 1のどちらかをとり、定数項は必ず1。

これを満たしていれば、係数は自由に選べる。

miのつくりかた ②

16
Primitive Polynomialの係数をXORで足し合わせる。

miのつくりかた ②

17
Primitive Polynomialの係数をXORで足し合わせる。

miの計算には m_(i-1) などが必要。
なのでm0やm1といった初期値は
なにか適当に決める必要があります。
→ Initial direction numbersと呼ばれる
これも 0<mi<2^i を満たす奇数。すべて 1も
可。
mi の計算例

18
Sobol sequenceの生成

19
Sobol sequenceの生成

次の演算をするだけ。

20
Gray code による改善

● Antonov and Saleevさんらによる改良

● Direction numbersからSobol sequenceを生成する際に、Gray codeを利用する

21
The (binary-reflected) Gray code

22
i (decimal) i (binary) i / 2 (binary) Gray (binary) Gray (decimal)
0 0000 0000 0000 0
1 0001 0000 0001 1
2 0010 0001 0011 3
3 0011 0001 0010 2
4 0100 0010 0110 6
5 0101 0010 0111 7
6 0110 0011 0101 5
7 0111 0011 0100 4
Gray codeは

同じ色で囲われた範囲の 

中で順序を入れ替える特性
をもつ。



また1 bitしか変化がない。 

The (binary-reflected) Gray code

23
i (decimal) i (binary) i / 2 (binary) Gray (binary) Gray (decimal)
0 0000 0000 0000 0
1 0001 0000 0001 1
2 0010 0001 0011 3
3 0011 0001 0010 2
4 0100 0010 0110 6
5 0101 0010 0111 7
6 0110 0011 0101 5
7 0111 0011 0100 4
1 bit 右にシフト
Gray codeは

同じ色で囲われた範囲の 

中で順序を入れ替える特性
をもつ。



また1 bitしか変化がない。 

The (binary-reflected) Gray code

24
i (decimal) i (binary) i / 2 (binary) Gray (binary) Gray (decimal)
0 0000 0000 0000 0
1 0001 0000 0001 1
2 0010 0001 0011 3
3 0011 0001 0010 2
4 0100 0010 0110 6
5 0101 0010 0111 7
6 0110 0011 0101 5
7 0111 0011 0100 4
i と i/2 のXOR
Gray codeは

同じ色で囲われた範囲の 

中で順序を入れ替える特性
をもつ。



また1 bitしか変化がない。 

Gray codeは

同じ色で囲われた範囲の 

中で順序を入れ替える特性
をもつ。



また1 bitしか変化がない。 

The (binary-reflected) Gray code

25
i (decimal) i (binary) i / 2 (binary) Gray (binary) Gray (decimal)
0 0000 0000 0000 0
1 0001 0000 0001 1
2 0010 0001 0011 3
3 0011 0001 0010 2
4 0100 0010 0110 6
5 0101 0010 0111 7
6 0110 0011 0101 5
7 0111 0011 0100 4
同じ色で囲われた範囲の中で
順序が入れ替わっている
Gray codeは

同じ色で囲われた範囲の 

中で順序を入れ替える特性
をもつ。



また1 bitしか変化がない。 

The (binary-reflected) Gray code

26
i (decimal) i (binary) i / 2 (binary) Gray (binary) Gray (decimal)
0 0000 0000 0000 0
1 0001 0000 0001 1
2 0010 0001 0011 3
3 0011 0001 0010 2
4 0100 0010 0110 6
5 0101 0010 0111 7
6 0110 0011 0101 5
7 0111 0011 0100 4
1 bitずつしか変化しない。
変化するbitの位置は、
rightmost zero-bit in the
binary representation of i
(太字にした箇所)
The (binary-reflected) Gray code

27
先程の性質を使うとXOR演算回数がN回→1回 

The (binary-reflected) Gray code

28
先程の性質を使うとXOR演算回数がN回→1回 

Computational Complexity
実際には O(1) になるわけでは
なく、rightmost zero-bitの特定
に O(log n) かかる。
余談: rightmost zero-bitの特定方法

● scipyの実装は改善の余地があったので修正してPR作成済み

○ https://github.com/scipy/scipy/pull/13514

● popcntというCPU命令を使うことで O(1) に近い高速化も可

○ ※ gccやclangのコンパイラ拡張として提供されているので、書き分けが必要。
また (v & -v) - 1の演算が直感的ではなく、可読性が悪いのでscipyへの導入
は諦めました (Goptunaには適用して高速化を確認)。

29
余談: rightmost zero-bitの特定方法

● scipyの実装は改善の余地があったので修正してPR作成済み

○ https://github.com/scipy/scipy/pull/13514

● popcntというCPU命令を使うことで O(1) に近い高速化も可

○ ※ gccやclangのコンパイラ拡張として提供されているので、書き分けが必要。
また (v & -v) - 1の演算が直感的ではなく、可読性が悪いのでscipyへの導入
は諦めました (Goptunaには適用して高速化を確認)。

30
x & (-x) 0101100 → 0000100
(x & (-x))-1 0101100 → 0000011
popcnt(...) 0101100 → 2
アルゴリズムの流れ (ふりかえり)

1. Direction numbersの生成

a. Primitive Polynomial (原始多項式) を用意

b. Initial direction numbers を用意

c. 多項式の各係数をXORで足し合わせる

2. Sobol sequenceの生成

a. 生成回数 i のbinary表現からGray codeを計算

b. Gray codeを導入すること、計算量を減らすことができる。

31
その他のテクニック

32
多次元のSobol sequenceの生成

● 2次元以上のSobol sequenceを生成する際には、各次元ごとにprimitive polynomial
とinitial direction numbersを用意して、先程と同じ手順を繰り返す。

● 注意点として、initial direction numbersは何でもいいといいつつ、多次元において
はよい値(low discrepancyになるもの)とそうでないものがある。

○ そのため良いinitial direction numbersを調べている人たちもいる。

33
Joe & Kuo

34
PyTorch, Scipyもここで配布されているprimitive polynomialsとinitial direction numbersを
埋め込んで使用する。21201次元まで対応

https://web.maths.unsw.edu.au/~fkuo/sobol/
Skipping initial points

Joe & Kuoさんらが公開しているノートにかかれていたテクニック



> Sobol’ sequence tends to perform better

> if an initial portion of the sequence is dropped:

> the number of points skipped is the largest

> power of 2 smaller than number of points to be

> used. However, we are less persuaded by such 

> recommendation ourselves.

35
まとめ

36
まとめ

● Sobol sequenceは、Low discrepancyという性質をもつ乱数のようなものがある。

○ 主な用途はOptimizationやNumerical Integration。

● 今回はアルゴリズムの解説論文を読んで紹介。

○ Theoretical underpinningsは触れられていないので、どうしてLow discrepancy
になるのかは、Sobolさんの論文を参照

37

Implementing sobol's quasirandom sequence generator

  • 1.
    ALGORITHM 659
 Implementing Sobol’s QuasirandomSequence Generator
 PAUL BRATLEY and BENNETT L. FOX
 ACM Transactions on Mathematical Software
 Vol. 14 (1988) 88;
 1
  • 2.
    読んできた論文
 ● Ilya M.Sobolさんという方が考えた Quasi-random Sequence Generatorの実装 方法について書かれた論文 ○ “Quasi” とついているように、通常の乱数とは少し違い、Low Discrepancyと いう性質を持つ。 ● Sobolさんの論文ではなく、アルゴリズムの実装にフォーカスして解説してくれてい る論文 (引用数: 961件) ○ Theoretical underpinningsにはあんまりちゃんと触れない。 ● 今回はGo言語で実装してGoptunaにも導入 2
  • 3.
  • 4.
    Low discrepancy
 4 ハイパーパラメーター最適化を例に考えると、 似たようなパラメーターは、 (経験的に)似たような 評価値になることが多いので、散らばってくれると 嬉しい。
  • 5.
    アルゴリズムの用途
 ● Optimization
 ○ BBOChallengeで自分たちも使用 
 ○ “We employ Quasi Monte Carlo sampling based on Sobol sequences instead of Latin hypercube sampling.” 
 ● Numerical Integration 
 ○ 疑似乱数(pseudorandom)を使用する Monte-Carlo Integrationと区別して ”Quasi-Monte Carlo Integration” と呼ばれ る
 5
  • 6.
  • 7.
  • 8.
  • 9.
    Direction numbers の形式と性質
 9 i1 2 3 4 5 6 mi 1 3 7 5 7 43 vi 0.1 0.11 0.111 0.0101 0.00111 0.101011
  • 10.
    Direction numbers の形式と性質
 10 i1 2 3 4 5 6 mi 1 3 7 5 7 43 vi 0.1 0.11 0.111 0.0101 0.00111 0.101011 こういう形式のdirection numbersを生成。 ※ 論文上はこのような浮動小数点数です が、実装の都合上プログラムでは uint32など で表現されています。
  • 11.
    Direction numbers の形式と性質
 11 i1 2 3 4 5 6 mi 1 3 7 5 7 43 vi 0.1 0.11 0.111 0.0101 0.00111 0.101011 direction numbersは、miという値から計 算されるので、miの作り方を知る必要が あります。
  • 12.
    Direction numbers の形式と性質
 12 i1 2 3 4 5 6 mi 1 3 7 5 7 43 vi 0.1 0.11 0.111 0.0101 0.00111 0.101011 miはこのような値を持ちます。
  • 13.
    Direction numbers の形式と性質
 13 i1 2 3 4 5 6 mi 1 3 7 5 7 43 vi 0.1 0.11 0.111 0.0101 0.00111 0.101011
  • 14.
  • 15.
    mi のつくりかた ①
 15 次のようなPrimitivePolynomial (原子多項式)をφ(2^d-1)/d個用意。
 各係数は0 or 1のどちらかをとり、定数項は必ず1。
 これを満たしていれば、係数は自由に選べる。

  • 16.
  • 17.
    miのつくりかた ②
 17 Primitive Polynomialの係数をXORで足し合わせる。
 miの計算にはm_(i-1) などが必要。 なのでm0やm1といった初期値は なにか適当に決める必要があります。 → Initial direction numbersと呼ばれる これも 0<mi<2^i を満たす奇数。すべて 1も 可。
  • 18.
  • 19.
  • 20.
  • 21.
    Gray code による改善
 ●Antonov and Saleevさんらによる改良
 ● Direction numbersからSobol sequenceを生成する際に、Gray codeを利用する
 21
  • 22.
    The (binary-reflected) Graycode
 22 i (decimal) i (binary) i / 2 (binary) Gray (binary) Gray (decimal) 0 0000 0000 0000 0 1 0001 0000 0001 1 2 0010 0001 0011 3 3 0011 0001 0010 2 4 0100 0010 0110 6 5 0101 0010 0111 7 6 0110 0011 0101 5 7 0111 0011 0100 4 Gray codeは
 同じ色で囲われた範囲の 
 中で順序を入れ替える特性 をもつ。
 
 また1 bitしか変化がない。 

  • 23.
    The (binary-reflected) Graycode
 23 i (decimal) i (binary) i / 2 (binary) Gray (binary) Gray (decimal) 0 0000 0000 0000 0 1 0001 0000 0001 1 2 0010 0001 0011 3 3 0011 0001 0010 2 4 0100 0010 0110 6 5 0101 0010 0111 7 6 0110 0011 0101 5 7 0111 0011 0100 4 1 bit 右にシフト Gray codeは
 同じ色で囲われた範囲の 
 中で順序を入れ替える特性 をもつ。
 
 また1 bitしか変化がない。 

  • 24.
    The (binary-reflected) Graycode
 24 i (decimal) i (binary) i / 2 (binary) Gray (binary) Gray (decimal) 0 0000 0000 0000 0 1 0001 0000 0001 1 2 0010 0001 0011 3 3 0011 0001 0010 2 4 0100 0010 0110 6 5 0101 0010 0111 7 6 0110 0011 0101 5 7 0111 0011 0100 4 i と i/2 のXOR Gray codeは
 同じ色で囲われた範囲の 
 中で順序を入れ替える特性 をもつ。
 
 また1 bitしか変化がない。 

  • 25.
    Gray codeは
 同じ色で囲われた範囲の 
 中で順序を入れ替える特性 をもつ。
 
 また1bitしか変化がない。 
 The (binary-reflected) Gray code
 25 i (decimal) i (binary) i / 2 (binary) Gray (binary) Gray (decimal) 0 0000 0000 0000 0 1 0001 0000 0001 1 2 0010 0001 0011 3 3 0011 0001 0010 2 4 0100 0010 0110 6 5 0101 0010 0111 7 6 0110 0011 0101 5 7 0111 0011 0100 4 同じ色で囲われた範囲の中で 順序が入れ替わっている
  • 26.
    Gray codeは
 同じ色で囲われた範囲の 
 中で順序を入れ替える特性 をもつ。
 
 また1bitしか変化がない。 
 The (binary-reflected) Gray code
 26 i (decimal) i (binary) i / 2 (binary) Gray (binary) Gray (decimal) 0 0000 0000 0000 0 1 0001 0000 0001 1 2 0010 0001 0011 3 3 0011 0001 0010 2 4 0100 0010 0110 6 5 0101 0010 0111 7 6 0110 0011 0101 5 7 0111 0011 0100 4 1 bitずつしか変化しない。 変化するbitの位置は、 rightmost zero-bit in the binary representation of i (太字にした箇所)
  • 27.
    The (binary-reflected) Graycode
 27 先程の性質を使うとXOR演算回数がN回→1回 

  • 28.
    The (binary-reflected) Graycode
 28 先程の性質を使うとXOR演算回数がN回→1回 
 Computational Complexity 実際には O(1) になるわけでは なく、rightmost zero-bitの特定 に O(log n) かかる。
  • 29.
    余談: rightmost zero-bitの特定方法
 ●scipyの実装は改善の余地があったので修正してPR作成済み
 ○ https://github.com/scipy/scipy/pull/13514
 ● popcntというCPU命令を使うことで O(1) に近い高速化も可
 ○ ※ gccやclangのコンパイラ拡張として提供されているので、書き分けが必要。 また (v & -v) - 1の演算が直感的ではなく、可読性が悪いのでscipyへの導入 は諦めました (Goptunaには適用して高速化を確認)。
 29
  • 30.
    余談: rightmost zero-bitの特定方法
 ●scipyの実装は改善の余地があったので修正してPR作成済み
 ○ https://github.com/scipy/scipy/pull/13514
 ● popcntというCPU命令を使うことで O(1) に近い高速化も可
 ○ ※ gccやclangのコンパイラ拡張として提供されているので、書き分けが必要。 また (v & -v) - 1の演算が直感的ではなく、可読性が悪いのでscipyへの導入 は諦めました (Goptunaには適用して高速化を確認)。
 30 x & (-x) 0101100 → 0000100 (x & (-x))-1 0101100 → 0000011 popcnt(...) 0101100 → 2
  • 31.
    アルゴリズムの流れ (ふりかえり)
 1. Directionnumbersの生成
 a. Primitive Polynomial (原始多項式) を用意
 b. Initial direction numbers を用意
 c. 多項式の各係数をXORで足し合わせる
 2. Sobol sequenceの生成
 a. 生成回数 i のbinary表現からGray codeを計算
 b. Gray codeを導入すること、計算量を減らすことができる。
 31
  • 32.
  • 33.
    多次元のSobol sequenceの生成
 ● 2次元以上のSobolsequenceを生成する際には、各次元ごとにprimitive polynomial とinitial direction numbersを用意して、先程と同じ手順を繰り返す。
 ● 注意点として、initial direction numbersは何でもいいといいつつ、多次元において はよい値(low discrepancyになるもの)とそうでないものがある。
 ○ そのため良いinitial direction numbersを調べている人たちもいる。
 33
  • 34.
    Joe & Kuo
 34 PyTorch,Scipyもここで配布されているprimitive polynomialsとinitial direction numbersを 埋め込んで使用する。21201次元まで対応
 https://web.maths.unsw.edu.au/~fkuo/sobol/
  • 35.
    Skipping initial points
 Joe& Kuoさんらが公開しているノートにかかれていたテクニック
 
 > Sobol’ sequence tends to perform better
 > if an initial portion of the sequence is dropped:
 > the number of points skipped is the largest
 > power of 2 smaller than number of points to be
 > used. However, we are less persuaded by such 
 > recommendation ourselves.
 35
  • 36.
  • 37.
    まとめ
 ● Sobol sequenceは、Lowdiscrepancyという性質をもつ乱数のようなものがある。
 ○ 主な用途はOptimizationやNumerical Integration。
 ● 今回はアルゴリズムの解説論文を読んで紹介。
 ○ Theoretical underpinningsは触れられていないので、どうしてLow discrepancy になるのかは、Sobolさんの論文を参照
 37