好買網 IT技術交易平臺

Sobol 序列並行化的實踐經驗

隨機數發生器並行化的常見策略

隨機數發生器的並行化通常有四種策略（文獻【2】）：

1. “隨機化”種子
2. 引數化隨機數發生器
3. 分塊策略
4. 蛙跳策略

“隨機化”種子

“隨機化”種子就是讓每個執行緒上執行的發生器使用不同的種子，這種策略常見於偽隨機數。不過，使用不同隨機數種子的發生器所產生的隨機數是否依然保持統計意義上的獨立性，這一點是未知的。

引數化隨機數發生器

某些隨機數發生器可以採用不同的引數配置，例如最簡單的線性同餘發生器，讓每個執行緒執行不同引數配置的發生器就可以實現並行化。不過，這種策略和“隨機化”種子有相同的安全隱患，無法確保產生的隨機數依然保持統計意義上的獨立性。另外，即使保持了獨立性，可用的引數配置組合可能數量十分有限，無法充分發揮硬體的並行優勢。

分塊策略

要確保產生的隨機數保持統計意義上的獨立性，最簡單的辦法就是使用一個發生器，在不同的執行緒上返回不同的數字。要做到這一點可以採用分塊生成的策略，下面舉個例子。

無論是偽隨機還是擬隨機，發生器產生的隨機數的順序是不變的。假設要在 4 個執行緒上並行發生 10000 個隨機數，分塊策略就是讓第 1 個執行緒返回第 1~2500 個隨機數；第 2 個執行緒返回第 2501~5000 個隨機數；第 3 個執行緒返回第 5001~7500 個隨機數；第 4 個執行緒返回第 7501~10000 個隨機數。

分塊策略要求發生器必須有一個“長距離跳轉”的功能，而且要保證計算成本儘量低。

蛙跳策略

想要使用一個發生器在不同的執行緒上返回不同的數字，還可以採用“蛙跳”的策略，下面舉個例子。

假設要在 4 個執行緒上並行發生 10000 個隨機數，蛙跳策略就是讓第 1 個執行緒返回第  4 k 4k 個隨機數；第 2 個執行緒返回第  4 k + 1 4k+1 個隨機數；第 3 個執行緒返回第  4 k + 2 4k+2 個隨機數；第 4 個執行緒返回第  4 k + 3 4k+3 個隨機數。

同分塊策略相比，每個執行緒上產生的數字是間隔開的，跳轉步長等於執行緒數。

因為要多次跳轉，蛙跳策略要求發生器必須有一個“極低成本”的“短距離跳轉”的功能。

Sobol 序列的原理和跳轉功能

Sobol 序列的產生依賴於一組事先確定“方向數”—— m m，以及“異或”運算。對於多維 Sobol 序列，每個維度有各自的方向數。

Sobol 序列的發生可以表示為非遞迴和遞迴兩種形式：

• 非遞迴形式： y n = g 1 m 1 ⊕ g 2 m 2 ⊕ g 3 m 3 ⊕ ⋯ yn=g1m1⊕g2m2⊕g3m3⊕⋯
• 遞迴形式： y n = y n − 1 ⊕ m f ( n − 1 ) , y 0 = 0 yn=yn−1⊕mf(n−1),y0=0

其中  ⊕ ⊕ 表示異或運算； g i gi 是  n n 的 Gray 碼的二進位制表示， n n 的 Gray 碼等於  n ⊕ ( n / 2 ) n⊕(n/2)，而  n ⊕ ( n / 2 ) = … g 3 g 2 g 1 n⊕(n/2)=…g3g2g1； f ( n ) f(n) 表示  n n 的二進位制編碼中最右側的 0 所在的位數。

若採用 32 位整數， x n = 2 − 32 y n xn=2−32yn， x n xn 就是服從單位均勻分佈的 Sobol 序列。

非遞迴形式表明 Sobol 序列有低成本的長距離跳轉功能，因為若採用 32 位整數，直接計算  y n yn 至多進行 32 次異或計算。計算 Sobol 序列的程式通常都會提供一個  skipTo 函式，實現長距離跳轉。

Sobol 序列並行化實踐

隨機數的並行化要求精準的操作執行緒，因此在最常見的“單機多核”情形下可以採用 OpenMP 框架。

分塊策略

假設已經具備了單執行緒上的 Sobol 序列發生器  SobolRsg，可以構造一個複合物件  OmpSobolRsg，內部包含  N N 個完全相同的  SobolRsg 物件。在集中呼叫之前，事先告知  OmpSobolRsg 未來要呼叫的次數  M M，將  N N 個  SobolRsg 物件分別跳轉到合適的位置，第  i i 個  SobolRsg 物件將專門在第  i i 個執行緒上呼叫。這樣就實現了 Sobol 序列的並行發生。

藉助 OpenMP 框架，上述設計很容易實現，實現細節不贅述。

蛙跳策略

在蛙跳策略下，需要對  SobolRsg 物件稍加改造，以便實現從  y n yn 直接到  y n + N yn+N 的遞迴。和分塊策略類似，需要構造一個複合物件  OmpSobolRsg，內部包含  N N 個完全相同的  SobolRsg 物件。無需事先告知未來要呼叫的次數，只需初始化階段將  N N 個  SobolRsg 物件分別跳轉，恰好錯開一步。第  i i 個  SobolRsg 物件將專門在第  i i 個執行緒上呼叫。這樣就實現了 Sobol 序列的並行發生。

藉助 OpenMP 框架，上述設計也很容易實現，實現細節不贅述。

蛙跳策略的計算量分析

回想遞迴形式，在單執行緒情況下要生成  K K 維的 Sobol 序列，需要進行  K K 次整數間的異或計算，以及  K K 次浮點數除法計算。蛙跳策略中的跳轉則會觸發若干次異或計算，當執行緒數量增多時，勢必會增加計算量，削弱並行化的效果。

減少異或計算的技巧

依據遞迴形式，跳轉  N N 步的話有如下遞迴表示式：

在原始情況下，跳轉  N N 步將發生  N N 次異或計算。但根據  f ( n ) f(n) 的定義， m f ( n − i ) mf(n−i) 中可能有不少重合的數字。

與此同時，注意到異或計算的幾個性質：

• a ⊕ b = b ⊕ a a⊕b=b⊕a
• a ⊕ ( b ⊕ c ) = ( a ⊕ b ) ⊕ c a⊕(b⊕c)=(a⊕b)⊕c
• a ⊕ a = 0 a⊕a=0
• a ⊕ 0 = a a⊕0=a
• 0 ⊕ 0 = 0 0⊕0=0

也就是說，偶數個  m m 可以歸併為 0，而 0 是可以忽略的。如果某個數字  m m 出現了偶數次，可以直接忽略；如果某個數字  m m 出現了奇數次，只保留一次就可以了。文獻【1】中就舉了這樣的例子，跳轉  2 p 2p 步只需 2 次異或計算。

儘管如此，同分塊策略相比蛙跳策略避免不了頻繁跳轉，也就避免不了多餘的異或計算，這可能會導致並行 Monte Carlo 模擬的加速效率大打折扣。

分塊策略不算缺點的缺點

和蛙跳策略相比，分塊策略避免了反覆跳轉，但要求事先告知要呼叫的次數，這一點和隨機數發生器通常的使用模式相背離。但是瑕不掩瑜，分塊策略可以獲得更好地加速比。筆者分別採用分塊和蛙跳策略並行化了 QuantLib 的 Monte Carlo 定價框架，在 8 核 CPU 上，分塊策略的加速比在 7 倍左右，而蛙跳策略的加速比在 5.4 倍左右。