快速计算:可以提高准确性而不会损失太多性能？
And*_*son 20 c# performance floating-accuracy exp
我正在尝试快速的Exp(x)函数,这个函数之前在这个回答中描述了一个提高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函数之间的商,并且从图中可以看出,相对差异似乎以几乎恒定的频率重复.
是否有可能利用这种规律性来进一步提高"快速exp"功能的准确性而不会显着降低计算速度,或者精度提高的计算开销是否会超过原始表达式的计算增益？
(作为旁注,我也尝试过在同一个SO问题中提出的替代方法之一,但这种方法在C#中似乎没有计算效率,至少在一般情况下并非如此.)
5月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等效于本文开头指定的Exp函数.ExpSeries8是我最初声称在.NET上效率不高的公式 ; 当它像Neil一样实现它实际上非常快.ExpSeries16是类似的公式,但有16个乘法而不是8 EXP1通过exp7是从下阿德里亚诺的回答的不同的功能.exp7的最终变体是检查x的符号的变体; 如果为负,则函数返回1/exp(-x).
不幸的是,Adriano列出的expN函数都不足以在我考虑的更广泛的负值范围内.Neil Coffey的系列扩展方法似乎更适合于"我的"值范围,尽管它与较大的负x一起过快地发散,特别是当使用"仅"8次乘法时.
Adr*_*tti 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()通过"scoofy"使用泰勒系列从tanh()"fuzzpilz" 的实现来计算的(无论他们是谁,我只是在我的代码上有这些引用).
- "fuzzpilz"大声笑.有些人对昵称有一种奇怪的感觉. (3认同)
- 原始泰勒系列近似由scoofy@inf.elte.hu在这里:http://www.musicdsp.org/showone.php？id = 222 - Upvoted,因为它是通过泰勒系列的简单解决方案,惊讶它尚未发布之前. (3认同)
如果有人想要复制问题中显示的相对误差函数,这里有一种使用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从尾数溢出到指数的时间完全一致.让我们通过丢弃成为指数的位(使其成为周期性)将数据分成二进制位,并仅保留高八位(使我们的查找表合理大小):
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);
相对误差降低到+/- .0006.当然,其他表格大小也是可能的(例如,64位条目的6位表格给出+/- .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];
}
加速与原始代码非常相似 - 对于我的计算机,编译为x86的速度提高约30%,x64的速度提高约3倍.使用mono on ideone,这是一个巨大的净亏损(但原始版本也是如此).
完整的源代码和测试用例: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());
}
}
}
- @BenVoigt您是否有C的实现？ (2认同)
泰勒级数近似(例如阿德里亚诺答案中的expX()函数)在零附近最精确,在-20或甚至-5时可能会出现巨大误差.如果输入具有已知范围,例如原始问题的-20到0,则可以使用小型查找表和一个额外的乘法来大大提高准确性.
诀窍是认识到exp()可以分为整数和小数部分.例如:
exp(-2.345) = exp(-2.0) * exp(-0.345)
小数部分将始终在-1和1之间,因此泰勒级数近似将非常准确.整数部分对于exp(-20)到exp(0)只有21个可能的值,因此它们可以存储在一个小的查找表中.
以下代码应满足精度要求，因为对于[-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;
}