快速近似兰伯特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-1ubuntu1 和 libc 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} $。