Rapide Exp calcul: possible pour améliorer la précision sans perdre trop de la performance?

Je suis en train de sortir de la fast Exp(x) la fonction qui, auparavant, a été décrit dans cette réponse à un DONC, la question sur l'amélioration de la vitesse de calcul en C#:

public static double Exp(double x)
{
  var tmp = (long)(1512775 * x + 1072632447);
  return BitConverter.Int64BitsToDouble(tmp << 32);
}

L'expression est l'utilisation de certaines virgule flottante IEEE "trucs et astuces" et est principalement conçu pour une utilisation dans des ensembles de neurones. La fonction est environ 5 fois plus rapide que le Math.Exp(x) fonction.

Malheureusement, la précision numérique est seulement -4% -- +2% par rapport à la régulière Math.Exp(x) fonction, idéalement, je voudrais avoir une précision à l'intérieur d'au moins la sous-gamme pour cent.

J'ai tracé le quotient entre l'approximatif et le régulier fonctions Exp, et comme on peut le voir dans le graphique de la différence relative semble être répété avec pratiquement constante de fréquence.

Rapide Exp calcul: possible pour améliorer la précision sans perdre trop de la performance?

Est-il possible de profiter de cette régularité pour améliorer la précision de la "fast exp" fonction " sans réduire la vitesse de calcul, ou serait le calcul de la surcharge d'une précision d'amélioration de l'emporter sur le calcul du gain de l'expression d'origine?

(Comme une note de côté, j'ai aussi essayé de d'un autre approches proposées dans la même DONC, la question, mais cette approche ne semble pas être efficaces de calcul en C#, du moins pas pour le cas général.)

MISE À JOUR LE 14 MAI

À la demande de @Adriano, j'ai effectué un très simple indice de référence. J'ai réalisé 10 millions de calculs à l'aide de chacune des autres exp fonctions pour les valeurs à virgule flottante dans l'intervalle [-100, 100]. Depuis la plage de valeurs, je suis intéressé à s'étend de -20 à 0 j'ai aussi explicitement la valeur de la fonction en x = -5. Voici les résultats:

      Math.Exp: 62.525 ms, exp(-5) = 0.00673794699908547
Empty function: 13.769 ms
     ExpNeural: 14.867 ms, exp(-5) = 0.00675211846828461
    ExpSeries8: 15.121 ms, exp(-5) = 0.00641270968867667
   ExpSeries16: 32.046 ms, exp(-5) = 0.00673666189488182
          exp1: 15.062 ms, exp(-5) = -12.3333325982094
          exp2: 15.090 ms, exp(-5) = 13.708332516253
          exp3: 16.251 ms, exp(-5) = -12.3333325982094
          exp4: 17.924 ms, exp(-5) = 728.368055056781
          exp5: 20.972 ms, exp(-5) = -6.13293614238501
          exp6: 24.212 ms, exp(-5) = 3.55518353166184
          exp7: 29.092 ms, exp(-5) = -1.8271053775984
      exp7 +/-: 38.482 ms, exp(-5) = 0.00695945286970704

ExpNeural est l'équivalent de la Exp fonction spécifiée dans le début de ce texte. ExpSeries8 est le formulation que j'ai à l'origine revendiquée n'était pas très efficace .NET; lors de la mise en œuvre exactement comme Neil il était en fait très rapide. ExpSeries16 est l'analogue de la formule, mais avec 16 multiplications au lieu de 8. exp1 par exp7 sont les différentes fonctions de Adriano réponse ci-dessous. La dernière variante de exp7 est une variante où le signe de x est cochée; si le résultat est négatif, la fonction renvoie 1/exp(-x) à la place.

Malheureusement, ni l'un ni l' expN fonctions énumérées par Adriano sont suffisantes dans la plus large de la valeur négative de la gamme que j'envisage. L'extension de la série approche par Neil Coffey semble être plus appropriée dans "ma" plage de valeurs, mais il est trop rapide divergentes avec les plus grandes négatif x, en particulier lors de l'utilisation de "seulement" 8 multiplications.

  • je suis curieux de connaître votre référence à des "ensembles de neurones." actuellement, je suis à la simulation d'un réseau de neurones à l'aide de C++ et face à la même exp goulot d'étranglement des performances que vous avez relevé. existe-il des papiers dans les neurosciences computationnelles qui ont proposé approximative exp fonctions qui sont très vite?
InformationsquelleAutor Anders Gustafsson | 2012-05-11
  1. 9

    Au cas où quelqu'un veut reproduire l'erreur relative de la fonction de la question, voici un moyen à l'aide de Matlab (le "rapide" de l'exposant n'est pas très rapide dans Matlab, mais il est précis):

    t = 1072632447+[0:ceil(1512775*pi)];
x = (t - 1072632447)/1512775;
ex = exp(x);
t = uint64(t);
import java.lang.Double;
et = arrayfun( @(n) java.lang.Double.longBitsToDouble(bitshift(n,32)), t );
plot(x, et./ex);

    Maintenant, la période de l'erreur exactement coïncide avec le moment où la valeur binaire de tmp déborde de la mantisse en l'exposant. Brisons nos données dans des bacs en jetant les bits qui font l'exposant (ce qui périodique), et de ne garder que le haut de huit bits restants (pour faire de notre table de recherche d'une taille raisonnable):

    index = bitshift(bitand(t,uint64(2^20-2^12)),-12) + 1;

    Maintenant nous calculons la moyenne d'ajustement requise:

    relerrfix = ex./et;
adjust = NaN(1,256);
for i=1:256; adjust(i) = mean(relerrfix(index == i)); end;
et2 = et .* adjust(index);

    L'erreur relative est réduite à +/- .0006. Bien sûr, d'autres tables de dimensions sont possibles (par exemple, un 6 bits de la table avec 64 entrées +/- .0025) et l'erreur est quasiment linéaire dans la taille de la table. L'interpolation linéaire entre les entrées de la table serait d'améliorer l'erreur encore plus loin, mais au détriment des performances. Puisque nous avons déjà rencontré la précision de l'objectif, il faut éviter toute nouvelle dégradation de leurs performances.

    À ce point, c'est trivial éditeur de compétences pour prendre les valeurs calculées par MatLab et de créer une table de recherche en C#. Pour chaque calcul, nous avons ajouter un masque de bits, la lecture de la table, et de double précision multiplier.

    static double FastExp(double x)
{
    var tmp = (long)(1512775 * x + 1072632447);
    int index = (int)(tmp >> 12) & 0xFF;
    return BitConverter.Int64BitsToDouble(tmp << 32) * ExpAdjustment[index];
}

    L'accélération est très semblable à l'original du code-pour mon ordinateur, c'est environ 30% plus rapide compilé que les x86 et environ 3x plus rapide pour x64. Avec mono sur ideone, c'est une importante perte nette (mais c'est l'original).

    Le code source complet et cas de test: http://ideone.com/UwNgx

    using System;
using System.Diagnostics;
namespace fastexponent
{
class Program
{
static double[] ExpAdjustment = new double[256] {
1.040389835,
1.039159306,
1.037945888,
1.036749401,
1.035569671,
1.034406528,
1.033259801,
1.032129324,
1.031014933,
1.029916467,
1.028833767,
1.027766676,
1.02671504,
1.025678708,
1.02465753,
1.023651359,
1.022660049,
1.021683458,
1.020721446,
1.019773873,
1.018840604,
1.017921503,
1.017016438,
1.016125279,
1.015247897,
1.014384165,
1.013533958,
1.012697153,
1.011873629,
1.011063266,
1.010265947,
1.009481555,
1.008709975,
1.007951096,
1.007204805,
1.006470993,
1.005749552,
1.005040376,
1.004343358,
1.003658397,
1.002985389,
1.002324233,
1.001674831,
1.001037085,
1.000410897,
0.999796173,
0.999192819,
0.998600742,
0.998019851,
0.997450055,
0.996891266,
0.996343396,
0.995806358,
0.995280068,
0.99476444,
0.994259393,
0.993764844,
0.993280711,
0.992806917,
0.992343381,
0.991890026,
0.991446776,
0.991013555,
0.990590289,
0.990176903,
0.989773325,
0.989379484,
0.988995309,
0.988620729,
0.988255677,
0.987900083,
0.987553882,
0.987217006,
0.98688939,
0.98657097,
0.986261682,
0.985961463,
0.985670251,
0.985387985,
0.985114604,
0.984850048,
0.984594259,
0.984347178,
0.984108748,
0.983878911,
0.983657613,
0.983444797,
0.983240409,
0.983044394,
0.982856701,
0.982677276,
0.982506066,
0.982343022,
0.982188091,
0.982041225,
0.981902373,
0.981771487,
0.981648519,
0.981533421,
0.981426146,
0.981326648,
0.98123488,
0.981150798,
0.981074356,
0.981005511,
0.980944219,
0.980890437,
0.980844122,
0.980805232,
0.980773726,
0.980749562,
0.9807327,
0.9807231,
0.980720722,
0.980725528,
0.980737478,
0.980756534,
0.98078266,
0.980815817,
0.980855968,
0.980903079,
0.980955475,
0.981017942,
0.981085714,
0.981160303,
0.981241675,
0.981329796,
0.981424634,
0.981526154,
0.981634325,
0.981749114,
0.981870489,
0.981998419,
0.982132873,
0.98227382,
0.982421229,
0.982575072,
0.982735318,
0.982901937,
0.983074902,
0.983254183,
0.983439752,
0.983631582,
0.983829644,
0.984033912,
0.984244358,
0.984460956,
0.984683681,
0.984912505,
0.985147403,
0.985388349,
0.98563532,
0.98588829,
0.986147234,
0.986412128,
0.986682949,
0.986959673,
0.987242277,
0.987530737,
0.987825031,
0.988125136,
0.98843103,
0.988742691,
0.989060098,
0.989383229,
0.989712063,
0.990046579,
0.990386756,
0.990732574,
0.991084012,
0.991441052,
0.991803672,
0.992171854,
0.992545578,
0.992924825,
0.993309578,
0.993699816,
0.994095522,
0.994496677,
0.994903265,
0.995315266,
0.995732665,
0.996155442,
0.996583582,
0.997017068,
0.997455883,
0.99790001,
0.998349434,
0.998804138,
0.999264107,
0.999729325,
1.000199776,
1.000675446,
1.001156319,
1.001642381,
1.002133617,
1.002630011,
1.003131551,
1.003638222,
1.00415001,
1.004666901,
1.005188881,
1.005715938,
1.006248058,
1.006785227,
1.007327434,
1.007874665,
1.008426907,
1.008984149,
1.009546377,
1.010113581,
1.010685747,
1.011262865,
1.011844922,
1.012431907,
1.013023808,
1.013620615,
1.014222317,
1.014828902,
1.01544036,
1.016056681,
1.016677853,
1.017303866,
1.017934711,
1.018570378,
1.019210855,
1.019856135,
1.020506206,
1.02116106,
1.021820687,
1.022485078,
1.023154224,
1.023828116,
1.024506745,
1.025190103,
1.02587818,
1.026570969,
1.027268461,
1.027970647,
1.02867752,
1.029389072,
1.030114973,
1.030826088,
1.03155163,
1.032281819,
1.03301665,
1.033756114,
1.034500204,
1.035248913,
1.036002235,
1.036760162,
1.037522688,
1.038289806,
1.039061509,
1.039837792,
1.040618648
};
static double FastExp(double x)
{
var tmp = (long)(1512775 * x + 1072632447);
int index = (int)(tmp >> 12) & 0xFF;
return BitConverter.Int64BitsToDouble(tmp << 32) * ExpAdjustment[index];
}
static void Main(string[] args)
{
double[] x = new double[1000000];
double[] ex = new double[x.Length];
double[] fx = new double[x.Length];
Random r = new Random();
for (int i = 0; i < x.Length; ++i)
x[i] = r.NextDouble() * 40;
Stopwatch sw = new Stopwatch();
sw.Start();
for (int j = 0; j < x.Length; ++j)
ex[j] = Math.Exp(x[j]);
sw.Stop();
double builtin = sw.Elapsed.TotalMilliseconds;
sw.Reset();
sw.Start();
for (int k = 0; k < x.Length; ++k)
fx[k] = FastExp(x[k]);
sw.Stop();
double custom = sw.Elapsed.TotalMilliseconds;
double min = 1, max = 1;
for (int m = 0; m < x.Length; ++m) {
double ratio = fx[m] / ex[m];
if (min > ratio) min = ratio;
if (max < ratio) max = ratio;
}
Console.WriteLine("minimum ratio = " + min.ToString() + ", maximum ratio = " + max.ToString() + ", speedup = " + (builtin / custom).ToString());
}
}
}
    • Travail fantastique, et une bonne explication! Merci beaucoup pour cette réponse, c'était juste le genre de progrès que j'avais espéré. Aviez-vous développé plus tôt, ou avez-vous mettre en œuvre à la suite de cette question?
    • Je suis totalement volé l'approche suggérée dans la question.
    • Après des essais à android NDK, il est plus lent que le système de std::exp(). Mais il est plus rapide dans le PC. (gist.github.com/maxint/0172c1dcd075d3589eeb)
    • Ce code a un comportement indéfini écrit sur elle. Utilisation memcpy pour votre type beaucoup les jeux de mots. De toute façon, selon que votre cible a virgule flottante, vous pouvez utiliser en simple précision pour la table de recherche. Nous parlons d'une erreur relative de .0006, sorte de double précision n'aide pas.
    • Avez-vous un C mise en œuvre de cette?
    • Il suffit de remplacer BitConverter fonctions par un memcpy, et déplacer le [] dans la définition de tableau. Le reste du code C# est valide C déjà.

    InformationsquelleAutor Ben Voigt
  2. 10

    Essayez les solutions suivantes (exp1 est plus rapide, exp7 est plus précis).

    Code

    public static double exp1(double x) { 
return (6+x*(6+x*(3+x)))*0.16666666f; 
}
public static double exp2(double x) {
return (24+x*(24+x*(12+x*(4+x))))*0.041666666f;
}
public static double exp3(double x) {
return (120+x*(120+x*(60+x*(20+x*(5+x)))))*0.0083333333f;
}
public static double exp4(double x) {
return 720+x*(720+x*(360+x*(120+x*(30+x*(6+x))))))*0.0013888888f;
}
public static double exp5(double x) {
return (5040+x*(5040+x*(2520+x*(840+x*(210+x*(42+x*(7+x)))))))*0.00019841269f;
}
public static double exp6(double x) {
return (40320+x*(40320+x*(20160+x*(6720+x*(1680+x*(336+x*(56+x*(8+x))))))))*2.4801587301e-5;
}
public static double exp7(double x) {
return (362880+x*(362880+x*(181440+x*(60480+x*(15120+x*(3024+x*(504+x*(72+x*(9+x)))))))))*2.75573192e-6;
}

    Précision

    De la fonction d'Erreur dans [-1...1] Erreur [3.14 3.14...] 
exp1 de 0,05% à 1,8% 8.8742 38.40% 
exp2 0.01 0.36% 4.8237 20.80% 
exp3 0.0016152 0.59% 2.28 9.80% 
exp4 0.0002263 de 0,0083% 0.9488 4.10% 
exp5 0.0000279 à 0,001% 0.3516 de 1,50% 
exp6 0.0000031 0.00011% 0.1172 de 0,50% 
exp7 0.0000003 0.000011% 0.0355 0.15%

    Crédits

    Ces implémentations de exp() ont été calculés par "scoofy" à l'aide de la série de Taylor d'une tanh() mise en œuvre de "fuzzpilz" (quels qu'ils soient, je viens d'avoir ces références sur mon code).

    • "fuzzpilz" LOL. Certaines personnes ont un drôle de sens de nick noms.
    • Origine de la série de Taylor de l'approximation par [email protected] ici: musicdsp.org/showone.php?id=222 - Upvoted que c'est une solution simple et efficace via le développement en série de Taylor, surpris qu'il n'ait pas été posté avant.
    • merci pour la référence, il est passé perdu depuis longtemps!
    • Merci beaucoup, Adriano, pour la fourniture de ces méthodes, et merci Mahmoud de fournir le lien vers le texte original. Je n'avais pas rencontré ces mises en avant. J'ai essayé exp7, et il est environ trois fois plus rapide que les Mathématiques.Exp. Pour x négatif il s'écarte assez rapidement, mais qui pourrait peut-être être pris en compte par le calcul de 1/exp(abs(x))? Je vais revoir mon propre code pour voir si l'intervalle des valeurs de x est suffisante pour l'utilisation de exp7, et si elle est, je suis plus que disposé à accepter votre réponse, Adriano.
    • veuillez mettre à jour votre question (ou poster une réponse) avec les résultats de votre test, je suis à l'aide de cette fonction (en fait exp5) dans un programme C++ donc je n'ai jamais essayé sa performance dans .NET. Je ne sais pas comment la performance de 1/exp(abs(x)) sera, peut-être pas si loin de Mathématiques.Exp. Si vous avez de l'utiliser dans un (par couche) de la boucle, vous pouvez envisager de mettre en œuvre cette fonction dans une bibliothèque C++/CLI, il vous serez capable d'utiliser inline assemblée avec des instructions SIMD sur la boucle (si cela ne vous force pas à restructurer le code de trop).
    • désolé de ne pas répondre plus tôt. Je mettrai à jour la question avec mes résultats de test dès que possible, et je vais le tester C/C++ exp pour la comparaison, à la même heure. J'ai été en utilisant le Parallèle fonctionnalité dans .NET 4, bien que dans ce cas particulier, je n'ai pas encore préparé le code pour la parallélisation. Ce que bien.
    • Merci, je suis assez intéressé. Je ne voulais pas dire la .NET Parallèle de la Bibliothèque (je suppose que pour le rendre parallèle cause quelques refactoring), mais quelque chose comme ceci (dans le cas où c'est possible): stackoverflow.com/a/412831/1207195
    • J'ai mis à jour ma question avec certains de référence timings maintenant. J'ai aussi fait un examen plus attentif de l'échelle de valeur qui est pertinent pour mon application, et les valeurs de l'ordre de 0 à négative des centaines. Les valeurs de -10 -15 -- auront un impact significatif sur le résultat final, donc la précision au moins de ces valeurs est pertinent. BTW, merci pour la suggestion sur le calcul en C++, je vois maintenant ce que tu veux dire. Il faudra une grande refactoring dans mon cas, mais peut-être que ça en vaut la peine; je vais lui donner une certaine pensée.
    • oui, il ne fonctionne bien que dans les +/- 3.14 plage mais merci pour la référence!!!
    • href="http://stackoverflow.com/a/14143184/1339280">La technique décrite dans ma réponse va grandement améliorer la précision sur une large gamme d'entrées.

    InformationsquelleAutor Adriano Repetti
  3. 9

    De la série de Taylor des approximations (comme le expX() fonctions dans Adriano répondre) sont plus précises proche de zéro et peut avoir d'énormes erreurs à -20 ou même à -5. Si l'entrée a une aire de répartition connue, comme de -20 à 0, comme la question d'origine, vous pouvez utiliser une petite table et un autre se multiplier afin d'améliorer considérablement la précision.

    L'astuce est de reconnaître que exp() peut être séparé en entier et de fractions. Par exemple:

    exp(-2.345) = exp(-2.0) * exp(-0.345)

    La partie fractionnaire sera toujours compris entre -1 et 1, donc un développement en série de Taylor rapprochement sera assez précis. La partie entière a seulement 21 valeurs possibles pour exp(-20) exp(0), de sorte que ceux-ci peuvent être stockés dans un petit look up table.

    InformationsquelleAutor shoelzer
  4. 5

    Le code suivant devrait répondre aux exigences de précision, comme pour les entrées dans [-87,88] les résultats ont erreur relative <= 1.73 e-3. Je ne sais pas C#, c'est donc le code en C, mais je suppose que la conversion devrait être assez simple.

    Je suppose que, puisque l'exigence d'exactitude est faible, l'utilisation de la précision de calcul est très bien. Un algorithme classique est utilisé, dans lequel le calcul de exp() est appliquée pour le calcul de exp2(). Après l'argument de conversion par multiplication par log2(e), exponentation par la partie fractionnaire est gérée à l'aide d'un minimax polynôme de degré 2, alors que exponentation par la partie intégrante de l'argument est effectuée par la manipulation directe de l'exposant partie de la norme IEEE-754 simple précision nombre.

    La volatilité de l'union facilite la ré-interprétation d'une séquence de bits comme un entier ou un flottant simple précision-nombre de point nécessaire pour l'exposant à la manipulation. Il ressemble à C# offre decidated ré-interprétation de fonctions pour ce, qui est beaucoup plus propre.

    Les deux potentiels problèmes de performances de la fonction floor() et float->int conversion. Traditionnellement, les deux ont été lents sur x86 en raison de la nécessité de gérer processeur de dynamique de l'état. Mais l'ESS (en particulier SSE 4.1) fournit des instructions qui permettent à ces opérations pour être rapide. Je ne sais pas si C# peuvent faire usage de ces instructions.

     /* max. rel. error <= 1.73e-3 on [-87,88] */
float fast_exp (float x)
{
volatile union {
float f;
unsigned int i;
} cvt;
/* exp(x) = 2^i * 2^f; i = floor (log2(e) * x), 0 <= f <= 1 */
float t = x * 1.442695041f;
float fi = floorf (t);
float f = t - fi;
int i = (int)fi;
cvt.f = (0.3371894346f * f + 0.657636276f) * f + 1.00172476f; /* compute 2^f */
cvt.i += (i << 23);                                          /* scale by 2^i */
return cvt.f;
}
    • Un grand merci pour un grand exemple et une bonne explication. Je vais essayer de convertir votre application en C# pour voir comment il performe par rapport à la Exp de la fonction. Je ne me souviens pas avoir vu cette solution n'importe où ailleurs, avez-vous trouvé cela comme un résultat de la question?
    • J'ai conçu et mis en œuvre des algorithmes pour diverses fonctions transcendantes à plusieurs reprises dans le passé. L'approche que j'ai choisi ci-dessus est très bien d'un jardin-variété de l'algorithme. Je n'ai créer un minimax approximation pour le polynôme spécifiquement en réponse à cette question. Il existe des outils pour cela, comme Mathematica, Maple, et d'autres; ils sont généralement basés sur des variantes de l'algorithme de Remez.
    • Notez qu'il n'est pas correct en C++ pour utiliser de l'union. Mais vous pouvez utiliser memcpy en C et C++, et l'optimiseur doit faire quelque chose de raisonnable, sans le briser, avec optimisation basée sur la stricte aliasing.
    • Avez-vous de l'ESS version (d'Exploitation sur __m128)? Je Vous Remercie.
    • Je n'ai pas du tout l'ESS code prêt à poster, et je n'ai pas programmé avec l'ESS intrinsèques depuis je suis passé au Gpu il y a dix ans.
    • J'ai un ESS version de ce code maintenant. Je suggère de poser une question pour une faible précision expf() SIMD-mise en œuvre et que je ne puis répondre.

    InformationsquelleAutor njuffa
  5. 3

    J'ai étudié la papier par Nicol Schraudolph où l'original C mise en œuvre de la fonction ci-dessus a été défini plus en détail maintenant. Il semble qu'il n'est probablement pas possible à la quasi-approuver l'exactitude de la exp calcul sans de graves impacts sur la performance. D'autre part, l'approximation n'est valable aussi pour les grandes amplitudes de x, jusqu'à +/- 700, ce qui est évidemment avantageux.

    L'implémentation de la fonction ci-dessus est réglé pour obtenir le minimum d'erreur quadratique moyenne. Schraudolph décrit la façon dont le terme additif dans le tmp expression peut être modifiée pour obtenir des rapprochement propriétés.

    "exp" >= exp for all x                      1072693248 -  (-1) = 1072693249
"exp" <= exp for all x                                 - 90253 = 1072602995
"exp" symmetric around exp                             - 45799 = 1072647449
Mimimum possible mean deviation                        - 68243 = 1072625005
Minimum possible root-mean-square deviation            - 60801 = 1072632447

    Il souligne aussi que, à un "microscopique" niveau approximatif "exp" fonction des expositions d'escalier en cas de comportement depuis 32 bits sont jetés dans la conversion de long à double. Cela signifie que la fonction est la pièce sage constant sur une très petite échelle, mais la fonction est au moins de ne jamais diminuer avec l'augmentation de la x.

    InformationsquelleAutor Anders Gustafsson