Fast Approximate Logarithms, Part I: The Basics

http://www.ebaytechblog.com/2015/05/01/fast-approximate-logarithms-part-i-the-basics/

by r4um • 11 years ago

70 46 11 years ago

sdenton42• 11 years ago

Awesome writeup.

I've occasionally thought that there's a lot of time to be saved in machine learning systems by making the floating point operations a bit more noisy. Use maybe eight bit (or even four bit) numbers for the matrix entries, with the philosophy that the overall degree of freedom in the system is enough to make up for this local rigidity, while you gain a lot in terms of both storage space and evaluation time.

Along these lines, I saw a great demo of a robot a while ago which used 256 binary devices for movement. Each of these devices is the diagonal of a quadrilateral; when the device is in it's 'on' state, the diagonal is long, and when it's in the 'off' state the diagonal is short. Chain a lot of these together, and you can get pretty good accuracy aiming for specific points in a two-dimensional space. The demo I saw was a three-dimensional version of the same concept, of course.

So likewise, we chain together lots of linear operations to build (say) a deep learning system. The weight space is very high-dimensional, and the decision space is the much lower dimensional classification/regression problem: one probably needs a lot less than 64 bits of accuracy in the weight space to approximate a point in the decision space.

benanne• 11 years ago

Indeed, 64 bits of accuracy is overkill for a lot of ML algorithms, where there is so much noise that the additional quantization noise due to low precision is negligible. Most deep neural nets are already being trained in single precision because that is what NVIDIA GPUs can do the fastest. There has been quite some research on reducing the bit depth further the last couple of years, see e.g. http://arxiv.org/abs/1502.02551 and http://arxiv.org/abs/1412.7024.

The current generation of GPUs already has limited support for half-precision operations, and the tools for using these operations in neural networks (dot product, convolution implementation) are starting to become available as well, which is awesome (see e.g. https://github.com/NervanaSystems/nervanagpu).

NVIDIA themselves have also noticed and better hardware support for low-precision arithmetic is coming in Pascal: http://techreport.com/news/27978/nvidia-pascal-to-feature-mi...

sdenton42• 11 years ago

Awesome, thanks for the pointers!

ot• 11 years ago

I also liked a lot this ICML paper [1], which is more theoretically principled than the two you reference (and weirdly they do not cite it).

[1] http://www.eecs.tufts.edu/~dsculley/papers/round-model-icml....

Houshalter8• 11 years ago

The number of bits does seem to affect neural network training and really kills it when you get really low (like 8 to 12 bits). The problem is that really small numbers are rounded down to zero, and then when summed together the result is really inaccurate.

Ideally you could do some kind of stochastic rounding, where the decimal value represents it's probability of being rounded up or down. On average the totals would still be the same. You could probably get the required number of bits really low then, perhaps even only one or two!

Something to look into when these algorithms are put into FPGAs and ASICs, which is starting to happen.

andor• 11 years ago

The problem is that really small numbers are rounded down to zero, and then when summed together the result is really inaccurate

Use log values instead?

TheLoneWolfling3• 11 years ago

Another interesting idea is ditching the format of a classic polynomial. Multiplication isn't the fastest thing ever.

I wonder what would happen if you ran the same minmax scoring scheme through an general instruction synthesizer for the same complexity cap?

me2you• 11 years ago

As in Math Toolkit for Real-Time Programming by Jack Crenshaw

paulpauper• 11 years ago

But why is this important? Can't computers already do this? Even my 1978 TI pocket calculator can do logs very fast. It's not like this is an important unsolved mathematics puzzle.

pjscott• 11 years ago

From the first sentence:

Performance profiling of some of the eBay code base showed the logarithm (log) function to be consuming more CPU than expected.

There are a lot of different types of "very fast", and these guys needed a faster one.

TheLoneWolfling3• 11 years ago

I wonder why eBay is using the log function?

rg2004• 11 years ago

If your bottleneck is log() and you're not memory-limited, why not use a lookup table?

tgbrter2• 11 years ago

This is faster than a look-up table at the cost of some accuracy. Also no cache misses.

walshemj• 11 years ago

which is how pre calculators you used to do it by hand using Log tables

Houshalter8• 11 years ago

Computing log tables is what the would-be first computer was intended to do.

minimax• 11 years ago

The author of this post, David Goldberg, is the same guy that wrote "What Every Computer Scientist Should Know About Floating-Point Arithmetic."

http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.ht...

bumbledraven• 11 years ago

Great find. There's an archive of his posts on the eBay Tech Blog. Looks like interesting reading.

http://www.ebaytechblog.com/author/dgoldberg/

boulos• 11 years ago

This is a pretty nice introduction, and makes me want to return to my side project where I used to play with these approximations but for SIMD math: https://github.com/boulos/syrah/blob/master/src/include/syra... .

I hope in a later post he mentions Sollya (http://sollya.gforge.inria.fr), I found it invaluable for playing around with different precision vs speed.

Houshalter8• 11 years ago

I tried evolving a good approximation with Eureqa: https://i.imgur.com/1CTTjYD.png?1

dalke4• 11 years ago

It looks like all of them (they use sqrt() and x^y) will be slower than calling log() directly.

Houshalter8• 11 years ago

Here is it without those operations: https://i.imgur.com/xhyNnU4.png?1

EDIT: Without sin/cos: https://i.imgur.com/mTwgKKO.png?

You could do better by splitting the function into different ranges and using a lookup table for a formula that fits each range the best. This is how many fast approximations are done. See how the residual error fits nicely into a few different curves? Eureqa doesn't support this though.

tagrun5• 11 years ago

From a practical point of view, in case you have Mathematica (which is a proprietary software unfortunately), you can do this easily for any function in a general way: http://reference.wolfram.com/language/FunctionApproximations...

(MiniMaxApproximation minimizes the relative maximum error)

For the case discussed in write up, it gives

MiniMaxApproximation[Log2[x + 1], {x, {1, 2}, 2, 0}]

log2(x+1) ~ 0.177392 + 0.943626 x - 0.120232 x^2

with a maximum error of -0.000786282.

In the 4th order (that is, (4,0)), it gives

0.0399284 + 1.27027 x - 0.391233 x^2 + 0.0909038 x^3 - 0.00986308 x^4

max error: -4.83031 10^-6.

However, if the division isn't too expensive, a (2,2) approximant also works for an even better precision:

(0.00383323 + 1.42412 x + 0.336161 x^2)/(1 + 0.704317 x + 0.0598015 x^2)

max error: -1.33902 10^-7

As for [0.75, 1.5) interval, I get

0.108715 + 1.05621 x - 0.165315 x^2

which gives a max error -0.000651425 and the max relative error is around 0.0006, which is again way better than what he gives in the write up, although the error is non zero at points 0.75 and 1.5 (which is not really a necessary constraint, but an artifact of the particular polynomial he chose). The overall plot for relative error in a wider range is smooth.

Bottom line: 1) Unless you are an expert in numerical methods, consult to a software/article/book by an expert (or to the expert him/herself if you can) in the field at some point. 2) While I appreciate his efforts and the write up, I wouldn't recommend using his polynomials in practice.

(BTW, one can use frexp instead of bit-twiddling to make the C code friendlier since not everyone is familiar with the bit layout of IEEE 754 numbers.)

stagger87• 11 years ago

I am assuming a vectorized version of the function will appear in a later update, which is why frexp is not used.

Range reduction -> polynomial -> scale is the process an expert/book will recommend for log functions.

tagrun5• 11 years ago

That is the general procedure for any function of course. The problem usually lies in implementing these steps, which requires an expert.

tagrun5• 11 years ago

I realized that I missed the part that he's talking about log2(x) rather than log(x+1), which changes the interval.

Since MiniMaxApproximation gets quirky when the function is 0 in the given interval (it minimizes the relative error, so a division by 0 occurs, unless you cancel it), I'm rather using the interval [0.25, 0.5] with t=1/4:

-3.64631 + 7.93796 x - 5.30406 x^2

This gives a maximum relative error of 0.00334047, which is still an order better than his.

As for the (2,2) approximant: (-5.83261 - 16.7012 x + 22.394 x^2)/(1 + 11.2638 x + 7.81124 x^2) max rel error: 1.69763 10^-6

I checked several other values of t, and the max rel error for (2,0) approximant is roughly the same. This leads me to believe that the reason he is getting 10 times larger error is because he is imposing that errors should vanish at certain points. Not only this increases the max rel error, but it also distorts the smoothness of the function in a bad way (the kinks in the plots are a bad sign in general).

gjm113• 11 years ago

Your polynomial, since it doesn't get it exactly right at 1/4 & 1/2, will have an enormous relative error when you use it for computing logs of numbers very close to 1.

Goldberg's initial function b ignores this issue, and a glance at the graph of its absolute errors will confirm that it's the minimax poly for absolute error or very close to it. His next candidate t half-fixes this, and again looking at the graph of the relative error makes it clear that it's right. (Note that avoiding intervals containing zero will indeed allow you to get better max relative errors -- keeping the relative error finite where the function is zero effectively reduces your polynomial degree by 1.) Then finally he fully fixes it with the polynomial g, solving the more interesting problem of minimizing the max relative error for the log function everywhere, and chooses a better functional form by shifting the interval he reduces to with the polynomial s.

(The fact that his earlier polynomials do the wrong thing is clearly not a mistake but a paedagogical strategy: start with something simple and incrementally improve it.)

My apologies for spelling all this out, but it really appears that you're accusing him of making elementary mistakes without actually understanding what he's doing.

"Imposing that errors should vanish at certain points" is not a mistake, it's a necessity if you're implementing log by argument reduction plus local approximation and you care about relative error in the resulting approx_log function, and you are getting better numbers than Goldberg by solving the wrong problem.

Your (2,2) rational approximant will be much more expensive to compute than the polynomial and won't yield the speed gains Goldberg is looking for. And, again, it will give you a log function with infinite relative error at 1.

[EDITED to fix a typo and clarify something.]

malkia• 11 years ago

I've learned to be extra careful with fast log approximation. I was the tools programmer at a game studio, and few years ago was given the task to speed things a bit. There was a process that was summing neighbouring pixels, doing some filtering/pre-processing then calculating a value through log, so I've profiled that this was the slowdown, and decided to provide a faster log version. After some googling, and testing it seems everything should've work.

Few days later, black mipmaps started to appear. Turns out, I was doing testing in a wrong way, where I was thinking I'm comparing new to old, where I was comparing old to old created by other machine (a bit different results due to be expected). All because I forgot to turn off caching.

And my approximation of log2 did not really worked for the most of the numbers involved in. It was using floats, not doubles to begin with, and producing a lot of either NaN's and infinities (hence the black colors).

Anyway it was reverted, and had to reconvert everything to get back to normal.

jgalt212• 11 years ago

Why is Ebay using the log function and why is it used so heavily that it causes a performance bottleneck for them? Is this an inherited code base from their Hunch purchase?

brandmeyer2• 11 years ago

Boost used to have a minimax approximation algorithm in its codebase. Unfortunately, nobody knows enough about how it works it to maintain it, and you have to look in much older releases of Boost to get it.

alnsn• 11 years ago

This one? #include <boost/math/tools/remez.hpp>

I used some bits from it to help me in building my own log approximation but I didn't have time to dig into it.

I ended up doing Chebyshev approximation over [1, 2) interval but I extended the interval slightly to make sure that an error term at 1.0 is minimal (you can also minimise an error term at 2.0 if you want some smoothness between [1, 2) and [2, 4) intervals).

It's easy to calculate new interval values for Chebyshev approximation and it shouldn't be a problem to modify Remez but I couldn't find any Remez implementation that would let me do what I needed.

Usually, interval ends in Remez algorithm are fixed, while internal points can be moved by the algorithm. I needed a way to to fix one or two internal points (1.0 and optionally 2.0) but let the algorithm move interval ends within certian limits.

If anyone knows of any tool that let me do it, I'd be happy to play with it and publish my code. My approach is a bit similar to [1] but I'd more likely use Intel assembly or DynASM.

[1] http://gallium.inria.fr/blog/fast-vectorizable-math-approx/

PS The author of [1] uses remez from e Sollya tool but it doesn't let me fix internal points. I looked at the implementation but it's more obscure than boost/remez.

UPDATE: it's not clear from the above what my objectives were. I wanted to to produce exact values for integer powers of 2 including log2(1). This approach also gives a small relative error around x=1.

brandmeyer2• 11 years ago

Yes, that one. It comes with a driver program using the Boost.Test framework that allows you to change the range and fix the values at the interval ends if that is an objective.

Their framework does not allow you to arbitrarily fix internal values, as per your spec however.

Fast Approximate Logarithms, Part I: The Basics

Related Stories