Fast Approximate Logarithms, Part II: Rounding Error

In Part I, I assumed all computations were done exactly. It’s time to bring in floating-point rounding error. There’s a key fact about rounding error that might seem surprising: If x \approx 1, then the computation of x-1 in floating-point has no rounding error. The surprise comes for those who’ve been taught to beware of rounding error when subtracting two similar quantities. For example, the computation in floating-point of x^2-1 might be far from the true value of x^2-1 when x \approx 1. That’s because x^2 can have rounding error, so the subtraction might cancel most of the good bits leaving mainly the incorrect (due to rounding) bits. But in x-1, neither x nor 1 have rounding error, so the computation of x-1 is exactly equal to x-1 when x \approx 1.

Here’s the explanation. Rounding error in x+y can happen in two ways. First if x > y and x has a different exponent from y, then x will have its fractional part shifted right to make the exponents match, and so x might drop some bits. Second, even if the exponents are the same, there may be rounding error if the addition of the fractional parts has a carry-out from the high order bit. In the case of x − 1, the exponents are the same if 1 ≤ x < 2. And if 1/2 ≤ x < 1, the larger number is 1 and it does not drop bits when shifted. So there is no rounding error in x − 1 if 1/2 < x < 2.

 
The rule of thumb is that when approximating a function f with f(a) = 0, severe rounding error can be reduced if the approximation is in terms of $x-a$. For us, \log 1 = 0 so a=1, and the rule suggests polynomials written in terms of x-1 rather than x. By the key fact above, there is no rounding error in x-1 when x \approx 1.

Let me apply that to two different forms of the quadratic polynomials used in Part I: the polynomial can be written in terms of x or x-1.

    \begin{eqnarray*} \log_2x & \approx & s_2 x^2 + s_1x + s_0 \\ \log_2x & \approx & \widehat{s}_2 (x-1)^2 + \widehat{s}_1 (x-1) + \widehat{s}_0 \end{eqnarray*}

If they are to be used on the interval [0.75, 1.5) and I want to minimize relative error, it is crucial that the polynomial be 0 when x=1, so they become

    \begin{eqnarray*} \log_2x & \approx & s_2 x^2 + s_1x + (-s_2 - s_1) \\ \log_2x & \approx & \widehat{s}_2 (x-1)^2 + \widehat{s}_1 (x-1) \end{eqnarray*}

The second equation has no constant term, so they both cost the same amount to evaluate, in that they involve the same number of additions and multiplications.

But one is much more accurate. You can see that empirically using an evaluation program (code below) that I will be using throughout to compare different approximations. I invoked the program as \mathit{eval\: 10\: 40960\: 1024} and \mathit{eval\: 10\: 10\: 4194304} and got the following:

SPACING IS 1/1024
using x    bits  5.5 at x=0.750000  2.06 nsecs nsec/bit=0.372  bits/nsec=2.69
using x-1  bits  5.5 at x=0.750000  2.10 nsecs nsec/bit=0.380  bits/nsec=2.63

SPACING IS 1/4194304 = 2^-22
using x    bits  1.7 at x=1-2.4e-07 2.08 nsecs nsec/bit=1.222  bits/nsec=0.82
using x-1  bits  5.5 at x=0.750000  2.12 nsecs nsec/bit=0.384  bits/nsec=2.61

When the approximating polynomials are evaluated at points spaced 1/1024 apart, they have similar performance. The accuracy of both is 5.5 bits, and the one using x-1 is slightly slower. But when they are evaluated at points spaced 2^{-22} apart, the polynomial using x has poor accuracy when x is slightly below 1. Specifically, the accuracy is only 1.7 bits when x \approx 1 - 2.4 \times 10^{-7}.

To see why, note that when x \approx 1, \widehat{s}_2 (x-1)^2 + \widehat{s}_1 (x-1) is summing two numbers that have rounding error, but are of different sizes, since (x-1)^2 \ll x-1. But s_2 x^2 + s_1x + s_0 is summing two numbers of similar size, since $s_0 = -s_1 - s_2$ and the sum of the first two terms is about s_1 + s_2. This is the bad case of subtracting two nearby numbers (cancellation), because they both have rounding error.

I suppose it is an arguable point whether full accuracy for all x is worth a time performance hit of about 2%. I will offer this argument: you can reason about your program if you know it has (in this case) 5.5 bits of accuracy on every input. You don’t want to spend a lot of time tracking down unexpectedly low accuracy in your code that came about because you used a log library function with poor precision on a small set of inputs.

Here’s some more information on the output of the evalution program displayed above. The first number is accuracy in bits measured in the usual way as -\log_2 E where E is the maximum relative error. Following \mathit{x=} is the value of x where the max error E occured. The execution time (e.g. 2.06 nsecs for the first line) is an estimate of the time it takes to do a single approximation to \log_2, including reducing the argument x to the interval [0.75, 1.5). The last two numbers are self explanatory.

Estimating execution time is tricky. For example on my MacBook, if the argument x must be brought into the cache, it will significantly affect the timings. That’s why the evaluation program brings \mathit{xarr} and \mathit{yarr} into the cache before beginning the timing runs.

For polynomials, using x-1 has almost the same cost and better accuracy, so there is a good argument that it is superior to using x. Things are not so clear when the approximation is a rational function rather than a polynomial. For example, q(x) =\frac{a(x-1)+b}{c(x-1)+d}. Because q(1) = 0, the numerator is actually a(x-1). And because you can multiply numerator and denominator by anything (as long as it’s the same anything), it further simplifies to q(x) = \frac{a(x-1)}{(x-1)+e}. This will have no floating-point cancellation, and will have good accuracy even when x \approx 1. But there’s a rewrite of this expression that is faster:

    \[ \frac{a(x-1)}{(x-1)+e} = a - \frac{ae}{(x-1) + e} \]

Unfortunately, this brings back cancellation, because when x \approx 1 there will be cancellation between a and the fraction. Because there’s cancellation anyway, you might as well make a further performance improvement eliminating the need to compute x-1, namely

    \[ \frac{a(x-1)}{(x-1)+e} = a - \frac{ae}{x + (e-1)} \]

Both sides have a division. In addition, the left hand side has a multiplication and 2 additions. The right hand side has no multiplications and 2 additions (ae is a constant and doesn’t involve a run-time multiplication, similarly for e-1). So there is one less multiplication, which should be faster. But at the cost of a rounding error problem when x \approx 1.

SPACING IS 1/1024
using x-1  bits  7.5 at x=0.750000  2.17 nsecs nsec/bit=0.289  bits/nsec=3.46
using x    bits  7.5 at x=0.750000  2.07 nsecs nsec/bit=0.275  bits/nsec=3.64

SPACING IS 1/4194304 = 2^-22
using x-1  bits  7.5 at x=0.750000  2.24 nsecs nsec/bit=0.298  bits/nsec=3.36
using x    bits  1.4 at x=1+2.4e-07 2.09 nsecs nsec/bit=1.522  bits/nsec=0.66

As expected, the rational function that has one less multiplication (the line marked using x) is faster, but has poor accuracy when x is near 1. There’s a simple idea for a fix. When |x-1| is small, use the Taylor series, \log_2 x \approx \log_2 e \left( (x-1) - \frac{1}{2} (x-1)^2 + \cdots\right). Using (\log_2 e) (x-1) is a subtraction and a multiplication, which is most likely cheaper than a division and two additions, What is the size cutoff? The error in the Taylor series is easy to compute: it is the next term in the series, (x-1)^2/2, so the relative error is about (x-1)/2. And I want to maintain an accuracy of 7.5 bits, or 2^{-7.5} \approx 0.0055. So the cutoff is 1 \pm \delta, where \delta/2 = 0.0055 or \delta = 0.011.

On my MacBook, the most efficient way to implement |x-1| < \delta appears to be x < 1 + \delta \small{\mathsf{\;AND\;}} x > 1 - \delta. In the evaluation program, most of the x are greater than 1, so only the first of the inequalities is executed. Despite this, adding the check still has a high cost, but no more accuracy than using x-1.

SPACING IS 1/1024
using x-1  bits  7.5 at x=0.750000  2.17 nsecs nsec/bit=0.289  bits/nsec=3.46
using x    bits  7.5 at x=0.750000  2.07 nsecs nsec/bit=0.275  bits/nsec=3.64
cutoff     bits  7.5 at x=0.750000  2.58 nsecs nsec/bit=0.343  bits/nsec=2.91

SPACING IS 1/4194304 = 2^-22
using x-1  bits  7.5 at x=0.750000  2.24 nsecs nsec/bit=0.298  bits/nsec=3.36
using x    bits  1.4 at x=1+2.4e-07 2.09 nsecs nsec/bit=1.522  bits/nsec=0.66
cutoff     bits  7.5 at x=0.989000  2.60 nsecs nsec/bit=0.347  bits/nsec=2.88

In Part I of this series, I noted that testing whether x was in the range 0.75 \leq x < 1.5 can be done with a bit operation rather than a floating-point one. The same idea could be used here. Instead of using the Taylor series when |x-1| < 0.011 or 0.989 \leq x \leq 1.011, use it in a slightly smaller interval

    \begin{alignat*}{3} & 0.992 && < x && < 1.007812 \\ & 1 - 2^{-7} && < x && < 1 + 2^{-7} \end{alignat*}

The latter can be converted to bit operations on $f$, the fraction part of x, as follows:

    \begin{eqnarray*} 1 \leq x < 1 + 2^{-7} & \small{\mathsf{\;\; OR \;\;}} & 1 - 2^{-7} < x < 1 \\ 1 \leq x < 2 \;\small{\mathsf{\; AND \;}} f < 2^{-7} & \small{\mathsf{\;\; OR \;\;}} & \frac{1}{2} \leq x < 1 \;\small{\mathsf{\; AND \;}} 1 - 2^{-6} < 2x-1 = f\\ \mathit{\; exp \;} = 0 \;\small{\mathsf{\; AND \;}} f < 2^{-7} & \small{\mathsf{\;\; OR \;\;}} & \mathit{\; exp \;} = -1 \;\small{\mathsf{\; AND \;}} f > 1 - 2^{-6} \end{eqnarray*}

As bit operations, this is

(exp == 0 && (f & 111111100...) == 0) 
   OR (exp = -1 && (f & 11111100...) == 111111000...)

When I tested this improvement (\mathit{fcutoff} in the table below). it was faster, but still slower than using x-1, at least on my MacBook.

SPACING IS 1/1024
using x-1  bits  7.5 at x=0.750000  2.17 nsecs nsec/bit=0.289  bits/nsec=3.46
using x    bits  7.5 at x=0.750000  2.07 nsecs nsec/bit=0.275  bits/nsec=3.64
cutoff     bits  7.5 at x=0.750000  2.58 nsecs nsec/bit=0.343  bits/nsec=2.91
fcutoff    bits  7.5 at x=0.750000  2.46 nsecs nsec/bit=0.327  bits/nsec=3.06

SPACING IS 1/4194304 = 2^-22
using x-1  bits  7.5 at x=0.750000  2.24 nsecs nsec/bit=0.298  bits/nsec=3.36
using x    bits  1.4 at x=1+2.4e-07 2.09 nsecs nsec/bit=1.522  bits/nsec=0.66
cutoff     bits  7.5 at x=0.989000  2.60 nsecs nsec/bit=0.347  bits/nsec=2.88
fcutoff    bits  7.5 at x=0.750001  2.44 nsecs nsec/bit=0.325  bits/nsec=3.08

Bottom line: having special case code when x \approx 1 appears to significantly underperform computing in terms of x-1.

In the first post, I recommended reducing to [0.75, 1.5) instead of [1, 2) because you get one extra degree of freedom, which in turns gives greater accuracy. Rounding error gives another reason for preferring [0.75, 1.5). When x = 1 - \delta, reduction to [1, 2) will have cancellation problems. Recall the g function that was optimal for the interval [1, 2), g(x) = -(1/3)(x-1)^2 + (4/3)(x-1). When x = 1 - \delta, x must be multiplied by two to move into [1, 2), and then to compensate, the result is g(2x) - 1. When x \approx 1, g(2x) \approx 1, and so you get cancellation. Below are the results of running the evaluation program on g. If there was no rounding error, g would be accurate to 3.7 bits. As you get closer to 1 (\delta \rightarrow 0) the accuracy drops.

SPACING IS 1/2^19
g    bits  3.5 at x=1-3.8e-06  2.05 nsecs nsec/bit=0.592  bits/nsec=1.69

SPACING IS 1/2^20
g    bits  2.9 at x=1-9.5e-07  2.05 nsecs nsec/bit=0.706  bits/nsec=1.42

SPACING IS 1/2^21
g    bits  2.9 at x=1-9.5e-07  2.05 nsecs nsec/bit=0.706  bits/nsec=1.42

SPACING IS 1/2^22
g    bits  1.7 at x=1-2.4e-07  2.05 nsecs nsec/bit=1.203  bits/nsec=0.83

The goal of this series of posts is to show that you can create logarithm routines that are much faster than the library versions and have a minimum guaranteed accuracy for all x. To do this requires paying attention to rounding error. Summarizing what I’ve said so far, my method for minimizing rounding error problems is to reduce x to the interval [0.75, 1.5) and write the approximating expression using x-1, for example a(x-1)^2 + b(x-1)). More generally, the approximating expression would be a polynomial

    \[ a_n(x-1)^n + a_{n-1}(x-1)^{n-1} + \cdots + a_1(x-1) \]

or rational function

    \[ \frac{ b_n(x-1)^n + b_{n-1}(x-1)^{n-1} + \cdots + b_1(x-1)}{ (x-1)^m + c_{m-1}(x-1)^{m-1} + \cdots + c_1(x-1) + c_0} \]

I close by giving the code for the evaluation program that was used to compare the time and accuracy of the different approximations:

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <math.h> 

/*
 *   Usage: eval [hi reps spacing]
 *       Evaluates an approximation to log2 in the interval [0.125, hi]
 *       For timing purposes, repeats the evaluation reps times.
 *       The evaluation is done on points spaced 1/spacing apart.	       
 */

int
main(argc, argv)
    char **argv;
{
    float x;
    struct timeval start, stop;
    float lo, hi, delta;
    int i, j, n, repetitions, one_over_delta;
    double xd;
    float *xarr, *lg2arr, *yarr;

    // parameters
    lo = 0.125;
    hi = 10.0;
    one_over_delta = 4194304.0; // 2^22
    repetitions = 1;
    if (argc > 1) {
        hi = atof(argv[1]);
        repetitions = atoi(argv[2]);
        one_over_delta = atoi(argv[3]);
    }
    delta = 1.0/one_over_delta;

    // setup
    n = ceil((hi - lo)/delta) + 1;
    xarr = (float *)malloc(n*sizeof(float));
    yarr = (float *)malloc(n*sizeof(float));
    lg2arr = (float *)malloc(n*sizeof(float));
    i = 0;
    for (xd = lo; xd <= hi; xd += delta) {
        x = xd;
        if (x == 1.0) // relative error would be infinity
            continue;
        xarr[i] = x;
        lg2arr[i++] = log2(x);
    }
    if (i >= n)  // assert (i < n)
        fprintf(stderr, "Help!!!\n");
    n = i;

    /* cache-in xarr[i], yarr[i] */
    yarr[0] = 0.0;
    for (i = 1; i < n; i++) {
        yarr[i] = xarr[i] + yarr[i-1];
    }
    fprintf(stderr, "cache-in: %f\n\n", yarr[n-1]); // to foil optimizer

    gettimeofday(&start, 0);
    for (j = 0; j < repetitions; j++) {
        for (i = 0; i < n; i++) {
            yarr[i] = approx_fn(xarr[i]);
        }
    }
    gettimeofday(&stop, 0);
    finish(&start, &stop, "name ", n, repetitions, xarr, yarr, lg2arr);

    exit(0);
}

// convert x to string, with special attention when x is near 1
char *format(float x) {
    static char buf[64];
    float y;

    if (fabs(x - 1) > 0.0001)
        sprintf(buf, "%f", x);
    else {
        y = x-1;
        if (y < 0)
            sprintf(buf, "1%.1e", y);
        else
            sprintf(buf, "1+%.1e", y);
    }
    return(buf);
}

void
finish(struct timeval *start, struct timeval *stop,
       char *str, int n, int repetitions, float *xarr, float *yarr, float *lg2arr)
{
    double elapsed; // nanosecs
    float max, rel;
    int maxi, i;
    double bits;

    elapsed = 1e9*(stop->tv_sec - start->tv_sec) + 1000.0*(stop->tv_usec - start->tv_usec);
    max = 0.0;
    for (i = 0; i < n; i++ ) {
        rel = fabs( (yarr[i] - lg2arr[i])/lg2arr[i]);
        if (rel > max) {
            max = rel;
            maxi = i;
        }
    }
    bits = -log2(max);
    elapsed = elapsed/(n*repetitions);
    printf("%s bits %4.1f at x=%s  %.2f nsecs nsec/bit=%.3f  bits/nsec=%.2f\n",
           str, bits, format(xarr[maxi]), elapsed, elapsed/bits, bits/elapsed);
}


Powered by QuickLaTeX

Mobile First – A Retrospective

Responsive_Web_DesignMost engineers would agree that simply having a mobile-first mindset is not enough to build a high-quality responsive website — we also need relevant mobile-first guidelines and principles. We recently explored the challenges of scaling a mobile-oriented site to tablet and desktop. What we found was not always pretty, but the lessons learned were valuable. In this article, we will explore those lessons, in hopes that we can all improve how we build across mobile, tablet, and desktop.

Building a fully scalable website requires a strong focus on code quality. Concepts such as modularity, encapsulation, and testability become extremely important as you move across domains. Whether we are scaling up to desktop or down to mobile, we need the code to stay consistent and maintainable. Every hacked, poorly planned, or rushed piece of code we might add reduces our ability to write elegant, scalable, responsive code.

Perhaps creating a responsive app is not high on your team’s priority list right now. But one day it will be — and the conversion time frame might be very tight when that day comes.

Ideally, all you need to do is add media query CSS and everything just works. But the only way that can happen is if the code readily adapts to responsive changes.

Below are some suggestions and fixes that will make conversion to responsive easier. Some are specific to responsive design while others are general good practices.

Media queries

Yes, we all know about media queries. How hard can they be? Sprinkle some on any page and you have a responsive website, right?

Using media queries on your pages is essential; they allow you to overwrite CSS values based on screen size. This technique might sound simple, but in a larger project it can quickly get out of hand. A few major problems can get in the way of using media queries properly:

  • Colliding media queries: It is easy to make the mistake of writing media queries that overwrite each other if you do not stick to a common pattern. We recommend using the same boilerplate throughout all projects, and have created one here.
  • Setting element styles from JS: This is a tempting, but inferior, approach to building responsive websites. When an element relies on JS logic to set its width, it is unable to properly use media queries. If the JS logic is setting width as an inline property, the width cannot be overwritten in CSS without using !important. In addition, you have to now maintain an ever-growing set of JS logic.
  • Media queries not at the bottom: If your queries are not loaded last, they will not override their intended targets. Every module might have its own CSS file, and the overall ordering might not place it at the bottom, which leads us to our next point.
  • CSS namespacing for encapsulation: If you are writing a module, its CSS selectors should be properly encapsulated via namespace. We recommend prefixing class names with the module name, such as navbar-parent. Following this pattern will prevent conflicts with other modules, and will ensure that media queries at the bottom of your module’s CSS file override their intended targets.
  • Too many CSS selectors: CSS specificity rules require media queries to use the same specificity in order to override. It is easy to get carried away in LESS, which allows you to nest CSS multiple levels deep. While it can be useful to go one or two levels deep for encapsulation, usually this is unnecessarily complicating your code. We recommend favoring namespacing over nested specifiers as it is cleaner and easier to maintain.
  • Using !important to override styles: Adding !important to your styles reduces maintainability. It is better to avoid relying on !important overrides and instead use CSS namespacing to prevent sharing between modules.

Responsive or adaptive?

Both responsive and adaptive web design techniques contain powerful tools, but it is important to understand the differences between the two. Responsive techniques usually include media queries, fluid grids, and CSS percentage values. Adaptive techniques, on the other hand, are focused more on JavaScript logic, and the adding or removing of features based on device detection or screen size.

So, which should you use? Responsive or adaptive? The answer depends on the feature you are trying to implement. It can be tempting to jump straight into applying adaptive techniques to your feature, but in many cases it may not be required. Worse, applying adaptive techniques can quickly over-complicate your design. An example of this that we saw in many places is the use of JavaScript logic to set CSS style attributes.

Use JavaScript sparingly

When styling your UI, JavaScript should be avoided whenever possible. Dynamic sizing, for example, is better done through media queries. For most UI designs, you will be deciding on layouts based on screen size, not on device type. Confusing the need for device detection with screen size can lead us to apply adaptive where responsive would be superior.

Rethink any design that requires CSS attributes to change based on device detection; in almost all cases it will be better to rely on screen size alone, via media queries. So, when should we use adaptive Javascript techniques?

When to use adaptive

Adaptive web design techniques are powerful, as they allow for selective loading of resources based on user agent or screen size. Logic that checks for desktop browsers, for example, can load high-resolution images instead of their mobile-optimized counterparts. Loading additional resources and features for larger screens can also be useful. Desktop browsers, for example, could show more functionality due to the increased screen size, browser capability, or bandwidth.

Ideally, additional resources will be lazy-loaded for their intended platforms. Lazily loading modules helps with site speed for mobile web, while still allowing for a full set of functionality for desktop and tablet web. This technique can be applied by checking the user agent on the client or server. If done on the server, only resources supported by the user’s platform should be returned. Alternatively, client-based lazy loading can use Ajax requests to load additional resources if they are supported. This effect can be achieved using client-side JavaScript, based on browser support or user agent. Client-side detection is generally preferred, as it allows feature detection based on actual browser functionality instead of potentially complicated user agent checks.

Simple flex grid example

A responsive flex grid doesn’t have to be complicated. In our live demo page, we show a simple implementation that creates a horizontally scrolling section of image containers. The images are centered, allowed to expand up to 100% of their container, and will maintain their original aspect ratio. In addition, the container height values are set to 100%, allowing us to adjust the height in the parent wrapper only, and keeping our media query overrides simple and easy to read.

The html and css source code use the concepts mentioned above. We plan to add more boilerplate patterns; please don’t hesitate to add your own as well. Pull requests are welcomed!

Best practices

We hope that the information above will come in handy when you are working on your next mobile-first web project. Below is a summary of what we mentioned above and other helpful tips.

On using media queries

  • Most responsive layout can and should be done with media queries. JS manipulation of CSS (maybe with the exception of adding/removing classes) should be avoided. Setting width in JS is not as maintainable or dynamic compared to CSS.
  • Use media query boilerplate to ensure you do not have contradicting media queries or have media queries that are always skipped.
  • Put media queries at the bottom. Media queries override CSS and should be the final overrides, whether page level or module level.
  • If your regular CSS rules have many selectors, your media query CSS rules will have to as well, due to CSS specificity rules. Use as few selectors as possible when defining CSS rules.

 On writing CSS

  • Use CSS classes, not CSS IDs, to avoid CSS specificity issues.
  • Use the fewest number of selectors possible to define your selector.
  • Reuse classes. If an element has the same look on different parts of the page, do not create two different classes. Make a generic class and reuse it.
  • Encapsulate your CSS selectors by using proper namespacing to prevent conflicts.

e.g., class=”module-name-parent”

  • It is very rare that you need to use !important. Before you use it, ask yourself whether you can instead add another class (parent or same level). And then ask yourself whether the rule you are trying to override has unnecessary selectors.

 On writing LESS

  • Use LESS nesting only where needed. Nesting is good for organization, but it is also a recipe for CSS specificity issues.
  • Check that you do not have a CSS rule that looks like this:
    #wrapper #body-content #content #left-side #text {
        border: 1px solid #000;
    }
  • Work with the design team and define LESS variables using good names. Then, use these LESS variables everywhere possible.
  • If you are using a set of CSS rules repeatedly, make it a LESS mixin.

 On adding element wrappers

  • Most dom structures are more complex than necessary.
  • Add a wrapper only when needed. Do not add a wrapper when proper CSS can do the same thing.
  • If you remove the wrapper and the layout does not change, you do not need it. Now, do a global search for this wrapper’s references (JS, CSS, rhtml, jsp, tag) and remove them.

On lazy loading

  • Add a placeholder to your component for lazy loading.
  • Lazy-loaded sections will start off empty, so make sure you reserve the correct amount of space for this behavior. Otherwise, you will see the page shift as modules load in.
  • Use media queries for the empty section so that it closely matches the filled size.

On cleaning up

  • If you are playing around with CSS to attempt a layout and it starts working, remember to remove the unnecessary CSS rules. Many of them are probably not needed anymore. Remove the unnecessary wrappers as well.

Image source:  http://upload.wikimedia.org/wikipedia/commons/e/e2/Responsive_Web_Design.png

Fast Approximate Logarithms, Part I: The Basics

Performance profiling of some of the eBay code base showed the logarithm (log) function to be consuming more CPU than expected. The implementation of log in modern math libraries is an ingenious wonder, efficiently computing a value that is correct down to the last bit. But for many applications, you’d be perfectly happy with an approximate logarithm, accurate to (say) 10 bits, especially if it were two or three times faster than the math library version. This post is the first of a series that examines approximate logarithms and the tradeoff between accuracy and speed. By the end of the series, you will see how to design fast log functions, and will discover that on a MacBook Pro, you can get a roughly 3x speedup if you are willing to settle for 7 bits. Alternatively, if you want more accuracy you can get 11 bits with a 2x speedup. The numbers might vary for your target machine, but I will present code so that you can determine the tradeoffs for your case.

You can find code for approximate logs on the web, but they rarely come with an evaluation of how they compare to the alternatives, or in what sense they might be optimal. That is the gap I’m trying to fill here. The first post in this series covers the basics, but even if you are familiar with this subject I think you will find some interesting nuggets. The second post considers rounding error, and the final post gives the code for a family of fast log functions.

A very common way to to compute log (meaning \ln = \log_e) is by using the formula \log x = (\log 2) \log_2 x \approx 0.6931472 \log_2x to reduce the problem to computing \log_2. The reason is that \log_2 for arbitrary x is easily reduced to the computation of \log_2 x for x in the interval [1, 2); details below. So for the rest of this series I will exclusively focus on computing \log_2. The red curve in the plot below shows \log_2 on [1, 2]. For comparison, I also plot the straight line y=x-1.

fig1_part1

If you’ve taken a calculus course, you know that \log x has a Taylor series about x=1 which is \log x = (x-1) - \frac{1}{2} (x-1)^2 + \cdots. Combining with \log_2 x = (\log_2 e)\log x gives (t for Taylor)

    \[ t(x) \approx -0.7213 x^2 + 2.8854 x - 2.1640 \]

How well does t(x) approximate \log_2 x?

fig2_part_1

The plot shows that the approximation is very good when x \approx 1, but is lousy for x near 2—so is a flop over the whole interval from 1 to 2. But there is a function that does very well over the whole interval: b(x) = -0.344845x^2 + 2.024658x -1.674873. I call it b for better. It is shown below in red (\log_2 in blue). The plot makes it look like a very good approximation.

fig3_part_1

A better way to see quality of the approximation is to plot the error \log_2 x - b(x). The largest errors are around x \approx 1.2 and x \approx 1.7.

fig4_part1

Now that I’ve shown you an example, let me get to the first main topic: how do you evaluate different approximations? The conventional answer is minimax. Minimax is very conservative—it only cares about the worst (max) error. It judges the approximation over the entire range 1 \leq x < 2 by its error on the worst point. As previously mentioned, in the example above the worst error occurs at x \approx 1.2, or perhaps at x \approx 1.7, since the two have very similar errors. The term minimax means you want to minimize the max error, in other words find the function with the minimum max error. The max error here is very close to 0.0050, and it is the smallest you can get with a quadratic polynomial. In other words, b solves the minimax problem.

Now, onto the first of the nuggets mentioned at the opening. One of the most basic facts about \log is that \log 1 = 0, whether it’s \log_2 or \ln or \log_{10}. This means there’s a big difference between ordinary and relative error when x \approx 1.

As an example, take x = 1.001. The error in b is quite small: \log_2 x - b(x) \approx -0.005. But most likely you care much more about relative error: (\log_2 x - b(x))/\log_2(x), which is huge, about -3.35. It’s relative error that tells you how many bits are correct. If \log_2x and b(x) agree to k bits, then (\log_2 x - b(x))/\log_2(x) is about 2^{-k}. Or putting it another way, if the relative error is \epsilon, then the approxmation is good to about -\log_2 \epsilon bits.

The function b that solved the minimax problem solved it for ordinary error. But it is a lousy choice for relative error. The reason is that its ordinary error is about -0.005 near x=1. As x \rightarrow 1 it follows that \log x \rightarrow 0, and so the relative error will be roughly -0.005/\log x \rightarrow \infty. But no problem. I can compute a minimax polynomial with respect to relative error; I’ll call it r(x) for relative. The following table compares the coefficients of the Taylor method t, minimax for ordinary error b and minimax for relative error r:

Rendered by QuickLaTeX.com

The coefficients of b and r are similar, at least compared to t, but r is a function that is always good to at least 5 bits, as the plot of relative error (below) shows.

fig5_part1

Here’s a justification for my claim that r is good to 5 bits. The max relative error for r occurs at x=1, x=1.5, and x=2. For example, at x=1.5

    \begin{eqnarray*} r(1.5) & = & 0.1001\textbf{1000}\ldots \\ \log_2(1.5) & = & 0.1001\textbf{0101}\ldots \\ \end{eqnarray*}

If you’re a nitpicker, you might question whether this is good to 5 bits as claimed. But if you round each expression to 5 bits, each is 0.10011_2.

 

Unfortunately, there’s a big problem we’ve overlooked. What happens outside the interval [1,2)? Floating-point numbers x are represented as f2^k with 1 \leq f < 2. This leads to the fact mentioned above: \log_2 x = \log_2 f2^k = \log_2 f + k. So you only need to compute \log_2 on the interval [1,2). When you use r(x) for 1 \leq x < 2 and reduce to this range for other x, you get

fig6_part1

The results are awful for x just below 1. After seeing this plot, you can easily figure out the problem. The relative error of r for x \approx 2 is about 0.02, and is almost the same as ordinary error (since the denominator is close to \log_2 2 = 1). Now take an x just below 1. Such an x is multiplied by 2 to move it into [1,2), and the approximation to \log_2 is r(2x) - 1, where the -1 compensates for changing x to 2x. The ordinary error (r(2x) - 1) - \log_2(x) \approx r(2) - 1 is still about 0.02. But \log_2 x is very small for x \approx 1, so the ordinary error of 0.02 is transformed to 0.02/\log_2 x, which is enormous. At the very least, a candidate for small relative error must satisfy r(2) = \log_2(2) = 1. But r(2) \approx 0.98. This can be fixed by finding the polynomial that solves the minimax problem for all x. The result is a polynomial g for global.

Rendered by QuickLaTeX.com

One surprise about g is that its coefficients appear to be simple rational numbers, suggesting there might be a simple proof that this polynomial is optimal. And there is an easy argument that it is locally optimal. Since g(x) = Cx2 + Bx + A must satisfy g(1) = 0 and g(2) = 1 it is of the form gC(x) = C(x-1)2 + (1−C)(x−1). When x > 1 the relative error is ε(x) = (gC(x)− log2(x)) ⁄ log2(x) and limx→ 1+ ε(x) = (1−C)log2 − 1. When x < 1 then ε(x) = (gC(2x) − 1− log2 (x)) ⁄ log2(x) and lim x→1ε(x) = 2(1+C)log2 − 1. The optimal gC has these two limits equal, that is (1−C)log2 − 1 = 2(1+C)log2 − 1, which has the solution C = −1/3.

 

fig7_part1

Globally (over all x), the blue curve g does dramatically better, but of course it comes at a cost. Its relative error is not as good as over the interval [1, 2). That’s because it’s required to satisfy g(2) = 1 in order to have a small relative error at x-2. The extra requirement reduces the degrees of freedom, and so g does less well on [1, 2].

fig8_part1

Finally, I come to the second nugget. The discussion so far suggests rethinking the basic strategy. Why reduce to the interval [1,2)? Any interval [t, 2t) will do. What about using [0.75, 1.5)? It is easy to reduce to this interval (as I show below), and it imposes only a single requirement: that g(1) = 0. This gives g an extra degree of freedom that can be used to do a better job of approximating \log_2. I call the function based on reduction [0.75, 1.5) s for shift, since the interval has been shifted.

Rendered by QuickLaTeX.com

fig9_part1

The result is a thing of beauty! The error of s is significantly less than the error in g. But you might wonder about the cost: isn’t it more expensive to reduce to [0.75, 1.5) instead of [1.0, 2.0)? The answer is that the cost is small. A floating-point number is represented as (1+f)2^e, with f stored in the right-most 23 bits. To reduce to [0.75, 1.5) requires knowing when f \geq \frac{1}{2}, and that is true exactly when the left-most of the 23 bits is one. In other words, it can be done with a simple bit check, not a floating-point operation.

Here is more detail. To reduce to [0.75, 1.5), I first need code to reduce to the interval [1,2). There are library routines for this of course. But since I’m doing this whole project for speed, I want to be sure I have an efficient reduction, so I write my own. That code combined with the further reduction to [0.75, 1.5) is below. Naturally, everything is written in single-precision floating-point. You can see that the extra cost of reducing to [0.75, 1.5) is a bit-wise operation to compute the value \mathsf{greater}, and a test to see if \mathsf{greater} is nonzero. Both are integer operations.

The code does not check that x > 0, much less check for infinities or NaNs. This may be appropriate for a fast version of log.

float fastlog2(float x)  // compute log2(x) by reducing x to [0.75, 1.5)
{
    // a*(x-1)^2 + b*(x-1) approximates log2(x) when 0.75 <= x < 1.5
    const float a =  -.6296735;
    const float b =   1.466967;
    float signif, fexp;
    int exp;
    float lg2;
    union { float f; unsigned int i; } ux1, ux2;
    int greater; // really a boolean 
    /* 
     * Assume IEEE representation, which is sgn(1):exp(8):frac(23)
     * representing (1+frac)*2^(exp-127)  Call 1+frac the significand
     */

    // get exponent
    ux1.f = x;
    exp = (ux1.i & 0x7F800000) >> 23; 
    // actual exponent is exp-127, will subtract 127 later

    greater = ux1.i & 0x00400000;  // true if signif > 1.5
    if (greater) {
        // signif >= 1.5 so need to divide by 2.  Accomplish this by 
        // stuffing exp = 126 which corresponds to an exponent of -1 
        ux2.i = (ux1.i & 0x007FFFFF) | 0x3f000000;
        signif = ux2.f;
        fexp = exp - 126;    // 126 instead of 127 compensates for division by 2
        signif = signif - 1.0;                    // <
        lg2 = fexp + a*signif*signif + b*signif;  // <
    } else {
        // get signif by stuffing exp = 127 which corresponds to an exponent of 0
        ux2.i = (ux1.i & 0x007FFFFF) | 0x3f800000;
        signif = ux2.f;
        fexp = exp - 127;
        signif = signif - 1.0;                    // <<--
        lg2 = fexp + a*signif*signif + b*signif;  // <<--
    }
    // lines marked <<-- are common code, but optimize better 
    //  when duplicated, at least when using gcc
    return(lg2);
}
You might worry that the conditional test if greater will slow things down. The test can be replaced with an array lookup. Instead of doing a bitwise | with 0x3f000000 in one branch and 0x3f800000 in the other, you can have a single branch that uses \mathsf{arr[greater]}. Similarly for the other difference in the branches, exp−126 versus exp−127. This was not faster on my MacBook.

 

In summary:

  • To study an approximation p(x), don’t plot \log_2 and p(x) directly, instead plot their difference.
  • The measure of goodness for p(x) is its maximum error.
  • The best p(x) is the one with the smallest max error (minimax).
  • For a function like \log_2 x, ordinary and relative error are quite different. The proper yardstick for the quality of an approximation to \log_2 x is the number of correct bits, which is relative error.
  • Computing \log requires reducing to an interval [t, 2t) but you don’t need t=1. There are advantages to picking t = 3/4 instead.

In the next post, I’ll examine rounding error, and how that affects good approximations.

Powered by QuickLaTeX