Fast logistic function

Question

I am looking for a way to implement a fast logistic function. The classic definition of logistic function is:

y(x) = 1 / (1 + (1/e^x)) where ^ is exponentiation.

or equally: y(x) = (e^x) / (e^x + 1)

However, my special logistic function has a base E instead of e, so:

y(x) = E^x / (E^x + 1)

E and x in my case are 32-bit integer, fixed point in base 2 of 16.16 type. E is so near as possible to the real constant e.

What is the fastest algorithm for such function ? I would prefer no lookup tables. Something like bit-manipulations should be perfect.

---

UPDATE:

Intuitively, I feel that there exists a very fast and simple method, based on Euler's formula:

e^(i*x) = cos(x) + i*sin(x)

where i is the imaginary unit.

Answer 1

Powers do not easily convert to bit shifts, because

E^x = 2^log2(E^x) = 2^(x*log2(E))

and this will give you a fractional number of bits to shift. You could calculate x*log2(E) and then decompose the power into a sum of separate bit shifts, but this is unlikely be as fast as a lookup table followed by a couple of iterations of Newton's method. Then again, this depends on how expensive your fixed-point arithmetic operations are.

Answer 2

Keep overflow and underflow in mind. The code exp(x)/(1 + exp(x)) will overflow for large positive values of x even though the result should be approximately 1. You can avoid this by using 1/(1 + exp(-x)). This will underflow to 0 for large negative values of x, but that may be OK depending on your context since the exact result is nearly zero in that case.

If your code is being called for large positive or negative values of x, you could save a little time by returning 0 or 1 for values that are going to be represented as 0 or 1 anyway and avoid calling exp. So your code could be something like

if (x > positive cutoff) return 1;
if (x < negative cutoff) return 0;
return 1 / (1 + exp(-x))

Whether that's actually more efficient depends on your specific environment and how often you get arguments past the cutoff values.

Answer 3

I'm going to talk about floats, not ints, because that is where the technology seems to be.

The standard way of computing functions like this is to use some special case logic so that you only have to represent the function within some range [a, b], approximate this by a rational function - one polynomial divided by another, and then patch up whatever you had to do to reduce the range. The source of exp(x) at http://www.netlib.org/fdlibm/e_exp.c seems to follow this pattern.

This gives you an approximation of exp(x) of the form a(x)/b(x). You actually want 1/(1+exp(-x)). You should be able to rearrange the implementation of a(x)/b(x) so as to get it to do b(-x)/(a(-x)+b(-x)), so that you have just one divide instruction, inside the rearranged exp(x), instead of one divide instruction inside it and one outside it.

This will save you something, depending on how much more expensive divide is on your machine - it might be noticeable if your inner loop really is 90% calls to the logistic function. The pattern of range reduction plus rational approximation is so firmly entrenched that I doubt that you will do much better without sacrificing a good deal of accuracy, although if you are resorting to integers you may be prepared to do that.

I dare say you could transfer this into the fixed point world if you worked hard enough. I'm afraid I'd be inclined to go back to linear interpolation between values in a table instead, that is, assuming that I couldn't find a way to get the logistic function out of the inner loop.