Notes de cours Python scientifique

Optimisation du code Python

« L'optimisation prématurée est la racine de tout Mal » , Donald Knuth

Ce chapitre traite des stratégies pour rendre le code Python plus rapide.

Commentez Donner une note  l'article (5)

Article lu   fois.

Les deux auteur et traducteur

Traducteur : Profil Pro

Liens sociaux

Viadeo Twitter Facebook Share on Google+   

I. Optimisation du flux de travail

  1. Faites-le fonctionner : écrivez le code d'une manière simple et lisible.
  2. Faites-le fonctionner de manière fiable : écrivez des scénarios de test automatisés, assurez-vous vraiment que votre algorithme soit correct et que, si celui-ci dysfonctionne, les tests captureront ce dysfonctionnement.
  3. Optimisez le code en définissant des cas d'utilisation simples pour identifier les goulots d'étranglement et accélérez-les, en recherchant un meilleur algorithme ou une meilleure implémentation. Gardez à l'esprit qu'un compromis doit être trouvé entre le profilage sur un exemple réaliste et la simplicité et la rapidité d'exécution du code. Pour un travail efficace, il est préférable de travailler avec des cycles de profilage d’environ 10 secondes.

II. Profilage du code Python

Pas d’optimisation sans mesure !

  • mesure : profilage, chronométrage ;
  • vous aurez des surprises : le code le plus rapide n’est pas toujours celui que vous pensez.

II-A. Timeit

Dans IPython, utilisez timeit pour mesurer la durée des opérations élémentaires :

 
Sélectionnez
In [1]: import numpy as npIn [2]: a = np.arange(1000)
In [3]: %timeit a ** 2
100000 loops, best of 3: 5.73 us per loop
In [4]: %timeit a ** 2.1
1000 loops, best of 3: 154 us per loop
In [5]: %timeit a * a
100000 loops, best of 3: 5.56 us per loop

Utilisez ceci pour guider votre choix entre les différentes stratégies

Pour les appels de longue durée, utilisez %time à la place de %timeit ; c'est moins précis, mais plus rapide.

II-B. Profileur

Utile si vous avez un gros programme à profiler, par exemple le fichier suivant :

 
Sélectionnez
#Pour que cet exemple fonctionne , vous aurez aussi besoin du fichier 'ica.py'

import numpy as np
from scipy import linalg
from ica import fastica
def test():
    data = np.random.random((5000, 100))
    u, s, v = linalg.svd(data)
    pca = np.dot(u[:, :10].T, data)
    results = fastica(pca.T, whiten=False)
if __name__ == '__main__':
    test()

Il s’agit d’une combinaison de deux techniques d’apprentissage non supervisées, l’analyse en composantes principales (PCA) et l’analyse en composantes indépendantes (ICA). L’ACP est une technique de réduction de dimensionnalité, c’est-à-dire un algorithme permettant d’expliquer la variance observée dans vos données en utilisant moins de dimensions. ACI est une technique de séparation de source permettant par exemple de mélanger plusieurs signaux enregistrés au moyen de plusieurs capteurs. Faire une ACP en premier, puis une ACI peut être utile si vous avez plus de capteurs que de signaux. Pour plus d'informations, voir : the FastICA example from scikits-learn.

Pour l'exécuter, vous devez également télécharger le module ICA. En IPython, nous pouvons chronométrer le script:

 
Sélectionnez
In [1]: %run -t demo.py
IPython CPU timings (estimated):
    User  :    14.3929 s.
    System:   0.256016 s.

et le profiler :

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
In [2]: %run -p demo.py

      916 function calls in 14.551 CPU seconds

Ordered by: internal time

ncalls  tottime  percall  cumtime  percall filename:lineno (function)
     1   14.457   14.457   14.479   14.479 decomp.py:849 (svd)
     1    0.054    0.054    0.054    0.054 {method 'random_sample' of 'mtrand.RandomState' objects}
     1    0.017    0.017    0.021    0.021 function_base.py:645 (asarray_chkfinite)
    54    0.011    0.000    0.011    0.000 {numpy.core._dotblas.dot}
     2    0.005    0.002    0.005    0.002 {method 'any' of 'numpy.ndarray' objects}
     6    0.001    0.000    0.001    0.000 ica.py:195 (gprime)
     6    0.001    0.000    0.001    0.000 ica.py:192 (g)
    14    0.001    0.000    0.001    0.000 {numpy.linalg.lapack_lite.dsyevd}
    19    0.001    0.000    0.001    0.000 twodim_base.py:204 (diag)
     1    0.001    0.001    0.008    0.008 ica.py:69 (_ica_par)
     1    0.001    0.001   14.551   14.551 {execfile}
   107    0.000    0.000    0.001    0.000 defmatrix.py:239 (__array_finalize__)
     7    0.000    0.000    0.004    0.001 ica.py:58 (_sym_decorrelation)
     7    0.000    0.000    0.002    0.000 linalg.py:841 (eigh)
   172    0.000    0.000    0.000    0.000 {isinstance}
     1    0.000    0.000   14.551   14.551 demo.py:1 (<module>)
    29    0.000    0.000    0.000    0.000 numeric.py:180 (asarray)
    35    0.000    0.000    0.000    0.000 defmatrix.py:193 (__new__)
    35    0.000    0.000    0.001    0.000 defmatrix.py:43 (asmatrix)
    21    0.000    0.000    0.001    0.000 defmatrix.py:287 (__mul__)
    41    0.000    0.000    0.000    0.000 {numpy.core.multiarray.zeros}
    28    0.000    0.000    0.000    0.000 {method 'transpose' of 'numpy.ndarray' objects}
     1    0.000    0.000    0.008    0.008 ica.py:97 (fastica)
     ...

Clairement, le SVD (dans decomp.py) est ce qui prend le plus de temps, autrement dit, c’est le goulot d’étranglement. Nous devons trouver un moyen d'accélérer cette étape ou de l'éviter (optimisation algorithmique). Passer du temps sur le reste du code est inutile.

II-C. Profilage en dehors d'IPython, en exécutant `` cProfile``

Un profilage similaire peut être effectué en dehors d’IPython, en appelant simplement les profileurs Python intégrés cProfile et profil.

 
Sélectionnez
1.
$  python -m cProfile -o demo.prof demo.py

L'utilisation du commutateur -o générera les résultats du profileur dans le fichier demo.prof à afficher avec un outil externe. Cela peut être utile si vous souhaitez traiter la sortie du profileur avec un outil de visualisation.

II-D. Profileur de ligne

Le profileur nous indique quelle fonction prend le plus de temps, mais pas où elle est appelée.

Pour cela, nous utilisons line_profiler : dans le fichier source, nous décorons quelques fonctions que nous voulons inspecter avec @profile (inutile de l'importer)

 
Sélectionnez
1.
2.
3.
4.
5.
6.
@profile
def test():
    data = np.random.random((5000, 100))
    u, s, v = linalg.svd(data)
    pca = np.dot(u[:, :10], data)
    results = fastica(pca.T, whiten=False)

Ensuite, nous exécutons le script à l'aide du programme kernprof.py , avec les commutateurs -l, --line-by-line et -v, --view pour utiliser le profileur ligne par ligne et afficher les résultats en plus de les enregistrer :

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
$ kernprof.py -l -v demo.py

Wrote profile results to demo.py.lprof
Timer unit: 1e-06 s

File: demo.py
Function: test at line 5
Total time: 14.2793 s

Line #      Hits         Time  Per Hit   % Time  Line Contents
=========== ============ ===== ========= ======= ==== ========
    5                                           @profile
    6                                           def test():
    7         1        19015  19015.0      0.1      data = np.random.random((5000, 100))
    8         1     14242163 14242163.0   99.7      u, s, v = linalg.svd(data)
    9         1        10282  10282.0      0.1      pca = np.dot(u[:10, :], data)
   10         1         7799   7799.0      0.1      results = fastica(pca.T, whiten=False)

Le SVD utilise tout le temps. Nous devons optimiser cette ligne.

III. Rendre le code plus rapide

Une fois les goulots d’étranglement identifiés, nous devons rendre plus rapide le code correspondant.

III-A. Optimisation algorithmique

La première chose à rechercher est l'optimisation algorithmique : existe-t-il des moyens de calculer moins ou mieux?

Pour une vue d'ensemble du problème, une bonne compréhension des calculs derrière l'algorithme peut être utile. Cependant, il n'est pas rare de trouver des changements simples, tels que le déplacement des calculs ou l'allocation de mémoire en dehors d'une boucle, qui apportent des gains importants.

III-A-1. Exemple du SVD

Dans les deux exemples ci-dessus, la décomposition en valeurs singulières (SVD) est ce qui prend le plus de temps. En effet, le coût de calcul de cet algorithme est approximativement égal à kitxmlcodeinlinelatexdvpn^3finkitxmlcodeinlinelatexdvpdans la taille de la matrice d'entrée.

Cependant, dans ces deux exemples, nous n'utilisons pas toute la sortie du SVD, mais seulement les premières lignes de son premier argument de retour. Si nous utilisons l'implémentation SVD de scipy, nous pouvons demander une version incomplète du SVD. Notez que les implémentations de l'algèbre linéaire dans scipy sont plus riches que celles de numpy et devraient être préférées.

 
Sélectionnez
In [3]: %timeit np.linalg.svd(data)
1 loops, best of 3: 14.5 s per loop

In [4]: from scipy import linalg

In [5]: %timeit linalg.svd(data)
1 loops, best of 3: 14.2 s per loop

In [6]: %timeit linalg.svd(data, full_matrices=False)
1 loops, best of 3: 295 ms per loop

In [7]: %timeit np.linalg.svd(data, full_matrices=False)
1 loops, best of 3: 293 ms per loop

Nous pouvons ensuite utiliser cette information pour optimiser le code précédent :

Les SVD vraiment incomplets, par exemple calculant uniquement les 10 premiers vecteurs propres, peuvent être calculés avec arpack, disponible dans scipy.sparse.linalg.eigsh.

Algèbre linéaire de calcul

Pour certains algorithmes, beaucoup de goulots d'étranglement seront des calculs d'algèbre linéaire. Dans ce cas, il est essentiel d’utiliser la bonne fonction pour résoudre le bon problème. Par exemple, un problème de valeur propre avec une matrice symétrique est plus facile à résoudre qu'avec une matrice générale. Par ailleurs, le plus souvent, vous pouvez éviter d’inverser une matrice et utiliser une opération moins coûteuse (et plus stable numériquement).

Connaissez votre algèbre linéaire de calcul. En cas de doute, explorez scipy.linalg et utilisez %timeit pour essayer différentes solutions sur vos données.

 
Sélectionnez
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
def test():
    data = np.random.random((5000, 100))
    u, s, v = linalg.svd(data, full_matrices=False)
    pca = np.dot(u[:, :10].T, data)
    results = fastica(pca.T, whiten=False)
In [1]: import demo

In [2]: %timeit demo.
demo.fastica   demo.np        demo.prof.pdf  demo.py        demo.pyc
demo.linalg    demo.prof      demo.prof.png  demo.py.lprof  demo.test

In [2]: %timeit demo.test()
ica.py:65: RuntimeWarning: invalid value encountered in sqrt
  W = (u * np.diag(1.0/np.sqrt(s)) * u.T) * W  # W = (W * W.T) ^{-1/2} * W
1 loops, best of 3: 17.5 s per loop

In [3]: import demo_opt

In [4]: %timeit demo_opt.test()
1 loops, best of 3: 208 ms per loop

IV. Écrire un code numérique plus rapide

Une discussion complète sur l'utilisation avancée de numpy se trouve au chapitre Advanced NumPy ou dans l'article The NumPy array: a structure for efficient numerical computation par van der Walt et al. Nous ne présentons ici que quelques astuces fréquemment rencontrées pour rendre le code plus rapide.

  • Vectoriser les boucles
    Trouvez des astuces pour éviter aux boucles d’utiliser les tableaux numpy. Pour cela, les tableaux de masques et d'index peuvent être utiles.
  • Utiliser « Broadcasting »
    Utilisez broadcasting pour effectuer des opérations sur des tableaux aussi petits que possible avant de les combiner.
  • Exécuter les opérations sur place
 
Sélectionnez
In [1]: a = np.zeros(1e7)

In [2]: %timeit global a ; a = 0*a
10 loops, best of 3: 111 ms per loop

In [3]: %timeit global a ; a *= 0
10 loops, best of 3: 48.4 ms per loop

Remarque : nous avons besoin de global a dans le timeit pour que cela fonctionne, comme il est assigné à a et le considère donc comme une variable locale.

  • Économiser la mémoire

Utilisez des vues, pas des copies

La copie de grands tableaux est aussi coûteuse que d'effectuer de simples opérations numériques :

 
Sélectionnez
In [1]: a = np.zeros(1e7)

In [2]: %timeit a.copy()
10 loops, best of 3: 124 ms per loop

In [3]: %timeit a + 1
10 loops, best of 3: 112 ms per loop
  • Se méfier des effets de cache

L'accès à la mémoire est moins coûteux lorsqu'il est groupé : l'accès à un grand tableau de manière continue est beaucoup plus rapide que l'accès aléatoire. Cela implique notamment que les petites étapes sont plus rapides (voir CPU cache effects) :

 
Sélectionnez
In [1]: c = np.zeros((1e4, 1e4), order='C')

In [2]: %timeit c.sum(axis=0)
1 loops, best of 3: 3.89 s per loop

In [3]: %timeit c.sum(axis=1)
1 loops, best of 3: 188 ms per loop

In [4]: c.strides
Out[4]: (80000, 8)

C’est la raison pour laquelle les commandes Fortran ou C peuvent faire toute la différence sur les opérations :

 
Sélectionnez
In [5]: a = np.random.rand(20, 2**18)

In [6]: b = np.random.rand(20, 2**18)

In [7]: %timeit np.dot(b, a.T)
1 loops, best of 3: 194 ms per loop

In [8]: c = np.ascontiguousarray(a.T)

In [9]: %timeit np.dot(b, c)
10 loops, best of 3: 84.2 ms per loop

Notez que copier les données pour contourner cet effet peut ne pas en valoir la peine :

 
Sélectionnez
In [10]: %timeit c = np.ascontiguousarray(a.T)
10 loops, best of 3: 106 ms per loop
  • Utiliser numexpr peut être utile pour optimiser automatiquement le code contre de tels effets.
  • Utiliser le code compilé
  • Le dernier recours, une fois que vous êtes certain que toutes les optimisations de haut niveau ont été explorées, consiste à transférer les points chauds, c'est-à-dire les quelques lignes ou fonctions pour lesquelles la plupart du temps est consacré au code compilé. Pour le code compilé, l’option recommandée consiste à utiliser Cython : il est facile de transformer le code Python existant en code compilé, et avec une bonne utilisation du support de numpy, on obtient un code efficace sur les tableaux numpy, par exemple en déroulant des boucles.

Pour tout ce qui précède : profilez et chronométrez vos choix. Ne basez pas votre optimisation sur des considérations théoriques.

IV-A. Liens additionnels

  • si vous avez besoin de profiler l'utilisation de la mémoire, vous pouvez essayer memory_profiler ;

  • si vous avez besoin de profiler dans des extensions C, vous pouvez utiliser gperftools de Python avec yep ;

  • si vous souhaitez suivre les performances de votre code dans le temps, c'est-à-dire lorsque vous effectuez de nouveaux commits dans votre référentiel, vous pouvez essayer asv ;

  • si vous avez besoin de visualisation interactive, pourquoi ne pas essayer RunSnakeRun ?

Vous avez aimé ce tutoriel ? Alors partagez-le en cliquant sur les boutons suivants : Viadeo Twitter Facebook Share on Google+   

  

Les sources présentées sur cette page sont libres de droits et vous pouvez les utiliser à votre convenance. Par contre, la page de présentation constitue une œuvre intellectuelle protégée par les droits d'auteur. Copyright © 2019 Gaël Varoquaux. Aucune reproduction, même partielle, ne peut être faite de ce site ni de l'ensemble de son contenu : textes, documents, images, etc. sans l'autorisation expresse de l'auteur. Sinon vous encourez selon la loi jusqu'à trois ans de prison et jusqu'à 300 000 € de dommages et intérêts.