2012-05-11 57 views
20

我想出來的快速EXP(x)中先前在this答案在C#描述的SO問題改善計算速度函數:快速計算:可以提高精度而不會損失太多性能?

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

的表達使用一些IEEE浮點「技巧」和主要用於神經系統。該功能比常規的Math.Exp(x)功能快大約5倍。

不幸的是,相對於常規的Math.Exp(x)函數,數值精度僅爲-4% - + 2%,理想情況下,我希望在至少百分比範圍內具有準確性。

我已經繪製了近似和正則Exp函數之間的商,並且可以在圖中看到相對差異似乎以幾乎恆定的頻率重複。

Quotient between fast and regular exp function

是否有可能利用這個規律性的優點是改進了「快EXP」函數的準確性進一步而不顯着降低運算速度,或將的精度的提高計算開銷大於計算增益原始表達?

(作爲一個方面說明,我也試圖在同一個SO問題提出one of the alternative的方法,但這種方法似乎並不在C#計算效率,至少對於一般的情況下)。

更新可能14

根據@Adriano的請求,我現在執行了一個非常簡單的基準測試。我已經使用替代exp函數爲浮點值在[-100,100]範圍內執行了1000萬次計算。由於我感興趣的範圍從-20到0,我也明確列出了x = -5的函數值。下面是結果:

 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等同於本文的開頭指定的精通功能。 ExpSeries8是我原來聲稱在.NET上效率不高的formulation;當它像Neil一樣實施時,它實際上非常快。 ExpSeries16是類似的公式,但有16個乘法而不是8個。exp1exp7是與以下Adriano的答案不同的功能。 exp7的最終變體是其中檢查了x的符號的變體;如果爲負,則函數返回1/exp(-x)

不幸的是,所列的expN功能在我正在考慮的更廣泛的負值範圍內都不足。 Neil Coffey的系列擴展方法似乎更適合於「我的」數值範圍,儘管它與更大的負值分離太快,尤其是在使用「僅」8次乘法時。

+0

我很好奇你對「神經系統」的提及。目前我正在模擬一個使用C++的神經網絡,並面臨着你面臨的同樣的'exp'性能瓶頸。在計算神經科學中有沒有提出近似'exp'功能非常快的論文? – dbliss

回答

7

如果有人想複製的問題顯示的相對誤差函數,這裏是用Matlab的方式(「快」的指數是不是在Matlab非常快,但它是正確的):

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

現在,錯誤的週期恰好與tmp的二進制值從尾數溢出到指數時一致。且讓我們打破數據轉換成箱通過丟棄成爲指數位(使其定期),且只保留高8其餘位(讓我們的查找表,一個合理的規模):

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

現在我們計算意味着所需的調整:

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

的相對誤差減小至+/- 0.0006。當然,其他表格大小也是可能的(例如,一個具有64個入口的6位表格給出+/- 0.0025),並且該誤差在表格大小中幾乎是線性的。表格條目之間的線性插值會進一步改善錯誤,但是會以犧牲性能爲代價。既然我們已經達到了準確性目標,那麼讓我們避免任何進一步的性能命中。

在這一點上,它是一些簡單的編輯技能,採取由MatLab計算的值,並在C#中創建一個查找表。對於每個計算,我們添加一個位掩碼,表查找和雙精度乘法。

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

的加速是非常類似的原代碼 - 爲我的電腦,這是約30%的速度編譯爲x86和3倍左右的速度爲64。隨着對ideone的單聲道,這是一個實質性的淨虧損(但是the original也是如此)。

完整的源代碼和測試用例: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()); 
     } 
    } 
} 
+0

神奇的工作,和一個很好的解釋!非常感謝提供這個答案,這只是我所希望的那種進步。如果您之前已經開發過這個功能,或者您是否因爲這個問題而實現了它? –

+0

@Anders:我完全偷了你在問題中提出的方法。 –

+0

在android NDK中測試之後,它比系統std :: exp()慢。但在個人電腦上速度更快。 (https://gist.github.com/maxint/0172c1dcd075d3589eeb) – maxint

10

請嘗試以下替代方案(exp1更快,exp7更精確)。

代碼

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

精密

 
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% 

現金
exp()這些實現已經從tanh()實施「fuzzpilz的計算由 「scoofy」 使用泰勒級數「(無論他們是誰,我只是在我的代碼上有這些引用)。

+2

「fuzzpilz」大聲笑。有些人對暱稱有陌生感。 –

+3

原始的泰勒系列近似由[email protected]在這裏:http://www.musicdsp.org/showone.php?id=222 - Upvoted,因爲它是通過泰勒系列的一個簡單的解決方案,驚訝它沒有之前發佈。 –

+0

@ MahmoudAl-Qudsi感謝您的參考,很久以前它已經消失了! –

4

我已經研究了Nicol Schraudolph的paper,其中現在更詳細地定義了上述函數的原始C實現。看起來似乎不可能大幅度批准計算的準確性,而不會嚴重影響性能。另一方面,近似值對於大量的x也是有效的,高達+/- 700,這當然是有利的。

上面的函數實現被調整以獲得最小的均方根誤差。 Schraudolph描述瞭如何改變表達式中的附加項以實現替代近似屬性。

"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 

他還指出,在「微觀」電平的近似「EXP」功能表現出階梯情況下行爲,因爲32位被在轉換丟棄從。這意味着該功能在很小的範圍內是分段恆定的,但功能至少永遠不會隨着增加而減小x

3

下面的代碼應該解決精度要求,對於[-87,88]中的輸入結果有相對誤差< = 1.73e-3。我不知道C#,所以這是C代碼,但轉換應該是簡單的故障處理。

我假設由於精度要求低,所以使用單精度計算很好。正在使用一種經典算法,其中將exp()的計算映射到exp2()的計算。在通過乘以log2(e)進行參數轉換之後,小數部分的指數用2級的最小多項式處理,而參數的整數部分的指數通過直接操作IEEE-754單元的指數部分 - 精確數字。

易失聯合促成指數操作所需的位模式重新解釋爲整數或單精度浮點數。它看起來像C#爲此提供了重新解釋函數,這更清晰。

兩個潛在的性能問題是floor()函數和float-> int轉換。傳統上,由於需要處理動態處理器狀態,所以x86上都很慢。但是SSE(特別是SSE 4.1)提供了允許這些操作快速的指令。我不知道C#是否可以使用這些說明。

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

非常感謝一個很好的例子和一個很好的解釋。我將嘗試在C#中轉換您的實現,以查看它與常規_Exp_函數相比的執行情況。我不記得在其他地方看到過這個解決方案,你是否因爲SO問題而提出這個問題? –

+0

我在過去曾多次設計/實現過各種超越函數的算法。我上面選擇的方法非常適合花園多樣化的算法。我確實爲多項式創建了自定義極大極小近似值來回應這個問題。有一些工具可用於Mathematica,Maple和其他工具;通常它們基於Remez算法的變體。 – njuffa

+0

請注意,在C++中使用union並不正確。但是你可以在C和C++中使用'memcpy',並且優化器應該做一些明智的事情,而不用基於嚴格別名的優化來破壞它。 –

9

泰勒級數近似(如Adriano's answerexpX()功能)接近零最準確的,可在-20℃或甚至-5有巨大的錯誤。如果輸入有一個已知的範圍,如-20至0,就像原始問題一樣,您可以使用小型查找表和一個附加乘法來大大提高準確性。

訣竅是認識到exp()可以分爲整數部分和小數部分。例如:

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

小數部分將始終在-1和1之間,所以泰勒級數近似將相當準確。整數部分exp(-20)到exp(0)只有21個可能的值,所以這些可以存儲在一個小查找表中。