Сколько FLOP нужно tanh?

Примечание. Этот ответ не относится к Python, но я не думаю, что что-то вроде tanh принципиально отличается от языка.

Tanh обычно реализуется путем определения верхней и нижней границ, для которых возвращаются 1 и -1 соответственно. Промежуточная часть аппроксимируется различными функциями следующим образом:

 Interval 0  x_small               x_medium               x_large 
  tanh(x) |  x  |  polynomial approx.  |  1-(2/(1+exp(2x)))  |  1

Существуют полиномы с точностью до одинарной точности с плавающей запятой, а также с двойной точностью. Этот алгоритм называется алгоритмом Коди-Уэйта.

Ссылаясь на это описание (дополнительную информацию о математика, например, как определить x_medium), рациональная форма Коди и Уэйта требует четырех умножений, трех сложений и одного деления с одинарной точностью и семи умножений, шести сложений и одного деления с двойной точностью.

Для отрицательного x вы можете вычислить |x| и перевернуть знак. Поэтому вам нужны сравнения, для которых находится интервал x, и оценивайте соответствующее приближение. Это в общей сложности:

1. Взяв абсолютное значение x
2. 3 сравнения за интервал
3. В зависимости от интервала и точности с плавающей запятой, от 0 до нескольких FLOPS для экспоненты, проверьте этот вопрос о том, как вычислить экспоненту.
4. Одно сравнение, чтобы решить, переворачивать ли знак.

Это отчет 1993 года, но я не думаю, что здесь что-то изменилось.

Если мы посмотрим на 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)). Однако нужно быть точным, потому что 1-exp(t) очень проблематично для малых t-значений из-за потери значимости, поэтому используется expm(x)=exp(x)-1 и вычисляется tanh(x)=-expm1(-2x)/(expm1(-2x)+2).

Таким образом, наихудший случай примерно в 2 раза превышает количество операций, необходимых для expm1, что является довольно сложной функцией. Вероятно, лучший способ — просто измерить время, необходимое для вычисления tanh(x), по сравнению со временем, необходимым для простого умножения двух чисел типа double.

Мои (небрежные) эксперименты на Intel-процессоре дали следующий результат, дающий примерное представление:

Таким образом, для очень маленьких и чисел> 22 затрат почти нет, для чисел до 0.1 мы платим 6 FLOPS, затем затраты возрастают примерно до 20 FLOPS за tanh-расчет.

Главный вывод: затраты на вычисление tanh(x) зависят от параметра x, а максимальные затраты находятся где-то между 10 и 100 FLOP.

Существует инструкция Intel под названием F2XM1, которая вычисляет 2^x-1 для -1.0<x<1.0, которую можно использовать для вычисление tanh, по крайней мере, для некоторого диапазона. Однако, если верить таблицам Агнера, затраты на эту операцию составляют около 60 FLOP.

Другая проблема заключается в векторизации — обычная glibc-реализация, насколько я понимаю, не векторизована. Поэтому, если ваша программа использует векторизацию и должна использовать невекторизованную реализацию tanh, это еще больше замедлит работу программы. Для этого в компиляторе Intel есть mkl-библиотека, которая векторизирует tanh среди других.

Как видно из таблиц, максимальные затраты составляют порядка 10 тактов на операцию (затраты на float-операцию около 1 такта).

Я предполагаю, что есть шанс, что вы можете выиграть несколько FLOP, используя параметр компилятора -ffast-math, что приводит к более быстрой, но менее точной программе (это вариант для 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...
}

Вопрос предполагает, что он задается в контексте машинного обучения, и поэтому основное внимание уделяется вычислениям с одинарной точностью и, в частности, использованию формата IEEE-754 binary32. Аскер заявляет, что графические процессоры NVIDIA представляют интерес. Я сосредоточусь на использовании этих графических процессоров с помощью CUDA, поскольку я не знаком с привязками Python для CUDA.

Говоря о FLOPS, существуют различные школы мысли о том, как их считать помимо простых сложений и умножений. Графические процессоры, например, вычисляют деления и квадратные корни в программном обеспечении. Менее двусмысленно идентифицировать инструкции с плавающей запятой и подсчитывать их, что я и сделаю здесь. Обратите внимание, что не все инструкции с плавающей запятой будут выполняться с одинаковой пропускной способностью, и это также может зависеть от архитектуры графического процессора. Некоторую важную информацию о пропускной способности инструкций можно найти в Руководстве по программированию CUDA.

Начиная с архитектуры Тьюринга (вычислительная способность 7.5), графические процессоры включают инструкцию MUFU.TANH для вычисления гиперболического тангенса одинарной точности с точностью около 16 бит. Функции одинарной точности, поддерживаемые многофункциональным блоком (MUFU), обычно вычисляются посредством квадратичной интерполяции в таблицах, хранящихся в ПЗУ. Насколько я могу судить, 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 мы можем реализовать встроенную функцию с очень близкими характеристиками с помощью алгебраического преобразования и использования машинных инструкций 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.EX2 и MUFU.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, показывает, что ветвь сохраняется. Это означает, что, как правило, во всех активных потоках некоторые из них будут следовать одной стороне ветви, а остальные — другой стороне ветви, что потребует последующей синхронизации. Таким образом, чтобы получить реалистичное представление о необходимых вычислительных затратах, нам нужно объединить количество инструкций с плавающей запятой для обоих путей выполнения. Получается тринадцать инструкций с плавающей запятой.

Границы ошибок в приведенных выше комментариях к коду были установлены исчерпывающими тестами со всеми возможными аргументами одинарной точности.