Mon problème est le suivant, j'ai un algorithme itératif tel qu'à chaque itération il faut effectuer plusieurs multiplications matrice-matrice dot ( A_i , B_i ), pour i = 1 ... k. Puisque ces multiplications sont effectuées avec le point de Numpy, je sais qu'ils appellent l'implémentation BLAS-3, qui est assez rapide. Le problème est que le nombre d'appels est énorme et que cela s'est avéré être un goulot d'étranglement dans mon programme. Je voudrais minimiser les frais généraux dus à tous ces appels en faisant moins de produits mais avec des matrices plus grandes.
Pour simplifier, considérez que toutes les matrices sont n x n (généralement n n'est pas grand, il est compris entre 1 et 1000). Une façon de contourner mon problème serait de considérer le diag de matrice diagonale de bloc ( A_i ) et d'exécuter le produit ci-dessous.
Ceci est juste un appel à la fonction dot mais maintenant le programme perd beaucoup de temps à effectuer des multiplications avec des zéros. Cette idée ne semble pas fonctionner mais elle donne le résultat [ A_1 B_1 , ..., A_k B_k ], c'est-à-dire tous les produits empilés dans une seule grande matrice .
Ma question est la suivante: existe-t-il un moyen de calculer [ A_1 B_1 , ..., A_k B_k ] avec un seul appel de fonction? Ou plus précisément, comment puis-je calculer ces produits plus rapidement que de faire une boucle de points Numpy?
3 Réponses :
Vous pouvez empiler les tableaux pour avoir la forme (k, n, n) et appeler numpy.matmul ou utiliser l'opérateur @ .
Par exemple,
In [32]: A = np.array([A0, A1, A2])
In [33]: A
Out[33]:
array([[[ 1, 2],
[ 3, 4]],
[[ 1, 2],
[-3, 5]],
[[ 4, 0],
[ 1, 1]]])
In [34]: B = np.array([B0, B1, B2])
In [35]: A @ B
Out[35]:
array([[[-5, 12],
[-9, 28]],
[[ 4, 3],
[-1, 2]],
[[-8, 36],
[-2, 10]]])
Ou, en utilisant @:
In [18]: A0 = np.array([[1, 2], [3, 4]])
In [19]: A1 = np.array([[1, 2], [-3, 5]])
In [20]: A2 = np.array([[4, 0], [1, 1]])
In [21]: B0 = np.array([[1, 4], [-3, 4]])
In [22]: B1 = np.array([[2, 1], [1, 1]])
In [23]: B2 = np.array([[-2, 9], [0, 1]])
In [24]: np.matmul([A0, A1, A2], [B0, B1, B2])
Out[24]:
array([[[-5, 12],
[-9, 28]],
[[ 4, 3],
[-1, 2]],
[[-8, 36],
[-2, 10]]])
merci, cela semble gérer tous les produits en un seul appel. Cette approche est-elle plus rapide?
Modifier
Pour les matrices nxn plus grandes (environ 20) un appel BLAS à partir du code compilé est plus rapide, pour les plus petits matrices Les noyaux Numba ou Cython personnalisés sont généralement plus rapides.
La méthode suivante génère des fonctions point personnalisées pour des formes d'entrée données. Avec cette méthode, il est également possible de bénéficier des optimisations liées au compilateur comme le déroulement de boucle, qui sont particulièrement importantes pour les petites matrices.
Il faut noter que la génération et la compilation d'un noyau prend environ. 1s, assurez-vous donc d'appeler le générateur uniquement si vous devez vraiment le faire.
Fonction du générateur
A=np.random.rand(1000,2,2)
B=np.random.rand(1000,2,2)
X=A@B
X2=np.einsum("xik,xkj->xij",A,B)
X3=dot_22(A,B) #avoid measurig compilation overhead
X4=dot(A,B) #avoid measurig compilation overhead
%timeit X=A@B
#262 µs ± 2.55 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit np.einsum("xik,xkj->xij",A,B,optimize=True)
#264 µs ± 3.22 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit X3=dot_22(A,B)
#5.68 µs ± 27.2 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit X4=dot(A,B)
#9.79 µs ± 61.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
Exemple d'utilisation
import numpy as np
import numba as nb
#Don't use this for larger submatrices
@nb.njit(fastmath=True,parallel=True)
def dot(A,B):
assert A.shape[0]==B.shape[0]
assert A.shape[2]==B.shape[1]
res=np.empty((A.shape[0],A.shape[1],B.shape[2]),dtype=A.dtype)
for ii in nb.prange(A.shape[0]):
for i in range(A.shape[1]):
for j in range(B.shape[2]):
acc=0.
for k in range(B.shape[1]):
acc+=A[ii,i,k]*B[ii,k,j]
res[ii,i,j]=acc
return res
@nb.njit(fastmath=True,parallel=True)
def dot_22(A,B):
assert A.shape[0]==B.shape[0]
assert A.shape[1]==2
assert A.shape[2]==2
assert B.shape[1]==2
assert B.shape[2]==2
res=np.empty((A.shape[0],A.shape[1],B.shape[2]),dtype=A.dtype)
for ii in nb.prange(A.shape[0]):
res[ii,0,0]=A[ii,0,0]*B[ii,0,0]+A[ii,0,1]*B[ii,1,0]
res[ii,0,1]=A[ii,0,0]*B[ii,0,1]+A[ii,0,1]*B[ii,1,1]
res[ii,1,0]=A[ii,1,0]*B[ii,0,0]+A[ii,1,1]*B[ii,1,0]
res[ii,1,1]=A[ii,1,0]*B[ii,0,1]+A[ii,1,1]*B[ii,1,1]
return res
Ancienne réponse
Une autre alternative, mais avec plus de travail à faire, serait de utiliser des implémentations BLAS spéciales, ce qui crée noyaux personnalisés pour de très petites matrices juste à temps et que d'appeler ces noyaux depuis C.
Exemple
A=np.random.rand(1000,2,2) B=np.random.rand(1000,2,2) dot22=gen_dot_nm(2,2,2) X=dot22(A,B) %timeit X3=dot22(A,B) #5.94 µs ± 21.3 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
def gen_dot_nm(x,y,z):
#small kernels
@nb.njit(fastmath=True,parallel=True)
def dot_numba(A,B):
"""
calculate dot product for (x,y)x(y,z)
"""
assert A.shape[0]==B.shape[0]
assert A.shape[2]==B.shape[1]
assert A.shape[1]==x
assert B.shape[1]==y
assert B.shape[2]==z
res=np.empty((A.shape[0],A.shape[1],B.shape[2]),dtype=A.dtype)
for ii in nb.prange(A.shape[0]):
for i in range(x):
for j in range(z):
acc=0.
for k in range(y):
acc+=A[ii,i,k]*B[ii,k,j]
res[ii,i,j]=acc
return res
#large kernels
@nb.njit(fastmath=True,parallel=True)
def dot_BLAS(A,B):
assert A.shape[0]==B.shape[0]
assert A.shape[2]==B.shape[1]
res=np.empty((A.shape[0],A.shape[1],B.shape[2]),dtype=A.dtype)
for ii in nb.prange(A.shape[0]):
res[ii]=np.dot(A[ii],B[ii])
return res
#At square matices above size 20
#calling BLAS is faster
if x>=20 or y>=20 or z>=20:
return dot_BLAS
else:
return dot_numba
J'ai une certaine expérience avec Numba, donc cette réponse est la bienvenue. Donnez-moi juste un peu de temps pour essayer ça. Une petite question: la ligne d'assert accélère-t-elle? Pourquoi? Merci!
@Integral Il peut donner des accélérations car il informe également le compilateur de la disposition exacte de la mémoire. Mais la principale raison est d'éviter les accidents. Il n'y a pas de vérification des limites, donc si vous sortez des limites, Python plantera simplement ou la fonction produira des absurdités (le même comportement qu'une routine C ou Cython avec la vérification des limites désactivée)
Je comprends que votre dot_22 est plus rapide car vous éliminez une boucle en écrivant chaque résultat explicitement. Vous pensez que cette approche fonctionne toujours pour un dot_nn générique? Je veux dire fonctionne dans le sens d'être plus rapide.
@Integral Oui, c'est possible. Avez-vous une large gamme de formes de tableaux différentes ou seulement quelques-unes?
Si vous ne voulez pas perdre de temps à multiplier les zéros, ce que vous voulez vraiment, ce sont des matrices clairsemées. Utilisation des matrices A et B de @WarrenWeckesser:
from scipy import sparse
sparse.block_diag((A0, A1, A2), format = "csr") @ np.concatenate((B0, B1, B2), axis = 0)
Out[]:
array([[-5, 12],
[-9, 28],
[ 4, 3],
[-1, 2],
[-8, 36],
[-2, 10]], dtype=int32)
Ceci est probablement une accélération pour les grandes matrices. Pour les plus petits, @ max9111 a probablement la bonne idée en utilisant numba.
D'autres ont remarqué que les
nappels àdotpeuvent être plus rapides qu'un appel avec un tableaunfois plus grand. Les frais généraux de gestion de la mémoire réduisent les gains de temps liés à la réduction du nombre d'itérations. Une autre façon de le dire, quelques «quelques» itérations sur une tâche complexe peuvent en fait être optimales. Dans ce cas, à moins qu'il n'y ait un code de niveau BLAS capable de gérer le bloc-point, votre itération d'origine pourrait bien être la plus rapide.@hpaulj merci pour ton commentaire. Pensez-vous à donner quelques références sur ce que d'autres ont remarqué? Je suis intéressé.
n x n allant de 1 à 1_000 est un domaine assez vaste. Avec des tableaux n (<20) assez petits, une fonction Cython ou Numba personnalisée est de loin plus rapide que l'opérateur @ (appelant une routine Blas. Pour n plus grand, la routine BLAS sera la plus rapide.