遅いぞ 5000倍の高速化 エンジニアリング。CUDAストリームで、GPUリソースの効率的な自動配分。クーパイ ライブラリで Python like なコーディング。

ショートストーリー: 「東京のプログラマとCUDAストリーム」

東京の高層ビル群の中、オフィスの片隅で一人のプログラマ、リュウイチが黙々とキーボードを叩いていた。彼は技術に対して熱い情熱を持ち、その日も新たなプロジェクトに取り組んでいた。テーマは、巨大なデータセットの高速解析。データサイエンスの世界で彼が目指しているのは、誰よりも早く、そして効率的に結果を出すことだった。

リュウイチはGPUの性能を最大限に引き出すための技術に興味を持っていた。彼が目をつけたのは「CUDAストリーム」だった。これは、GPUリソースを効率的に自動配分し、複数のタスクを並列に処理する技術だった。彼は、この技術を使ってデータ解析の速度を飛躍的に向上させることを目指していた。

その日のプロジェクトは、東京全域の気象データをリアルタイムで解析し、異常気象を予測するというものだった。データ量は膨大で、通常の処理では数時間かかるものだった。リュウイチはストリームを使うことで、複数のデータ解析プロセス（温度解析、湿度解析、風速解析）を並列に実行することにした。

リュウイチはプログラムをセットアップし、ストリームを起動した。ストリームはGPUリソースを効率的に割り当て、各プロセスを最適に並列化した。彼はコンソールに映る進捗を見守りながら、ストリームの力を感じた。ストリームはまるで彼の背後で羽ばたく風のように、データをスムーズに処理していった。

「やった、予測が出た！」リュウイチは興奮しながら叫んだ。予測結果は驚くほど早く出揃い、その精度も非常に高かった。彼はストリームがどれだけ効率的にGPUを使っているかを実感し、満足げに笑みを浮かべた。

その夜、東京の夜景を見下ろしながら、リュウイチは未来の可能性について考えた。ストリームのような先端技術があれば、どんなに複雑な問題でも解決できると彼は信じていた。彼の中には、新たなプロジェクトへの情熱が湧き上がっていた。

「次は何を解析しようか？」リュウイチは静かに呟いた。彼の目は次なる挑戦を見据えて輝いていた。ストリームの力を使って、東京のプログラマはまた新たな冒険に挑もうとしていた。

CuPyを使用して、pythonコードをGPU計算に対応させる例を以下に示します。CuPyはNumPy互換のインターフェースを提供し、NVIDIAのCUDAを利用して高速な数値計算を行うためのライブラリです。

グラフ上の5本の放物線の各点を 全て 同時 並列で計算しています、

実行結果。5本の各放物線の計算を行い、その処理時間と結果をプロットする。

python loop 処理時間: 1.631545 秒

CuPyのカーネルモジュール (起動したスレッド数: 5000) 処理時間: 0.000185 秒

CuPyでベクトル化された計算 処理時間: 0.000754 秒

CuPyでベクトル化された計算 のストリーム最適化。 処理時間: 0.000350 秒

(CUDAストリームは、GPUリソースの効率的な配分を自動で行い、待ち時間の最小化と並列処理の最適化を実現します。)

3つの方法は、それぞれ異なるアプローチで計算を行い、特にGPUの活用方法に違いがあります。それぞれの方法について、以下に簡単に説明します。

CUDAカーネル（CuPyでのカーネル使用）:

概要: カーネル関数を使用してGPUの各スレッドで並列に計算を行う方法。
特長: 各スレッドに対して明示的にインデックスを指定し、細かく制御可能。高度な並列処理を実現できる。
利点: スレッドレベルでの細かい制御が可能で、大量のデータを並列処理するのに適している。
欠点: カーネルの設計とスレッド管理が複雑で、コードの書き方に注意が必要。

通常のPythonループ:

概要: 通常のPythonのforループを使用して計算を逐次的に行う方法。
特長: シンプルで理解しやすいが、GPUの並列処理能力は活かされない。
利点: コードが簡単で、デバッグがしやすい。
欠点: 並列処理が行われないため、計算時間が長くなる可能性が高い。

CuPyのベクトル化:

概要: CuPyのベクトル操作を利用して、配列全体に対する計算を一度に行う方法。
特長: 配列全体に対するベクトル演算を利用することで、GPUの並列計算能力を活用。
利点: コードが簡潔で、CuPyの自動的な並列処理によって高速な計算が可能。
欠点: スレッドインデックスなどの細かい制御ができないため、カスタムカーネルほどの柔軟性はない。

このコードは、Pythonのループを使用して5本の各放物線の計算を行い、その結果をプロットするものです。

import cupy as cp
import numpy as np
import time
import matplotlib.pyplot as plt

# データの準備
width = 1000
num_parabolas = 5
N = width * num_parabolas
a = cp.array([0.5 * i for i in range(1, num_parabolas + 1)], dtype=cp.float32)  # 各放物線の係数
b = cp.zeros(N, dtype=cp.float32)

# GPU計算
start_time = time.time()
for y in range(num_parabolas):
    for x in range(width):
        t = x / float(width - 1)  # 正規化
        a_param = a[y]  # 放物線の係数
        b[y * width + x] = a_param * t * t  # 放物線の計算
end_time = time.time()

# 処理時間の表示
print("処理時間: {:.6f} 秒".format(end_time - start_time))

# 結果の転送 (必要に応じて)
b_cpu = cp.asnumpy(b)

# 結果のプロット
x = np.linspace(0, 1, width)
for i in range(num_parabolas):
    plt.plot(x, b_cpu[i * width:(i + 1) * width], label=f'Parabola {i+1}')

plt.title('5 Parabolas with Different Parameters')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

!pip install cupy

import cupy as cp
import numpy as np
import time
import matplotlib.pyplot as plt

# CUDAカーネルの定義 (CuPy形式)
kernel_code = """
extern "C" __global__ void compute_parabolas(const float *a, float *b, int width) {
    int x = threadIdx.x;
    int y = blockIdx.x; // 5本の放物線を表す　２次元のスレッドインデックスでスレッドを起動。

    if (x < width) {
        float t = x / float(width - 1); // 正規化
        float a_param = a[y]; // 放物線の係数
        b[y * width + x] = a_param * t * t; // 放物線の計算
    }
}
"""

# データの準備
width = 1000
num_parabolas = 5
N = width * num_parabolas
a = np.array([0.5 * i for i in range(1, num_parabolas + 1)], dtype=np.float32)  # 各放物線の係数
b = np.zeros(N, dtype=np.float32)

# CuPyのカーネルモジュールのコンパイル
mod = cp.RawModule(code=kernel_code, backend='nvcc')
compute_parabolas = mod.get_function("compute_parabolas")

# データの転送
a_gpu = cp.array(a)
b_gpu = cp.zeros(N, dtype=cp.float32)

# スレッドとブロックの設定
block_size_x = width
grid_size_y = num_parabolas

# CuPyのカーネルの呼び出し
start_time = time.time()
compute_parabolas(grid=(grid_size_y, 1, 1), block=(block_size_x, 1, 1), args=(a_gpu, b_gpu, np.int32(width)))
cp.cuda.Stream.null.synchronize()
end_time = time.time()

# 結果の転送
b = cp.asnumpy(b_gpu)

# 処理時間とスレッド数の表示
print("処理時間: {:.6f} 秒".format(end_time - start_time))
print("ブロックサイズ: ({}, 1, 1)".format(block_size_x))
print("グリッドサイズ: ({}, 1)".format(grid_size_y))
print("起動したスレッド数: {}".format(grid_size_y * block_size_x))

# 結果のプロット
x = np.linspace(0, 1, width)
for i in range(num_parabolas):
    plt.plot(x, b[i * width:(i + 1) * width], label=f'Parabola {i+1}')

plt.title('5 Parabolas with Different Parameters')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

コードの説明
CUDAカーネル定義: カーネルコード内で、Xはスレッドのインデックス(threadIdx.x)であり、各放物線に沿った計算を担当します。Yはブロックのインデックス(blockIdx.x)であり、異なる放物線を表します。
処理時間の計測: timeモジュールを使用してカーネルの実行前後の時間を計測し、処理時間を表示しています。
スレッド数の計算: 使用したブロックサイズとグリッドサイズを用いて、起動したスレッド数を計算し、表示しています。
このコードは、XとYのスレッドインデックスを使用して各放物線の計算を行い、結果をプロットするものです。

#!pip install cupy

import cupy as cp
import numpy as np
import time
import matplotlib.pyplot as plt

# データの準備
width = 1000
num_parabolas = 5
N = width * num_parabolas
a = cp.array([0.5 * i for i in range(1, num_parabolas + 1)], dtype=cp.float32)  # 各放物線の係数
b = cp.zeros((num_parabolas, width), dtype=cp.float32)  # 各放物線の結果

# ベクトル化された計算

t = cp.linspace(0, 1, width)  # 正規化されたx座標

start_time = time.time()
b = a[:, cp.newaxis] * t**2   # 各放物線の計算　この１行でOK。簡単でシンプルで速い。
end_time = time.time()

# 処理時間の表示
print("処理時間: {:.6f} 秒".format(end_time - start_time))

# 結果の転送 (必要に応じて)
b_cpu = cp.asnumpy(b)

# 結果のプロット
x = np.linspace(0, 1, width)
for i in range(num_parabolas):
    plt.plot(x, b_cpu[i], label=f'Parabola {i+1}')

plt.title('5 Parabolas with Different Parameters')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

重要なポイント CuPyではベクトル化された計算を行うために、配列全体に対して一度に計算を行うことができます。

t = cp.linspace(0, 1, width): これにより、すべてのスレッドで同時に計算するためのベクトルが作成されます。
b = a[:, cp.newaxis] * t2: 各放物線の係数aとt2のベクトルを一度に計算することで、GPUの並列性を活かしています。

import cupy as cp
import numpy as np
import time
import matplotlib.pyplot as plt

# データの準備
width = 1000
num_parabolas = 5
a = cp.array([0.5 * i for i in range(1, num_parabolas + 1)], dtype=cp.float32)  # 各放物線の係数
t = cp.linspace(0, 1, width)  # 正規化されたx座標

# ストリームの作成　で高度に自動で最適化。
stream = cp.cuda.Stream()

# ベクトル化された計算
with stream:
    start_time = time.time()
    b = a[:, cp.newaxis] * t**2  # 各放物線の計算
    stream.synchronize()
    end_time = time.time()

# 処理時間の表示
print("処理時間: {:.6f} 秒".format(end_time - start_time))

# 結果の転送 (必要に応じて)
b_cpu = cp.asnumpy(b)

# 結果のプロット
x = np.linspace(0, 1, width)
for i in range(num_parabolas):
    plt.plot(x, b_cpu[i], label=f'Parabola {i+1}')

plt.title('5 Parabolas with Different Parameters')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

説明
ストリームの作成: cp.cuda.Stream()を使用して、新しいCUDAストリームを作成します。これにより、同じストリームに追加された操作は順番に実行され、異なるストリームに追加された操作は並列に実行される可能性があります。

ストリーム内での計算: with stream:のブロック内で行われる計算は、指定したストリームで実行されます。このブロック内のすべての操作は、streamにより順番に処理されます。

同期と計測: ストリーム内での計算が完了するまで待機するためにstream.synchronize()を使用します。これにより、GPUでの計算が完全に終了した時点での時間を測定します。

参考。

CuPyはNumpy互換のインターフェースを持つため、Numpyと同様の関数が多く提供されています。CuPyは、GPU上で高速な計算を実行するために設計された配列操作ライブラリであり、多くの科学技術計算やデータ解析に使用される関数が含まれています。以下は、CuPyに準備されている代表的な関数のリストです。

基本的な配列操作
cupy.array: 配列の作成
cupy.zeros: 0で初期化された配列の作成
cupy.ones: 1で初期化された配列の作成
cupy.empty: 未初期化の配列の作成
cupy.full: 指定された値で初期化された配列の作成
cupy.arange: 等差数列の配列の作成
cupy.linspace: 等間隔の数値を持つ配列の作成
cupy.eye: 単位行列の作成
配列の形状変更
cupy.reshape: 配列の形状を変更
cupy.transpose: 配列の転置
cupy.ravel: 配列のフラット化
cupy.concatenate: 配列の結合
cupy.split: 配列の分割
数学関数
cupy.add: 配列同士の加算
cupy.subtract: 配列同士の減算
cupy.multiply: 配列同士の乗算
cupy.divide: 配列同士の除算
cupy.power: 累乗
cupy.exp: 指数関数
cupy.log: 自然対数
cupy.sin: 正弦
cupy.cos: 余弦
cupy.tan: 正接
cupy.sqrt: 平方根
cupy.dot: ドット積
cupy.sum: 配列要素の総和
cupy.mean: 平均
cupy.var: 分散
cupy.std: 標準偏差
統計関数
cupy.min: 最小値
cupy.max: 最大値
cupy.argmin: 最小値のインデックス
cupy.argmax: 最大値のインデックス
cupy.median: 中央値
cupy.percentile: パーセンタイル
線形代数
cupy.linalg.inv: 逆行列
cupy.linalg.det: 行列式
cupy.linalg.eig: 固有値・固有ベクトル
cupy.linalg.svd: 特異値分解
cupy.linalg.norm: ノルム
cupy.matmul: 行列積
ランダム数生成
cupy.random.rand: 一様分布に従うランダム数
cupy.random.randn: 標準正規分布に従うランダム数
cupy.random.randint: 整数のランダム数
cupy.random.choice: 要素のランダム選択
FFT（高速フーリエ変換）
cupy.fft.fft: 1次元の高速フーリエ変換
cupy.fft.ifft: 1次元の逆高速フーリエ変換
cupy.fft.fftn: n次元の高速フーリエ変換
cupy.fft.ifftn: n次元の逆高速フーリエ変換
その他
cupy.sort: 配列のソート
cupy.argsort: インデックスのソート
cupy.unique: ユニークな要素の抽出
cupy.where: 条件を満たす要素のインデックス

ベクトル化とは、配列全体に対する要素単位の操作を、一度にまとめて効率的に行うことです。これにより、ループを使用することなく、高速で簡潔なコードを書くことができます。CuPy（およびNumPy）では、配列全体に対して同時に演算を行うことが可能であり、これがベクトル化の基本的なアイデアです。

以下は、先ほどのCuPyを使ったコードのベクトル化について詳しく説明します。

元のPythonループを使った実装

# Pythonでのループを使った例
b = np.zeros(N, dtype=np.float32)

for i in range(num_parabolas):
    a_param = a[i]
    for x in range(width):
        t = x / (width - 1)  # 正規化
        b[i * width + x] = a_param * t * t  # 放物線の計算

この例では、外側のループが異なる放物線を扱い、内側のループが各放物線の各点を計算しています。各点の計算は個別に行われており、計算が逐次的に行われています。

ベクトル化されたCuPyの実装

import cupy as cp

# データの準備
width = 1000
num_parabolas = 5
N = width * num_parabolas
a = cp.array([0.5 * i for i in range(1, num_parabolas + 1)], dtype=cp.float32)
b = cp.zeros(N, dtype=cp.float32)

# ベクトル化された操作
x = cp.arange(width, dtype=cp.float32) / (width - 1)  # 正規化した x の配列
x = x[cp.newaxis, :]  # 2次元に拡張して放物線の係数とブロードキャスト可能にする
a = a[:, cp.newaxis]  # 2次元に拡張して x の配列とブロードキャスト可能にする

# ベクトル化計算: 全ての放物線のすべての x の値を同時に計算
b = (a * x**2).ravel()

ベクトル化の詳細

データの準備:

xは、width個の等間隔な値の配列です。これを放物線の計算に使用するため、0から1に正規化しています。
次元の拡張とブロードキャスト:

x = x[cp.newaxis, :]では、1次元のx配列を2次元に拡張しています。これにより、a配列とのブロードキャストが可能になります。
a = a[:, cp.newaxis]でも同様に、1次元のa配列を2次元に拡張しています。

ブロードキャスト:

CuPyのブロードキャスト機能により、aとxの配列を自動的に拡張してサイズを揃え、要素ごとに計算を行います。この例では、aは各放物線の係数を持ち、xは各点の位置を持ちます。ブロードキャストにより、全ての放物線に対して同時に計算が行われます。

計算:

b = (a * x**2).ravel()により、全ての放物線の全ての点に対して計算を行い、結果をフラット化してbに保存します。

ベクトル化の利点
パフォーマンス: ループを回すことなく、CuPyの内部で最適化された処理により、計算が高速化されます。
シンプルさ: ベクトル化により、コードが簡潔になり、可読性が向上します。
このように、CuPyを使用したベクトル化によって、複数のデータに対する同時並行的な計算が簡単に実現できます。

配列生成の高速化:

CuPyを使って配列を生成する場合、例えばcp.arangeやcp.zerosなどの関数を使用すると、これらの配列はGPUメモリ上に直接作成されます。これにより、CPUからGPUへのデータ転送が不要になり、処理が高速化されます。

二次元配列とブロードキャスト:

配列の次元を拡張することで、複数の配列に対する一括操作が可能です。例えば、1次元の配列を2次元に拡張し、他の配列とブロードキャストすることで、配列の形を自動的に揃えて計算を行います。これにより、多次元データの操作が容易になります。

並列計算:

GPUには多数の計算ユニット（コア）があり、CuPyを使った操作では、これらのコアが並列に動作します。これは特に、要素ごとの計算や行列演算など、同じ操作を複数のデータに対して行う場合に効果的です。たとえば、上記の放物線の計算では、各スレッドが異なるデータポイントの計算を同時に行います。

CUDAストリームは、GPUの並列処理能力を最大限に活用するための手段であり、異なる計算タスクを同時に実行するために利用されます。具体的には、以下のような機能を提供します：

主なポイント
待ちの最小化: ストリームを使用することで、計算タスクが互いにブロックされることなく、非同期に実行されるため、GPUのリソースが無駄になる時間が減少します。

リソースの最適化: ストリームはGPUのコアやメモリ帯域幅を効率的に使用するため、利用可能なリソースが最大限に活用されます。これにより、複数の計算が効率的に進行し、全体の処理時間が短縮されます。

使いやすさ: ストリームを活用することで、開発者は複雑な並列処理の詳細を手動で管理する必要がなくなります。ストリームの管理はCUDAによって自動化されます。

これらの機能により、CUDAストリームを使用することで、GPUの計算能力を最大限に引き出し、計算プロセスを効率的に並列化することができます。