Notes de cours Python scientifique

Dans les deux exemples précédents, la SVD décomposition en valeurs singulières est ce qui nous prend le plus de temps.

En effet, le coût de traitement de cet algorithme est environ n³ de la taille de la matrice fournie en entrée.

Toutefois, dans ces deux exemples, nous n'utilisons pas toutes les informations fournies par la SVD, mais seulement les premières lignes du premier argument renvoyé. Si nous utilisons la méthode SVD de scypy, nous pouvons demander une version incomplète de la SVD. Notez que les méthodes d'algèbre linéaire en scipy sont plus complètes que celles en numpy et devraient leur être préférées.

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 donc utiliser ces résultats pour optimiser le code précédent :

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

Des décompositions de valeurs singulières incomplètes sur des réels (ex. : calculer seulement les 10 premiers vecteurs propres) peuvent être réalisées avec arpack, disponible dans scipy.sparse.li, alg.eigsh.

Algèbre linéaire informatique

Pour certains algorithmes, beaucoup de goulots d'étranglement seront des calculs d'algèbre linéaire. Dans ce cas, la solution est d'utiliser la bonne fonction au bon moment. Par exemple, un problème de valeur propre est plus facile à résoudre avec une matrice symétrique plutôt qu'avec une matrice générale. De plus, bien souvent, vous pouvez éviter d'inverser une matrice et utiliser une opération moins gourmande.

Maîtrisez l'algèbre linéaire informatique. En cas de doute, explorez scipy.linalg et utilisez %timeit pour tester d'autres approches avec vos données.

Vous pouvez trouver une analyse complète de l'utilisation avancée de numpy au chapitre Advanced Numpy, ou dans l'article The NumPy array : a structure for efficient numerical computation de van derWalt et al. Ici nous aborderons seulement quelques astuces fréquemment utilisées pour accélérer votre code.

• Méthodes de vectorisation pour les boucles

Trouvez des astuces pour éviter l'utilisation des tableaux numpy dans vos boucles. Pour ce faire, les masques et les tableaux d'indices pourront vous être utiles,

• Méthode de broadcast

Utilisez la méthode de broadcast pour réduire autant que possible les opérations sur les tableaux avant de les combiner.

• Les instructions « en place »

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

• Soyez économe en mémoire, utilisez des vues pas des copies

Copier de grands tableaux est aussi gourmand que de leur appliquer des opérations numériques simples :

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

• Attention aux effets de cache

Les accès mémoire sont plus légers lorsqu'ils sont groupés : accéder à un tableau de manière continue est plus rapide que de manière inopinée. Cela implique surtout que de plus petits traitements sont plus rapides (voir Les effets de cache CPU) :

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)

Voilà pourquoi les commandes Fortran ou C pourront faire une grande différence sur vos traitements :

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 qu'il ne sera probablement pas efficace de copier les données pour contourner ce problème.

In [10]: %timeit c = np.ascontiguousarray(a.T)
10 loops, best of 3: 106 ms per loop

L'utilisation de numexpr peut être utile pour optimiser automatiquement votre code afin de parer à ce genre d'effets

• Utiliser du code compilé

Le dernier recours, une fois que vous êtes sûr d'avoir exploré toutes les autres méthodes d'optimisations, est de transformer les points problématiques (les quelques lignes ou fonctions qui prennent le plus de temps) en code compilé. Pour ce dernier, la meilleure option est d'utiliser Cython : il est facile de transformer du code Python sortant en code compilé, et une bonne utilisation de l'assistance numpy génère un code performant pour les tableaux numpy, en déroulant les boucles par exemple.