http://minosys.hateblo.jp/entry/2020/07/09/070407
ゼロから作る Deep Learning❸を読んだ感想を追加しました。
気ままな技術者生活から人生について考える
http://minosys.hateblo.jp/entry/2020/07/09/070407
ゼロから作る Deep Learning❸を読んだ感想を追加しました。
Anaconda の numpy がいつの間にかデフォルト mkl になっていたので、改めて numpy 対 cupy 対決をしてみました。
結果を示すと以下の様になりました。
np(1000000):0:00:16.926656 sec
sum:250000055.953125
cp(1000000):0:00:01.184336 sec
sum:250006660.0
np(100000):0:00:01.684878 sec
sum:24998295.525390625
cp(100000):0:00:00.144552 sec
sum:24996460.0
np(10000):0:00:00.175715 sec
sum:2499707.635986328
cp(10000):0:00:00.143209 sec
sum:2499635.8
1万回ループでは殆ど差異はなく、10万回以上で10倍程度の差になっています。
CPUが今となっては非力な点があり、コア数が多く、周波数が高くAVX512 が使える 9XXX 番台だともっと高い数値が期待できると思います。
ちなみに計算時間の殆どは乱数計算が占めています。複雑な並列計算はGPU圧勝という感じでした。
最後に今回のスクリプトを掲載します。
# -*- coding: utf-8 -*-
# numpy(INTEL mkl edition) vs CuPy
import numpy as np
import cupy as cp
import datetime
# initiation...
cp.cuda.set_allocator(cp.cuda.MemoryPool().malloc)
xp = None
size=100000
loop = 1000
def func(r1, r2):
return xp.dot(r1, r2)
def prepare(func, size, rem, el):
st = datetime.datetime.now()
r1 = xp.random.rand(size).astype(np.float32)
r2 = xp.random.rand(size).astype(np.float32)
rem += func(r1, r2)
el += datetime.datetime.now() - st
return (rem, el)
def measure(amble, size, func, count):
rem = 0.0
el = datetime.timedelta()
for r in range(count):
(rem, el) = prepare(func, size, rem, el)
print (amble + ":" + str(el) + " sec")
print ("sum:" + str(rem))
sizes = (1000000, 100000, 10000)
engines = (('np', np), ('cp', cp))
for s in sizes:
for e in engines:
xp = e[1]
amble = e[0] + "(" + str(s) + ")"
measure(amble, s, func, loop)
「ゼロから作るディープラーニング2」の書評を追加しました。
PWM, ADC, OUTポート、I2C, SIO と次々に切っていっても
電池を消費してしまうため、試しに SHT31 だけを電池に
つないで放置してみた。
2時間経過しても電源電圧に変化なしだったので、
これは TWELite 側の問題だなと確信し、デバイスを別の
ものに変えてみた。
3時間経過後、電源電圧は降下していていなかったので、
おそらく予想が当たっていたのではないだろうか。
1000要素のベクトルのドット積を求めるプログラムで、
numpy, cupy, numba の実行速度を比較してみました。
numba 64bit: 0:00:00.001167 numba 32bit: 0:00:00.000692 cupy 64bit: 0:00:00.000074 cupy 32bit: 0:00:00.000060 numpy 64bit: 0:00:00.000006 numpy 32bit: 0:00:00.000003
となり、numpy 圧勝。cupy が若干遅くなっているのはメイン
メモリとグラフィックメモリの相互移動が発生しているためと
推察します。
numba が著しく遅いのは和の計算に atomic add を使っている
ためでしょう。
(バイナリツリーを用いると劇的に速くなりますが、
面倒なので実装せず。)
最後にテストに用いたプログラムを載せておきます。
# -*- coding: utf-8 -*- from numba import cuda import numpy as np import cupy as cp import datetime @cuda.jit def dotProduct(c, a, b): i = cuda.threadIdx.x + cuda.blockIdx.x * cuda.blockDim.x cuda.atomic.add(c, 0, a[i] * b[i]) N = 1000 a = np.random.rand(N).astype(np.float64) b = np.random.rand(N).astype(np.float64) d = np.zeros(1, dtype=np.float64) a32 = np.random.rand(N).astype(np.float32) b32 = np.random.rand(N).astype(np.float32) d32 = np.zeros(1, dtype=np.float32) ac = cp.asarray(a) bc = cp.asarray(b) dc = cp.asarray(d) ac32 = cp.asarray(a32) bc32 = cp.asarray(b32) dc32 = cp.asarray(d32) threadN = 256 blockN = (N + threadN - 1) // threadN # burn-in dotProduct[blockN, threadN](d, a, b) dotProduct[blockN, threadN](d32, a32, b32) dc = cp.dot(ac, bc) dc32 = cp.dot(ac32, bc32) d_gold = np.dot(a, b) d32 = np.dot(a32, b32) d = np.zeros(1, dtype=np.float64) st = datetime.datetime.now() dotProduct[blockN, threadN](d, a, b) ed = datetime.datetime.now() print("numba 64bit:", ed - st) d32 = np.zeros(1, dtype=np.float32) st = datetime.datetime.now() dotProduct[blockN, threadN](d32, a32, b32) ed = datetime.datetime.now() print("numba 32bit:", ed - st) st = datetime.datetime.now() dc = cp.dot(ac, bc) ed = datetime.datetime.now() print("cupy 64bit:", ed - st) st = datetime.datetime.now() dc32 = cp.dot(ac32, bc32) ed = datetime.datetime.now() print("cupy 32bit:", ed - st) st = datetime.datetime.now() d_gold = np.dot(a, b) ed = datetime.datetime.now() print("numpy 64bit:", ed - st) st = datetime.datetime.now() d_gold32 = np.dot(a32, b32) ed = datetime.datetime.now() print("numpy 32bit:", ed - st)
書評を追加しました。
CUDA サンプルにもある nbody (多体問題) を cupy で解いてみた。
うまくはまれば namba を使うよりも少ない行数で記述できる。
実行速度は namba とほぼ変わらなかった。
初速は銀河面上を円運動するように与えているため、シミュレーション
時間を長くしてもあまり面白い絵になっていませんので、
適宜工夫してみて下さい。
ソースコードは以下の通り。
# -*- coding: utf-8 -*- import datetime import cupy as cp import numpy as np import matplotlib.pyplot as plt # 質点数 N = 1024 REP = 5000 AREAX = 1e3 AREAY = 1e3 FLIP = 0 # 物理定数 cp.cuda.set_allocator(cp.cuda.MemoryPool().malloc) MASS = cp.float32(1.0) M = cp.ones(N, dtype=cp.float32) * MASS G = cp.float32(1.0) EPSILON = cp.float32(1e-6) DT = cp.float32(0.01) # 分布、初速を決定する DAX = cp.zeros(N, dtype=cp.float32) DAY = cp.zeros(N, dtype=cp.float32) def advance(x0, y0, xx1, xx2, yy1, yy2, M, EPSILON): diffx = x0[xx2] - x0[xx1] diffy = y0[yy2] - y0[yy1] r = cp.sqrt(diffx * diffx + diffy * diffy) r3 = r * r * r + EPSILON dax = cp.sum(M * diffx / r3, axis = 1) day = cp.sum(M * diffy / r3, axis = 1) return dax, day def update(x0, y0, vx0, vy0, DAX, DAY, G, DT): vx1 = vx0 + DAX * G * DT vy1 = vy0 + DAY * G * DT x1 = x0 + vx1 * DT y1 = y0 + vy1 * DT return x1, y1, vx1, vy1 def eccentric(x): return (cp.exp(2 - 8 * x) - cp.exp(2)) / (cp.exp(2 - 8) - cp.exp(2)) def init_dist(): global AREAX, AREAY, MASS, G RR = eccentric(cp.random.rand(N)) + 1e-7 THETA = cp.random.rand(N) XS = RR * cp.cos(THETA * 2 * cp.pi) YS = RR * cp.sin(THETA * 2 * cp.pi) VAREA = cp.ones(N, dtype=cp.float32) * cp.sqrt(G/N) / RR RR = cp.random.rand(N) / 4 THETA = cp.random.rand(N) VSX = RR * cp.cos(THETA * 2 * cp.pi) VSY = RR * cp.sin(THETA * 2 * cp.pi) r = cp.sqrt(XS * XS + YS * YS) VSX = VAREA * (-YS + VSX * (1 - r)) VSY = VAREA * (XS + VSY * (1 - r)) XS *= AREAX YS *= AREAY # 重心を求める x0 = cp.sum(XS) / N y0 = cp.sum(YS) / N vx0 = cp.sum(VSX) / N vy0 = cp.sum(VSY) / N # 重心座標系に移行する XS -= x0 YS -= y0 VSX -= vx0 VSY -= vy0 return XS, YS, VSX, VSY # 分布、初速を決定する XS, YS, VSX, VSY = init_dist() X = [XS, cp.zeros(N, dtype=cp.float32)] Y = [YS, cp.zeros(N, dtype=cp.float32)] VX = [VSX, cp.zeros(N, dtype=cp.float32)] VY = [VSY, cp.zeros(N, dtype=cp.float32)] print("intialized") # REP 時間繰り返す x0s = range(N) y0s = range(N) xx1, xx2 = np.meshgrid(x0s, x0s) yy1, yy2 = np.meshgrid(y0s, y0s) st = datetime.datetime.now() for t in range(REP): DAX, DAY = advance(X[FLIP], Y[FLIP], xx1, xx2, yy1, yy2, M, EPSILON) X[1-FLIP], Y[1-FLIP], VX[1-FLIP], VY[1-FLIP] = \ update(X[FLIP], Y[FLIP], VX[FLIP], VY[FLIP], \ DAX, DAY, G, DT) if t % (REP // 50) == 0: print('*', end = '', flush = True) # 結果をグラフで表示 if t == 0: plt.plot(cp.asnumpy(X[FLIP]), cp.asnumpy(Y[FLIP]), 'x', color='red') FLIP = 1 - FLIP ed = datetime.datetime.now() print("\nelapsed time:", ed - st) x = cp.asnumpy(X[FLIP]) y = cp.asnumpy(Y[FLIP]) plt.plot(x, y, '.', color='green') plt.show()
cupy の中の人の投稿によると、cupy.cuda.set_allocator() で
メモリプールを使い回すようにすると高速になると書いてある。
実際そのとおりなのだが、投稿にある test() を2回連続で
回すと、アロケータをセットした方が圧倒的に遅くなった。
ソースコードを見ていないが、どこかに実行結果をキャッシュしている?
numpy: 0:00:00.004890 sec cupy: 0:00:00.000179 sec cupy(memory pool): 0:00:00.000456 sec
なお、記事では触れていないが、いわいる「バーンイン」を
行わないで GPU 回りのメソッドを呼ぶと初期化設定のため、
最初の1回だけは非常に時間がかかることに注意。
使用したテストコードは以下の通り。
# -*- coding: utf-8 -*- import numpy as np import cupy as cp import datetime def test(xp): x = xp.arange(1000000).reshape(1000, -1) return x.T * 2 def test_with_time(title, xp): t1 = datetime.datetime.now() test(xp) t2 = datetime.datetime.now() print(title, t2 - t1, "sec") # バーンイン #cp.arange(1) test(cp) test_with_time("numpy:", np) test_with_time("cupy:", cp) cp.cuda.set_allocator(cp.cuda.MemoryPool().malloc) test_with_time("cupy(memory pool):", cp)
書評を追加しました。