SSDとは何なのか？ここで言うところのSSDはSolid State Drive(ソリッドステートドライブ)のことではないのであしからず。ここで言うSSDとはSum of Squared Differences(ssd)のことで、 SSD between 2 Numpy arraysと言えば、2つのナンピー(ナンパイが正解だが日本人ならナンピーだろ)配列間の差分二乗和を意味する。今回はこのSSDのnumpyでの表記法とパイキューダ(pycudaはパイクーダと言うらしいが日本人なら当然パイキューダだろ)との速度比較をやる。

SSD Numpyバージョン¶

先ずはこのサイトとこのサイトを参考にしてnumpy版のSSDのコードを書く。

import numpy as np

a = np.random.randn(5)
b = np.random.randn(5)
print(a)
print(b)
np.sum((a-b)**2)

[ 0.50086932  0.36577004  0.02242422 -0.03547456  0.10152153]
[ 1.6365283  -0.37077582 -0.04609223  2.0547539   1.23728546]

7.495930359813565

np.sum((b-a)**2)

7.495930359813565

a = np.asarray([[1,2,4],
       [3,1,2]])
b = np.asarray([[2,1,1],
       [3,2,3],
       [4,1,2],
       [2,2,1],])

np.sum((b-a)**2)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-20-064a97fa313d> in <module>()
      6        [2,2,1],])
      7 
----> 8 np.sum((b-a)**2)

ValueError: operands could not be broadcast together with shapes (4,3) (2,3)

a = np.asarray([[1,2,4],
       [3,1,2]])
b = np.asarray([[2,1,1],
       [3,2,3],
       [4,1,2],
       [2,2,1],])

np.sum((b[:,None] - a)**2)

a = np.asarray([[1,2,4],
       [3,1,2]])
b = np.asarray([[2,1,1],
       [3,2,3],
       [4,1,2],
       [2,2,1],])

np.square(b[:,None] - a).sum(axis=2).T

array([[11,  5, 14, 10],
       [ 2,  2,  1,  3]])

np.square(b[:,None] - a).sum(axis=1).T

array([[ 2,  4, 10,  2],
       [ 1,  1,  1,  1],
       [10,  2,  4, 10]])

np.square(b[:,None] - a).sum(axis=0).T

array([[15,  3],
       [ 2,  2],
       [23,  3]])

c = np.square(b[:,None] - a).sum(axis=0).T
np.sum(c)

PyCUDA vs. Numpy¶

このサイトからコードを拝借して速度比較をする。

import pycuda.gpuarray as gpuarray
import pycuda.reduction as reduction
import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule
import numpy
import time

krnl = reduction.ReductionKernel(numpy.float32, neutral="0",
        reduce_expr="a+b", map_expr="x[i]*y[i]",
        arguments="float *x, float *y")
ssd = reduction.ReductionKernel(numpy.float32, neutral="0",
        reduce_expr="a+b", map_expr="(x[i] - y[i])*(x[i] - y[i])",
        arguments="float *x, float *y")

for i in range(5):
    a = numpy.random.randn(140**i)
    b = numpy.random.randn(140**i)

    a_gpu = gpuarray.to_gpu(a.astype(numpy.float32))
    b_gpu = gpuarray.to_gpu(b.astype(numpy.float32))

    start = time.time()
    numpy_dot = numpy.dot(a,b)
    numpy_ssd = numpy.sum((a - b) ** 2)
    end = time.time()
    dt = end - start

    print ("CPU time", dt)
    print ("numpy_dot", numpy_dot)
    print ("numpy_euclid", numpy_ssd)

    start = time.time()
    my_dot_prod = krnl(a_gpu, b_gpu).get()
    my_euclid = ssd(a_gpu, b_gpu).get()
    end = time.time()

    dt1 = end - start
    print (("GPU time", dt1))
    print (("my dot product", my_dot_prod))
    print (("my euclid", my_euclid))
    print('Speed difference: {:0.1f}x'.format(dt / dt1))
    print (("\n"))

CPU time 5.2928924560546875e-05
numpy_dot -0.7079061558881492
numpy_euclid 4.496492793847235
('GPU time', 0.0003001689910888672)
('my dot product', array(-0.7079062, dtype=float32))
('my euclid', array(4.4964924, dtype=float32))
Speed difference: 0.2x


CPU time 2.4318695068359375e-05
numpy_dot -12.117062199760833
numpy_euclid 258.74412697331815
('GPU time', 0.00030350685119628906)
('my dot product', array(-12.1170635, dtype=float32))
('my euclid', array(258.74414, dtype=float32))
Speed difference: 0.1x


CPU time 9.465217590332031e-05
numpy_dot 71.99858292120105
numpy_euclid 39102.06162600858
('GPU time', 0.008191585540771484)
('my dot product', array(71.998604, dtype=float32))
('my euclid', array(39102.062, dtype=float32))
Speed difference: 0.0x


CPU time 0.014576911926269531
numpy_dot -873.7404364995434
numpy_euclid 5484635.181583515
('GPU time', 0.001306295394897461)
('my dot product', array(-873.7409, dtype=float32))
('my euclid', array(5484635., dtype=float32))
Speed difference: 11.2x


CPU time 1.6815850734710693
numpy_dot 2721.003087284951
numpy_euclid 768284025.8217536
('GPU time', 0.03971982002258301)
('my dot product', array(2721.0151, dtype=float32))
('my euclid', array(7.68284e+08, dtype=float32))
Speed difference: 42.3x

配列が小さいとCPUが圧倒的に速いが、配列がでかくなるに従ってGPUの独壇場になる。しかし、巨大なアレイになると途端にメモリ不足を起こすので、やはりメインメモリは64G、GPUメモリは16Gは欲しいところだ。