numba を使ってGPU計算してみよう

numba について

numba のインストール

numba は次のようにインストールします。

code: bash

$ pip install numba

Continuum Analytics, Incは、Numbaの拡張バージョンであるNumbaProをリリースしていあす。

このバージョンは、学術研究機関では無料で利用することができます(アカデミックフリー)。

開発者がNumPyとうまく統合できる最適化されたコードを迅速に作成できるようにする

プレミアム機能が追加されています。

Numba のCUDAツールキットの検察順序

Numba では CUDAツールキットを利用してGPU用のコードを生成します。

Numbaは、CUDA Toolkitのインストールを次の順序で検索します。

Conda コマンドでインストールしたcudatoolkitパッケージ

ユーザーはCUDA_HOMEを使用して特定のバージョンを選択できる

Numpy ufunc

Numpyユニバーサル関数またはufuncは、要素ごとにnumpy配列を操作する関数です。 たとえば、numpy配列の二乗を取る場合、ufuncは、結果の配列を返す前に、各要素の二乗を計算します。

code: ipython

...: import numpy as np

...: x = np.arange(10)

...: x**2

...:

ほとんどの数学関数は、numpyのufuncとして利用することができます。 たとえば、numpy配列のすべての要素をべき乗するには、次のようにします。

code: ipython

array([1.00000000e+00, 2.71828183e+00, 7.38905610e+00, 2.00855369e+01,

5.45981500e+01, 1.48413159e+02, 4.03428793e+02, 1.09663316e+03,

2.98095799e+03, 8.10308393e+03])

ほとんどのufuncはコンパイルされたCコードで実装されているため、Pythonでコードするよりもはるかに高速に計算されます。

次のコードは、大きな配列を Pure Pythonとnumpyの両方で各要素の対数を計算したものです。

code: python

...: import numpy as np

...: import math

...: x = np.arange(int(1e6))

...: # %timeit np.sqrt(x)

...:

280 ms ± 24.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

7.46 ms ± 221 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

numpy ufuncは約40倍高速であることがわかります。 それでも、ufuncでは、計算は並列化されません。平方根は、配列内の各要素ごとにCPU上で順次計算されます。

GPUは、CPUとは異なり、多数の操作を同時に実行できます。

次にGPUのufuncを作成して、多くの要素で同時に計算を実行する方法を説明します。

GPU ufunc

GPUのufuncは、多数の要素の平方根を計算することができます。

まず、0から1000万の範囲のポイントのサンプルを作成します。 GPUは、少数のビットでエンコードされた数値を使用すると、より効率的になります。 そして、多くの場合、非常に高い精度は必要ありません。 したがって、float32番号のサンプルを作成します（デフォルトはfloat64です）。

code: python

import numpy as np

npoints = int(1e7)

a = np.arange(npoints, dtype=np.float32)

https://gyazo.com/d1fe2336c40a17de2dbc8bf67255aac2

https://gyazo.com/ba28593b573cb813f98d61c79ceb93a1

戻り値は関数として考えることができ、入力値を引数に与えているわけです。 また、これらの型データに適合している必要があることに注意してください。

code: python

TypeError: No matching version. GPU ufunc requires array arguments to have the exact types. This behaves like regular ufunc with casting='no'.

ここで、性能差を確認してみましょう。

https://gyazo.com/ad95477470d3941b589003e9694263b8

あれ！？　CPUバージョンの方が４倍以上速い結果となっていますね。

これには単純な理由があります。 GPUで実行しているときは内部で次のことが発生しています。

結果の配列はホストシステムに転送する

計算が単純すぎる場合や、データ転送が完了するまで長時間待機する場合では、高速並列処理のためにデータをデバイスに送信しても効果が期待できません。つまり、処理時間のほとんどはデータ転送に費やされ、GPUのパワーが活かしきれないわけです。

複雑な計算をしてみる

もっと複雑な計算をさせてみましょう。

https://gyazo.com/4daf360d871ee3e6be935834ada293f9

２次元座標を極座標に変換すると便利なことがよくあります。

https://gyazo.com/174fc6f004d4f65274c2a8b1e2d40a08

$ r = \sqrt{ x^2 + y^2 }

https://gyazo.com/1029d607774bbb493204b67729297a78

次に、GPUで同じ計算を実行してみましょう。

今度の入力値は2つになります。

code: ipython

import math

def gpu_arctan2(y, x):

theta = math.atan2(y,x)

return theta

code: python

一見問題ないように見えますが、実はこれはエラーになります。

次のようにスライスを指示した結果を連続した配列に変換します。

code: python

https://gyazo.com/dd1f57cec8428b5b24b14a66457d3d0b

1,000万ポイントに増やして計算してみましょう。

https://gyazo.com/3600428b57772e7c1d3d468dee5077f3

もう一度、実行速度を比較してみましょう。

https://gyazo.com/d937844fbd9d8b674485409c58264c39

今度は、期待どおりにGPUバージョンがNumPyと比較して9.6倍、Pure Python バージョンとでは176倍速い結果となっています。

わかりやすく言い換えると、PurePython で１週間かかる計算がGPUでは1時間で終わるということです。

Generalized ufuncs をGPUで使用する

通常のufuncでは、計算は入力配列の各要素で実行され、スカラーを返します。

Generalized ufuncs（gufuncs）では、計算で入力配列のサブ配列を処理し、さまざまな次元の配列を返すことができます。

前述の２次元座標から極座標への変換のような処理は、gufuncs を使用すると簡単になります。

code: python

from numba import guvectorize

'(i)->(i)',

target='cuda')

def gpu_polar(vec, out):

補足説明：シグニチャ(signature)

シグニチャ（signature）は、もともとは署名の意味ですが、関数や変数、配列などのエンティティーの名前と型、その並び順組み合わせた情報のことを言います。

https://gyazo.com/8d2b64506197c837d3269a764cd65d41

https://gyazo.com/f24eb108872ddbb0e4af95b045a68a2a

gufuncsがどのように機能するかをよく理解するために、2D配列の各行の値の平均を計算するgufuncsを作ってみましょう。

code: python

'(n)->()',

target='cuda')

def gpu_average(array, out):

acc = 0

for val in array:

acc += val

２D配列を作っておきます。

https://gyazo.com/862783d8393c2110b8bfa336e681dba1

https://gyazo.com/57dd06ba3e88666d2519e358a32bb09c

デバイス関数

GPUでufuncまたはgufuncのいずれかの単一の関数を実行してきましたが、すべての処理を単一の関数としてコードする必要はありません。

ヘルパー関数をコンパイルすることも可能で、これらの関数をデバイス関数と呼びます。デバイス関数をGPUで使用するとコードをモジュール化することができます。

前述の2D配列の値の平均を計算するgufuncをもう一度見てみましょう。

code: python

from numba import cuda

@cuda.jit(device=True)

def add(array):

acc = 0

for val in array:

acc += val

return acc

'(n)->()',

target='cuda')

def gpu_average_2(array, out):

'(n)->()',

target='cuda')

def gpu_sum(array, out):

https://gyazo.com/0008b09eca7654951a70dcf89473596b

複雑なアルゴリズムを実装する場合、コードの重複を避けることができるのでデバイス関数は非常に便利です。

デバイス配列によるメモリ管理

一番最初のufuncの例でわかるように、GPUを使うといつもパフォーマンスが向上するとは限りません。

実際、GPUの計算能力を使用するためには、その前後にデータの転送が必要です。

データをメモリからデバイスに転送する

計算結果はデバイスからメモリに転送する

パフォーマンスを向上させるポイントは、ホストシステムのメモリとデバイス間のデータ転送を最小限に抑えることです。これは、デバイス配列(DeviceArray)を使用して実行できます。

画像処理用のニューラルネットワークを実装することを考えてみます。ネットワークの隠れ層は、次のことを行う必要がある場合があります。

画像のグレースケール値を正規化する

それらの重みを計算する

アクティベーション関数を適用する

これらの3つのタスクはそれぞれ、GPU上で並行して実行できます。

まず最初に、単純なnumpyを使用してCPUで処理してみます。

code: python

n = 100_000_000

# 0 〜 255 のランダムな値

greyscales = np.floor(np.random.uniform(0, 256, n).astype(np.float32))

# 0.5を中心とした幅が0.1の正規分布に従うランダムな重み

weights = np.random.normal(.5, .1, n).astype(np.float32)

def normalize(grayscales):

return grayscales / 255

def weigh(values, weights):

return values * weights

def activate(values):

return ( np.exp(values) - np.exp(-values) ) / \

( np.exp(values) + np.exp(-values) )

https://gyazo.com/41a8777b4fc22718f6d2f8e754b53ae7

次に、GPU用にこのアルゴリズムの並列バージョンを実装します。

code: python

def gpu_normalize(x):

return x / 255

def gpu_weigh(x, w):

return x * w

def gpu_activate(x):

return ( math.exp(x) - math.exp(-x) ) / ( math.exp(x) + math.exp(-x) )

https://gyazo.com/99b17a2b452b1f691c5b89a3967ab72e

そこそこ良い結果が得られました。

しかし、ホストとGPUの間での次のデータ転送の時間があるはずです。

実際に必要なこものは、次の転送だけのはずです。

https://gyazo.com/0f9c0784fb38670e44976d7bc48e39ec

不要なデータ転送をなくすことで処理速度が向上したのがわかります。

次のようにGPUとの間の転送を完全に制御することもできます。

code: python

from numba import cuda

# 入力データをGPUに転送

greyscales_gpu = cuda.to_device(greyscales)

weights_gpu = cuda.to_device(weights)

# GPUで中間配列と出力配列を作成

normalized_gpu = cuda.device_array(shape=(n,),

dtype=np.float32)

weighted_gpu = cuda.device_array(shape=(n,),

dtype=np.float32)

activated_gpu = cuda.device_array(shape=(n,),

dtype=np.float32)

https://gyazo.com/3e2c56abf2ef47a347d54ed5cc9abde1

この結果には、入力データと出力データに必要な転送時間が含まれていないため、かなり意図的な計測結果ですが、データ転送を制御する方法を示すために例示しています。

しかし、ここで定義したデバイス配列のいずれかを再利用したいときは、まだデバイス上にあるためすぐ再利用することができます。

結果を取得する方法は次のとおりです。

https://gyazo.com/a7dc12df639b39df5604120327a1eeb5

Numba.Cudaの詳細

CUDAカーネルの作成

CUDAには、CPUのプログラミングに使用される従来のシーケンシャルモデルとは異なる実行モデルがあります。 CUDAでは、作成したコードは一度に複数のスレッド（多くの場合、数百または数千）によって実行されます。ソリューションは、グリッド、ブロック、およびスレッドのスレッド階層を定義することによってモデル化されます。

NumbaのCUDAサポートは、このスレッドの階層を宣言や、管理するための機能を公開します。施設は、NVidiaのCUDA C言語(nvcc)で公開されているものとほぼ同じです。

Numbaは、グローバルデバイスメモリ（GPUに接続されている大きくて比較的低速のオフチップメモリ​​）、オンチップ共有メモリ、ローカルメモリの3種類のGPUメモリも公開しています。最も単純なアルゴリズムを除くすべての場合、帯域幅の要件と競合を最小限に抑えるために、メモリの使用方法とアクセス方法を慎重に検討することが重要です。

カーネル宣言

カーネル関数は、CPUコードから呼び出されることを意図したGPU関数で、次の2つの基本的な特徴があります。

カーネルは明示的に値を返すことはできない

すべての結果データは、関数に渡される配列に書き込む必要があります（スカラーを計算する場合は、おそらく1要素の配列を渡します）。

カーネルは、呼び出されたときにスレッド階層を明示的に宣言する

つまり、スレッドブロックの数とブロックあたりのスレッドの数です。コンパイルされたカーネルは、異なるブロックサイズまたはグリッドサイズで複数回呼び出すことができます。

補足説明：

新しいCUDAデバイスは、デバイス側のカーネル起動をサポートしています。

この機能は動的並列処理と呼ばれますが、現在のNumbaではまだサポートしていません

一見すると、Numbaを使用してCUDAカーネルを作成することは、CPU用のJIT関数を作成することと非常によく似ています。

code: python

@cuda.jit

def increment_by_one(an_array):

# ...

カーネル呼び出し

カーネルは通常、次の方法で起動されます。

code: python

threadsperblock = 32

blockspergrid = (an_array.size +（threadsperblock-1)) // threadsperblock

ここで2つのステップがあることがわかります。

カーネルを実行し、入力配列（および必要に応じて個別の出力配列）を渡します。 デフォルトでは、カーネルの実行は同期的です。カーネルの実行が終了し、データが同期されたときに関数が戻ります。

ブロックサイズの選択

カーネルに必要なスレッドの数を宣言するときに、2レベルの階層を持つのは不思議に思えるかもしれません。多くの場合、ブロックサイズ（つまり、ブロックあたりのスレッド数）は重要です。

ソフトウェア側では、ブロックサイズによって、共有メモリの特定の領域を共有するスレッドの数が決まります。

ハードウェア側では、ブロックサイズは実行ユニットを完全に占有するのに十分な大きさである必要があります。

多次元ブロックとグリッド

スレッドの配置

カーネルを実行する場合、カーネル関数のコードはすべてのスレッドによって1回実行されます。したがって、どの配列要素が原因であるかを知るために、どのスレッドにあるかを知る必要があります

1つの方法は、スレッドがグリッドとブロック内の位置を決定し、対応する配列の位置を手動で計算することです。

code: python

@cuda.jit

def increment_by_one(an_array):

# 1DブロックのスレッドID

tx = cuda.threadIdx.x

# 1DグリッドのブロックID

ty = cuda.blockIdx.x

# ブロック幅、つまりブロックあたりのスレッド数

bw = cuda.blockDim.x

# 配列内のフラット化されたインデックスを計算

pos = tx + ty * bw

if pos < an_array.size: # 配列の境界を確認

numba.cuda.threadIdx

numba.cuda.blockDim

カーネルをインスタンス化するときに宣言された、スレッドのブロックの形状。この値は、異なるブロックに属している（つまり、各ブロックが「フル」である）場合でも、特定のカーネル内のすべてのスレッドで同じです。

numba.cuda.blockIdx

numba.cuda.gridDim

ブロックのグリッドの形状、つまり、カーネルをインスタンス化するときに宣言された、このカーネル呼び出しによって起動されたブロックの総数。

絶対位置

単純なアルゴリズムは、上記の例に示されているのと同じ方法で常にスレッドインデックスを使用する傾向があります。 Numbaは、このような計算を自動化するための追加機能を提供します。

numba.cuda.grid(ndim)

numba.cuda.gridsize(ndim)

これらの関数を使用すると、インクリメントの例は次のようになります。

code: python

@cuda.jit

def increment_by_one(an_array):

pos = cuda.grid(1)

if pos < an_array.size:

２次元配列とスレッドのグリッドの同じ例は次のようになります。

code: pyton

@cuda.jit

def increment_a_2D_array(an_array):

x, y = cuda.grid(2)

カーネルをインスタンス化するときのグリッド計算は、引き続き手動で実行する必要があることに注意してください。

code: python

threadsperblock = (16, 16)

blockspergrid = (blockspergrid_x, blockspergrid_y)

メモリ管理

データ転送

NumbaはNumPyアレイをデバイスに自動的に転送できますが、カーネルの終了時に常にデバイスのメモリをホストに転送することによってのみ、保守的に転送できます。 読み取り専用配列の不要な転送を回避するには、次のAPIを使用して転送を手動で制御できます。

numba.cuda.device_array(shape, dtype=np.float, strides=None, order='C', stream=0)

numba.cuda.device_array_like(ary, stream=0)

numba.cuda.to_device(obj, stream=0, copy=True, to=None)

numpy配列をホストからGPUに転送する

code: python

ary = np.arange(10)

d_ary = cuda.to_device(ary)

転送をストリームにキューに入れる

code: python

d_ary = cuda.to_device(ary, stream=stream)

結果をGPUからホストに転送する

code: python

hary = d_ary.copy_to_host()

結果をGPUっからホストの既存のアレイにコピーする

code: python

ary = np.empty(shape=d_ary.shape, dtype=d_ary.dtype)

d_ary.copy_to_host(ary)

転送をストリームにキューに入れる

code: python

hary = d_ary.copy_to_host(stream=stream)

Numbaは、デバイスアレイに加えて、cudaアレイインターフェイスを実装する任意のオブジェクトを使用できます。 これらのオブジェクトは、次のAPIを使用してGPUバッファーのビューを作成することにより、手動でNumbaデバイス配列に変換することもできます。

numba.cuda.as_cuda_array(obj)

numba.cuda.is_cuda_array(obj)

インターフェイスの有効性を検証しません。

DeviceArrayクラス

class numba.cuda.cudadrv.devicearray.DeviceNDArray(shape, strides, dtype, stream=0, writeback=None, gpu_data=None)

オンGPU配列タイプ

copy_to_host(self, ary=None, stream=0)

CUDAストリームが指定されている場合、転送は指定されたストリームの一部として非同期で行われます。 それ以外の場合、転送は同期的です。コピーが終了した後、関数は戻ります。

常にホスト配列を返します。

code: python

import numpy as np

from numba import cuda

arr = np.arange(1000)

d_arr = cuda.to_device(arr)

result_array = d_arr.copy_to_host()

is_c_contiguous(self)

配列がC形式配列の場合、trueを返します。

is_f_contiguous(self)

配列がFortran形配列に隣接している場合は、trueを返します。

ravel(self, order='C', stream=0)

reshape(self, *newshape, **kws)

code: python

d_arr = d_arr.reshape(20, 50, order='F')

ピンメモリ（Pinnedメモリ）

ページングされないメモリは，page-locked memoryやpinned memoryなどと呼ばれます。

numba.cuda.pinned( *arylist)

numba.cuda.pinned_array(shape, dtype=np.float, strides=None, order='C')

ストリーム

numba.cuda.stream()

デバイスのコマンドキューを表すCUDAストリームを作成します。

numba.cuda.default_stream()

デフォルトのCUDAストリームを取得します。一般に、CUDAセマンティクスでは、デフォルトストリームは、使用されているCUDA APIに応じて、レガシーデフォルトストリームまたはスレッドごとのデフォルトストリームのいずれかになります。 Numbaでは、レガシーデフォルトストリームのAPIが常に使用されていますが、スレッドごとのデフォルトストリームにAPIを使用するオプションが将来提供される可能性があります。

numba.cuda.legacy_default_stream()

従来のデフォルトのCUDAストリームを取得します。

numba.cuda.per_thread_default_stream()

スレッドごとのデフォルトのCUDAストリームを取得します。

numba.cuda.external_stream(ptr)

Numbaの外部に割り当てられたストリームのNumbaストリームオブジェクトを作成します。

CUDAストリームには次の方法があります。

class numba.cuda.cudadrv.driver.Stream(context, handle, finalizer, external=False)

auto_synchronize(self)

このストリーム内のすべてのコマンドが実行されるのを待機し、コンテキストの終了時に保留中のメモリ転送をコミットするコンテキストマネージャー。

synchronize(self)

このストリーム内のすべてのコマンドが実行されるのを待ちます。これにより、保留中のメモリ転送がコミットされます。

共有メモリとスレッドの同期

必要に応じて、限られた量の共有メモリをデバイスに割り当てて、データへのアクセスを高速化できます。そのメモリは、特定のブロックに属するすべてのスレッド間で共有され（つまり、読み取りと書き込みの両方）、通常のデバイスメモリよりもアクセス時間が速くなります。また、スレッドが特定のソリューションで連携できるようにします。手動で管理されるデータキャッシュと考えることができます。

従来の動的メモリ管理とは異なり、メモリはカーネルの期間中に1回割り当てられます。

numba.cuda.shared.array(shape, type)

指定された形状とタイプの共有配列をデバイスに割り当てます。この関数は、デバイス上で呼び出す必要があります（つまり、カーネルまたはデバイス関数から）。 shapeは、配列の次元を表す整数または整数のタプルのいずれかであり、単純な定数式である必要があります。 typeは、配列に格納する必要のある要素のNumbaタイプです。

返された配列のようなオブジェクトは、通常のデバイス配列と同じように読み書きできます（インデックス作成など）。

numba.cuda.syncthreads()

同じスレッドブロック内のすべてのスレッドを同期します。この関数は、従来のマルチスレッドプログラミングのバリアと同じパターンを実装します。この関数は、ブロック内のすべてのスレッドがそれを呼び出すまで待機し、その時点ですべての呼び出し元に制御を返します。

ローカルメモリ

ローカルメモリは、各スレッド専用のメモリ領域です。ローカルメモリを使用すると、スカラーローカル変数が十分でない場合にスクラッチパッド領域を割り当てるのに役立ちます。従来の動的メモリ管理とは異なり、メモリはカーネルの期間中に1回割り当てられます。

numba.cuda.local.array(shape, type)

指定された形状とタイプのローカル配列をデバイスに割り当てます。 shapeは、配列の次元を表す整数または整数のタプルのいずれかであり、単純な定数式である必要があります。 typeは、配列に格納する必要のある要素のNumbaタイプです。配列は現在のスレッドに対してプライベートです。配列のようなオブジェクトが返されます。これは、任意の標準配列と同様に読み取りおよび書き込みが可能です（インデックス付けなど）。

コンスタントメモリ

コンスタントメモリは、読み取り専用、キャッシュ、オフチップのメモリ領域であり、すべてのスレッドからアクセスでき、ホストによって割り当てられます。定数メモリに配列を作成する方法は、次を使用することです。

numba.cuda.const.array_like(arr)

配列のようなarrに基づいて、定数メモリ内の配列を割り当ててアクセスできるようにします。

割り当て解除の動作

このセクションでは、Numbaの内部メモリ管理の割り当て解除動作について説明します。外部メモリ管理プラグインが使用されている場合、割り当て解除の動作が異なる場合があります。この動作は設定で解除することもできます。

すべてのCUDAリソースの割り当て解除は、コンテキストごとに追跡されます。デバイスメモリへの最後の参照が削除されると、基になるメモリの割り当てが解除されるようにスケジュールされます。割り当て解除はすぐには発生しません。保留中の割り当て解除のキューに追加されます。この設計には2つの利点があります。

リソース割り当て解除APIにより、デバイスが同期される場合があります。したがって、非同期実行を中断します。割り当て解除を延期すると、パフォーマンスが重要なコードセクションの遅延を回避できます。

一部の割り当て解除エラーにより、残りのすべての割り当て解除が失敗する場合があります。割り当て解除エラーが続くと、CUDAドライバーレベルで重大なエラーが発生する可能性があります。場合によっては、これはCUDAドライバーのセグメンテーション違反を意味する可能性があります。最悪の場合、これによりシステムGUIがフリーズし、システムをリセットした場合にのみ回復する可能性があります。割り当て解除中にエラーが発生すると、残りの保留中の割り当て解除はキャンセルされます。割り当て解除エラーが報告されます。プロセスが終了すると、CUDAドライバーは終了したプロセスによって割り当てられたすべてのリソースを解放できます。

次のイベントが発生するとすぐに、割り当て解除キューが自動的にフラッシュされます。

メモリ不足エラーのため、割り当てに失敗しました。すべての割り当て解除をフラッシュした後、割り当てが再試行されます。

コードセクションが終了するまで、リソースの割り当て解除を延期することが望ましい場合があります。ほとんどの場合、ユーザーは割り当て解除による暗黙的な同期を避けたいと考えています。これは、次のコンテキストマネージャーを使用して実行できます。

numba.cuda.defer_cleanup()

メモリの割り当て解除を一時的に無効にします。これを使用して、リソースの割り当て解除が非同期実行を中断するのを防ぎます。

code: python

with defer_cleanup():

# すべてのクリーンアップはここで延期される

do_speed_critical_code()

# ここでクリーンアップを実行

このコンテキストマネージャーはネストすることができます。

デバイス関数の記述

CUDAデバイス関数は、カーネルまたは別のデバイス関数によって、デバイス内からのみ呼び出すことができます。

code: python

from numba import cuda

@cuda.jit(device=True)

def a_device_function(a, b):

return a + b

カーネル関数とは異なり、デバイス関数は通常の関数と同じように値を返すことができます。

参考