Fast Exp-Berechnung: möglich, um die Genauigkeit zu verbessern, ohne zu viel Leistung zu verlieren?

Question

Ich versuche, die schnell Exp (x) Funktion aus, die zuvor in this Antwort auf eine Frage SO auf die Verbesserung der Rechengeschwindigkeit in C# beschrieben wurde:

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

Der Ausdruck wird mit einigen IEEE Floating-Point „Tricks“ und ist hauptsächlich für den Einsatz in neuralen Sets vorgesehen. Die Funktion ist ca. 5 mal schneller als die normale Math.Exp(x) Funktion.

Leider ist die numerische Genauigkeit nur -4% - + 2% relativ zur regulären Math.Exp(x) Funktion, idealerweise hätte ich gerne eine Genauigkeit innerhalb von mindestens dem Sub-Prozent-Bereich.

Ich habe den Quotienten zwischen den approximierten und den regulären Exp-Funktionen aufgetragen, und wie in der Grafik zu sehen ist, scheint der relative Unterschied mit praktisch konstanter Frequenz wiederholt zu werden.

Ist es möglich, die Vorteile dieser Regelmäßigkeit, um die Richtigkeit der „schnellen exp“ Funktion weiter ohne wesentliche Verringerung der Rechengeschwindigkeit, oder würde der Rechenaufwand einer Genauigkeit Verbesserung schwerer wiegen als die rechnerische Gewinn zu verbessern des ursprünglichen Ausdrucks?

(Als Randbemerkung, habe ich versucht, auch one of the alternative in der gleichen Frage SO vorgeschlagenen Ansätze, aber dieser Ansatz scheint nicht rechnerisch in C# effizient zu sein, zumindest nicht für den allgemeinen Fall.)

UPDATE MAI 14

Auf Anfrage von @Adriano, habe ich jetzt einen sehr einfachen Benchmark durchgeführt. Ich habe 10 Millionen Berechnungen mit jeder der exp Funktionen für Fließkommawerte im Bereich [-100, 100] durchgeführt. Da der Bereich der Werte, die mich interessieren, von -20 bis 0 reicht, habe ich auch explizit den Funktionswert bei x = -5 aufgelistet. Hier sind die Ergebnisse:

 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 entspricht der Exp Funktion am Anfang dieses Textes angegeben. ExpSeries8 ist die formulation, die ich ursprünglich behauptete, war nicht sehr effizient auf .NET; Bei der Implementierung wie Neil war es tatsächlich sehr schnell. ExpSeries16 ist die analoge Formel aber mit 16 Multiplikationen statt 8. exp1 bis exp7 sind die verschiedenen Funktionen von Adriano Antwort unten. Die letzte Variante exp7 ist eine Variante, bei der das Vorzeichen x geprüft wird; Wenn negativ, gibt die Funktion stattdessen 1/exp(-x) zurück.

Leider sind keine der expN Funktionen, die von Adriano aufgeführt werden, in dem breiteren negativen Wertebereich, den ich in Betracht ziehe, ausreichend. Der Serienexpansionsansatz von Neil Coffey scheint in "meinem" Wertebereich besser geeignet zu sein, obwohl er bei größeren negativen Multiplikationen zu schnell divergiert, insbesondere wenn "nur" 8 Multiplikationen verwendet werden.

Answer 1

Falls jemand will, um die relative Fehlerfunktion in der Frage gezeigt replizieren, ist hier eine Möglichkeit, Matlab (der „schnelle“ Exponent ist nicht sehr schnell in Matlab, aber es ist genau):

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

Nun stimmt die Fehlerdauer exakt überein, wenn der Binärwert tmp von der Mantisse in den Exponenten überläuft. Lassen Sie sich unsere Daten in Bins durchbrechen, indem die Bits zu verwerfen, dass der Exponent wird (es periodisch zu machen), und nur die hohe verbleibenden acht Bits zu halten (zu unserer Lookup-Tabelle eine vernünftige Größe machen):

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

Jetzt berechnen wir die mittlere erforderliche Anpassung:

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

Der relative Fehler wird auf +/- 0,0006 verringert. Natürlich sind auch andere Tabellengrößen möglich (zum Beispiel ergibt eine 6-Bit-Tabelle mit 64 Einträgen +/- 0,0025) und der Fehler ist in der Tabellengröße fast linear. Lineare Interpolation zwischen Tabelleneinträgen würde den Fehler noch weiter verbessern, jedoch auf Kosten der Leistung. Da wir das Genauigkeitsziel bereits erreicht haben, vermeiden wir weitere Leistungseinbußen.

An diesem Punkt ist es einige triviale Editor-Fähigkeiten, die von MatLab berechneten Werte zu nehmen und eine Nachschlagetabelle in C# zu erstellen. Für jede Berechnung fügen wir eine Bitmaske, eine Tabellensuche und eine Multiplikation mit doppelter Genauigkeit hinzu.

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

Der Speedup ist sehr ähnlich dem Original-Code - für den Computer, das ist ca. 30% schneller als x86 kompiliert und etwa 3x so schnell für x64. Mit Mono auf Ideon ist es ein beträchtlicher Nettoverlust (aber auch the original).

komplette Quellcode und Testfall: 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()); 
     } 
    } 
} 

Answer 2

Versuchen Sie folgende Alternativen (exp1 ist schneller, exp7 ist genauer).

-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; 
} 

Precision

 
Function  Error in [-1...1]    Error in [3.14...3.14] 

exp1   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  0.0083%   0.9488   4.10% 
exp5   0.0000279  0.001%   0.3516   1.50% 
exp6   0.0000031  0.00011%  0.1172   0.50% 
exp7   0.0000003  0.000011%  0.0355   0.15% 

Credits
Diese Implementierungen von exp() wurden von "scoofy" mit Taylor-Reihe von einer tanh() Umsetzung von „fuzzpilz berechnet "(Wer auch immer sie sind Ich hatte gerade diese Referenzen auf meinem Code).

Answer 3

Ich habe die paper von Nicol Schraudolph untersucht, wo die ursprüngliche C-Implementierung der oben genannten Funktion im Detail jetzt definiert wurde. Es scheint, dass es wahrscheinlich nicht möglich ist, im Wesentlichen die Genauigkeit der exp Berechnung zu genehmigen, ohne die Leistung ernsthaft zu beeinträchtigen. Andererseits gilt die Approximation auch für große Beträge von x bis zu +/- 700, was natürlich vorteilhaft ist.

Die obige Funktionsimplementierung ist darauf abgestimmt, einen minimalen mittleren quadratischen Fehler zu erhalten. Schraudolph beschreibt, wie der additive Term in dem Ausdruck geändert werden kann, um alternative Approximationseigenschaften zu erhalten.

"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 

Er weist auch darauf hin, dass bei einer „mikroskopisch“ -Niveau der ungefähren „exp“ -Funktion Stufen-Verhalten zeigt, da 32 Bits in der Umwandlung von langen zu Doppeln verworfen. Dies bedeutet, dass die Funktion stückweise in einem sehr kleinen Maßstab konstant ist, aber die Funktion wird zumindest nie mit zunehmendem x abnehmen.

Answer 4

Der folgende Code sollte die Genauigkeitsanforderungen erfüllen, da für die Eingaben in [-87,88] die Ergebnisse einen relativen Fehler haben < = 1.73e-3. Ich kenne C# nicht, also ist das C-Code, aber die Konvertierung sollte einfach sein.

Ich gehe davon aus, dass, da die Genauigkeitsanforderung niedrig ist, die Verwendung der Single-Precision-Berechnung in Ordnung ist. Ein klassischer Algorithmus wird verwendet, in dem die Berechnung von exp() auf die Berechnung von exp2() abgebildet wird. Nach der Konvertierung des Arguments durch Multiplikation mit log2 (e) wird die Potenzierung durch den Bruchteil mit einem Minimax-Polynom Grad 2 behandelt, während die Exponentiation durch den ganzzahligen Teil des Arguments durch direkte Manipulation des Exponententeils des IEEE-754-Singles erfolgt -präzise Nummer.

Die volatile union erleichtert die Neuinterpretation eines Bitmusters als Ganzzahl oder Gleitkommazahl mit einfacher Genauigkeit, die für die Exponentenmanipulation benötigt wird. Es sieht so aus, als ob C# dafür entschiedene Neuinterpretationsfunktionen bietet, die viel sauberer sind.

Die beiden möglichen Leistungsprobleme sind die Funktion floor() und float-> int. Traditionell waren beide auf x86 wegen der Notwendigkeit, den dynamischen Prozessorzustand zu handhaben, langsam. Aber SSE (insbesondere SSE 4.1) stellt Anweisungen bereit, die diese Operationen erlauben, schnell zu sein. Ich weiß nicht, ob C# diese Anweisungen verwenden kann.

/* 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; 
} 

Answer 5

Taylorreihe Annäherungen (wie die expX() Funktionen in Adriano's answer) sind die genauesten in der Nähe von Null und große Fehler bei -20 oder sogar -5 haben. Wenn die Eingabe einen bekannten Bereich hat, z. B. -20 bis 0 wie die ursprüngliche Frage, können Sie eine kleine Nachschlagetabelle und eine zusätzliche Multiplikation verwenden, um die Genauigkeit erheblich zu verbessern.

Der Trick besteht darin zu erkennen, dass exp() in ganze und gebrochene Teile getrennt werden kann. Zum Beispiel:

Der Bruchteil wird immer zwischen -1 und 1 liegen, so dass eine Taylor-Näherung ziemlich genau ist. Der ganzzahlige Teil hat nur 21 mögliche Werte für exp (-20) bis exp (0), so dass diese in einer kleinen Nachschlagetabelle gespeichert werden können.