Bonjour les Zér0s!

Aujourd'hui j'ai un petit problème de performances avec Python, et plus particulièrement avec Numpy. Je dois effectuer 15 produits matriciels à la suite (jusque là aucun problème), mais avec des matrices de taille plutôt imposante: (5278, 49013) x (49013, 5278).

Le second problème, et non des moindre, c'est que ces matrices ne sont pas des matrices creuses, donc impossible d'optimiser le calcul avec le module sparse fourni par Scipy.

Après quelques rapides benchmark, je trouve des chiffres plutôt énormissimes:

import numpy.random as random
import numpy as np

# Liste de matrices de différentes tailles
L = [random.randint(0, 10000, size=(i, 5000)) for i in [2**j for j in range(1, 12)]]

np.dot(L[0], L[0].T)   #25.4 µs
np.dot(L[1], L[1].T)   #96.2 µs
np.dot(L[2], L[2].T)   #379 µs
np.dot(L[3], L[3].T)   #1.52 ms
np.dot(L[4], L[4].T)   #6.1 ms
np.dot(L[5], L[5].T)   #24.7 ms
np.dot(L[6], L[6].T)   #138 ms
np.dot(L[7], L[7].T)   #530 ms
np.dot(L[8], L[8].T)   #2.12 s
np.dot(L[9], L[9].T)   #8.71 s
np.dot(L[10], L[10].T) #34.6 s

Si on regarde un peu, a chaque ligne la dimension de la matrice double, et le temps quadruple (à peu près). En voyant la taille de mes matrices, et en faisant un petit calcul avec ma conjecture précédente:

L[10] fait (2048, 5000) donc...

• L[11] fait (4096, 5000) et s'exécute en 34.6*(4^1) = 138.4s = 2m18s
• ....
• L[15] fait (65536, 5000) et s'exécute en 34.6*(4^5) = 35430.4s = 590m30s soit presque 10h *_*

10h pour un calcul matriciel, sachant que je dois en réaliser 15... N'y as-t'il pas un moyen d'optimiser tout ça?

Non je crains que tu vas devoir trouver autre chose pour résoudre ton problème. J'ai trouvé cet article sur Stackoverflow où l'auteur du message suggère de trouver des méthodes pour rendre ta matrice creuse et donc pouvoir optimiser les calculs.

C'est vrai que c'est une belle façon de procéder!

Par contre, je me suis rendu compte (trop tard :p) que mon point de repère était en fait une machine de guerre... En fait je lançais l'algorithme en Matlab, et le calcul matriciel ne durait jamais plus de 15s, mais j'avais 'oublié' que derrière ce n'était pas mon PC qui faisait tourner Matlab, mais un monstre Du coup je n'ai pas encore pu tester (pas de droit d'admin sur la machine, pas d'admin dans les parages, et pas de numpy installé), mais je pense que mon problème est résolu parce que l'opération devrait durer beaucoup moins longtemps sur une machine digne de ce nom (et pas sur le vieux Centrino 2 que j'ai ><)

Merci pour ta réponse, le lien était très intéressant

En fait Matlab est quand même bien plus rapide, même avec mon petit PC... Voici les tests que j'ai réalisé sur le même PC (mais sous Windows XP, les tests précédents de Python et Numpy étaient sous Ubuntu Mate 15.04):

Code lancé (Matlab):

for k=1:13
    Y = rand(2^k, 5000);
    disp(size(Y))
    tic;
    Y'*Y;
    toc;
end

Output:

           2        5000

Elapsed time is 0.496342 seconds.
           4        5000

Elapsed time is 0.454630 seconds.
           8        5000

Elapsed time is 0.455890 seconds.
          16        5000

Elapsed time is 0.460292 seconds.
          32        5000

Elapsed time is 0.489033 seconds.
          64        5000

Elapsed time is 0.565867 seconds.
         128        5000

Elapsed time is 0.711033 seconds.
         256        5000

Elapsed time is 1.068826 seconds.
         512        5000

Elapsed time is 1.696833 seconds.
        1024        5000

Elapsed time is 2.983439 seconds.
        2048        5000

Elapsed time is 5.631079 seconds.
        4096        5000

Elapsed time is 10.886500 seconds.
        8192        5000

Elapsed time is 22.431605 seconds.

Il y a quand même une sacrée différence de temps, et en plus de ça une différence qui se creuse énormément! Comment ça se fait?

Voici de quoi te divertir : Linear Algebra on Large Arrays

Est-ce que je peux savoir, si tu t'en souviens encore, les mots clés qui ont servis à ta recherche? Parce que je n'avait pas trouvé cette page durant mes recherches alors que c'était exactement ce que je cherchais

Donc, voici un récapitulatif de ce que j'ai pu trouver, pour les éventuelles personnes qui se frotteront au même problème que moi:

I. Rangement des tableaux dans Numpy:

Vous pouvez aller vous renseigner sur:

• La doc Scipy
• Le site donnée par mon voisin du dessus, section Linear Algebra on Large Arrays.

Pour les anglophobes, voici un petit résumé de ce qui est dit dans ces liens (surtout le second lien en fait :p):

Numpy stocke ses tableaux par défaut dans un style C, c'est à dire qu'il stocke les lignes d'un tableau dans des lieux de mémoire adjacents. C'est à dire que tous les éléments d'une même (et unique) ligne seront stockés sur un même bloc de mémoire, sans 'trous' au milieu.

Il est possible de changer ce comportement pour stocker un tableau donné en style Fortran, qui est l'analogue du C: le style Fortran stocke les colonnes dans des cases mémoire adjacentes.

Ce qu'il faut en plus savoir, c'est qu'accéder successivement à des zones mémoires adjacentes est beaucoup plus rapide que d'accéder à deux cases mémoires aléatoires (voir cache miss: résumé et plus complet).

Pour connaître le rangement de vos tableaux:

>>> A.flags
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False

Si on transpose A, alors le rangement de la matrice A va changer:

>>> A.T.flags
  C_CONTIGUOUS : False
  F_CONTIGUOUS : True
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False

II. Produit matriciel Numpy:

Maintenant qu'on sait qu'on peut ordonner nos tableaux de façons différentes, voyons les performances:

import numpy.random as r

A = r.randint(0, 10000, size=(1000, 1000))
# A est en C_CONTIGUOUS
B = A.T
# B est en F_CONTIGUOUS

np.dot(A, A.T)   #1.61s
np.dot(B, B.T)   #7.08s

Déjà ça fais un énorme différence. Maintenant il est possible que vous n'ayez pas des temps très bons, essayez:

>>> np.dot.__module__
'numpy.core._dotblas'   # Bon
'numpy.core.mutliarray' # Pas bon

ou encore:

>>> np.__config__.show()
# Résultat chez moi:
blas_opt_info:
    library_dirs = ['/usr/lib']
    libraries = ['openblas']
    language = f77
openblas_info:
    library_dirs = ['/usr/lib']
    libraries = ['openblas']
    language = f77
openblas_lapack_info:
    library_dirs = ['/usr/lib']
    libraries = ['openblas']
    language = f77
blas_mkl_info:
  NOT AVAILABLE
lapack_opt_info:
    library_dirs = ['/usr/lib']
    libraries = ['openblas']
    language = f77

Le mieux étant de passer par les méthodes de calcul bas niveau réalisées en Fortran par des gens (très très très très) bons:

import scipy.linalg.blas as b # Méthodes en Fortran
import numpy.random as r

A = r.randint(0, 10000, size=(1000, 1000))
# A est en C_CONTIGUOUS
B = A.T
# B est en F_CONTIGUOUS

# Double précision
b.dgemm(alpha=1.O, a=A, b=A, trans_a=True)   #146ms
b.dgemm(alpha=1.0, a=B, b=B, trans_a=True)   #154ms
# Simple précision
b.sgemm(alpha=1.O, a=A, b=A, trans_a=True)   #76.6ms
b.sgemm(alpha=1.0, a=B, b=B, trans_a=True)   #71.9ms

Attention, si vous n'avez pas installé BLAS/LAPACK avant Scipy (compilées avec un compilateur Fortran), alors votre Scipy ne pourra pas utiliser ces fonctions.

Toutes ces améliorations se ressentent aussi énormément sur une matrice plus grande:

import scipy.linalg.blas as b # Méthodes en Fortran
import numpy.random as r

A = r.randint(0, 10000, size=(5000, 5000))
# A est en C_CONTIGUOUS
B = A.T
# B est en F_CONTIGUOUS

# Double précision
b.dgemm(alpha=1.0, a=A, b=A, trans_a=True)   #16.4s
b.dgemm(alpha=1.0, a=B, b=B, trans_a=True)   #16s

Normalement la méthode était censée marcher plus vite avec B qui était déjà aligné en Fortran, mais ça ne semble pas le cas.

-
Edité par Anonyme 20 mai 2015 à 11:38:01