这里Matrix multiplication using hdf5我使用hdf5(pytables)进行大矩阵乘法,但我感到惊讶,因为使用hdf5比使用普通的numpy.dot和在RAM中存储矩阵的速度更快,这种行为的原因是什么?
也许在python中有一些更快的矩阵乘法函数,因为我仍然使用numpy.dot进行小块矩阵乘法。
下面是一些代码:
假设矩阵可以放在RAM中:测试矩阵10×1000×1000。
使用默认的numpy(我认为没有BLAS lib)。 普通的numpy数组在RAM中:时间9.48
如果A,B在RAM中,C在磁盘上:时间1.48
如果A,B,C在磁盘上:时间372.25
如果我使用numpy和MKL,结果是:0.15,0.45,43.5。
结果看起来很合理,但我仍然不明白为什么在第一种情况下块乘法更快(当我们把A,B存储在RAM中时)。
n_row=1000
n_col=1000
n_batch=10
def test_plain_numpy():
A=np.random.rand(n_row,n_col)# float by default?
B=np.random.rand(n_col,n_row)
t0= time.time()
res= np.dot(A,B)
print (time.time()-t0)
#A,B in RAM, C on disk
def test_hdf5_ram():
rows = n_row
cols = n_col
batches = n_batch
#using numpy array
A=np.random.rand(n_row,n_col)
B=np.random.rand(n_col,n_row)
#settings for all hdf5 files
atom = tables.Float32Atom() #if store uint8 less memory?
filters = tables.Filters(complevel=9, complib='blosc') # tune parameters
Nchunk = 128 # ?
chunkshape = (Nchunk, Nchunk)
chunk_multiple = 1
block_size = chunk_multiple * Nchunk
#using hdf5
fileName_C = 'CArray_C.h5'
shape = (A.shape[0], B.shape[1])
h5f_C = tables.open_file(fileName_C, 'w')
C = h5f_C.create_carray(h5f_C.root, 'CArray', atom, shape, chunkshape=chunkshape, filters=filters)
sz= block_size
t0= time.time()
for i in range(0, A.shape[0], sz):
for j in range(0, B.shape[1], sz):
for k in range(0, A.shape[1], sz):
C[i:i+sz,j:j+sz] += np.dot(A[i:i+sz,k:k+sz],B[k:k+sz,j:j+sz])
print (time.time()-t0)
h5f_C.close()
def test_hdf5_disk():
rows = n_row
cols = n_col
batches = n_batch
#settings for all hdf5 files
atom = tables.Float32Atom() #if store uint8 less memory?
filters = tables.Filters(complevel=9, complib='blosc') # tune parameters
Nchunk = 128 # ?
chunkshape = (Nchunk, Nchunk)
chunk_multiple = 1
block_size = chunk_multiple * Nchunk
fileName_A = 'carray_A.h5'
shape_A = (n_row*n_batch, n_col) # predefined size
h5f_A = tables.open_file(fileName_A, 'w')
A = h5f_A.create_carray(h5f_A.root, 'CArray', atom, shape_A, chunkshape=chunkshape, filters=filters)
for i in range(batches):
data = np.random.rand(n_row, n_col)
A[i*n_row:(i+1)*n_row]= data[:]
rows = n_col
cols = n_row
batches = n_batch
fileName_B = 'carray_B.h5'
shape_B = (rows, cols*batches) # predefined size
h5f_B = tables.open_file(fileName_B, 'w')
B = h5f_B.create_carray(h5f_B.root, 'CArray', atom, shape_B, chunkshape=chunkshape, filters=filters)
sz= rows/batches
for i in range(batches):
data = np.random.rand(sz, cols*batches)
B[i*sz:(i+1)*sz]= data[:]
fileName_C = 'CArray_C.h5'
shape = (A.shape[0], B.shape[1])
h5f_C = tables.open_file(fileName_C, 'w')
C = h5f_C.create_carray(h5f_C.root, 'CArray', atom, shape, chunkshape=chunkshape, filters=filters)
sz= block_size
t0= time.time()
for i in range(0, A.shape[0], sz):
for j in range(0, B.shape[1], sz):
for k in range(0, A.shape[1], sz):
C[i:i+sz,j:j+sz] += np.dot(A[i:i+sz,k:k+sz],B[k:k+sz,j:j+sz])
print (time.time()-t0)
h5f_A.close()
h5f_B.close()
h5f_C.close()
当
np.dot
发送到BLAS时float32
、float64
、complex32
或complex64
,并且否则,它默认使用自己的、缓慢的矩阵乘法例程。
检查BLAS链接的说明如下here。简而言之,检查NumPy安装中是否有文件
_dotblas.so
或类似文件。如果有的话,请检查它链接到哪个BLAS库;参考BLAS速度慢,ATLAS速度快,OpenBLAS和特定于供应商的版本(如Intel MKL)甚至更快。注意多线程BLAS实现,因为它们是Python的multiprocessing
的don't play nicely。接下来,通过检查数组的
flags
来检查数据对齐。在1.7.2之前的NumPy版本中,np.dot
的两个参数都应该是C顺序的。在NumPy>;=1.7.2中,这不再像Fortran数组的特殊情况那样重要。如果您的NumPy没有链接到BLAS,可以(简单地)重新安装它,或者(硬)使用SciPy中的BLAS
gemm
(广义矩阵乘法)函数:这看起来很简单,但几乎没有任何错误检查,所以你必须真正知道你在做什么。
相关问题 更多 >
编程相关推荐