7 ms·

Fast Approximate Logarithms, Part I: The Basics

70 points by r4um 11y ago

sdenton4 11y 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 11y 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 http://arxiv.org/abs/1502.02551 and http://arxiv.org/abs/1412.7024 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 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-mixed-precision-mode-up-to-32gb-of-ram http://techreport.com/news/27978/nvidia-pascal-to-feature-mi...
  - sdenton4 11y ago
    Awesome, thanks for the pointers!
  - ot 11y 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.pdf http://www.eecs.tufts.edu/~dsculley/papers/round-model-icml....
- Houshalter 11y 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 11y 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?
    - Houshalter 11y ago
      Using the log domain is cool but log summation still has the same problem. Logs just transform numbers into a different range, they don't magically add more bits of precision. That is adding a very small number in the log domain to another much larger number in the log domain, will just round it down.
deleted 11y ago
[deleted]
TheLoneWolfling 11y 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 11y ago
As in Math Toolkit for Real-Time Programming by Jack Crenshaw
paulpauper 11y 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 11y 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.
  - TheLoneWolfling 11y ago
    I wonder why eBay is using the log function?
    - pmiller2 11y ago
      I'm more curious why they use it so extensively that tripling the speed of it is significant to them.
      - TheLoneWolfling 11y ago
        That's a much better way to phrase that, thanks.
    - smitherfield 11y ago
      Probably has to do with the weighting of search results, reviews, etc.
rg2004 11y ago
If your bottleneck is log() and you're not memory-limited, why not use a lookup table?
- tgbrter 11y ago
  This is faster than a look-up table at the cost of some accuracy. Also no cache misses.
- walshemj 11y ago
  which is how pre calculators you used to do it by hand using Log tables
  - Houshalter 11y ago
    Computing log tables is what the would-be first computer was intended to do.
minimax 11y 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.html http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.ht...
- bumbledraven 11y 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/ http://www.ebaytechblog.com/author/dgoldberg/
boulos 11y 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/syrah/FixedVectorMath.h https://github.com/boulos/syrah/blob/master/src/include/syra... . I hope in a later post he mentions Sollya (http://sollya.gforge.inria.fr http://sollya.gforge.inria.fr), I found it invaluable for playing around with different precision vs speed.
Houshalter 11y ago
I tried evolving a good approximation with Eureqa: https://i.imgur.com/1CTTjYD.png?1 https://i.imgur.com/1CTTjYD.png?1
- dalke 11y ago
  It looks like all of them (they use sqrt() and x^y) will be slower than calling log() directly.
  - Houshalter 11y ago
    Here is it without those operations: https://i.imgur.com/xhyNnU4.png?1 https://i.imgur.com/xhyNnU4.png?1 EDIT: Without sin/cos: https://i.imgur.com/mTwgKKO.png 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.
    - dalke 11y ago
      The original posting ended with a second-order polynomial, a few bitwise operations, and an if statement. All of the ones you just posted are more computationally difficult. The non-trig ones, for example, include a division, and the best fit is a higher-order polynomial. There's no reason to believe it will be as fast to evaluate, so the question would be, is it more accurate and is the accuracy worth the additional time? It would be ill-advised to use one of those if it ended up being both less accurate than log() and slower. Concerning lookup tables, as tgbrter answered in this thread, they are likely slower. It still needs the scaling, and also needs two value lookups followed by an interpolation. The cost of doing the interpolation will be about the same as evaluating the second order polynomial, but the pipeline will be waiting on the two memory lookups. Assuming everything in is L2 cache, that might not be a problem. But depending on the problem, it may be in L3 or (shudder) main memory, which is eons away. So as tgbrter also pointed out, you may end up with variable calculation time depending on the cache.
jgalt212 11y 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?
malkia 11y 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.
tagrun 11y 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/tutorial/FunctionApproximations.html 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 11y 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.
  - tagrun 11y ago
    That is the general procedure for any function of course. The problem usually lies in implementing these steps, which requires an expert.
- 11y ago
brandmeyer 11y 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 11y 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/ 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.
  - brandmeyer 11y 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.

70 points by r4um 11y 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 11y 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 http://arxiv.org/abs/1502.02551 and http://arxiv.org/abs/1412.7024 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 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-mixed-precision-mode-up-to-32gb-of-ram http://techreport.com/news/27978/nvidia-pascal-to-feature-mi...
- sdenton4 11y ago
  Awesome, thanks for the pointers!
- ot 11y 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.pdf http://www.eecs.tufts.edu/~dsculley/papers/round-model-icml....
Houshalter 11y 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 11y 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?
  - Houshalter 11y ago
    Using the log domain is cool but log summation still has the same problem. Logs just transform numbers into a different range, they don't magically add more bits of precision. That is adding a very small number in the log domain to another much larger number in the log domain, will just round it down.

deleted 11y ago

[deleted]

TheLoneWolfling 11y 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 11y ago

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

paulpauper 11y 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 11y 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.
- TheLoneWolfling 11y ago
  I wonder why eBay is using the log function?
  - pmiller2 11y ago
    I'm more curious why they use it so extensively that tripling the speed of it is significant to them.
    - TheLoneWolfling 11y ago
      That's a much better way to phrase that, thanks.
  - smitherfield 11y ago
    Probably has to do with the weighting of search results, reviews, etc.

rg2004 11y ago

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

tgbrter 11y ago
This is faster than a look-up table at the cost of some accuracy. Also no cache misses.
walshemj 11y ago
which is how pre calculators you used to do it by hand using Log tables
- Houshalter 11y ago
  Computing log tables is what the would-be first computer was intended to do.

minimax 11y 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.html http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.ht...

bumbledraven 11y 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/ http://www.ebaytechblog.com/author/dgoldberg/

boulos 11y 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/syrah/FixedVectorMath.h https://github.com/boulos/syrah/blob/master/src/include/syra... . I hope in a later post he mentions Sollya (http://sollya.gforge.inria.fr http://sollya.gforge.inria.fr), I found it invaluable for playing around with different precision vs speed.

Houshalter 11y ago

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

dalke 11y ago
It looks like all of them (they use sqrt() and x^y) will be slower than calling log() directly.
- Houshalter 11y ago
  Here is it without those operations: https://i.imgur.com/xhyNnU4.png?1 https://i.imgur.com/xhyNnU4.png?1 EDIT: Without sin/cos: https://i.imgur.com/mTwgKKO.png 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.
  - dalke 11y ago
    The original posting ended with a second-order polynomial, a few bitwise operations, and an if statement. All of the ones you just posted are more computationally difficult. The non-trig ones, for example, include a division, and the best fit is a higher-order polynomial. There's no reason to believe it will be as fast to evaluate, so the question would be, is it more accurate and is the accuracy worth the additional time? It would be ill-advised to use one of those if it ended up being both less accurate than log() and slower. Concerning lookup tables, as tgbrter answered in this thread, they are likely slower. It still needs the scaling, and also needs two value lookups followed by an interpolation. The cost of doing the interpolation will be about the same as evaluating the second order polynomial, but the pipeline will be waiting on the two memory lookups. Assuming everything in is L2 cache, that might not be a problem. But depending on the problem, it may be in L3 or (shudder) main memory, which is eons away. So as tgbrter also pointed out, you may end up with variable calculation time depending on the cache.

jgalt212 11y 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?

malkia 11y 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.

tagrun 11y 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/tutorial/FunctionApproximations.html 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 11y 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.
- tagrun 11y ago
  That is the general procedure for any function of course. The problem usually lies in implementing these steps, which requires an expert.
11y ago

brandmeyer 11y 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 11y 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/ 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.
- brandmeyer 11y 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.