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。 SSE指令集包含以下指标功能 比较测试,当与 ``浮点数和'' 指令最终计算出无分支三元运算符。

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

无分支兰伯特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} $。

1条评论:

  1. 谢谢,这很有用。我在机器学习上下文中遇到了Lambert W的一种用法,尽管在W(exp(x))会更有用的情况下也是如此。

    就像您暗示的那样,对大x乘以y = exp(x)仅用于计算W(y)可能会产生溢出的风险,此后才有效地再次获取(接近)对数。

    不过,您的初始化方法很容易适应W(exp(x)),加油!

    回复删除