6
$\begingroup$

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?

$\endgroup$
  • 5
    $\begingroup$ Read this. $\endgroup$ – horchler Apr 23 '16 at 13:10
  • $\begingroup$ matlab also has syms so you could use syms x; exp(log(x))-x and this gives 0 $\endgroup$ – Nasser Apr 23 '16 at 14:22
  • $\begingroup$ Whenever you see an error of 1e-15 it is bound to be roundoff errors of floating point representation. A 64-bit "double precession" floating point number is good to decimal 15 digits about. $\endgroup$ – ja72 May 6 '16 at 22:08
10
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ This is correct. The virtual representation of any number has limited precision (matlab uses double precision by default) and manipulations (e.g. exp(log(x))) are computed with finite precision. After subtracting x, there is a loss of precision. $\endgroup$ – Charles Apr 22 '16 at 20:58
  • 1
    $\begingroup$ I was about to say "happy round-off to you". $\endgroup$ – EngrStudent - Reinstate Monica Apr 23 '16 at 2:23
7
$\begingroup$

[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: the plane of nummerical errors

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');

enter image description here

$\endgroup$
  • 1
    $\begingroup$ Interesting! I guess this is really rooted in floating point arithmetics, as gnuplot gives the same result both on AMD and Intel CPUs. $\endgroup$ – Joce May 13 '16 at 11:34
2
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

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

syms x
y = exp(log(x))-x
% y = 0
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.