Faster computation of $ke^{-(x - h)^2}$

Question

The question is quite simple; almost every computer language today provides the $\exp(x)$ function in their standard library to compute expressions like $ke^{-(x - h)^2}.$ However, I would like to know whether this function is the fastest way to compute the above expression. In other words, is there some way to compute $ke^{-(x - h)^2}$ faster than $\exp(x)$ in standard libraries while keeping the result very accurate?

I would like to specify that Taylor series will not work for my application, nor will any other polynomial approximations.

In the end, I wish to compute a series of these expressions, $$\sum^n_{i = 0} k_ie^{-(x - h_i)^2}$$ for some integral $n \geq 0.$

Answer 1

There are several reasons why your standard library implementation of exp is slower than you might hope:

• It should to be accurate to 0.5ulp (or 1ulp at worst), no matter what number you give it. This means that special care may need to be taken, say, in the subnormal range.
• It needs to work for all floating point numbers, and return a reasonable answer for all of them. Assuming that you're on an IEEE-754 platform, this includes NaN, Inf, etc. The standard explains what the "correct" answer is for all of those cases, and the library function needs to get it right.
• It may not be specialised for your specific instruction set (i.e. your CPU may have some advanced floating-point instructions that libc doesn't use by default).

This is what it means to be general-purpose.

Chances are that your standard library has had a lot of effort put into it, and you cannot beat it and still be general-purpose. To get more performance, you need to relax some restrictions. So we're going to need some details. How much accuracy is enough? Can you assume anything about the range of inputs?

I've gone through a typical implementation elsewhere (see this previous answer), and that should help you get an idea of how these things can be implemented.

Answer 2

TL;DR: Using the exp function of your library is likely the fastest way to compute $ke^{-(x-h)^{2}}$, but it's usually worthwhile to do some limited experiments to double-check that assumption.

Generally speaking, when programming for a reasonably mature platform, standard math libraries tend to be highly optimized by both domain experts and microarchitecture experts. While not all functions may be optimized equally well, the most commonly used math functions (exp, log, sin, cos, atan{2}) tend to be the most heavily optimized.

I am assuming you have already profiled your code to establish that the multiple calls to exp are a bottleneck in your code, and that you have double-checked your algorithm(s) to minimize the calls to this function. I further assume that you have already established that you cannot perform the computation at lower precision (say, float instead of double), which results in a significant performance increase on most platforms.

Are you using the latest compiler and libraries available for your platform? Performance improvements are incorporated all the time, so recent tool chains with their associated libraries tend to offer the highest performance. Are you targeting the compiler's code generation to the architecture that most closely reflects your processor's architecture? Newer processors tend to add performance enhancing hardware, such as fused multiply-add (FMA) units and wider SIMD operations and compiler often need to be instructed to use them via compiler flags, e.g. -march=core-avx2.

Also, make sure you are maxing out compiler optimizations. Some advanced optimizations may require adding compiler switches by hand as they are not subsumed under -O3. Examples could be auto-vectorization, whole-program optimization (by use of an optimizing linker), or profile-guided optimizations. Your math library may offer multiple levels of performance / accuracy trade-offs. For example, Intel's MKL provides three modes: high accuracy (maximum error < 1 ulp), lower accuracy (maximum eror < 4 ulp), enhanced performance. The lower the accuracy requirement, the higher the performance.

Note that the overall numerical error in the evaluation of the expression will very likely be dominated by the error in the exp argument magnified through exponentiation. Depending on the magnitude of the argument, a 1 ulp error in the input can well turn into a 1000 ulp error in the output. In light of that, the exp function itself does not have to be extremely accurate.

Standard math library functions need to follow the relevant language specification exactly, which includes overhead for the handling of special cases and detection of errors. Standards may also mandate certain accuracy requirements. If your use case allows the elimination of special case handling and a reduction in accuracy, you could try to roll your own function, like the exemplary C implementation below, which requires hardware support for FMA. It is usually a good idea to use the tool chain specific attributes to force inlining of any custom function to eliminate function call overhead and improve instruction scheduling flexibility.

#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>

double uint64_as_double (uint64_t a)
{
    double r;
    memcpy (&r, &a, sizeof r);
    return r;
}

uint64_t double_as_uint64 (double a)
{
    uint64_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

/* Compute exponential function e**x. Maximum error found in testing: < 0.9 ulp */
double my_exp (double a)
{
    const double ln2_hi = 6.9314718055829871e-01;
    const double ln2_lo = 1.6465949582897082e-12;
    const double l2e = 1.4426950408889634; // log2(e)
    const double cvt = 6755399441055744.0; // 3 * 2**51
    double f, j, p, r;
    uint64_t i;

    // exp(a) = exp2(i) * exp(f); i = rint (a / log(2))
    j = fma (l2e, a, cvt);
    i = double_as_uint64 (j);
    j = j - cvt;
    f = fma (j, -ln2_hi, a);
    f = fma (j, -ln2_lo, f);
    // approximate p = exp(f) on interval [-log(2)/2, +log(2)/2]
    p =            2.5022018235176802e-8;  // 0x1.ade0000000000p-26
    p = fma (p, f, 2.7630903491116071e-7); // 0x1.28af3fcaa8f70p-22
    p = fma (p, f, 2.7557514543681978e-6); // 0x1.71dee62382584p-19
    p = fma (p, f, 2.4801491039342422e-5); // 0x1.a01997c8b03e6p-16
    p = fma (p, f, 1.9841269589067952e-4); // 0x1.a01a01475dae0p-13
    p = fma (p, f, 1.3888888945916467e-3); // 0x1.6c16c1852b7d7p-10
    p = fma (p, f, 8.3333333334557717e-3); // 0x1.11111111224c6p-7
    p = fma (p, f, 4.1666666666519782e-2); // 0x1.55555555502a5p-5
    p = fma (p, f, 1.6666666666666477e-1); // 0x1.5555555555511p-3
    p = fma (p, f, 5.0000000000000122e-1); // 0x1.000000000000bp-1
    p = fma (p, f, 1.0000000000000000e+0); // 0x1.0000000000000p+0
    p = fma (p, f, 1.0000000000000000e+0); // 0x1.0000000000000p+0

    // exp(a) = 2**i * exp(f);
    uint64_t ri = (double_as_uint64 (p) + (i << 52));
    r = uint64_as_double (ri);

    // handle special cases
    double fa = fabs (a);
    if (! (fa < 708.0)) { // |a| >= 708 requires double scaling
        i = (a > 0.0) ? 0ULL : 0x8030000000000000ULL;
        r = uint64_as_double (0x7fe0000000000000ULL + i);
        r = r * uint64_as_double (ri - i - 0x3ff0000000000000ULL);
        if (! (fa < 746.0)) { // |a| >= 746 severe overflow / underflow
            r = (a > 0.0) ? INFINITY : 0.0;
            if (isnan (a)) {
                r = a + a;
            }
        }
    }
    return r;
}

Answer 3

See njuffa’s answer first.

You can improve on it if x is often close to h. Njuffa’s answer finds the result for cases where the result is between sqrt(0.5) and sqrt(2) and scaling by a power of 2. If abs(h-x) is small enough that the exponent is less than sqrt(1/2) then you can remove all the scaling code; if x is even closer to h then you can use a lower degree polynomial. If you evaluate this term gazillion times with the same k, then k can be incorporated in the polynomial.

And on typical processors you can evaluate a polynomial if high degree with less latency then using the Horner scheme by evaluating even and odd powers in parallel.