Sobol 序列並行化的實踐經驗
Sobol 序列並行化的實踐經驗
隨機數發生器並行化的常見策略
隨機數發生器的並行化通常有四種策略（文獻【2】）：
- 「隨機化」種子
- 參數化隨機數發生器
- 分塊策略
- 蛙跳策略
「隨機化」種子
「隨機化」種子就是讓每個執行緒上運行的發生器使用不同的種子，這種策略常見於偽隨機數。不過，使用不同隨機數種子的發生器所產生的隨機數是否依然保持統計意義上的獨立性，這一點是未知的。
參數化隨機數發生器
某些隨機數發生器可以採用不同的參數配置，例如最簡單的線性同餘發生器，讓每個執行緒運行不同參數配置的發生器就可以實現並行化。不過，這種策略和「隨機化」種子有相同的安全隱患，無法確保產生的隨機數依然保持統計意義上的獨立性。另外，即使保持了獨立性，可用的參數配置組合可能數量十分有限，無法充分發揮硬體的並行優勢。
分塊策略
要確保產生的隨機數保持統計意義上的獨立性，最簡單的辦法就是使用一個發生器，在不同的執行緒上返回不同的數字。要做到這一點可以採用分塊生成的策略，下面舉個例子。
無論是偽隨機還是擬隨機，發生器產生的隨機數的順序是不變的。假設要在 4 個執行緒上並行發生 10000 個隨機數，分塊策略就是讓第 1 個執行緒返回第 1~2500 個隨機數；第 2 個執行緒返回第 2501~5000 個隨機數；第 3 個執行緒返回第 5001~7500 個隨機數；第 4 個執行緒返回第 7501~10000 個隨機數。
分塊策略要求發生器必須有一個「長距離跳轉」的功能，而且要保證計算成本盡量低。
蛙跳策略
想要使用一個發生器在不同的執行緒上返回不同的數字，還可以採用「蛙跳」的策略，下面舉個例子。
假設要在 4 個執行緒上並行發生 10000 個隨機數，蛙跳策略就是讓第 1 個執行緒返回第 \(4k\) 個隨機數；第 2 個執行緒返回第 \(4k + 1\) 個隨機數；第 3 個執行緒返回第 \(4k + 2\) 個隨機數；第 4 個執行緒返回第 \(4k + 3\) 個隨機數。
同分塊策略相比，每個執行緒上產生的數字是間隔開的，跳轉步長等於執行緒數。
因為要多次跳轉，蛙跳策略要求發生器必須有一個「極低成本」的「短距離跳轉」的功能。
Sobol 序列的原理和跳轉功能
Sobol 序列的產生依賴於一組事先確定「方向數」——\(m\)，以及「異或」運算。對於多維 Sobol 序列，每個維度有各自的方向數。
Sobol 序列的發生可以表示為非遞歸和遞歸兩種形式：
- 非遞歸形式：\(y_n = g_1 m_1 \oplus g_2 m_2 \oplus g_3 m_3 \oplus \cdots\)
- 遞歸形式：\(y_n = y_{n-1} \oplus m_{f(n-1)}, y_0 = 0\)
其中 \(\oplus\) 表示異或運算；\(g_i\) 是 \(n\) 的 Gray 碼的二進位表示，\(n\) 的 Gray 碼等於 \(n \oplus (n/2)\)，而 \(n \oplus (n/2) = \dots g_3 g_2 g_1\)；\(f(n)\) 表示 \(n\) 的二進位編碼中最右側的 0 所在的位數。
若採用 32 位整數，\(x_n = 2^{-32} y_n\)，\(x_n\) 就是服從單位均勻分布的 Sobol 序列。
非遞歸形式表明 Sobol 序列有低成本的長距離跳轉功能，因為若採用 32 位整數，直接計算 \(y_n\) 至多進行 32 次異或計算。計算 Sobol 序列的程式通常都會提供一個 skipTo 函數，實現長距離跳轉。
Sobol 序列並行化實踐
隨機數的並行化要求精準的操作執行緒，因此在最常見的「單機多核」情形下可以採用 OpenMP 框架。
分塊策略
假設已經具備了單執行緒上的 Sobol 序列發生器 SobolRsg，可以構造一個複合對象 OmpSobolRsg，內部包含 \(N\) 個完全相同的 SobolRsg 對象。在集中調用之前，事先告知 OmpSobolRsg 未來要調用的次數 \(M\)，將 \(N\) 個 SobolRsg 對象分別跳轉到合適的位置，第 \(i\) 個 SobolRsg 對象將專門在第 \(i\) 個執行緒上調用。這樣就實現了 Sobol 序列的並行發生。
藉助 OpenMP 框架，上述設計很容易實現，實現細節不贅述。
蛙跳策略
在蛙跳策略下，需要對 SobolRsg 對象稍加改造，以便實現從 \(y_{n}\) 直接到 \(y_{n+N}\) 的遞歸。和分塊策略類似，需要構造一個複合對象 OmpSobolRsg，內部包含 \(N\) 個完全相同的 SobolRsg 對象。無需事先告知未來要調用的次數，只需初始化階段將 \(N\) 個 SobolRsg 對象分別跳轉，恰好錯開一步。第 \(i\) 個 SobolRsg 對象將專門在第 \(i\) 個執行緒上調用。這樣就實現了 Sobol 序列的並行發生。
藉助 OpenMP 框架，上述設計也很容易實現，實現細節不贅述。
蛙跳策略的計算量分析
回想遞歸形式，在單執行緒情況下要生成 \(K\) 維的 Sobol 序列，需要進行 \(K\) 次整數間的異或計算，以及 \(K\) 次浮點數除法計算。蛙跳策略中的跳轉則會觸發若干次異或計算，當執行緒數量增多時，勢必會增加計算量，削弱並行化的效果。
減少異或計算的技巧
依據遞歸形式，跳轉 \(N\) 步的話有如下遞歸表達式：
\]
在原始情況下，跳轉 \(N\) 步將發生 \(N\) 次異或計算。但根據 \(f(n)\) 的定義，\(m_{f(n-i)}\) 中可能有不少重合的數字。
與此同時，注意到異或計算的幾個性質：
- \(a \oplus b = b \oplus a\)
- \(a \oplus (b \oplus c) = (a \oplus b) \oplus c\)
- \(a \oplus a = 0\)
- \(a \oplus 0 = a\)
- \(0 \oplus 0 = 0\)
也就是說，偶數個 \(m\) 可以歸併為 0，而 0 是可以忽略的。如果某個數字 \(m\) 出現了偶數次，可以直接忽略；如果某個數字 \(m\) 出現了奇數次，只保留一次就可以了。文獻【1】中就舉了這樣的例子，跳轉 \(2^p\) 步只需 2 次異或計算。
儘管如此，同分塊策略相比蛙跳策略避免不了頻繁跳轉，也就避免不了多餘的異或計算，這可能會導致並行 Monte Carlo 模擬的加速效率大打折扣。
分塊策略不算缺點的缺點
和蛙跳策略相比，分塊策略避免了反覆跳轉，但要求事先告知要調用的次數，這一點和隨機數發生器通常的使用模式相背離。但是瑕不掩瑜，分塊策略可以獲得更好地加速比。筆者分別採用分塊和蛙跳策略並行化了 QuantLib 的 Monte Carlo 定價框架，在 8 核 CPU 上，分塊策略的加速比在 7 倍左右，而蛙跳策略的加速比在 5.4 倍左右。
圖中可以看到，當執行緒數較少時，兩種策略難分伯仲；但當執行緒數較多時，蛙跳策略的加速效率明顯下降，這是頻繁跳轉造成的拖累。
參考文獻
- 【1】Parallelisation Techniques for Random Number Generators
- 【2】Tina’s Random Number Generator Library
