显示带有标签的帖子 快速逼近. 显示所有帖子
显示带有标签的帖子 快速逼近. 显示所有帖子

2011年7月17日,星期日

快速近似兰伯特W

毕竟,这只是个人用途的问题, 兰伯特的W 功能 一旦 在机器学习上下文中,它以$ W(\ exp(x))-x $出现,最好直接近似。尽管如此,生成这些快速的近似函数使娱乐变得有趣,并且与例如进一步发展我的计算机游戏能力相比,对某处某人有用的可能性要高得多。

无分支的生活方式

Lambert W在一个方面是典型的函数:有一个很好的渐近逼近可以用来初始化一个 户主法,但仅适用于部分域。特别是对于大$ x $,\ [
W(x)\ approx \ log(x)-\ log \ log(x)+ \ frac {\ log \ log(x)} {\ log(x)},
\] 很棒,因为 对数便宜。不幸的是,对于$ x \ leq 1 $,无法评估,对于$ 1< x <2 $它给出了非常差的结果。一种自然的处理方式是为$ x使用不同的初始化策略< 2$.
if (x < 2)
  { 
    // alternate initialization
  }
else
  {
    // use asymptotic approximation
  }

// householder steps here
这种简单的方法存在两个问题。在标量上下文中,这可以 挫败流水线架构 现代CPU。在向量上下文中,这甚至更成问题,因为组件可能会落入条件的不同分支中。

该怎么办?事实证明像这样
a = (x < 2) ? b : c
看起来像有条件的,但不是必须的,因为它们可以重写为
a = f (x < 2) * b + (1 - f (x < 2)) * c
$ f $是一个指示符函数,根据参数的真值返回0或1。 上证所 指令集包含以下指标功能 比较测试 ,当与 ``浮点数和'' 指令最终计算出无分支三元运算符。

底线是,如果计算条件的两个分支,则可以使确定性执行具有确定性,并且在足够简单的情况下,可以直接硬件支持快速执行此操作。

无分支兰伯特W

因此,这里的主要思想是为Householderer步骤进行另一种初始化,以便可以对无输入方式进行计算,因为对于大的输入使用渐近逼近。因此,我寻找\\形式的近似值
W(x)\大约a + \ log(c x + d)-\ log \ log(c x + d)+ \ frac {\ log \ log(c x + d)} {\ log(c x + d)},
\] 对于大$ x $,$ a = 0 $,$ c = 1 $和$ d = 0 $。我通过Mathematica找到了$ a $,$ c $,$ d $的值和$ x $的临界值。 (好奇的可以查看 Mathematica笔记本)。矢量版本最终看起来像
// WARNING: this code has been updated.  Do not use this version.
// Instead get the latest version from http://code.google.com/p/fastapprox

static inline v4sf
vfastlambertw (v4sf x)
{
  static const v4sf threshold = v4sfl (2.26445f);

  v4sf under = _mm_cmplt_ps (x, threshold);
  v4sf c = _mm_or_ps (_mm_and_ps (under, v4sfl (1.546865557f)),
                      _mm_andnot_ps (under, v4sfl (1.0f)));
  v4sf d = _mm_and_ps (under, v4sfl (2.250366841f));
  v4sf a = _mm_and_ps (under, v4sfl (-0.737769969f));

  v4sf logterm = vfastlog (c * x + d);
  v4sf loglogterm = vfastlog (logterm);

  v4sf w = a + logterm - loglogterm + loglogterm / logterm;
  v4sf expw = vfastexp (w);
  v4sf z = w * expw;
  v4sf p = x + z;

  return (v4sfl (2.0f) * x + w * (v4sfl (4.0f) * x + w * p)) /
         (v4sfl (2.0f) * expw + p * (v4sfl (2.0f) + w));
}
您可以从 快速约 项目。

时间和准确性

时序测试是通过编译完成的 -O3 -finline功能-fast-math,在运行64位Ubuntu Lucid(所以 gcc 4:4.4.3-1ubuntu1libc 2.11.1-0ubuntu7.6)。我还测量了以$(\ frac {1} {2} U(-1 / e,1)+ \ frac {1} {2} U(0,100))$分布的$ x $的平均相对准确度,即,从两个均匀分布中抽取50-50。将精度与牛顿方法的20次迭代进行比较,当$ x时,初始点为0<5 $,否则渐近逼近。我还测试了 gsl实施 精度更高,但速度明显慢得多。
 Fastlambertw  average relative error = 5.26867e-05
fastlambertw max relative error (at 2.48955e-06) = 0.0631815
fasterlambertw average relative error = 0.00798678
fasterlambertw max relative error (at -0.00122776) = 0.926378
vfastlambertw average relative error = 5.42952e-05
vfastlambertw max relative error (at -2.78399e-06) = 0.0661513
vfasterlambertw average relative error = 0.00783347
vfasterlambertw max relative error (at -0.00125244) = 0.926431
gsl_sf_lambert_W0 average relative error = 5.90309e-09
gsl_sf_lambert_W0 max relative error (at -0.36782) = 6.67586e-07
fastlambertw million calls per second = 21.2236
fasterlambertw million calls per second = 53.4428
vfastlambertw million calls per second = 21.6723
vfasterlambertw million calls per second = 56.0154
gsl_sf_lambert_W0 million calls per second = 2.7433
这些平均精度在最小域中隐藏了相对较差的性能。就在$ x = -e ^ {-1} $,这是该域的最小值, Fastlambertw 近似差(-0.938,而正确答案是-1;因此相对误差为$ 6 \%$);但是在$ x = -e ^ {-1} + \ frac {1} {100} $时,相对误差降为$ 2 \乘以10 ^ {-4} $。

2011年7月6日,星期三

快速逼近集合

我将各种快速(矢量化)近似值捆绑在一起 在Google代码上进行项目.

CISC的复仇

因此,就像我们现在所听到的那样,我的童年Commodore 64时钟速度的千倍增长已不再重复。像其他科学计算社区一样,机器学习社区正认真地拥抱多核和集群计算,以便进一步扩展,这是一件好事。

但是,这项练习告诉我的是,计算机的速度越来越快:不是以时钟速度而言,而是每个时钟滴答完成了多少工作。芯片设计人员不仅会发布新的令人敬畏的指令,而且还会与每个芯片版本一起使用,以提高其效率(例如,更少的延迟,停顿)。如果您使用的是最新版本的编译器,例如gcc 4,则某些收益会自动累积 自动向量化 但是gcc 3没有。这就是为什么我``仅''看到与之相比又增加了20%的原因之一 LDA 实现的向量化 在Vowpal Wabbit中:编译器已经对明显的东西进行了矢量化处理(例如,点积,循环累加计数器),但还不足以识别我的digamma,lgamma和指数近似值也可以进行矢量化处理。如果您使用gcc 3编译Vowpal Wabbit,则可能要付出速度损失。

好吧,那又如何呢?要达到PB级,还是需要一个集群,对吗?那么,存在很多问题,在笔记本电脑上拥有几GB的样本数据对于快速实验很有用。因此,它有助于以最快的速度实现在笔记本电脑上运行的学习算法。要快速实现,请注意指令集!

2011年6月26日,星期日

更快的LDA:第二部分

我能够提高其他效率 Vowpal兔子 的LDA实施 通过 向量化 。基本思想是找到重复的昂贵操作的循环并通过 SIMD 同时执行多个昂贵的操作。特别是我正在利用 上证所 指令系统。

向量化范例

对于我的简单功能 近似对数和指数,转换为SSE是一项最直接的机械任务(如此机械,实际上,我想知道为什么编译器无法为我完成此任务?)。这是一个示例:首先,原始的快速以2为底的对数,
// WARNING: this code has been updated.  Do not use this version.
// Instead get the latest version from http://code.google.com/p/fastapprox

inline float
fastlog2 (float x)
{
  union { float f; uint32_t i; } vx = { x };
  union { uint32_t i; float f; } mx = { (vx.i & 0x007FFFFF) | (0x7e << 23) };
  float y = vx.i;
  y *= 1.0f / (1 << 23);

  return
    y - 124.22544637f - 1.498030302f * mx.f - 1.72587999f / (0.3520887068f + mx.f);
}
接下来是SSE版本
// WARNING: this code has been updated.  Do not use this version.
// Instead get the latest version from http://code.google.com/p/fastapprox

typedef float v4sf __attribute__ ((__vector_size__ (16)));
typedef int v4si __attribute__ ((__vector_size__ (16)));
#define v4sf_dup_literal(x) ((v4sf) { x, x, x, x })
#define v4si_dup_literal(x) ((v4si) { x, x, x, x })
#define v4sf_from_v4si __builtin_ia32_cvtdq2ps

inline v4sf
vfastlog2 (v4sf x)
{
  union { v4sf f; v4si i; } vx = { x };
  union { v4si i; v4sf f; } mx = {
    (vx.i & v4si_dup_literal (0x007FFFFF)) | v4si_dup_literal ((0x7e << 23))
  };
  v4sf y = v4sf_from_v4si (vx.i);
  y *= v4sf_dup_literal ((1.0f / (1 << 23)));
  
  const v4sf c_124_22544637 = v4sf_dup_literal (124.22544637f);
  const v4sf c_1_498030302 = v4sf_dup_literal (1.498030302f);
  const v4sf c_1_725877999 = v4sf_dup_literal (1.72587999f);
  const v4sf c_0_3520087068 = v4sf_dup_literal (0.3520887068f);
  
  return 
    y - c_124_22544637 - c_1_498030302 * mx.f - c_1_725877999 / (c_0_3520087068 + mx.f);
}
vfastlog2 运行速度与 fastlog2 (实际上有时会稍快一些,具体取决于 fastlog2 并使用x87运算进行编译),并计算几乎相同的内容(有时差异很小,比近似误差小得多)。但是,它可以一次执行4次计算。

有一个陷阱:它总是一次执行4次计算,因此,如果数组的长度不是16个字节的倍数,则最后一次计算将从内存的其他位读取(可能导致分段错误)。此外,输入必须对齐16字节,否则会出现分段错误。有人告诉我,在超高性能应用程序中,开发人员总是将阵列排列为16字节对齐,并且是16字节长的倍数。但是,没有简单的方法来安排此操作,因此我不得不进行一些防御性的编码。我使用的模式是
// WARNING: this code has been updated.  Do not use this version.
// Instead get the latest version from http://code.google.com/p/fastapprox

typedef float v4sf_aligned __attribute__ ((__vector_size__ (16))) __attribute__ ((aligned (16)));

float f (float);
v4sf vectorized_f (v4sf);

inline void
conservative_map (float* x,
                  size_t n)
{
  unsigned int i;

  // handle non-aligned prefix
  for (i = 0; i < n && (((uintptr_t) (x + i)) % 16) > 0; ++i)
    { 
      x[i] = f (x[i]);
    }

  // the good stuff
  for (; i + 4 < n; i += 4)
    { 
      * (v4sf_aligned *) (x + i) = vectorized_f (* (v4sf_aligned*) (x + i));
    }

  // handle remainder
  for (; i < n; ++i)
    { 
      x[i] = f (x[i]);
    }
}
对于大约100个主题的地图,这在最坏的情况下会导致25个操作之外的6个附加操作,因此平均可能会有10%的罚款。

结果

我在与我的基准相同的基准上进行了一些其他计时测试 以前的帖子 。 \ [
\ begin {array} {c | c}
\ mbox {代码版本}&\ mbox {总挂钟} \\ \ hline
\ mbox {原始}&\ mbox {79秒} \\
\ mbox {近似值,无SSE}&\ mbox {49秒} \\
\ mbox {近似值,含SSE}&\ mbox {39秒}
\ end {array}
\]
近似加速比矢量化加速大,这很好,因为即使对于不支持SSE的环境,近似加速始终可用。

2011年6月24日,星期五

更快的LDA

现在,我已经掌握了 字节重新解释的邪恶力量 在整数和IEEE浮点数之间,我一直在寻找加速的东西。最近发布的 快速实施LDA 由Smola等等提醒我看Vowpal Wabbit的 LDA 实施 。 的 官方网站 该文件充满了看起来昂贵的先验功能,可以使用近似版本来加快速度,希望不会破坏学习算法。

更快的Log Gamma和Digamma

除了许多我已经有了快速近似的指数和对数调用之外,还有一些对 对数伽玛 地伽玛 。幸运的是, 斯特林近似 是对数,因此我可以利用我的快速对数近似值 以前的博客文章.

我发现对于Log Gamma和Digamma而言,两次应用移位公式可以在精度和速度之间取得最佳折衷。
// WARNING: this code has been updated.  Do not use this version.
// Instead get the latest version from http://code.google.com/p/fastapprox

inline float
fastlgamma (float x)
{
  float logterm = fastlog (x * (1.0f + x) * (2.0f + x));
  float xp3 = 3.0f + x;

  return
    -2.081061466f - x + 0.0833333f / xp3 - logterm + (2.5f + x) * fastlog (xp3);
}

inline float
fastdigamma (float x)
{
  float twopx = 2.0f + x;
  float logterm = fastlog (twopx);

  return - (1.0f + 2.0f * x) / (x * (1.0f + x))
         - (13.0f + 6.0f * x) / (12.0f * twopx * twopx)
         + logterm;
}
时间和准确性
时序测试是通过编译完成的 -O3 -finline功能-fast-math,在运行64位Ubuntu Lucid(所以 gcc 4:4.4.3-1ubuntu1libc 2.11.1-0ubuntu7.6)。我还测量了$ x \ in [1/100,10] $的平均相对准确度(使用 lgammaf 来自libc作为Log Gamma的黄金标准,并且 boost :: math :: digamma 作为Digamma的黄金标准)。
fastlgamma relative accuracy = 0.00045967
fastdigamma relative accuracy = 0.000420604
fastlgamma million calls per second = 60.259
lgammaf million calls per second = 21.4703
boost::math::lgamma million calls per second = 1.8951
fastdigamma million calls per second = 66.0639
boost::math::digamma million calls per second = 18.9672
好吧,Log Gamma的速度大约是3倍(但是为什么boost :: math :: lgamma这么慢?)却略好于Digamma的3倍。

LDA 时序

所以我将这些近似值放入 官方网站 并使用 wiki1K 来自的文件 测试/训练集/ 目录;但是我将文件本身连接了100次以减慢速度。
原始运行
% (cd test; time ../vw --lda 100 --lda_alpha 0.01 --lda_rho 0.01 --lda_D 1000 -b 13 --minibatch 128 train-sets/wiki1Kx100.dat)
your learning rate is too high, setting it to 1
using no cache
Reading from train-sets/wiki1Kx100.dat
num sources = 1
Num weight bits = 13
learning rate = 1
initial_t = 1
power_t = 0.5
learning_rate set to 1
average    since       example  example    current  current  current
loss       last        counter   weight      label  predict features
10.296875  10.296875         3      3.0    unknown   0.0000       38
10.437156  10.577436         6      6.0    unknown   0.0000       14
10.347227  10.239314        11     11.0    unknown   0.0000       32
10.498633  10.650038        22     22.0    unknown   0.0000        2
10.495566  10.492500        44     44.0    unknown   0.0000      166
10.469184  10.442189        87     87.0    unknown   0.0000       29
10.068007  9.666831        174    174.0    unknown   0.0000       17
9.477440   8.886873        348    348.0    unknown   0.0000        2
9.020482   8.563524        696    696.0    unknown   0.0000      143
8.482232   7.943982       1392   1392.0    unknown   0.0000       31
8.106654   7.731076       2784   2784.0    unknown   0.0000       21
7.850269   7.593883       5568   5568.0    unknown   0.0000       25
7.707427   7.564560      11135  11135.0    unknown   0.0000       39
7.627417   7.547399      22269  22269.0    unknown   0.0000       61
7.583820   7.540222      44537  44537.0    unknown   0.0000        5
7.558824   7.533827      89073  89073.0    unknown   0.0000      457

finished run
number of examples = 0
weighted example sum = 0
weighted label sum = 0
average loss = -nan
best constant = -nan
total feature number = 0
../vw --lda 100 --lda_alpha 0.01 --lda_rho 0.01 --lda_D 1000 -b 13 --minibatc  69.06s user 13.60s system 101% cpu 1:21.59 total
近似运行
% (cd test; time ../vw --lda 100 --lda_alpha 0.01 --lda_rho 0.01 --lda_D 1000 -b 13 --minibatch 128 train-sets/wiki1Kx100.dat)               your learning rate is too high, setting it to 1
using no cache
Reading from train-sets/wiki1Kx100.dat
num sources = 1
Num weight bits = 13
learning rate = 1
initial_t = 1
power_t = 0.5
learning_rate set to 1
average    since       example  example    current  current  current
loss       last        counter   weight      label  predict features
10.297077  10.297077         3      3.0    unknown   0.0000       38
10.437259  10.577440         6      6.0    unknown   0.0000       14
10.347348  10.239455        11     11.0    unknown   0.0000       32
10.498796  10.650243        22     22.0    unknown   0.0000        2
10.495748  10.492700        44     44.0    unknown   0.0000      166
10.469374  10.442388        87     87.0    unknown   0.0000       29
10.068179  9.666985        174    174.0    unknown   0.0000       17
9.477574   8.886968        348    348.0    unknown   0.0000        2
9.020435   8.563297        696    696.0    unknown   0.0000      143
8.482178   7.943921       1392   1392.0    unknown   0.0000       31
8.106636   7.731093       2784   2784.0    unknown   0.0000       21
7.849980   7.593324       5568   5568.0    unknown   0.0000       25
7.707124   7.564242      11135  11135.0    unknown   0.0000       39
7.627207   7.547283      22269  22269.0    unknown   0.0000       61
7.583952   7.540696      44537  44537.0    unknown   0.0000        5
7.559260   7.534566      89073  89073.0    unknown   0.0000      457

finished run
number of examples = 0
weighted example sum = 0
weighted label sum = 0
average loss = -nan
best constant = -nan
total feature number = 0
../vw --lda 100 --lda_alpha 0.01 --lda_rho 0.01 --lda_D 1000 -b 13 --minibatc  41.97s user 7.69s system 102% cpu 48.622 total
因此,从81秒降低到49秒,速度提高了大约40%。但这有效吗?好吧,我让马特(Matt)担任最终裁判(我给他寄了一封修改过的 官方网站 (用于测试),但初步迹象表明,使用昂贵功能的近似版本不会破坏学习过程。

下一步:矢量化

肯定会有更多的速度,因为 官方网站 加载像
for (size_t i = 0; i<global.lda; i++)
    {
      sum += gamma[i];
      gamma[i] = mydigamma(gamma[i]);
    }
只是为了尖叫 向量化 。 敬请关注!

2011年6月20日,星期一

快速近似对数,指数,幂和逆根

受启发 伊恩·斯蒂芬森(Ian Stephenson) ,我发现了一个更快更准确的$ \ log_2 $近似值,这反过来又提高了我对$ p $根的反向估计。所产生的技巧通常对于计算能力$ x ^ p $运作良好,这也是计算$ L_p $范数的一部分。还有一个快速的近似对数和指数,它们在广义上相关,因为这些运算在实践中被广泛地用于 表示概率.

更新资料 :这里的代码示例已过期;小问题已修复,您应该从 Google代码存储库 。 ( 更新资料 2:Google代码已失效,但是 版本存在于github

快速近似对数

快速的\ log_2 $近似值意味着对自然对数或编译时已知的任何其他基数的快速近似值,因为差是一个乘法。因此,考虑一个浮点​​数$ x $和值$ I_x $,它们将$ x $的字节表示形式解释为整数\ [
\ begin {aligned}
x&=(1 + m_x)2 ^ {e_x},\\
\ log_2(x)&= e_x + \ log_2(1 + m_x),\\
I_x&=(e_x + B)L + m_x L,\\
I_x / L-B&= e_x + m_x,\\
\ log_2(x)&= I_x / L-B + \ log_2(1 + m_x)-m_x。
\ end {aligned}
\] 因此,需要一种从$ x $中提取尾数并近似为[0,1)$中的$ z \ g_2(1 + z)-z $的尾数的方法。我采用了有理近似值,这似乎是准确性和速度的最佳结合点,
// WARNING: this code has been updated.  Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.

inline float
fastlog2 (float x)
{
  union { float f; uint32_t i; } vx = { x };
  union { uint32_t i; float f; } mx = { (vx.i & 0x007FFFFF) | (0x7e << 23) };
  float y = vx.i;
  y *= 1.0 / (1 << 23);

  return
    y - 124.22544637f - 1.498030302f * mx.f - 1.72587999f / (0.3520887068f + mx.f);
}

inline float
fastlog (float x)
{
  return 0.69314718f *  fastlog2  (x);
}
如果您想要更快,更粗糙,则$ g(z)$的平均值为$ \ frac {-2 + \ log(8)} {\ log(4)} $
// WARNING: this code has been updated.  Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.

inline float
fasterlog2 (float x)
{
  union { float f; uint32_t i; } vx = { x };
  float y = vx.i;
  y *= 1.0 / (1 << 23);
  return y - 126.94269504f;
}

inline float
fasterlog (float x)
{
  return 0.69314718f * fasterlog2 (x);
}
时间和准确性
时序测试是通过编译完成的 -O3 -finline功能-fast-math,在运行64位Ubuntu Lucid(所以 gcc 4:4.4.3-1ubuntu1libc 2.11.1-0ubuntu7.6)。我还测量了$ x \ in [1/100,10] $中的平均相对精度。
  fastlog2   relative accuracy = 2.09352e-05
fastlog relative accuracy = 2.09348e-05
fasterlog2 relative accuracy = 0.0130367
fasterlog relative accuracy = 0.0130367
fastlog2 million calls per second = 160.141
fastlog million calls per second = 143.552
fasterlog2 million calls per second = 218.345
fasterlog million calls per second = 210.435
log2f million calls per second = 40.8511

快速近似指数

同样,2是最方便的基础。这里$ 2 ^ x $是我要近似的值,$ I_ {2 ^ x} $是与将$ 2 ^ x $的字节模式解释为整数相关的值。 \ [
\ begin {aligned}
2 ^ x&= 2 ^ {x-\ lfloor x \ rfloor} 2 ^ {\ lfloor x \ rfloor} \\
&=(1 + 2 ^ {x-\ lfloor x \ rfloor}-1)2 ^ {\ lfloor x \ rfloor},\\
I_ {2 ^ x}&=(\ lfloor x \ rfloor + B)L +(2 ^ {x-\ lfloor x \ rfloor}-1)L \\
&= L(x + B-1 + 2 ^ {x-\ lfloor x \ rfloor}-x + \ lfloor x \ rfloor)。
\ end {aligned}
\] 校正项的格式为$ g(z)= 2 ^ z-z $ for $ z \ in [0,1)$,并且我又一次从有理函数逼近中得到了最大的收益。这导致
// WARNING: this code has been updated.  Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.

inline float
fastpow2 (float p)
{
  union { float f; uint32_t i; } vp = { p };
  int sign = (vp.i >> 31);
  int w = p;
  float z = p - w + sign;
  union { uint32_t i; float f; } v = { (1 << 23) * (p + 121.2740838f + 27.7280233f / (4.84252568f - z) - 1.49012907f * z) };
  return v.f;
}

inline float
fastexp (float p)
{
  return fastpow2 (1.442695040f * p);
}
如果您想要更快,更粗糙,则$ g(z)$的平均值为$ \ frac {2-\ log(2)} {2 \ log(2)} $
// WARNING: this code has been updated.  Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.

inline float
fasterpow2 (float p)
{
  union { uint32_t i; float f; } v = { (1 << 23) * (p + 126.94269504f) };
  return v.f;
}

inline float
fasterexp (float p)
{
  return fasterpow2 (1.442695040f * p);
}
时间和准确性
时序测试是通过编译完成的 -O3 -finline功能-fast-math,在运行64位Ubuntu Lucid(所以 gcc 4:4.4.3-1ubuntu1libc 2.11.1-0ubuntu7.6)。我还在[1/20,20] $中测量了$ p的平均相对准确度。当我测量相对准确度时,我尝试了每个$ p $和每个$ p $逆根(即$ -1 / p $)。
fastpow2 relative accuracy (positive p) = 1.58868e-05
fastexp relative accuracy (positive p) = 1.60712e-05
fasterpow2 relative accuracy (positive p) = 0.0152579
fasterexp relative accuracy (positive p) = 0.0152574
fastpow2 relative accuracy (inverse root p) = 1.43517e-05
fastexp relative accuracy (inverse root p) = 1.7255e-05
fasterpow2 relative accuracy (inverse root p) = 0.013501
fasterexp relative accuracy (inverse root p) = 0.0111832
fastpow2 million calls per second = 153.561
fastexp million calls per second = 143.311
fasterpow2 million calls per second = 215.006
fasterexp million calls per second = 214.44
expf million calls per second = 4.16527

快速近似功率

快速近似对数和快速近似指数的组合产生快速近似幂。
// WARNING: this code has been updated.  Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.

inline float
fastpow (float x,
         float p)
{
  return fastpow2 (p *  fastlog2  (x));
}
时间和准确性
时序测试是通过编译完成的 -O3 -finline功能-fast-math,在运行64位Ubuntu Lucid(所以 gcc 4:4.4.3-1ubuntu1libc 2.11.1-0ubuntu7.6)。我还测量了$(x,p)\ [1/200,5] \ times [1/40,10] $中的平均相对准确度。当我测量相对准确度时,我尝试了每个$ p $和每个$ p $逆根(即$ -1 / p $)。
fastpow relative accuracy (positive p) = 0.000165618
fastpow relative accuracy (inverse root p) = 0.00011997
fastpow million calls per second = 58.8533
powf million calls per second = 8.44577

快速近似逆根

通过反根,我的意思是为$ p \ geq 1 $解决$ y = x ^ {-1 / p} $。这与快速近似功率不一样吗?是的,但是在这种情况下,以下方法相当准确且速度更快。初始近似值基于修改后的 反平方根hack ,\ [
I_y =-\ frac {1} {p} I_x + \ frac {p + 1} {p}(B-\ sigma)L。
\] 使用近似精炼 哈雷法 在$ f(y)= \ log(y)+ \ frac {1} {p} \ log(x)$,\ [
\ begin {aligned}
y_ {n + 1}&= y_n \ frac {1-\ delta} {1 + \ delta} \ approx y_n(1-\ delta)^ 2 \\
\ delta&= \ frac {\ log(2)} {2} \ left(\ log_2(y)+ \ frac {1} {p} \ log_2(x)\ right)
\ end {aligned}
\] $ \ log_2 $是使用上面的快速近似对数计算的。总共变成
// WARNING: this code has been updated.  Do not use this version.
// Instead, see http://code.google.com/p/fastapprox for the latest version.

inline float
fasterinvproot (float x,
                float p)
{
  union { float f; uint32_t i; } v = { x };
  unsigned int R = 0x5F375A84U;   // R = (3/2) (B - sigma) L
  unsigned int pR = ((2 * (p + 1)) / (3 * p)) * R;
  unsigned int sub = v.i / p;
  v.i = pR - sub;
  return v.f;
}

inline float
fastinvproot (float x,
              float p)
{
  float y = fasterinvproot (x, p);
  float log2x =  fastlog2  (x);
  float log2y =  fastlog2  (y);
  float err = 1.0 - 0.34657359f * (log2x / p + log2y);
  return y * err * err;
}
注意: fastinvproot 是改进(更快,更准确)的 fast_inverse_p_one 呈现在我的 以前的博客文章.
时间和准确性
时序测试是通过编译完成的 -O3 -finline功能-fast-math,在运行64位Ubuntu Lucid(所以 gcc 4:4.4.3-1ubuntu1libc 2.11.1-0ubuntu7.6)。我还测量了$(x,p)\ [1/200,5] \ times [1/40,10] $中的平均相对准确度。当我测量相对准确度时,我尝试了每个$ p $和每个$ p $逆根(即$ -1 / p $)。
  fastinvproot   relative accuracy (positive p) = 0.000727901
fastinvproot relative accuracy (inverse root p) = 0.00300208
fastinvproot million calls per second = 82.7019
所以这有点令人困惑,但是这里的``正p''情况意味着要取一个反$ p $的根,而``反根p''的意思实际上是试图让这个例程陷入计算$ y = x ^ p $。因此,这些数字表示(出于我不完全理解的原因) fastinvproot 查找逆根比使用一般幂函数更准确。

买者自负

除了精度较低外,这些例程还不能处理浮点数 特殊价值 像标准例程一样。另外,如果您不适当地向这些例程提供负输入,则会发生奇怪的事情。

2011年6月16日,星期四

快速逆根

最近,我想知道如何快速近似一个数字的$ p $根,\ [
y = x ^ {-\ frac {1} {p}},
\] 使用时出现 $ L_p $归一化.

当$ p = 2 $时 反平方根hack 适用,因此我在Wikipedia上对$ p \ neq 2 $修改了派生结果,导致\ [
I_y =-\ frac {1} {p} I_x + \ frac {p + 1} {p}(B-\ sigma)L,
\] 导致了以下代码,
// WARNING: this code has been updated.  Do not use this version.
// Instead get the latest version from http://code.google.com/p/fastapprox

float
fast_inverse_p (float x,
                float p)
{
  union { float f; uint32_t i; } v = { x };

  unsigned int R = 0x5F375A84U;   // R = (3/2) (B - sigma) L
  unsigned int pR = ((2 * (p + 1)) / (3 * p)) * R;
  unsigned int sub = v.i / p;

  v.i = pR - sub;

  return v.f;
}
这个例程做的很合理,但是缺乏 牛顿法 反平方根hack具有的抛光步骤选项。由于要避免调用powf,因此我将牛顿方法应用于$ f(y)= \ log(y)+ \ frac {1} {p} \ log(x)$,这样会产生迭代\ [
y_ {n + 1} = y_n \ left(1-\ log(2)\ left(\ frac {1} {p} \ log_2(x)+ \ log_2(y)\ right)\ right)。
\] 的想法是将其与$ \ log_2 $的廉价近似值一起使用,以利用浮点表示形式。我尝试的第一个近似是一个有理近似,得出了以下结论。
// WARNING: this code has been updated.  Do not use this version.
// Instead get the latest version from http://code.google.com/p/fastapprox

float
fast_inverse_p_one (float x,
                    float p)
{
  float y = fast_inverse_p (x, p);
  union { float f; uint32_t i; } vx = { x };
  union { float f; uint32_t i; } vy = { y };

  // equivalent to: mx = frexp (x, &ex), ex += 126
  int ex = ((vx.i >> 23) & 0xFF);
  union { uint32_t i; float f; } mx = { (vx.i & 0x007FFFFF) | (0x7e << 23) };
  int ey = ((vy.i >> 23) & 0xFF);
  union { uint32_t i; float f; } my = { (vy.i & 0x007FFFFF) | (0x7e << 23) };

  float log2x = ex - 123.52964741f - 4.29802897f / (0.73958333f + mx.f);
  float log2y = ey - 123.52964741f - 4.29802897f / (0.73958333f + my.f);

  return y * (1.0f - 0.69314718f * (log2x / p + log2y));
}
那有些精确,有些贵。为了获得更高的精度,我采用了多项式逼近法和两次牛顿迭代法。
// WARNING: this code has been updated.  Do not use this version.
// Instead get the latest version from http://code.google.com/p/fastapprox

float
fast_inverse_p_two (float x,
                    float p)
{
  float y = fast_inverse_p (x, p);
  union { float f; uint32_t i; } vx = { x };
  union { float f; uint32_t i; } vy = { y };

  int ex = ((vx.i >> 23) & 0xFF);
  union { uint32_t i; float f; } mx = { (vx.i & 0x007FFFFF) | (0x7e << 23) };
  int ey = ((vy.i >> 23) & 0xFF);
  union { uint32_t i; float f; } my = { (vy.i & 0x007FFFFF) | (0x7e << 23) };

  float powmx = mx.f;
  float log2x = ex - 129.78533457081f + 10.07763081376f * powmx;
  powmx *= mx.f; log2x -= 13.90982855008f * powmx;
  powmx *= mx.f; log2x += 12.64852864324f * powmx;
  powmx *= mx.f; log2x -= 6.39532347130f * powmx;
  powmx *= mx.f; log2x += 1.364335673877f * powmx;
  log2x /= p;

  float powmy = my.f;
  float log2y = ey - 129.78533457081f + 10.07763081376f * powmy;
  powmy *= my.f; log2y -= 13.90982855008f * powmy;
  powmy *= my.f; log2y += 12.64852864324f * powmy;
  powmy *= my.f; log2y -= 6.39532347130f * powmy;
  powmy *= my.f; log2y += 1.364335673877f * powmy;

  union { float f; uint32_t i; } vz = { y * (1.0f - 0.69314718f * (log2x + log2y)) };
  int ez = ((vz.i >> 23) & 0xFF);
  union { uint32_t i; float f; } mz = { (vz.i & 0x007FFFFF) | (0x7e << 23) };

  float powmz = mz.f;
  float log2z = ez - 129.78533457081f + 10.07763081376f * powmz;
  powmz *= mz.f; log2z -= 13.90982855008f * powmz;
  powmz *= mz.f; log2z += 12.64852864324f * powmz;
  powmz *= mz.f; log2z -= 6.39532347130f * powmz;
  powmz *= mz.f; log2z += 1.364335673877f * powmz;

  return vz.f * (1.0f - 0.69314718f * (log2x + log2z));
}
这更准确,更昂贵。

这是一些示例值以及一些相对详尽的输入范围内的平均相对误差估计。
x       p       fast0   fast1   fast2   powf
7.000   0.870   0.10889 0.10663 0.10681 0.10681
7.000   2.488   0.45789 0.45704 0.45744 0.45744
7.000   4.106   0.63713 0.62284 0.62256 0.62256
7.000   5.724   0.73334 0.71151 0.71181 0.71181
7.000   7.342   0.78715 0.76660 0.76718 0.76718
7.000   8.960   0.82151 0.80411 0.80480 0.80479
69.000  0.870   0.00749 0.00769 0.00770 0.00770
69.000  2.488   0.18674 0.18229 0.18236 0.18236
69.000  4.106   0.36593 0.35650 0.35659 0.35658
69.000  5.724   0.47131 0.47700 0.47726 0.47726
69.000  7.342   0.56050 0.56189 0.56176 0.56176
69.000  8.960   0.63580 0.62382 0.62341 0.62341
211.000 0.870   0.00217 0.00213 0.00213 0.00213
211.000 2.488   0.11642 0.11628 0.11636 0.11636
211.000 4.106   0.27032 0.27150 0.27161 0.27161
211.000 5.724   0.40273 0.39225 0.39260 0.39260
211.000 7.342   0.47678 0.48220 0.48243 0.48243
211.000 8.960   0.54817 0.55025 0.55031 0.55030
err0 = 0.021138
err1 = 0.000680451
err2 = 7.20003e-06
我做了一个时序测试,用 -O3 -finline功能-fast-math,在运行Ubuntu Lucid(所以 gcc 4:4.4.3-1ubuntu1libc 2.11.1-0ubuntu7.6 )。
fast_inverse_p million calls per second = 206.677
fast_inverse_p_one million calls per second = 77.9861
fast_inverse_p_two million calls per second = 26.531
powf million calls per second = 9.04992
巧合的是,每次精度提高都会使吞吐量损失大约三倍。