Calcul rapide d'une gaussienne.

Bonjour,

je calcule une valeur gaussienne de la manière suivante :

Code:

---

1 2 3 4 5

private double gaussian(double x, double y, double sig) { return Math.exp(-(x*x+y*y)/(2*sig*sig)) / (2.0*Math.PI*sig*sig) ; }

---

Seulement cette fonction est appelées des millions de fois, donc je souhaiterai savoir s'il n'y a pas un moyen plus rapide de faire le calcul. Bon à vrai dire je serai surpris que oui.

J'avais pensé à :
- pré-calculer des valeurs, mais les paramètres x, y et sig sont TRES variables.
- faire des approximations car si x et y sont grands, le résultat est quasi nul.

Merci par avance...


PS : quelle est la complexité de la fonction Math.exp() ?
En gros, est elle très gourmande en terme de CPU ?

La fonction est calculée par le FPU du processeur et doit prendre le temps de quelques multiplications. Mon PC en calcule 20 millions en environ 2sec.

On peut approximer exp() en utilisant des hack sur la représentation en flottant, si les valeurs sont assez petites.

Code:

---

1 2 3 4

public static double exp(double val) { long hack = (long) (1512775 * val) + 1072632447; return Double.longBitsToDouble(hack << 32); }

---

On peut aussi utiliser interpoler une table de valeurs précalculées.

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14

static double[] EXPNEGTABLE = new double[1000]; static { for(int i=0;i<1000;i++) EXPNEGTABLE[i]=Math.exp((double)-i/32); }   public static double exp(double val) { if (val>0) throw new ArithmeticException(); int i=(int)(-val*32); if (i>=999) return 0.0; double a=(-val*32)-i; return EXPNEGTABLE[i]+a*(EXPNEGTABLE[i+1]-EXPNEGTABLE[i]); }

---

Hormis faire des approximations, je ne pense pas qu'on puisse accélérer les calculs. La méthode Math.exp() en Java est un appel natif, alors ca va être difficile de faire mieux.

Salut!
Tout d'abord, il est inutile de calculer deux fois 2.0*sig*sig. Si ça n'est pas assez efficace, pourquoi ne pas passer à un langage plus efficace, comme le C, le Fortran ou même Matlab?
Jean-Marc Blanc

Pour ma curiosité personnelle, combien de temps mets ton ordi à calculer une exponentielle ? Et quel système utilises-tu ?

Bonsoir,

Citation:

---

Envoyé par pseudocode

On peut approximer exp() en utilisant des hack sur la représentation en flottant, si les valeurs sont assez petites.

Code:

---

1 2 3 4

public static double exp(double val) { long hack = (long) (1512775 * val) + 1072632447; return Double.longBitsToDouble(hack << 32); }

---

On peut aussi utiliser interpoler une table de valeurs précalculées.

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14

static double[] EXPNEGTABLE = new double[1000]; static { for(int i=0;i<1000;i++) EXPNEGTABLE[i]=Math.exp((double)-i/32); }   public static double exp(double val) { if (val>0) throw new ArithmeticException(); int i=(int)(-val*32); if (i>=999) return 0.0; double a=(-val*32)-i; return EXPNEGTABLE[i]+a*(EXPNEGTABLE[i+1]-EXPNEGTABLE[i]); }

---

Hormis faire des approximations, je ne pense pas qu'on puisse accélérer les calculs. La méthode Math.exp() en Java est un appel natif, alors ca va être difficile de faire mieux.

---

- 1 - je ne suis pas fan du hack.
- 2 - en gros il me faut pré-calculer les valeurs de l'exponentielle pour gagner du temps. Mais est ce réellement significatif d'après toi ? Comme tu le dis toi même, c'est du natif derrière.
Deuxième question : pourquoi diviser par 32 et faire ensuite ce calcul bizarre ?

Citation:

---

Envoyé par FR119492

Salut!
Tout d'abord, il est inutile de calculer deux fois 2.0*sig*sig. Si ça n'est pas assez efficace, pourquoi ne pas passer à un langage plus efficace, comme le C, le Fortran ou même Matlab?

---

J'avais déjà corrigé le double calcul.
Je souhaite optimiser ce calcul qui tourne pour une application en java.

Citation:

---

Envoyé par Nebulix

Pour ma curiosité personnelle, combien de temps mets ton ordi à calculer une exponentielle ? Et quel système utilises-tu ?

---

Je n'ai pas fait de test :-(
Je suis sous Mac OS X.

Citation:

---

Envoyé par ToTo13

- 2 - en gros il me faut pré-calculer les valeurs de l'exponentielle pour gagner du temps. Mais est ce réellement significatif d'après toi ? Comme tu le dis toi même, c'est du natif derrière.
Deuxième question : pourquoi diviser par 32 et faire ensuite ce calcul bizarre ?

---

Avec la table de précalcul tu gagnes un facteur 10 à 15 sur ma machine (Java 6-u21 64bits) . C'est pas négligeable tout de même. :D

La division par 32, c'est juste la fréquence d'échantillonnage de la gaussienne.

table[0] = exp(-0)
table[1] = exp(-1/32)
table[2] = exp(-2/32)
...
table[1000] = exp(-1000/32)

Citation:

---

Envoyé par pseudocode

Avec la table de précalcul tu gagnes un facteur 10 à 15 sur ma machine (Java 6-u21 64bits) . C'est pas négligeable tout de même. :D

La division par 32, c'est juste la fréquence d'échantillonnage de la gaussienne.

---

Merci et merci.
Je garde le facteur d'échantillonnage qui m'intéresse.
Pour le pas je prendrai quelque chose d'un peu plus précis.

Je viens de faire des tests et cette indexation préalable m'apporte un gain de 20% sur la vitesse globale de mon application. Donc c'est super !!!

Juste une dernière question : je ne comprends pas ça :

Code:

---

1 2 3

double a=(-val*32)-i; return EXPNEGTABLE[i]+a*(EXPNEGTABLE[i+1]-EXPNEGTABLE[i]);

---

Citation:

---

Envoyé par ToTo13

Merci et merci.
Je garde le facteur d'échantillonnage qui m'intéresse.
Pour le pas je prendrai quelque chose d'un peu plus précis.

---

Ok. Déjà avec 1/32, l'erreur est au max de 10^-4.

Citation:

---

Envoyé par ToTo13

Je viens de faire des tests et cette indexation préalable m'apporte un gain de 20% sur la vitesse globale de mon application. Donc c'est super !!!

Juste une dernière question : je ne comprends pas ça :

Code:

---

1 2	return EXPNEGTABLE[i]+a*(EXPNEGTABLE[i+1]-EXPNEGTABLE[i]);

---

---

bah... heu. C'est une surpuissante formule d'interpolation linéaire. :aie:

C'est peut être plus compréhensible comme cela:

Code:

---

return (1-a)*EXPNEGTABLE[i] + a*EXPNEGTABLE[i+1];

---

Il se fait tard et mes capacités (limitées d'origine) sont diminuées.

Pour ma part j'ai fait cela et ça semble très bien fonctionner :

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14

private float[] exp = new float[10000] ; ... for (int i=0 ; i < exp.length ; i++) exp[i] = (float)Math.exp(-(double)i/1000.0) ; ... private float gaussian(float x, float y, float sig) { float s = 2.0f*sig*sig ; //return (float) Math.exp(-(x*x+y*y)/s) / (pi*s) ; int n = (int)((x*x+y*y)/s * 1000.0f) ; if ( n >= exp.length ) return 0.0f ; return exp[n] / (pi*s) ; }

---

Citation:

---

Envoyé par pseudocode

bah... heu. C'est une surpuissante formule d'interpolation linéaire. :aie:
C'est peut être plus compréhensible comme cela:

Code:

---

return (1-a)*EXPNEGTABLE[i] + a*EXPNEGTABLE[i+1];

---

---

Ok... je viens enfin de comprendre (que je ferai mieux d'aller dormir) !
Effectivement, comme je partais du principe que l'erreur produite par le pas était négligeable. Je n'avais pas pensé à interpoler la valeur pour la corriger.
Comme d'habitude, j'ai été plus bourin et j'ai joué sur le pas !

Citation:

---

Envoyé par ToTo13

Ok... je viens enfin de comprendre (que je ferai mieux d'aller dormir) !
Effectivement, comme je partais du principe que l'erreur produite par le pas était négligeable. Je n'avais pas pensé à interpoler la valeur pour la corriger.
Comme d'habitude, j'ai été plus bourin et j'ai joué sur le pas !

---

Ah, oui. J'aurais dû préciser que les valeurs précalculées étaient ensuite interpolées. J'ai tellement l'habitude de faire ainsi que je ne le précise plus. :P

Citation:

---

Envoyé par ToTo13

Je n'ai pas fait de test :-(
Je suis sous Mac OS X.

---

Citation:

---

Envoyé par pseudocode

Avec la table de précalcul tu gagnes un facteur 10 à 15 sur ma machine (Java 6-u21 64bits)

---

Citation:

---

Envoyé par ToTo13

Je viens de faire des tests et cette indexation préalable m'apporte un gain de 20% sur la vitesse globale de mon application.

---

Là je ne comprends plus rien. Votre machine est câblée pour calculer une exponentielle en 100 nanosecondes. Vous réinventez la roue ? et vous la réinventez carrée ?

Citation:

---

Envoyé par Nebulix

Là je ne comprends plus rien. Votre machine est câblée pour calculer une exponentielle en 100 nanosecondes. Vous réinventez la roue ? et vous la réinventez carrée ?

---

Ah ? Ton processeur calcule les exponentielles ? C'est quoi l'opcode assembleur ?

Citation:

---

Envoyé par pseudocode

Ah ? Ton processeur calcule les exponentielles ? C'est quoi l'opcode assembleur ?

---

Par raison de commodité, je ne programme plus beaucoup en assembleur, Delphi (entre autres) est plus confortable avec des pertes de performance acceptables. J'ai retrouvé une vieille doc du 8087.(il y a peut-être mieux depuis)
L'instruction calcule 2^x-1, elle s'appelle F2XM1 et son opcode est 9B D9 F0.
Le calcul se fait avec 19 chiffres significatifs.
Voici le bout de code qui m'a permis d'évaluer la vitesse : 20 sec environ pour 200 millions d'appels

Code:

---

1 2 3 4 5 6 7 8 9 10 11

procedure TForm1.Button1Click(Sender: TObject); var i,j,n:integer;x,y:double; n1,n2:tdatetime; begin y:=1.33; n:=1000*1000*10; n1:=now; for j:=1 to n do for i:=1 to 20 do x:=exp(y*i ); n1:=now-n1; caption:=datetimetostr(n1); end;

---

Je suis vraiment très curieux de savoir combien tu mesures sur ton système

F2XM1 calcule 2 puissance (x-1) en virgule flottante.

Le code x86 habituel pour calculer une exponentiel (= le code de la fonction exp() dans la librairie mathématique standard) est:

Code:

---

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

.386 _TEXT segment use32 para public 'CODE' public _exp   _exp proc near assume cs:_TEXT push ebp mov ebp,esp sub esp,8 ; Allocate temporary space fld qword ptr [ebp+8] ; Load real from stack fldl2e ; Load log base 2(e) fmulp st(1),st ; Multiply x * log base 2(e) fst st(1) ; Push result frndint ; Round to integer fsub st(1),st ; Subtract fxch ; Exchange st, st(1) f2xm1 ; Compute 2 to the (x - 1) fld1 ; Load real number 1 fadd ; 2 to the x fscale ; Scale by power of 2 fstp st(1) ; Set new stack top and pop fst qword ptr [ebp-8] ; Throw away scale factor mov esp,ebp ; Deallocate temporary space pop ebp ret _exp endp   _TEXT ends end

---

Comme on peut le voir il y a toute une suite d'opérations en virgule flottantes a réaliser, dont certaines assez couteuses.

Il n'empêche que tout çà ne prend que 100 ns. Est-ce vraiment trop long pour Toto ou y a-t-il un autre problème ?

Citation:

---

Envoyé par Nebulix

Il n'empêche que tout çà ne prend que 100 ns. Est-ce vraiment trop long pour Toto ou y a-t-il un autre problème ?

---

Je cherche à avoir le code le plus rapide possible !
Un pré-calcul permet un accès direct à la valeur et évite des calculs forcément plus coûteux.