Why is $\exp(\ln(x))-x\neq0$ in floating point arithmetic?

Question

Analytically, the expression

$$\exp(\ln(x))-x \enspace,$$

should give 0.

However, in Matlab, it does not.

x = linspace(1, 10, 10);
exp(log(x)) - x;

for $x \in [1, 10]$, the answer is

0   
0   
-4.44089209850063e-16   
0   
0   
0   
0   
-1.77635683940025e-15   
1.77635683940025e-15    
1.77635683940025e-15

What is happening and how can I fix this?

Answer 1

The numeric precision is not perfect. You get rounding errors during your computation. When working with floats, don't check if they are = 0, but check if their absolute distance to 0 is smaller than some epsilon.

Answer 2

[EDIT] An alternate view: 64-bit floating numbers represent a discrete set $S$. For a function $f$ to be exactly invertible, it should be a bijection from $S$ to $S$. Suppose we are interested in a fast growing function like $\exp$. At one point, $\exp x > x$. If $M$ is the maximum element from $S$, then $f(M)\ge M$. The strict version $f(M)> M$ is not allowed, because there is no higher number that $M$ in $S$. Morally, a monotonic function acting like a bijection on a finite set is a permutation, and has very limited degrees of freedom, it is likely to be limited to the identity, or or a counter-identity which reverts the order of elements in the set. Outside these trivial examples, monotonic functions cannot be bijections.

In practice, it is due to finite-precision numerical computations (for instance 64-bit floating point), that especially affect functions involving irrationals numbers. Integer or sum-of-powers-of-two calculations may be less sensitive. If you dare to plot the numerical error, you get something like:

with Matlab code

x=linspace(1,10000,10000);
y = exp(log(x))-x;
plot(x,y,'.');
xlabel('x');
ylabel('y');

So you can see that the maximum error spreads with the value. If you want to do something, you can check and threshold a relative error, especially if your variable spread over a large scale.

It get even weirder, if you just swap the functions, say: $\log(\exp(x))-x$, the error is in this case numerically zero.

x=linspace(1,700,10000);
y = log(exp(x))-x;
plot(x,y,'.');
xlabel('x');
ylabel('y');

Answer 3

The basic issue is that generally $\exp(\ln(x)) \neq x$ when using a fixed-precision floating-point format. This is true independent of the accuracy with which the exponential and logarithm functions are being computed, due to the pidgeonhole principle.

Without loss of generality, this effect is most easily demonstrated by using a binary logarithm log2 in conjunction with one of the IEEE-754 binary floating-point formats. With such a format, each binade comprises an identical number of encodings. If we compute log2 of two consecutive binades $[4,8)$ and $[8,16)$ all results will all fall into the single binade $[2,4)$. Since there are twice as many source encodings than target encodings are available, we must have plenty of collisions by the pidgeonhole principle.

All logarithms are contracting operations. This means that when using a fixed-precision floating-point computation, multiple inputs can map to the same output. Information is lost due to this compression, and as a consequence, when we expand the numbers in $[2,4)$ back to $[4,16)$ with exp2, for many $x$ we cannot get back the starting value even if both functions use a correctly rounded implementation.

Reversing the order of operations, that is, use of an expanding operation followed by the contracting inverse operation, we can sometimes guarantee that every operand $x$ round-trips successfully, as long as there is no intermediate overflow or underflow, and assuming correctly-rounded operations. I do not know (off-hand) whether this guarantees $\ln(\exp(x)) = x$, but the graph in Laurent Duval's answer seems to suggest it. This paper shows that the square root of the square of a floating-point number $x$ is exactly $\lvert x \lvert$, under the restrictions stated.

Answer 4

If you have the symbolic toolbox you can fix the problem by

syms x
y = exp(log(x))-x
% y = 0