J'ai implémenté une fonction pour appliquer np.dot aux blocs qui sont explicitement lus dans la mémoire centrale à partir de la matrice mappée en mémoire :
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
J'ai ensuite procédé à une analyse comparative de mes blockwise_dot à la fonction normale np.dot appliquée directement à un tableau mappé en mémoire (voir ci-dessous pour le script d'évaluation comparative). J'utilise numpy 1.9.0.dev-205598b lié à OpenBLAS v0.2.9.rc1 (compilé à partir des sources). La machine est un ordinateur portable quad-core fonctionnant sous Ubuntu 13.10, avec 8 Go de RAM et un SSD, et j'ai désactivé le fichier swap.
Résultats
Comme l'a prédit @Bi Rico, le temps nécessaire pour calculer le produit scalaire est magnifique. O(n) par rapport aux dimensions de A . Fonctionnement sur des blocs en cache de A permet d'améliorer considérablement les performances par rapport à l'appel de la fonction normale np.dot sur l'ensemble du tableau de la mémoire :
![enter image description here]()
Il est étonnamment insensible à la taille des blocs traités - il y a très peu de différence entre le temps pris pour traiter le tableau en blocs de 1 Go, 2 Go ou 4 Go. J'en conclus que quelle que soit la mise en cache np.memmap que les tableaux implémentent nativement, il semble être très sous-optimal pour le calcul des produits scalaires.
Autres questions
C'est toujours un peu pénible de devoir mettre en œuvre manuellement cette stratégie de mise en cache, puisque mon code devra probablement être exécuté sur des machines ayant des quantités différentes de mémoire physique, et potentiellement des systèmes d'exploitation différents. Pour cette raison, je suis toujours intéressé de savoir s'il existe des moyens de contrôler le comportement de la mise en cache des tableaux mappés en mémoire afin d'améliorer les performances des applications suivantes np.dot .
J'ai remarqué un comportement étrange en matière de gestion de la mémoire lors de l'exécution des tests de référence. np.dot sur l'ensemble de A Je n'ai jamais vu la taille de l'ensemble résident de mon processus Python dépasser environ 3,8 Go, même si j'ai environ 7,5 Go de RAM libre. Cela m'amène à penser qu'il existe une limite imposée à la quantité de mémoire physique qu'un processus Python peut utiliser. np.memmap J'avais précédemment supposé qu'il utiliserait toute la RAM que le système d'exploitation lui permettrait d'utiliser. Dans mon cas, il pourrait être très bénéfique de pouvoir augmenter cette limite.
Est-ce que quelqu'un a des informations supplémentaires sur le comportement de la mise en cache de np.memmap des tableaux qui permettraient d'expliquer cela ?
Benchmarking 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)
