Как получить код быстрее, чем numpy.точка для умножения матрицы?
Здесь умножение матриц с помощью hdf5 я использую hdf5 (pytables) для большого умножения матриц, но я был удивлен, потому что с помощью hdf5 он работает еще быстрее, чем с помощью простого numpy.расставьте точки и храните матрицы в оперативной памяти, в чем причина такого поведения?
И, возможно, есть какая-то более быстрая функция для умножения матриц в python, потому что я все еще использую numpy.точка для умножения матрицы малого блока.
Вот некоторый код:
Предположим, что матрицы могут поместиться в ОЗУ: тест на матрица 10*1000 x 1000.
Используя по умолчанию numpy (я думаю, что нет BLAS lib). Простые массивы numpy находятся в оперативной памяти: time 9.48
Если A, B в оперативной памяти, C на диске: время 1.48
Если A, B, C на диске: время 372.25
Если я использую numpy с MKL, то результаты будут: 0.15, 0.45, 43.5.
Результаты выглядят разумными,но я все еще не понимаю, почему в 1-м случае умножение блоков происходит быстрее(когда мы храним A, B в оперативной памяти).
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()
1 ответ:
np.dot
депеши в Блас КогдаВ противном случае он по умолчанию использует свою собственную медленную процедуру умножения матрицы.
- NumPy был скомпилирован для использования BLAS,
- реализация BLAS доступна во время выполнения,
- ваши данные имеют один из типов dtypes
float32
,float64
,complex32
илиcomplex64
, и- данные соответствующим образом выровнены в памяти.
Проверка вашей связи BLAS описаназдесь . Короче говоря, проверьте, есть ли файл
_dotblas.so
или аналогичный в вашей установке NumPy. Если есть, проверьте, с какой библиотекой BLAS он связан; эталонный BLAS медленный, ATLAS быстрый, OpenBLAS и версии для конкретных поставщиков, такие как Intel MKL, еще быстрее. Будьте осторожны с многопоточными реализациями BLAS, поскольку они не играют хорошо С Pythonmultiprocessing
.Затем проверьте выравнивание данных, проверив
flags
ваших массивов. В версиях NumPy до 1.7.2 оба аргумента вnp.dot
должны быть С-приказал. В NumPy >= 1.7.2 это больше не имеет значения, так как были введены специальные случаи для массивов Fortran.>>> X = np.random.randn(10, 4) >>> Y = np.random.randn(7, 4).T >>> X.flags C_CONTIGUOUS : True F_CONTIGUOUS : False OWNDATA : True WRITEABLE : True ALIGNED : True UPDATEIFCOPY : False >>> Y.flags C_CONTIGUOUS : False F_CONTIGUOUS : True OWNDATA : False WRITEABLE : True ALIGNED : True UPDATEIFCOPY : False
Если ваш NumPy не связан с BLAS, либо (легко) переустановите его, либо (трудно) используйте функцию BLAS
gemm
(generalized matrix multiply) из SciPy:>>> from scipy.linalg import get_blas_funcs >>> gemm = get_blas_funcs("gemm", [X, Y]) >>> np.all(gemm(1, X, Y) == np.dot(X, Y)) True
Это выглядит легко, но это вряд ли делает проверку ошибок, так что вы должны действительно знать, что вы делаете.