2012-05-11 55 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个可能的值,所以这些可以存储在一个小查找表中。