GPU加速04:將CUDA應用於金融領域，使用Python Numba加速B-S期權估值模型

2019 年 12 月 26 日
筆記

很多領域尤其是機器學習場景對GPU計算力高度依賴，所幸一些成熟的軟體或框架已經對GPU調用做了封裝，使用者無需使用CUDA重寫一遍，但仍需要對GPU計算的基本原理有所了解。對於一些無法調用框架的場景，當數據量增大時，非常有必要進行GPU優化。量化金融是一個非常好的應用GPU並行編程的領域。

本文為NVIDIA GPU計算加速系列的第四篇，主要基於前三篇文章的內容，以金融領域期權估值案例來進行實戰練習。前三篇文章為：

AI時代人人都應該了解的GPU知識：主要介紹了CPU與GPU的區別、GPU架構、CUDA軟體棧簡介。
超詳細Python Cuda零基礎入門教程：主要介紹了CUDA核函數，Thread、Block和Grid概念，記憶體分配，並使用Python Numba進行簡單的並行計算。
讓Cuda程式如虎添翼的優化技巧：主要從並行度和記憶體控制兩個方向介紹了多流和共享記憶體兩個優化技術。

閱讀完以上文章後，相信讀者已經對NVIDIA GPU編程有了初步的認識，這篇文章將談談如何將GPU編程應用到實際問題上，並使用Python Numba給出具體的B-S模型實現。

2017年於瑞士

應用場景

我在本系列開篇就曾提到目前GPU的應用場景非常廣泛：金融建模、自動駕駛、智慧機器人、新材料發現、神經科學、醫學影像…不同學科一般都有相應的軟體，比如分子動力學模擬軟體AMBER 16在NVIDIA 的GPU上的運行速度比僅使用CPU的系統快15倍；金融領域則需要使用GPU加速的機器學習來對各類金融產品做分析和預測。

GPU計算加速使用最廣泛的領域要數機器學習和深度學習了。各行各業（包括金融量化）都可以將本領域的問題轉化為機器學習問題。幸運的是，一些大神幫我們做了封裝，做成了框架供我們直接調用，省去了自己編寫機器學習演算法的時間，並且對GPU的支援非常強，無論是Google家的TensorFlow，還是Facebook家的PyTorch，以及XGBoost都對NVIDIA GPU有了非常不錯的支援，程式設計師幾乎不需要了解太多CUDA編程技術。雖然機器學習工程師不需要了解編程知識，但還是需要了解一些GPU的基礎知識（詳見本系列第一篇文章），非常有助於其深度學習項目的落地。關於TensorFlow等框架如何調用GPU，大家可先參考這些框架各自的官方文檔。

還有很多問題是與具體場景高度相關的，並不能直接用這些框架和庫，需要編程人員針對具體問題來編程。例如量化金融領域常常使用蒙特卡洛模擬，而CUDA對蒙特卡洛模擬也有非常好的支援，當數據量增大時，CUDA的優勢非常明顯。本文以金融領域著名的Black-Scholes模型為案例來展示如何使用Python Numba進行CUDA並行加速。B-S模型為Python Numba官方提供的樣常式序，我在原來基礎上做了一些簡單修改。

Black-Scholes模型簡介

Black-Scholes模型，簡稱B-S模型，是一種對金融產品估價的數學模型。該模型由Fischer Black和Myron Scholes提出，後由Robert Merton完善，這幾位學者憑藉該模型榮獲1997年榮獲諾貝爾經濟學獎。金融主要是在研究現在的錢與未來的錢的價值問題，B-S模型就是一種對期權產品初始價格和交割價格估值的方法。模型的公式如下。

B-S模型

使用上面這個公式，給定期權現價S、交割價格K，期權時間t，可以計算出看漲期權（Call Option）和看跌期權（Put Option）的價格。我本人並不是金融科班出身，就不在此班門弄斧解釋這個模型的金融含義了。對於程式設計師來說，一個重要的能力就是不需要對業務有太深入理解，也能使用程式碼實現需求。

B-S模型的Python實現

這裡我隨機生成了一組數據，包括期權現價S、交割價格K和期權時間t，數據維度分別為1000、100000, 1000000, 4000000。分別用"Python + Numpy"和"CUDA"方式實現，在高性能的Intel E5-2690 v4 CPU和Telsa V100 PCI-E版上運行，運行耗時如下圖所示。數據量越小，Python和Numpy在CPU上運行的程式越有優勢，隨著數據量增大，CPU程式耗時急速上升，GPU並行計算的優勢凸顯。當數據量為400萬時，CUDA程式可以獲得30+倍速度提升！試想，如果一個程式之前需要在CPU上跑一天，改成CUDA並行計算後，可能只需要一個小時，這是何等程度的生產力提升啊！

耗時對比

在實現B-S模型時，編寫了一個正態分布的累計概率分布函數（Cumulative Distribution Function）：cnd。關於概率密度函數和累計概率分布函數我這裡不做贅述，本科的概率論課程都會涉及，網路上也有很多詳細介紹。我隨機初始化了一些數據，並保存在了numpy向量中。注意，在CPU上使用numpy時，盡量不要用for對數組中每個數據處理，而要使用numpy的向量化函數。對於CPU程式來說，numpy向量盡量使用numpy.log()、numpy.sqrt()、numpy.where()等函數，因為numpy在CPU上做了大量針對向量的計算優化。但是對於標量，numpy可能比math庫慢。還需要注意的是，Numba的CUDA有可能不支援部分numpy向量操作。其他CPU的Python加速技巧，我會在後續文章中分享。

Python Numba庫支援的Numpy特性：https://numba.pydata.org/numba-doc/dev/reference/numpysupported.html

整個程式如下：

import numpy as np  import math  from time import time  from numba import cuda  from numba import jit  import matplotlib  # 使用無顯示器的伺服器進行計算時，需加上下面這行，否則matplot報錯  matplotlib.use('Agg')  import matplotlib.pyplot as plt    def timeit(func ,number_of_iterations, input_args):        # 計時函數      start = time()      for i in range(number_of_iterations):          func(*input_args)      stop = time()      return stop - start    def randfloat(rand_var, low, high):        # 隨機函數      return (1.0 - rand_var) * low + rand_var * high      RISKFREE = 0.02    VOLATILITY = 0.30    def cnd(d):      # 正態分布累計概率分布函數      A1 = 0.31938153      A2 = -0.356563782      A3 = 1.781477937      A4 = -1.821255978      A5 = 1.330274429      RSQRT2PI = 0.39894228040143267793994605993438      K = 1.0 / (1.0 + 0.2316419 * np.abs(d))      ret_val = (RSQRT2PI * np.exp(-0.5 * d * d) *                 (K * (A1 + K * (A2 + K * (A3 + K * (A4 + K * A5))))))      return np.where(d > 0, 1.0 - ret_val, ret_val)      # 上面的Numpy函數比下面使用循環效率高很多      # for i in range(len(d)):      #     if d[i] > 0:      #         ret_val[i] = 1.0 - ret_val[i]      # return ret_val    def black_scholes(stockPrice, optionStrike, optionYears, riskFree, volatility):        # Python + Numpy 實現B-S模型      S = stockPrice      K = optionStrike      T = optionYears      r = riskFree      V = volatility      sqrtT = np.sqrt(T)      d1 = (np.log(S / K) + (r + 0.5 * V * V) * T) / (V * sqrtT)      d2 = d1 - V * sqrtT      cndd1 = cnd(d1)      cndd2 = cnd(d2)      expRT = np.exp(- r * T)      callResult = S * cndd1 - K * expRT * cndd2      putResult = K * expRT * (1.0 - cndd2) - S * (1.0 - cndd1)      return callResult, putResult    @cuda.jit(device=True)    def cnd_cuda(d):      # 正態分布累計概率分布函數      # CUDA設備端函數      A1 = 0.31938153      A2 = -0.356563782      A3 = 1.781477937      A4 = -1.821255978      A5 = 1.330274429      RSQRT2PI = 0.39894228040143267793994605993438      K = 1.0 / (1.0 + 0.2316419 * math.fabs(d))      ret_val = (RSQRT2PI * math.exp(-0.5 * d * d) *                 (K * (A1 + K * (A2 + K * (A3 + K * (A4 + K * A5))))))      if d > 0:          ret_val = 1.0 - ret_val      return ret_val    @cuda.jit    def black_scholes_cuda_kernel(callResult, putResult, S, K,                           T, r, V):      i = cuda.threadIdx.x + cuda.blockIdx.x * cuda.blockDim.x      if i >= S.shape[0]:          return      sqrtT = math.sqrt(T[i])      d1 = (math.log(S[i] / K[i]) + (r + 0.5 * V * V) * T[i]) / (V * sqrtT)      d2 = d1 - V * sqrtT      cndd1 = cnd_cuda(d1)      cndd2 = cnd_cuda(d2)      expRT = math.exp((-1. * r) * T[i])      callResult[i] = (S[i] * cndd1 - K[i] * expRT * cndd2)      putResult[i] = (K[i] * expRT * (1.0 - cndd2) - S[i] * (1.0 - cndd1))    def black_scholes_cuda(stockPrice, optionStrike,                            optionYears, riskFree, volatility):      # CUDA實現B-S模型      blockdim = 1024      griddim = int(math.ceil(float(len(stockPrice))/blockdim))      device_callResult = cuda.device_array_like(stockPrice)      device_putResult = cuda.device_array_like(stockPrice)      device_stockPrice = cuda.to_device(stockPrice)      device_optionStrike = cuda.to_device(optionStrike)      device_optionYears = cuda.to_device(optionYears)      black_scholes_cuda_kernel[griddim, blockdim](              device_callResult, device_putResult, device_stockPrice, device_optionStrike,              device_optionYears, riskFree, volatility)      callResult = device_callResult.copy_to_host()      putResult= device_putResult.copy_to_host()      cuda.synchronize()      return callResult, putResult      def main():        pure_time_list = []      cuda_time_list = []      dtype = np.float32      data_size = [10, 1000, 100000, 1000000, 4000000]      for OPT_N in data_size:          print("data size :" + str(OPT_N))          stockPrice = randfloat(np.random.random(OPT_N), 5.0, 30.0).astype(dtype)          optionStrike = randfloat(np.random.random(OPT_N), 1.0, 100.0).astype(dtype)          optionYears = randfloat(np.random.random(OPT_N), 0.25, 10.0).astype(dtype)          input_args=(stockPrice, optionStrike, optionYears, RISKFREE, VOLATILITY)          pure_duration = timeit(black_scholes, 20, input_args)          pure_time_list.append(pure_duration)          cuda_duration = timeit(black_scholes_cuda, 20, input_args)          cuda_time_list.append(cuda_duration)      print(pure_time_list)      print(cuda_time_list)      plt.plot(pure_time_list[1:], label='pure python')      plt.plot(cuda_time_list[1:], label='cuda')      plt.legend(loc='upper right')      plt.xticks(np.arange(5), ('1000', '100000', '1000000', '4000000'))      #設置坐標軸名稱      plt.ylabel('duration (second)')      plt.xlabel('option number')      plt.savefig("benchmark.png")      plt.show()    if __name__ == "__main__":        main()

程式都是非常基礎的CUDA使用技巧，在我的第二篇文章中都有提到，並沒有使用太多優化技巧。其中，cnd_cuda函數使用了@cuda.jit(device=True)修飾，表示這個函數只是GPU端做計算的設備函數。

小結

很多領域尤其是機器學習場景對GPU計算力高度依賴，所幸一些成熟的軟體或框架已經對GPU調用做了封裝，使用者無需使用CUDA重寫一遍，但仍需要對GPU計算的基本原理有所了解。對於一些無法調用框架的場景，當數據量增大時，非常有必要進行GPU優化。量化金融中經常使用蒙特卡洛模擬和機器學習等技術，是一個非常好的應用GPU並行編程的領域。