Efficacité calcul itératif

a marqué ce sujet comme résolu.
Auteur du sujet

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 :)

+0 -0

Cette réponse a aidé l’auteur du sujet

Salut,

Difficile à dire sans voir le code ni les équations que tu résous. Les problèmes de perf peuvent être dues à plein de choses avant d’être un problème de langage. Ça peut être la faute à un algo de résolution pourri, une mauvaise utilisation de numpy, une mauvaise utilisation de Python, etc.

Pour répondre un peu plus directement, il y a peu de chance que tu puisses faire mieux que les développeurs de numpy sur la partie array elle-même, mais si ton bottleneck est du à une partie en Python ça peut valoir le coup soit de la coder plus intelligemment pour exploiter numpy d’une meilleure façon ou bien d’utiliser un langage plus bas niveau. Mais les autres pistes sont prioritaires.

Sinon, tu as essayé dask, numba ou autre outil similaire ?

I don’t mind that you think slowly, but I do mind that you are publishing faster. — W. Pauli

+1 -0
Auteur du sujet

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

+0 -0

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…

I don’t mind that you think slowly, but I do mind that you are publishing faster. — W. Pauli

+3 -0

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. :)

Il y a bien des façons de passer à l’acte. Se taire en est une. Attribué à Jean-Bertrand Pontalis

+0 -0

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.

+0 -0
Auteur du sujet

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 …
+0 -0

Cette réponse a aidé l’auteur du sujet

Je ne suis pas sûr pour Python, mais cela m’évoque des détails de micro-optimisation relativement aux C&C++.

pow() y est une plaie qui coûte extrêmement cher. pow(x, constexpr int) est cependant optimisé par les compilos contemporains pour quelques valeurs. Et donc il y a deux choses qui me surprennent dans le code python.

  • x * y**(-1) mais… pourquoi faire compliqué? Si le pow n’est pas optimisé, il sera plus lent que x/y. Ni clang, ni gcc ne savent l’optimiser. Si numba ne sait pas le faire non plus, ce sera moins efficace que la simple division

  • a * x**4 + b * x**3 + c * x**2 +d, en temps normal, le truc qui marche bien pour l’évaluation de polynômes, c’est la méthode de Horner. Au lieu d’employer la l’écriture naturelle, on emploie ((((a * x) + b) * x) + c) *x ** 2 + d.

    Cela présente deux avantages: un: une meilleure précision; deux: de meilleurs temps de calculs. De nouveau, je ne sais pas si numba fait l’effort d’interpréter la formule pour employer la méthode de Horner.

    Maintenant, ce n’est pas exactement une évaluation de polynôme qui est réalisée, mais un truc plus compliqué avec des complexes qui mélange partie réelle et norme. A voir si cela peut être transformé et si cela améliorerait les perfs ou non.

Tout ça me fait penser… il n’y a pas plus efficace que np.abs(cplx)**2 pour calculer la norme euclidienne d’un complexe sans déclencher le calcul de racine carrée? (on l’a en C++, quid de Python? J’ai l’impression que cela demande un peu d’huile de coude)

NB: pas besoin du calcul de racine carrée pour la comparaison à epsilon. Mais de nouveau, il faut évaluer si numba sait interpréter le code pour le simplifier ou s’il compile bêtement ce que l’on exprime sans chercher à le simplifier.

+1 -0
Auteur du sujet

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

Pour l’expression du polynôme, une simplification est possible si zz est réel (on aurait (z+1)(z3)3(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 zz réel est assez rentable finalement (~5–10% d’économie, ce qui est pas mal)

Édité par Holosmos

+0 -0

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 zzˉ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 ?

+0 -0
Auteur du sujet

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 zzˉ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.

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 ?

Aabu

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

+0 -0

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.

+0 -0

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 !

+0 -0
Auteur du sujet

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

Édité par Holosmos

+0 -0

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

goldman_trace < 0., ça donne quoi ?

Il y a bien des façons de passer à l’acte. Se taire en est une. Attribué à Jean-Bertrand Pontalis

+0 -0
Auteur du sujet

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 <3

Édité par Holosmos

+0 -0

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.

I don’t mind that you think slowly, but I do mind that you are publishing faster. — W. Pauli

+0 -0
Vous devez être connecté pour pouvoir poster un message.
Connexion

Pas encore inscrit ?

Créez un compte en une minute pour profiter pleinement de toutes les fonctionnalités de Zeste de Savoir. Ici, tout est gratuit et sans publicité.
Créer un compte