Szybkie obliczanie Exp: czy można poprawić dokładność bez utraty zbyt dużej wydajności?
Wypróbowuję funkcję fast Exp (x), która wcześniej była opisana w this odpowiedź na pytanie SO dotyczące poprawy szybkości obliczeń w C#:
public static double Exp(double x)
{
var tmp = (long)(1512775 * x + 1072632447);
return BitConverter.Int64BitsToDouble(tmp << 32);
}
Wyrażenie używa pewnych "sztuczek" zmiennoprzecinkowych IEEE i jest przeznaczone przede wszystkim do użycia w zestawach neuronowych. Funkcja jest około 5 razy szybsza niż zwykła funkcja Math.Exp(x).
Niestety, dokładność liczbowa wynosi tylko -4% -- +2% w stosunku do zwykłej funkcji Math.Exp(x), najlepiej chciałbym mieć dokładność w zakresie co najmniej poniżej procentu.
Wykreśliłem iloraz pomiędzy przybliżoną i regularną funkcją Exp i jak widać na wykresie, względna różnica wydaje się powtarzać z praktycznie stałą częstotliwością.
Czy możliwe jest wykorzystanie tej regularności do dalszej poprawy dokładności funkcji "fast exp" bez znacznego zmniejszenia prędkości obliczeniowej, czy też narzut obliczeniowy poprawa dokładności przewyższa zysk obliczeniowy oryginalnego wyrażenia?
(na marginesie, próbowałem również jednego z alternatywnych metod zaproponowanych w tym samym pytaniu SO, ale to podejście nie wydaje się być wydajne obliczeniowo w C#, przynajmniej nie w ogólnym przypadku.)
AKTUALIZACJA 14 MAJA
Na prośbę @Adriano, wykonałem teraz bardzo prosty benchmark. Wykonałem 10 milionów obliczeń przy użyciu każdego z alternatywne exp funkcje dla wartości zmiennoprzecinkowych z zakresu [-100, 100]. Ponieważ zakres wartości interesują mnie zakresy od -20 do 0, również jawnie wymieniłem wartość funkcji na X = -5. Oto wyniki:
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 {[30] } jest odpowiednikiem Exp funkcji określonej na początku tego tekstu. ExpSeries8 to sformułowanie , które pierwotnie twierdziłem, że nie było zbyt wydajne na. Net; podczas implementacji dokładnie tak jak Neil, było bardzo szybko. ExpSeries16 jest analogicznym wzorem, ale z 16 mnożeniami zamiast 8. exp1 do exp7 są różne funkcje z odpowiedzi Adriano poniżej. Ostateczny wariant exp7 jest wariantem, w którym sprawdzany jest znak x ; jeśli jest ujemny, funkcja zwraca 1/exp(-x).
Niestety, żadna z funkcji expN wymienionych przez Adriano nie jest wystarczająca w szerszy zakres wartości ujemnych rozważam. Podejście do rozszerzania serii przez Neila Coffey'a wydaje się bardziej odpowiednie w "moim" przedziale wartości, chociaż zbyt szybko odbiegają od większego ujemnegox , zwłaszcza gdy używa się "tylko" 8 mnożeń.
5 answers
W przypadku, gdy ktoś chce replikować względną funkcję błędu pokazaną w pytaniu, oto sposób za pomocą Matlaba ("szybki" wykładnik nie jest bardzo szybki w Matlabie, ale jest dokładny):
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);
Teraz okres błędu dokładnie pokrywa się z okresem, w którym wartość binarna tmp przepełnia się z mantysy do wykładnika. Podzielmy nasze dane na pojemniki, odrzucając bity, które stają się wykładnikiem( czyniąc je okresowymi), i zachowując tylko wysokie osiem pozostałych bitów (aby nasze {
index = bitshift(bitand(t,uint64(2^20-2^12)),-12) + 1;
Teraz obliczamy średnią wymaganą korektę:
relerrfix = ex./et;
adjust = NaN(1,256);
for i=1:256; adjust(i) = mean(relerrfix(index == i)); end;
et2 = et .* adjust(index);
Błąd względny jest zmniejszany do+/ -.0006. Oczywiście możliwe są również inne rozmiary tabel (na przykład 6-bitowa tabela z 64 wpisami daje +/ -.0025) i błąd jest prawie liniowy w wielkości tabeli. Interpolacja liniowa między wpisami w tabeli jeszcze bardziej poprawiłaby błąd, ale kosztem wydajności. Skoro już osiągnęliśmy cel dokładności, unikajmy wszelkich kolejne hity.
W tym momencie jest kilka trywialnych umiejętności edytorskich, aby wziąć wartości obliczone przez MatLab i utworzyć tabelę wyszukiwania w C#. Dla każdego obliczenia dodajemy maskę bitową, przeszukiwanie tabeli i mnożenie podwójnej precyzji.
static double FastExp(double x)
{
var tmp = (long)(1512775 * x + 1072632447);
int index = (int)(tmp >> 12) & 0xFF;
return BitConverter.Int64BitsToDouble(tmp << 32) * ExpAdjustment[index];
}
Speedup jest bardzo podobny do oryginalnego kodu-dla mojego komputera jest to o 30% szybsze skompilowane jako x86 I O 3x tak szybkie dla x64. Z Mono na ideone, to spora strata netto (ale tak jest oryginał ).
Kompletny kod źródłowy i testcase: 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());
}
}
}
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2014-10-08 00:57:33
Spróbuj podążać za alternatywami (exp1 jest szybsze, exp7 jest bardziej precyzyjne).
Kod
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;
}
Precyzja
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%
Kredyty
Te implementacje exp() zostały obliczone przez" scoofy "używając szeregów Taylora z tanh() implementacji "fuzzpilz" (kimkolwiek oni są, właśnie miałem te odniesienia na moim kodzie).
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2013-01-03 16:24:42
Przybliżenia szeregów Taylora (takie jak funkcje expX() w odpowiedzi Adriano ) są najdokładniejsze w pobliżu zera i mogą mieć ogromne błędy przy -20 lub nawet -5. Jeśli Dane wejściowe mają znany zakres, taki jak -20 do 0, jak w przypadku pierwotnego pytania, możesz użyć małej tabeli wyszukiwania i jednego dodatkowego pomnożenia, aby znacznie poprawić dokładność.
Sztuczka polega na rozpoznaniu, że EXP() można rozdzielić na części całkowite i ułamkowe. Na przykład:
exp(-2.345) = exp(-2.0) * exp(-0.345)
Część ułamkowa będzie zawsze będzie między -1 a 1, więc przybliżenie szeregów Taylora będzie dość dokładne. Część całkowita ma tylko 21 możliwych wartości dla exp (-20) do exp (0), więc można je zapisać w małej tabeli wyszukiwania.
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2017-05-23 10:32:39
Przestudiowałem Artykuł autorstwa Nicola Schraudolpha, gdzie oryginalna implementacja C powyższej funkcji została teraz bardziej szczegółowo zdefiniowana. Wydaje się, że prawdopodobnie nie jest możliwe znaczne zatwierdzenie dokładności obliczeń exp bez poważnego wpływu na wydajność. Z drugiej strony, przybliżenie jest ważne Również dla dużych magnitudo x , do + / - 700, co jest oczywiście korzystne.
Powyższa implementacja funkcji jest dostrojony, aby uzyskać minimalny średni błąd kwadratowy pierwiastka. Schraudolph opisuje, w jaki sposób termin addytywny w wyrażeniu tmp może być zmieniony w celu uzyskania alternatywnych właściwości aproksymacji.
"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
Zwraca również uwagę, że na poziomie" mikroskopowym "przybliżona funkcja" exp " wykazuje zachowanie stair-case, ponieważ 32 bity są odrzucane w konwersji z long do double. Oznacza to, że funkcja jest stałą w bardzo małej skali, ale funkcja jest przynajmniej nigdy nie maleje wraz ze wzrostem x .
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2012-06-12 18:41:25
Poniższy kod powinien uwzględniać wymagania dotyczące dokładności, ponieważ dla danych wejściowych w [-87,88] wyniki mają błąd względny
Zakładam, że ponieważ wymóg dokładności jest niski, użycie obliczeń z pojedynczą precyzją jest w porządku. Używany jest klasyczny algorytm, w którym obliczenia exp() są odwzorowywane na obliczenia exp2(). Po konwersji argumentu przez mnożenie przez log2 (e), wykładnik przez część ułamkową jest obsługiwany za pomocą wielomianu minimax stopnia 2, podczas gdy wykładnik przez integralną część argumentu jest wykonywany przez bezpośrednią manipulację częścią wykładniczą liczby pojedynczej precyzji IEEE-754.
Związek zmiennoprzecinkowy umożliwia ponowną interpretację wzorca bitowego jako liczby całkowitej lub liczby zmiennoprzecinkowej o pojedynczej precyzji, potrzebnej do manipulacji wykładnikiem. Wygląda na to, że C# oferuje dekodowane funkcje reinterpretacji dla tego, co jest dużo czystsze.
Dwa potencjalne problemy z wydajnością to funkcja float () i konwersja float- > int. Tradycyjnie oba były wolne na x86 ze względu na konieczność obsługi dynamicznego stanu procesora. Jednak SSE (w szczególności SSE 4.1) udostępnia instrukcje, które pozwalają na szybkie wykonywanie tych operacji. Nie wiem, czy C# może korzystać z tych instrukcji.
/* 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;
}
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2012-05-29 02:06:57