tanh需要多少FLOPs?

12
我希望计算LeNet-5 (文章)每一层需要多少次浮点运算(FLOPs)。有些论文提供了其他架构总体的FLOPs (如1, 2, 3), 然而这些论文没有详细说明如何计算FLOPs,我也不知道非线性激活函数需要多少FLOPs。例如,计算tanh(x)需要多少FLOPs?

我猜这将是特定于实现和硬件的。然而,我主要想知道数量级。我们是在谈论10 FLOPs吗?100 FLOPs?1000 FLOPs?所以请任选一个架构/实现来回答您的问题。(尽管我会欣赏那些接近“常见”设置的答案,比如Intel i5 / Nvidia GPU / Tensorflow)

可能有一个TensorFlow的解决方案:https://github.com/tensorflow/tensorflow/issues/899 - Martin Thoma
在早期的MatLab中,它有一个名为“flops”的函数,可以告诉您它执行了多少个操作。这是非常有用的,可以初步估计算法在C实现中的实时性能。由于MatLab现在使用了大量的外部代码(例如FFTW而不是FFT.m),所以它不再具备这个功能了。 - bazza
3个回答

11
如果我们看一下glibc实现的tanh(x),我们可以看到:
  1. 对于x值大于22.0且双精度的情况下,可以安全地假定tanh(x)为1.0,因此几乎没有任何成本。
  2. 对于非常小的x(比如x < 2 ^(-55)),另一种便宜的近似方法是可能的:tanh(x)=x(1+x),因此只需要两个浮点运算。
  3. 对于中间的值,可以重写tanh(x)=(1-exp(-2x))/(1+exp(-2x))。然而,必须要准确,因为由于显著性损失,对于小的t值,1-exp(t)非常有问题,因此使用expm(x)=exp(x)-1,并计算tanh(x)=-expm1(-2x)/(expm1(-2x)+2)。
基本上,最坏情况下所需操作次数大约是 `expm1` 所需次数的2倍,而这是一个相当复杂的函数。最好的方法可能就是测量计算 `tanh(x)` 所需的时间,与两个双精度数乘法所需的时间进行比较。
我在 Intel 处理器上进行了(草率的)实验,得出了以下结果,可以给出一个大致的想法:

enter image description here

对于非常小的数字和大于22的数字,几乎没有任何成本,对于小于0.1的数字,我们支付6 FLOPS,然后成本上升到每个tanh计算约20 FLOPS。
关键是:计算tanh(x)的成本取决于参数x,最大成本在10至100 FLOPs之间。

有一种英特尔指令叫做F2XM1,用于计算-1.0<x<1.0时的2^x-1,可用于计算tanh,至少在某些范围内。然而,如果Agner's tables是可信的,这个操作的成本约为60 FLOPs。


另一个问题是向量化 - 据我所知,普通的glibc实现没有进行向量化。因此,如果您的程序使用了向量化,并且必须使用未向量化的tanh实现,则会使程序变得更慢。为此,英特尔编译器具有mkl库,其中对其他函数一样,对tanh进行了向量化

从表中可以看出,最大成本约为每个操作10个时钟周期(一个浮点操作的成本约为1个时钟周期)。


我猜你可以通过使用-ffast-math编译器选项赢得一些FLOPS,这会导致程序更快但不太精确(这是Cuda或c / c ++的选项,不确定是否可用于python / numpy)。


这段C++代码用于生成图表数据(使用g++ -std=c++11 -O2编译)。它的目的不是给出精确数字,而是让人们对成本有第一印象。

#include <chrono>
#include <iostream>
#include <vector>
#include <math.h>

int main(){
   const std::vector<double> starts={1e-30, 1e-18, 1e-16, 1e-10, 1e-5, 1e-2, 1e-1, 0.5, 0.7, 0.9, 1.0, 2.0, 10, 20, 23, 100,1e3, 1e4};
   const double FACTOR=1.0+1e-11;
   const size_t ITER=100000000; 


   //warm-up:
   double res=1.0;
      for(size_t i=0;i<4*ITER;i++){
      res*=FACTOR;
   }
   //overhead:
   auto begin = std::chrono::high_resolution_clock::now();
   for(size_t i=0;i<ITER;i++){
      res*=FACTOR;
   }
   auto end = std::chrono::high_resolution_clock::now();
   auto overhead=std::chrono::duration_cast<std::chrono::nanoseconds>(end-begin).count(); 
   //std::cout<<"overhead: "<<overhead<<"\n";


   //experiments:
   for(auto start : starts){
       begin=std::chrono::high_resolution_clock::now();
       for(size_t i=0;i<ITER;i++){
           res*=tanh(start);
           start*=FACTOR;
       }
       auto end = std::chrono::high_resolution_clock::now();
       auto time_needed=std::chrono::duration_cast<std::chrono::nanoseconds>(end-begin).count();
       std::cout<<start<<" "<<time_needed/overhead<<"\n"; 
   }

   //overhead check:
   begin = std::chrono::high_resolution_clock::now();
   for(size_t i=0;i<ITER;i++){
      res*=FACTOR;
   }
   end = std::chrono::high_resolution_clock::now();
   auto overhead_new=std::chrono::duration_cast<std::chrono::nanoseconds>(end-begin).count(); 
   std::cerr<<"overhead check: "<<overhead/overhead_new<<"\n";
   std::cerr<<res;//don't optimize anything out...
}

你是怎么得到这个图形的? - Martin Thoma
1
@MartinThoma,我添加了代码,但是正如我所说的那样,它只能给你一个大概的印象。 - ead

8
注意:本答案不特定于Python,但我认为像tanh这样的东西在各种语言中基本上是相同的。
通常通过定义上限和下限来实现tanh,分别返回1和-1。其中间部分使用不同的函数进行近似,如下所示:
 Interval 0  x_small               x_medium               x_large 
  tanh(x) |  x  |  polynomial approx.  |  1-(2/(1+exp(2x)))  |  1

存在一些多项式可以精确地计算单精度浮点数和双精度浮点数,这个算法被称为Cody-Waite算法。
引用this description(您也可以在那里找到有关数学的更多信息,例如如何确定x_medium),Cody和Waite的有理形式需要4次乘法、3次加法和1次除法来进行单精度运算,并需要7次乘法、6次加法和1次除法来进行双精度运算。
对于负数x,您可以计算|x|并翻转符号。因此,您需要比较x所在的区间,并评估相应的近似值。总共包括:
1. 取x的绝对值; 2. 3次比较区间; 3. 根据区间和浮点精度,0到几个FLOPS用于指数计算,请参见this question以了解如何计算指数; 4. 一个比较来决定是否翻转符号。

现在,这是一份1993年的报告,但我认为这里并没有太多改变。


我认为可能存在一个tanh汇编指令(适用于x86 / nvidia GPU)。例如,我不太确定如何理解这个英特尔页面或者nvidia的配置支持意味着什么。 - Martin Thoma
但是如果有一个:一个汇编指令是否意味着1个FLOP? - Martin Thoma
1
确实存在一个问题,就是如何定义浮点运算。它可以是一条单独的CPU指令,用于处理浮点数据(维基百科是这样定义的),也可以是由基准测试定义的某种操作。例如,融合乘加法最终被放入X64中,对芯片每秒执行多少浮点指令并没有什么影响,但极大地提高了FFT基准测试的性能。顺便说一句,英特尔推迟将FMA纳入x64的AVX,以保持Itanium(一直具有这种指令)的相关性。 - bazza
他刚刚通过Facebook告诉我,x86没有tanh指令,并提供https://en.wikipedia.org/wiki/X86_instruction_listings作为参考来源。 - Martin Thoma

1
该问题表明它是在机器学习的背景下提出的,因此重点放在单精度计算上,特别是使用IEEE-754 binary32格式。提问者表示NVIDIA GPU是一个感兴趣的平台。我将专注于使用CUDA的这些GPU,因为我不熟悉CUDA的Python绑定。
当谈到FLOPS时,有各种思路来计算它们,除了简单的加法和乘法。例如,GPU在软件中计算除法和平方根。更少模棱两可的方法是识别浮点指令并计数,这是我在这里要做的。请注意,并非所有浮点指令都具有相同的吞吐量,这也可能取决于GPU架构。有关指令吞吐量的一些相关信息可以在CUDA编程指南中找到。

从图灵架构(计算能力7.5)开始,GPU包括指令MUFU.TANH,用于计算具有约16位精度的单精度双曲正切。多功能单元(MUFU)支持的单精度函数通常通过在存储在ROM中的表中进行二次插值来计算。据我所知,MUFU.TANH在虚拟汇编语言PTX的级别上公开,但(截至CUDA 11.2)不作为设备函数内部函数。

但是,鉴于该功能在PTX级别上公开,我们可以轻松地使用一行内联汇编代码创建自己的内部函数:

// Compute hyperbolic tangent for >= sm75. maxulperr = 133.95290, maxrelerr = 1.1126e-5
__forceinline__ __device__ float __tanhf (float a)
{
    asm ("tanh.approx.f32 %0,%1; \n\t" : "=f"(a) : "f"(a));
    return a;
}

在计算能力小于7.5的旧GPU架构上,我们可以通过代数变换和使用机器指令MUFU.EX2和MUFU.RCP来实现具有非常相似特征的内在函数,分别计算基于2的指数和倒数。对于数量级较小的参数,我们可以使用tanh(x) = x,并通过实验确定两种逼近之间的良好切换点。
// like copysignf(); when first argument is known to be positive
__forceinline__ __device__ float copysignf_pos (float a, float b)
{
    return __int_as_float (__float_as_int (a) | (__float_as_int (b) & 0x80000000));
}

// Compute hyperbolic tangent for < sm_75. maxulperr = 108.82848, maxrelerr = 9.3450e-6
__forceinline__ __device__ float __tanhf (float a)
{
    const float L2E = 1.442695041f;
    float e, r, s, t, d;
    s = fabsf (a);
    t = -L2E * 2.0f * s;
    asm ("ex2.approx.ftz.f32 %0,%1;\n\t" : "=f"(e) : "f"(t));
    d = e + 1.0f;
    asm ("rcp.approx.ftz.f32 %0,%1;\n\t" : "=f"(r) : "f"(d));
    r = fmaf (e, -r, r);
    if (s < 4.997253418e-3f) r = a;
    if (!isnan (a)) r = copysignf_pos (r, a);
    return r;
}

使用CUDA 11.2编译此代码,针对目标sm_70,然后使用cuobjdump --dump-sass反汇编二进制文件,显示八个浮点指令。我们还可以看到生成的机器代码(SASS)是无分支的。
如果我们需要具有完全单精度精度的双曲正切函数,可以在参数小的情况下使用最小化多项式逼近,同时在参数较大的情况下使用代数变换和机器指令MUFU.EX2MUFU.RCP。超过一定的参数范围,结果将为±1。
// Compute hyperbolic tangent. maxulperr = 1.81484, maxrelerr = 1.9547e-7
__forceinline__ __device__ float my_tanhf (float a)
{
    const float L2E = 1.442695041f;
    float p, s, t, r;
    t = fabsf (a);
    if (t >= 307.0f/512.0f) { // 0.599609375
        r = L2E * 2.0f * t;
        asm ("ex2.approx.ftz.f32 %0,%1;\n\t" : "=f"(r) : "f"(r));
        r = 1.0f + r;
        asm ("rcp.approx.ftz.f32 %0,%1;\n\t" : "=f"(r) : "f"(r));
        r = fmaf (r, -2.0f, 1.0f);
        if (t >= 9.03125f) r = 1.0f;
        r = copysignf_pos (r, a);
    } else {
        s = a * a;
        p =              1.57394409e-2f;  //  0x1.01e000p-6
        p = fmaf (p, s, -5.23025580e-2f); // -0x1.ac766ap-5
        p = fmaf (p, s,  1.33152470e-1f); //  0x1.10b23ep-3
        p = fmaf (p, s, -3.33327681e-1f); // -0x1.5553dap-2
        p = fmaf (p, s, 0.0f);
        r = fmaf (p, a, a);
    }
    return r;
}

这段代码包含一个数据相关的分支。通过 CUDA 11.2 生成的 sm75 目标机器码可以看出该分支被保留。这意味着在一般情况下,所有活动线程中的一些将遵循分支的一侧,而其余线程将遵循分支的另一侧,需要后续同步。因此,为了获得所需的计算工作量的实际概念,我们需要结合两个执行路径的浮点指令计数。这总共有十三条浮点指令。
上面代码注释中的误差界限是通过针对所有可能的单精度参数进行详尽测试而建立的。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接