Schnelle Exp-Berechnung: Ist es möglich, die Genauigkeit zu verbessern, ohne zu viel Leistung zu verlieren?

Question

Ich probiere gerade die schnelle Funktion Exp(x) aus, die zuvor in diese Antwort auf eine SO-Frage zur Verbesserung der Rechengeschwindigkeit in C#:

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

Der Ausdruck verwendet einige IEEE-Gleitkomma-"Tricks" und ist in erster Linie für die Verwendung in neuronalen Gruppen gedacht. Die Funktion ist etwa 5 Mal schneller als die reguläre Math.Exp(x) Funktion.

Leider liegt die numerische Genauigkeit nur bei -4% bis +2% im Vergleich zur normalen Math.Exp(x) Funktion möchte ich im Idealfall eine Genauigkeit im unteren Prozentbereich erreichen.

Ich habe den Quotienten zwischen der approximativen und der regulären Exp-Funktion aufgetragen, und wie in der Grafik zu sehen ist, scheint sich der relative Unterschied mit praktisch konstanter Häufigkeit zu wiederholen.

Ist es möglich, diese Regelmäßigkeit zu nutzen, um die Genauigkeit der Funktion "fast exp" weiter zu verbessern, ohne die Berechnungsgeschwindigkeit wesentlich zu verringern, oder würde der Rechenaufwand einer Genauigkeitsverbesserung den rechnerischen Gewinn des ursprünglichen Ausdrucks überwiegen?

(Nebenbei bemerkt, habe ich auch versucht eine der Alternativen Ansätze, die in der gleichen SO-Frage vorgeschlagen wurden, aber dieser Ansatz scheint in C# nicht rechnerisch effizient zu sein, zumindest nicht für den allgemeinen Fall).

UPDATE 14. MAI

Auf Anfrage von @Adriano habe ich nun einen sehr einfachen Benchmark durchgeführt. Ich habe 10 Millionen Berechnungen mit jeder der folgenden Alternativen durchgeführt exp Funktionen für Fließkommazahlen im Bereich [-100, 100]. Da der Wertebereich, an dem ich interessiert bin, von -20 bis 0 reicht, habe ich auch den Funktionswert bei x = -5 explizit aufgeführt. 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 ist gleichbedeutend mit dem Exp Funktion, die am Anfang dieses Textes beschrieben wird. ExpSerie8 es el Formulierung die ich ursprünglich als nicht sehr effizient unter .NET bezeichnet hatte; bei der Implementierung genau wie Neil war sie tatsächlich sehr schnell. ExpSerie16 ist die analoge Formel, aber mit 16 Multiplikationen anstelle von 8. exp1 über exp7 sind die verschiedenen Funktionen aus Adriano's Antwort unten. Die letzte Variante von exp7 ist eine Variante, bei der das Vorzeichen von x wird geprüft; bei negativem Ergebnis gibt die Funktion 1/exp(-x) stattdessen.

Leider ist keiner der beiden expN Funktionen, die von Adriano sind in dem von mir in Betracht gezogenen breiteren negativen Wertebereich ausreichend. Der Ansatz der Reihenentwicklung von Neil Coffey scheint in "meinem" Wertebereich besser geeignet zu sein, obwohl er bei größeren negativen Werten zu schnell abweicht x , insbesondere wenn "nur" 8 Multiplikationen verwendet werden.

Answer 1

Taylor-Reihen-Approximationen (wie die expX() Funktionen in Adriano's Antwort ) sind in der Nähe von Null am genauesten und können bei -20 oder sogar -5 große Fehler aufweisen. Wenn die Eingabe einen bekannten Bereich hat, z. B. -20 bis 0 wie in der ursprünglichen 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 einen ganzzahligen und einen gebrochenen Teil aufgeteilt werden kann. Zum Beispiel:

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

Der Bruchteil wird immer zwischen -1 und 1 liegen, so dass eine Taylor-Reihen-Approximation 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.

Answer 2

Versuchen Sie folgende Alternativen ( exp1 ist schneller, exp7 ist präziser).

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äzision

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%

Kredite
Diese Implementierungen von exp() wurden von "scoofy" anhand von Taylor-Reihen aus einer tanh() Implementierung von "fuzzpilz" (wer auch immer das ist, ich hatte gerade diese Referenzen in meinem Code).

Answer 3

Für den Fall, dass jemand die in der Frage gezeigte relative Fehlerfunktion nachbilden möchte, hier ein Weg mit Matlab (der "schnelle" Exponent ist in Matlab nicht sehr schnell, aber er 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);

Der Zeitraum des Fehlers fällt nun genau mit dem Zeitpunkt zusammen, zu dem der Binärwert von tmp von der Mantisse in den Exponenten überläuft. Wir teilen unsere Daten in Bins auf, indem wir die Bits, die zum Exponenten werden, verwerfen (wodurch sie periodisch werden) und nur die hohen acht verbleibenden Bits behalten (um unsere Nachschlagetabelle auf eine vernünftige Größe zu bringen):

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

Nun berechnen wir den mittleren Anpassungsbedarf:

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 +/- .0006 verringert. Natürlich sind auch andere Tabellengrößen möglich (z. B. ergibt eine 6-Bit-Tabelle mit 64 Einträgen +/- 0,0025), und der Fehler ist fast linear zur Tabellengröße. Eine lineare Interpolation zwischen den Tabelleneinträgen würde den Fehler noch weiter verbessern, allerdings auf Kosten der Leistung. Da wir das Ziel der Genauigkeit bereits erreicht haben, sollten wir weitere Leistungseinbußen vermeiden.

An dieser Stelle ist es mit ein paar trivialen Editor-Fähigkeiten möglich, die von MatLab berechneten Werte zu nehmen und eine Nachschlagetabelle in C# zu erstellen. Für jede Berechnung fügen wir eine Bitmaske, eine Tabellennachschlagefunktion 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 Geschwindigkeitszuwachs ist dem des Originalcodes sehr ähnlich - auf meinem Computer ist dies etwa 30% schneller als x86 kompiliert und etwa 3x so schnell für x64. Mit Mono auf Ideone ist es ein erheblicher Nettoverlust (aber das ist das Original ).

Vollständiger 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 4

Der folgende Code sollte die Genauigkeitsanforderungen erfüllen, da die Ergebnisse für Eingaben in [-87,88] einen relativen Fehler <= 1,73e-3 aufweisen. Da ich C# nicht kenne, handelt es sich um C-Code, aber ich nehme an, dass die Konvertierung recht einfach sein sollte.

Ich gehe davon aus, dass aufgrund der geringen Genauigkeitsanforderungen die Verwendung von Berechnungen mit einfacher Genauigkeit in Ordnung ist. Es wird ein klassischer Algorithmus verwendet, bei dem die Berechnung von exp() auf die Berechnung von exp2() abgebildet wird. Nach der Umwandlung des Arguments durch Multiplikation mit log2(e) wird die Potenzierung des gebrochenen Teils mit einem Minimax-Polynom vom Grad 2 durchgeführt, während die Potenzierung des ganzzahligen Teils des Arguments durch direkte Manipulation des Exponententeils der IEEE-754-Zahl mit einfacher Genauigkeit erfolgt.

Die flüchtige Vereinigung erleichtert die Neuinterpretation eines Bitmusters entweder als Ganzzahl oder als Gleitkommazahl mit einfacher Genauigkeit, die für die Manipulation des Exponenten benötigt wird. Es sieht so aus, als ob C# hierfür dezidierte Uminterpretationsfunktionen anbietet, die viel sauberer sind.

Die beiden potenziellen Leistungsprobleme sind die Funktion floor() und die Konvertierung float->int. Traditionell waren beide auf x86 langsam, da sie einen dynamischen Prozessorstatus verarbeiten mussten. Aber SSE (insbesondere SSE 4.1) bietet Anweisungen, die diese Operationen schnell machen. Ich weiß nicht, ob C# diese Befehle nutzen 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

Ich habe die Papier von Nicol Schraudolph, in der die ursprüngliche C-Implementierung der obigen Funktion nun genauer definiert wurde. Es scheint, dass es wahrscheinlich nicht möglich ist, die Genauigkeit der Funktion wesentlich zu verbessern. exp Berechnungen, ohne die Leistung stark zu beeinträchtigen. Andererseits gilt die Näherung auch für große Größenordnungen von x bis zu +/- 700, was natürlich von Vorteil ist.

Die obige Funktionsimplementierung ist so abgestimmt, dass der mittlere quadratische Fehler (Root Mean Square Error) minimal ist. Schraudolph beschreibt, wie der additive Term in der tmp Ausdruck kann geändert werden, um alternative Näherungseigenschaften zu erreichen.

"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 auf "mikroskopischer" Ebene die ungefähre "exp"-Funktion ein treppenförmiges Verhalten aufweist, da 32 Bits bei der Umwandlung von lang a doppelt . Dies bedeutet, dass die Funktion auf einer sehr kleinen Skala stückweise konstant ist, aber die Funktion ist zumindest niemals abnehmend mit zunehmendem x .