3

It is well known that when x is close to 0, computing exp(x) - 1 introduces significant cancelation errors. As such, we have expm1 implemented in c99 and python.

My question is how expm1 avoids cancellation error? Can anyone give me a general idea without too many nitty-gritty implementation details?

user1559897
  • 1,807
  • 2
    I don't know how it's implemented, but the idea is likely this. If $x$ is close to zero then we can use a Taylor series expansion $e^x - 1 = x + x^2/2 + \ldots$. Note that the usual routine for $e^x$ might also use a Taylor series, but for very small $x$ we have $1+x \simeq 1$ so the precision will be lost as soon as we add $1$ (e.g. we can easily represent $10^{-100}$, but $1+10^{-100} = 1$). – Winther Oct 06 '19 at 20:01
  • You can do all kinds of things if you are not constrained to first put $\exp(x)$ into a floating-point variable with only $53$ binary bits of precision, and then subtract $1.$ You could simulate a $1000$-bit floating-point register, subtract $1$ from the result there, and then round to $53$ significant bits in order to get a return value. I'm reasonably sure that is not how it is done, but the point is we're not in a good position to guess. – David K Oct 06 '19 at 20:08
  • There is mention of multiple implementations in http://www.math.utah.edu/~beebe/reports/expm1.pdf -- what you're looking for may or may not be in that list. – David K Oct 06 '19 at 20:09
  • @Winther If you do the taylor series this way, when x is negative, the even terms are positive and odd terms are negative. When adding them up, there could be significant cancellation errors. It cant be implemented using this taylor series. – user1559897 Oct 06 '19 at 21:18
  • I'm talking about how to implement it for small $x$. When $x$ is small then $x^2$ is tiny so there is no big problem with that. It's not like the positive terms and the negative terms are of the same size which seems to be what you are thinking of. – Winther Oct 06 '19 at 21:20
  • @Winther well, the positive and negative terms have to be the same size. When x is smaill like -1e-10, exp(x) - 1 is 0. – user1559897 Oct 06 '19 at 21:22
  • @user1559897 If $x = 10^{-10}$ then $e^x -1 \approx 10^{-10}$. The corrections to this first order approximation is $10^{-20}$. Thus just using just the first term gives the result to a fractional accuracy of $10^{-10}$, i.e. $0.0000000001$ while the correct result being something like $0.000000000100000000001\ldots$. – Winther Oct 06 '19 at 21:25
  • @Winther Do you agree that if we were to calculate exp(x) when x is very negative like -100, we cant use the Taylor series exp(x) = 1 + x + x^2/2 ? – user1559897 Oct 06 '19 at 21:29
  • @user1559897 As I wrote many times, the comments above is for small $x$ because you wrote "It is well known that when x is close to 0" (and this is the region where we loose a lot of precision). Very negative $x$ is another case, however in this case it's simply $-1$ and no I would not use an alternating Taylor series to compute $e^{-100}$. – Winther Oct 06 '19 at 21:33
  • I think your confusion seems to be that you think that this problem is caused by cancellation errors (like if we were to compute $e^{-100}$ using a high order Taylor series). That is not the problem. The main use of this function is for $x$ close to $0$ and the problem is the issue I pointed out in my answer below. – Winther Oct 06 '19 at 21:39

2 Answers2

3

The reason for making the expm1 function is to avoid precision loss when computing $e^x-1$ when $x$ is close to $0$. For large $x$ or very negative $x$ the expm1 function will give the same result as computing the usual exp and then subtracting $1$ so we will only focus on $x\approx 0$ where differences are found. But I can mention that a common way of computing exp (and a similar thing works for for expm1) for large $|x|$ is to write $e^{x} = (e^{x/2})^2 = ((e^{x/4})^2)^2$ and so on (this trick is explained more in the other answer). This allows us to compute exp for a much smaller argument, for which a Taylor series-like approach works well, and then do a few multiplications to get back the value we want.

Now to explain what happens close to $x = 0$ lets use the simplest way of computing $e^x$ for very small $x$, namely ${\rm exp}(x) = 1 + x$ as an example. Likewise our implementation of expm1 will be ${\rm expm1}(x) = x$. This is not far away from how an actual implementation would do it for the smallest $x$.

To understand why expm1$(x)$ is superior to exp$(x)-1$ close to $x = 0$ let's quickly review how floating point numbers works. They take the form $a\cdot 2^b$ where we have a certain number of bits to represent $a$ and $b$ (here $a$ is some number between $1$ and $2$ and $b$ some number between $-N$ and $N$). For the case of double precision numbers $a$ has a precision of $10^{-16}$ and $N\sim 1023$. This is very simplified, but it's all we need here.

The smallest number we can represent above $1$ will be $\approx 1 + 10^{-16}$ (the exponent $b=0$ so the precision is the precision of $a$). On the other hand the smallest number next to zero we can represent is $2^{-1023} \approx 10^{-308}$ (here $a=1$ and the precision is determined by $b$). This is the key point to understand: numbers we can represent aren't spaced evenly (like $0,\epsilon,2\epsilon,3\epsilon,\ldots, 1-\epsilon,1,1+\epsilon, \ldots$) but the size of the spacing between two numbers we can represent depends on which number you are looking at and it's much smaller close to $0$ than close to $1$.

Because of this if $x$ is less than $10^{-16}$ then $1+x$ simply evaluates to $1$ and $(1+x) - 1$ evaluates to $0$. If we want to compute $\frac{e^x-1}{x}$ using the $e^x$ function this will therefore give the very wrong result $\frac{0}{x} = 0$. If we on the other hand use the expm1 function then we have $\frac{(e^x-1)}{x} = \frac{x}{x} = 1$ which is correct (of course if $|x|<10^{-308}$ then this will also give $0$ as we cannot represent smaller numbers).

This was just a simple example, but it explains the main problem: we get huge cancellation errors from adding $1$ to something small and then later subtracting $1$ when computing $e^x-1$. To read more about how it's actually implemented see this note and the Google Go implementation (these links were given in comments by David K and Bungo).

Winther
  • 24,478
  • As I said in the comment, if you do the taylor series this way, when x is negative, the even terms are positive and odd terms are negative. When adding them up, there could be significant cancellation errors. It cant be implemented using this taylor series. – user1559897 Oct 06 '19 at 21:20
  • 1
    @user1559897: That is why you reduce the argument as $x=k\ln(2)+u$, where $\ln2$ is coded with a much higher precision to avoid rounding errors, and then only have to compute $\exp(u)$ in only a small interval around $0$. Then to make the computation symmetric about $x=0$ use $\exp(x)=\frac{\exp(x/2)}{\exp(-x/2)}=\frac{1+\tanh(x/2)}{1-\tanh(x/2)}$ in some way and use a polynomial approximation of $xg(x^2)=\tanh(x/2)$ or $h(x^2)=x\coth(x/2)$ to get a finite formula. – Lutz Lehmann Oct 07 '19 at 07:40
2

Around $x=0$ the exponential is computed as some version of $1+xp(x)$ or $\frac{1+xg(x^2)}{1-xg(x^2)}$ or similar with some polynomials that approximate the exact term behind them in an uniform fashion on some interval.

For other values the logarithm laws $e^x=(e^{x/2^m})^{2^m}$ or $e^x=2^k3^me^{x-k\ln2-m\ln3}$ are used to reduce the argument to some small value.

Then it is easy to modify those approximating expressions to return the expm1(x) value, in the mentioned cases $xp(x)$ or $\frac{2xg(x^2)}{1-xg(x^2)}$.

Lutz Lehmann
  • 126,666