Multi-thread entier de multiplication de matrice dans NumPy/SciPy
faire quelque chose comme
import numpy as np
a = np.random.rand(10**4, 10**4)
b = np.dot(a, a)
utilise plusieurs cœurs, et il fonctionne bien.
les éléments de a sont des flotteurs 64 bits (ou 32 bits dans les plateformes 32 bits?), et j'aimerais multiplier entier 8 bits tableaux. Essayer la suivante, cependant:
a = np.random.randint(2, size=(n, n)).astype(np.int8)
fait en sorte que le produit dot n'utilise pas plusieurs cœurs, et fonctionne donc ~1000x plus lentement sur mon PC.
array: np.random.randint(2, size=shape).astype(dtype)
dtype shape %time (average)
float32 (2000, 2000) 62.5 ms
float32 (3000, 3000) 219 ms
float32 (4000, 4000) 328 ms
float32 (10000, 10000) 4.09 s
int8 (2000, 2000) 13 seconds
int8 (3000, 3000) 3min 26s
int8 (4000, 4000) 12min 20s
int8 (10000, 10000) It didn't finish in 6 hours
float16 (2000, 2000) 2min 25s
float16 (3000, 3000) Not tested
float16 (4000, 4000) Not tested
float16 (10000, 10000) Not tested
je comprends que NumPy utilise BLAS, qui ne supporte pas les entiers, mais si j'utilise les emballages scipy BLAS, c'est à dire.
import scipy.linalg.blas as blas
a = np.random.randint(2, size=(n, n)).astype(np.int8)
b = blas.sgemm(alpha=1.0, a=a, b=a)
le calcul est multi-threaded. Maintenant, blas.sgemm fonctionne avec exactement le même timing que np.dot pour les flotteurs 32, mais pour les non-flotteurs il convertit tout en float32 et les flotteurs de sortie, ce qui est quelque chose np.dot ne fonctionne pas. (En outre, b est maintenant dans F_CONTIGUOUS ordre, qui est un moindre problème).
donc, si je veux faire la multiplication de matrice entière, je dois faire l'un des suivants:
- utilisez Numpy's douloureusement lent
np.dotet soyez heureux que je puisse garder les entiers 8-bit. - utilisez
sgemmde SciPy et utilisez la mémoire 4x. - utilisez
np.float16de Numpy et utilisez seulement la mémoire 2x, avec la mise en garde quenp.dotest beaucoup plus lent sur les réseaux float16 que sur float32 tableaux, plus que int8. - trouve une bibliothèque optimisée pour la multiplication de matrices entières multi-threadées (en fait, Mathematica fait cela, mais je préférerais une solution Python), supportant idéalement des tableaux 1-bit, bien que les tableaux 8-bit soit également très bien... (Je suis en fait visant à faire la multiplication des matrices sur le champ fini Z / 2Z, et je sais que je peux le faire avec Sage , qui est tout à fait pythonique, mais, encore une fois, est il y a quelque chose de strictement Python?)
puis-je suivre l'option 4? Cette bibliothèque existe pas?
clause de non-responsabilité: en fait, J'exécute num Py + MKL, mais j'ai essayé un test similaire sur vanilly num Py, avec des résultats similaires.
2 réponses
- Option 5 - Roll a custom solution: Partition du produit de matrice dans quelques sous-produits et effectuer ceux-ci en parallèle. Cela peut être relativement facile à implémenter avec les modules Python standard. Les sous-produits sont calculés avec
numpy.dot, qui libère le verrou de l'interpréteur global. Ainsi, il est possible d'utiliser threads qui sont relativement légers et peuvent accéder aux tableaux à partir du fil principal pour l'efficacité de la mémoire.
mise en Œuvre:
import numpy as np
from numpy.testing import assert_array_equal
import threading
from time import time
def blockshaped(arr, nrows, ncols):
"""
Return an array of shape (nrows, ncols, n, m) where
n * nrows, m * ncols = arr.shape.
This should be a view of the original array.
"""
h, w = arr.shape
n, m = h // nrows, w // ncols
return arr.reshape(nrows, n, ncols, m).swapaxes(1, 2)
def do_dot(a, b, out):
#np.dot(a, b, out) # does not work. maybe because out is not C-contiguous?
out[:] = np.dot(a, b) # less efficient because the output is stored in a temporary array?
def pardot(a, b, nblocks, mblocks, dot_func=do_dot):
"""
Return the matrix product a * b.
The product is split into nblocks * mblocks partitions that are performed
in parallel threads.
"""
n_jobs = nblocks * mblocks
print('running {} jobs in parallel'.format(n_jobs))
out = np.empty((a.shape[0], b.shape[1]), dtype=a.dtype)
out_blocks = blockshaped(out, nblocks, mblocks)
a_blocks = blockshaped(a, nblocks, 1)
b_blocks = blockshaped(b, 1, mblocks)
threads = []
for i in range(nblocks):
for j in range(mblocks):
th = threading.Thread(target=dot_func,
args=(a_blocks[i, 0, :, :],
b_blocks[0, j, :, :],
out_blocks[i, j, :, :]))
th.start()
threads.append(th)
for th in threads:
th.join()
return out
if __name__ == '__main__':
a = np.ones((4, 3), dtype=int)
b = np.arange(18, dtype=int).reshape(3, 6)
assert_array_equal(pardot(a, b, 2, 2), np.dot(a, b))
a = np.random.randn(1500, 1500).astype(int)
start = time()
pardot(a, a, 2, 4)
time_par = time() - start
print('pardot: {:.2f} seconds taken'.format(time_par))
start = time()
np.dot(a, a)
time_dot = time() - start
print('np.dot: {:.2f} seconds taken'.format(time_dot))
avec cette implémentation j'obtiens une accélération d'environ x4, qui est le nombre physique de cœurs dans ma machine:
running 8 jobs in parallel
pardot: 5.45 seconds taken
np.dot: 22.30 seconds taken
Pourquoi est-il plus rapide d'effectuer une multiplication flottante par matrice flottante que par int? " explique pourquoi les entiers sont si lents: tout d'abord, les CPU ont des pipelines à haut débit en virgule flottante. Deuxièmement, BLAS n'a pas de type entier.
solution: Convertir les matrices en valeurs float32 obtient de grandes vitesses. Comment est 90X speedup sur un MacBook Pro 2015? (En utilisant float64 est la moitié de bon.)
import numpy as np
import time
def timeit(callable):
start = time.time()
callable()
end = time.time()
return end - start
a = np.random.random_integers(0, 9, size=(1000, 1000)).astype(np.int8)
timeit(lambda: a.dot(a)) # ≈0.9 sec
timeit(lambda: a.astype(np.float32).dot(a.astype(np.float32)).astype(np.int8) ) # ≈0.01 sec