Ответ 1
Я реализовал функцию для применения np.dot к блокам, которые явно считываются в память ядра из массива, сопоставленного с памятью:
import numpy as np
def _block_slices(dim_size, block_size):
"""Generator that yields slice objects for indexing into
sequential blocks of an array along a particular axis
"""
count = 0
while True:
yield slice(count, count + block_size, 1)
count += block_size
if count > dim_size:
raise StopIteration
def blockwise_dot(A, B, max_elements=int(2**27), out=None):
"""
Computes the dot product of two matrices in a block-wise fashion.
Only blocks of `A` with a maximum size of `max_elements` will be
processed simultaneously.
"""
m, n = A.shape
n1, o = B.shape
if n1 != n:
raise ValueError('matrices are not aligned')
if A.flags.f_contiguous:
# prioritize processing as many columns of A as possible
max_cols = max(1, max_elements / m)
max_rows = max_elements / max_cols
else:
# prioritize processing as many rows of A as possible
max_rows = max(1, max_elements / n)
max_cols = max_elements / max_rows
if out is None:
out = np.empty((m, o), dtype=np.result_type(A, B))
elif out.shape != (m, o):
raise ValueError('output array has incorrect dimensions')
for mm in _block_slices(m, max_rows):
out[mm, :] = 0
for nn in _block_slices(n, max_cols):
A_block = A[mm, nn].copy() # copy to force a read
out[mm, :] += np.dot(A_block, B[nn, :])
del A_block
return out
Затем я сделал некоторый бенчмаркинг, чтобы сравнить мою функцию blockwise_dot с нормальной функцией np.dot, примененной непосредственно к массиву с отображением памяти (см. ниже для бенчмаркинга script). Я использую numpy 1.9.0.dev-205598b, связанный с OpenBLAS v0.2.9.rc1 (скомпилирован из источника). Машина представляет собой четырехъядерный ноутбук под управлением Ubuntu 13.10 с 8 ГБ оперативной памяти и SSD, и я отключил файл подкачки.
Результаты
Как предсказал @Bi Rico, время, затраченное на вычисление точечного произведения, красиво O (n) относительно размеров A. Работа с кэшированными блоками A дает огромное улучшение производительности только при вызове нормальной np.dot функции во всем массиве с отображением памяти:
Это удивительно нечувствительно к размеру обрабатываемых блоков - очень мало различий между временем, затрачиваемым на обработку массива в блоках 1 ГБ, 2 ГБ или 4 ГБ. Я пришел к выводу, что все кэширование массивов np.memmap изначально реализуется, оно кажется очень субоптимальным для вычисления точечных продуктов.
Дополнительные вопросы
По-прежнему очень сложно вручную реализовать эту стратегию кэширования, поскольку мой код, вероятно, придется запускать на компьютерах с разным объемом физической памяти и потенциально разных операционных системах. По этой причине меня все еще интересует, существуют ли способы управления кэшированием поведения массивов с отображением памяти, чтобы повысить производительность np.dot.
Я заметил некоторое нечетное поведение обработки памяти, когда я запускал тесты - когда я вызывал np.dot во всем A, я никогда не видел, чтобы размер резидентного набора моего процесса Python превышал примерно 3,8 ГБ, хотя у меня есть около 7,5 ГБ оперативной памяти. Это заставляет меня подозревать, что существует ограничение на количество физической памяти, которой может занимать массив np.memmap - я ранее предполагал, что он будет использовать любую оперативную память, которую ОС позволяет ей захватывать. В моем случае было бы очень полезно увеличить этот предел.
Есть ли у кого-нибудь дальнейшее понимание поведения кеширования массивов np.memmap, которые помогут объяснить это?
Бенчмаркинг script
def generate_random_mmarray(shape, fp, max_elements):
A = np.memmap(fp, dtype=np.float32, mode='w+', shape=shape)
max_rows = max(1, max_elements / shape[1])
max_cols = max_elements / max_rows
for rr in _block_slices(shape[0], max_rows):
for cc in _block_slices(shape[1], max_cols):
A[rr, cc] = np.random.randn(*A[rr, cc].shape)
return A
def run_bench(n_gigabytes=np.array([16]), max_block_gigabytes=6, reps=3,
fpath='temp_array'):
"""
time C = A * B, where A is a big (n, n) memory-mapped array, and B and C are
(n, o) arrays resident in core memory
"""
standard_times = []
blockwise_times = []
differences = []
nbytes = n_gigabytes * 2 ** 30
o = 64
# float32 elements
max_elements = int((max_block_gigabytes * 2 ** 30) / 4)
for nb in nbytes:
# float32 elements
n = int(np.sqrt(nb / 4))
with open(fpath, 'w+') as f:
A = generate_random_mmarray((n, n), f, (max_elements / 2))
B = np.random.randn(n, o).astype(np.float32)
print "\n" + "-"*60
print "A: %s\t(%i bytes)" %(A.shape, A.nbytes)
print "B: %s\t\t(%i bytes)" %(B.shape, B.nbytes)
best = np.inf
for _ in xrange(reps):
tic = time.time()
res1 = np.dot(A, B)
t = time.time() - tic
best = min(best, t)
print "Normal dot:\t%imin %.2fsec" %divmod(best, 60)
standard_times.append(best)
best = np.inf
for _ in xrange(reps):
tic = time.time()
res2 = blockwise_dot(A, B, max_elements=max_elements)
t = time.time() - tic
best = min(best, t)
print "Block-wise dot:\t%imin %.2fsec" %divmod(best, 60)
blockwise_times.append(best)
diff = np.linalg.norm(res1 - res2)
print "L2 norm of difference:\t%g" %diff
differences.append(diff)
del A, B
del res1, res2
os.remove(fpath)
return (np.array(standard_times), np.array(blockwise_times),
np.array(differences))
if __name__ == '__main__':
n = np.logspace(2,5,4,base=2)
standard_times, blockwise_times, differences = run_bench(
n_gigabytes=n,
max_block_gigabytes=4)
np.savez('bench_results', standard_times=standard_times,
blockwise_times=blockwise_times, differences=differences)