Efficacité calcul itératif

Hello,

Pour ma recherche, j’ai codé un projet de simulation en Python avec Numpy, et je me pose des questions sur l’efficacité. J’envisage de réécrire des parties en C++ mais je me demande si c’est bien utile.

La partie qui m’intéresse consiste en un bout de code qui itère un numpy.array de taille 3x3 avec des numpy.clongdouble une centaine de fois. De plus, j’ai des threads différents qui font la même chose pour des matrices différentes. Je trouve la parallélisation en Python pas tiptop (parce que le partage de mémoire est casse pieds), et je me demande si en C++ ça ne marcherait pas mieux.

Je sais que numpy est connu pour avoir de bonnes performances sur de grands data set, ce qui n’est pas vraiment le cas ici ; et a tendance à pêcher sur le calcul itératif (ce qui est le cas ici). Mais je sais aussi qu’il est codé en dur par du C++/Fortran et je ne sais pas vraiment ce qu’une implémentation en C++ pourrait améliorer.

Et puis une autre problématique est quant à la portabilité et réutilisation. Les matheux sont rarement de bons codeurs, et il est quasi certain que la partie C++ sera trop difficile pour eux à lire. Donc je chercherai un bon facteur d’amélioration pour que ça vaille le coût.

Avez-vous un avis expérimenté sur le sujet ? Merci

Je ne connaissais pas Numba, je vais faire des tests. Dask me parait pas trop adapté car mes data sets sont assez petits (matrice 3x3 essentiellement).

C’est vraiment np.dot qui me parait être le facteur limitant. Le reste ce sont des tests conditionnels peu coûteux en calcul

Encore une fois, sans algo et sans faire un profilage du code, il est souvent difficile de se faire une bonne idée de ce qui coince et comment l’améliorer…

J’essayerai demain de poster un code suffisamment court et représentatif… ça demande un peu de travail

Un truc que j’ai fait et qui marche bien dans certains cas, c’est de coder une fonction en Fortran, qui va utiliser directement des arrays numpy. Ça se fait très bien avec f2py (fourni avec numpy). Lorsque tu as une fonction qui concentre les difficultés de calcul (des boucles en python, typiquement), ça permet de ne réécrire qu’elle, tout en utilisant partout des arrays numpy sans problème.

L’interfaçage C++/python est parfois un peu complexe, je crois.

Après, comme le dit adri1, ça dépend du point bloquant. Et c’est pas sûr que ce soit Python. J’attends le code court et représentatif.

L’interfaçage C++/python est parfois un peu complexe, je crois.

Il faut juste savoir où aller piocher. Il y a des frameworks de diverses générations pour aider à interfacer.

pybind11 est celui qui m’a le plus séduit de ceux que j’ai regardés, mais c’est loin d’être le seul. Cf https://wiki.python.org/moin/IntegratingPythonWithOtherLanguages

Dans tous les cas en aveugle on ne peut rien dire. Sur du numérique, matriciel & co, Python se débrouille déjà très bien avec numpy — en amateur, il n’est pas aussi simple que ça de faire mieux en C, C++ ou Fortran. numba est une excellente solution de contournement pour gratter facilement des perfs. Et si vraiment ça ne suffit pas ben… il reste, C, C++ et Fortran pour écrire les briques bas niveau, mais … il faut un minimum s’investir et avoir des notions d’optim.

Hello,

Tout d’abord j’ai essayé numba, et il donne de bons résultats.

C’est assez difficile pour moi de vous donner une idée fidèle de ce qu’il se passe sans pouvoir donner plusieurs centaines de lignes de code.

J’ai fait un code très élémentaire mais qui représente un peu ce qu’il se passe.

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from os import system
from time import time
import numpy as np

from numba import jit


C_DTYPE = np.cdouble
R_DTYPE = np.double

DEFAULT_MATRIX = np.array([
[1., 1., -1/2.],
[0., 1., 1.],
[0., 1., 1.]
], dtype=np.dtype(C_DTYPE))

ITERATIONS_NUMBER = 500

EPSILON = 1e-10


@jit(nopython=True, cache=True, fastmath=True)
def goldman_trace(matrix):

    z = matrix[0][0] + matrix[1][1] + matrix[2][2]
    z2 = np.abs(z)**2
    return z2**2 - 8.*(z**3).real + 18.*z2 - 27.

# Calcule l'orbite et cherche convergence
@jit(nopython=True, cache=True, fastmath=True)
def iterate(matrix=DEFAULT_MATRIX):

    point = np.array([C_DTYPE(-0.1), C_DTYPE(0.1), C_DTYPE(1.)])

    if goldman_trace(matrix) < 1e-10:
        # verifie que la dynamique est hyperbolique
        return (False,point)

    last_point = point

    for _ in range(ITERATIONS_NUMBER):

        last_point = point

        point = np.dot(matrix, point)

        z_inverse = point[2]**(-1)

        point = point * z_inverse

        if (np.abs(point - last_point) < EPSILON).all():

            return (True, point)

    return (False,point)

@jit(nopython=True, cache=True)
def compute():

    for _ in range(1000):

        for _ in range(1000):

            iterate()

    return 0


t = time()
compute()
print(time() - t)

Quelques remarques :

• En réalité, bien évidemment, la matrice n’est pas toujours la même et est calculée à la volée juste avant qu’elle soit nécessaire. C’est un peu couteux, mais indépendamment de moi (provient d’ailleurs).
• Plusieurs tests ont été supprimés dans iterate() mais ne sont pas couteux.
• Je me demande si numba dans cet exemple n’exploite pas intelligemment le fait que la matrice soit toujours la même …. en pratique numba m’a donné un facteur 5.

Quelques soucis avec numba:

• j’aurais aimé l’employer à plusieurs endroits … mais il ne manipule pas les dictionnaires ;
• numba ne prends pas les clongdouble de numpy, ce qui est assez dommage …

Merci pour ces remarques, que j’ai prises en compte.

Pour l’expression du polynôme, une simplification est possible si $z$ est réel (on aurait $(z+1)(z-3)^3$). Mais ce cas n’est pas forcément très courant et je crains que le gain soit perdu avec la conditionnelle systématique …

edit : le nonsense de la division en deux étapes provient du fait que j’ai une conditionnelle entre temps pour la certification du calcul (mais j’ai appliqué ton conseil sur **(-1))

edit2 : après essais et mesures, il semblerait que la simplification pour $z$ réel est assez rentable finalement (~5–10% d’économie, ce qui est pas mal)

Tu dis que les conditionnelles sont peu coûteuses, mais il faut voir aussi que selon les cas, ça peut avoir des effets assez importants sur la performance à cause de la manière dont les processeurs modernes fonctionnent (cache, prédiction de branche, etc.). Je ne sais pas à quel point c’est impactant sur ton code et pour ce genre d’environnement Python ceci dit. Quand c’est possible, certaines personnes font du branch-free spécialement pour ça (y compris dans des langages performants), quitte à utiliser des astuces de calcul qu’on sait qu’elles sont compilées sans instructions de saut.

Je ne connais rien à la compilation à la volée en soi, mais je vois des choses qui pourraient améliorer la performance et qui font rester en Python pur. Voilà quelques trucs qui pourraient être testées pour voir si le gain est significatif (je ne suis pas sûr du comportement pour chaque, il y a peut être des choses éliminables d’emblée) :

• Est-ce que numpy.trace ne serait pas envisageable et plus rapide que la calcul explicite dans goldman_trace ?
• Qu’est-ce qu’on gagne en faisant $z \bar z$ au lieu de la norme au carré ?
• en fonction de comment Python se débrouille la liste utilisée par np.array([...]) dans iterate pourrait être potentiellement recrée à chaque fois, ce qui est coûteux.
• L’appel de petites fonctions peut être coûteux en Python, parce que l’appel n’est à ma connaissance pas inliné. Si ta fonction de trace est véritablement petite et peu utilisée par ailleurs, ça peut valoir le coup de l'inliner (vu qu’elle semble appelée jusqu’à 500 millions de fois ici).

j’ai une conditionnelle entre temps pour la certification du calcul

Qu’est-ce que tu cherches à certifier exactement, si ce n’est pas indiscret ?

Je me suis posé les mêmes questions que toi :

• Est-ce que numpy.trace ne serait pas envisageable et plus rapide que la calcul explicite dans goldman_trace ?

Visiblement, numba ne prend pas en charge cette méthode

• Qu’est-ce qu’on gagne en faisant $z \bar z$ au lieu de la norme au carré ?

J’ai essayé, avec le type de numba, mais je n’ai pas réussi …

• en fonction de comment Python se débrouille la liste utilisée par np.array([...]) dans iterate pourrait être potentiellement recrée à chaque fois, ce qui est coûteux.

j’ai corrigé ça

• L’appel de petites fonctions peut être coûteux en Python, parce que l’appel n’est à ma connaissance pas inliné. Si ta fonction de trace est véritablement petite et peu utilisée par ailleurs, ça peut valoir le coup de l'inliner (vu qu’elle semble appelée jusqu’à 500 millions de fois ici).

Ah, ça je savais pas, je vais tester.

Je calcule des ensembles limites, qui sont des attracteurs dynamiques. On cherche à certifier que les points obtenus par convergence en sont bien … Je fais ça en vérifiant qu’ils restent bornés et que la division ne change pas trop le module

Le code sera in fine opensource si ça t’intéresse

Salut,

Une petite idée comme ça, mais il faudrait que quelqu’un confirme : ce genre de calculs devrait bien tourner sur GPU (multiplication de matrices notamment), et les librairies python qui permettent d’utiliser le GPU sont bien développées grâce à nos amis les data scientists. Ca vaudrait peut être le coup de tenter si tu as un GPU sous la main.

Salut,

J’ai survolé mais si jamais tu veux optimiser un calcul en intégrant un autre langage, Rust peut aussi être éventuellement une bonne option. C’est sûrement plus simple à interfacer, plus abordable pour des développeurs novices, et plus straighforward pour du multithreading, que du C++.

Exemple ici, intégration d’un calcul multithreadé dans du Python : https://doc.rust-lang.org/1.5.0/book/rust-inside-other-languages.html

Bonne journée !

J’ai un bug assez absurde que je n’arrive pas à comprendre

Le test goldman_trace < 1e-6 est beaucoup plus rapide que goldman_trace < 0 ou même np.signbit(goldman_trace)

Quelqu’un peut m’expliquer ce qu’il se passe ?

edit: d’ailleurs @Aabu, np.trace fonctionne bien, j’ai du mal comprendre une erreur plus tot dans la journée

Au risque de dire une bêtise : conversion implicite ?

goldman_trace < 0., ça donne quoi ?

Ralenti aussi

Je remonte mon sujet, car ce bug commence à me chiffonner.

J’ai essayé de reproduire la chose avec timeit mais sans succès : 0. et 1e-6 donnent les mêmes résultats.

Cependant, avec numpy seul, alors que np.signbit devrait être plus rapide (normalement c’est juste un accès en mémoire, non?) il met en fait deux fois plus de temps …

% python3 -m timeit -s "import numpy as np" "np.longdouble(0.1)>0."
500000 loops, best of 5: 872 nsec per loop
% python3 -m timeit -s "import numpy as np" "np.signbit(np.longdouble(0.1))"
100000 loops, best of 5: 1.99 usec per loop

J’aimerais bien pourtant avoir accès à un test très rapide de signe. Je n’arrive pas à comprendre comme numba joue la chose, si quelqu’un sait, ce serait top

Si ça se trouve, la différence est juste due au look up dynamique de ce qu’est np.signbit, alors que longdouble.__gt__ doit être inline. Faudrait regarder le bytecode pour être sûr.

Comment tu as trouvé longdouble.__gt__ ?

a > b en Python, c’est du sucre syntaxique pour a.__gt__(b).