There are frequent posts on the MATLAB newsgroup as well as lots of questions posed to Technical Support about the floating point accuracy of MATLAB. Many people think they have found a bug when some seemingly simple arithmetic doesn't give the intuitive answer. The issue that arises has to do with the finite number of bits that a computational device uses to store and operate on values.
Contents
Calculator Demo
Have you ever impatiently waited for an elevator and continued pushing the UP button in the hopes that the elevator will arrive faster? If so, you may also have impatiently pushed buttons repetitively on a calculator (assuming you know what a calculator is!). It's not just MATLAB that computes values in a way that seems to defy our expectations sometimes.
clear
format longHere's what I tried. I put the number 7 into my calculator. I think it calculates in single precision, but I'm not sure. It shows 7 digits after the decimal point when I take the square root, and doing it several times, I get
sqrts7inSingle = single([7 2.6457513 1.6265765 1.275373])
sqrts7inSingle = 7.0000000 2.6457512 1.6265765 1.2753730
Next, I take the final number and square it 3 times, yielding
squaringSqrt7inSingle = single([1.275373 1.6265762 2.6457501 6.99999935])
squaringSqrt7inSingle = 1.2753730 1.6265762 2.6457500 6.9999995
Compare the difference between the original and final values.
sqrts7inSingle(1)-squaringSqrt7inSingle(end)
ans = 4.7683716e-007
This number turns out to be very interesting in this context. Compare it to the distance between a single precision value of 7 and its next closest value, also know as the floating-point relative accuracy.
eps(single(7))
ans = 4.7683716e-007
We get the same value! So, we can demonstrate the effects of finite storage size for values on a handheld calculator. Let's try the equivalent in MATLAB now.
Accuracy in MATLAB
Let's do a similar experiment in MATLAB using single precision so we have a similar situation to the one when I used the calculator. In this case, I am going to keep taking the square root until the value reaches 1.
num = single(7); count = 0; tol = eps('single'); dnum = sqrt(num); while abs(dnum-1)>=tol count = count+1; dnum = sqrt(dnum); end count
count =
23
So, it took 23 iterations before we reached the number 1. But now we're in a funny situation. There is no way to square exactly 1 and ever reach our original value of 7.
Using Double Precision
Let's repeat the same experiment from the calculator in MATLAB using double precision.
num = 7; sqrtd = sqrt(num); sqrtd(2) = sqrt(sqrtd(1)); sqrtd(3) = sqrt(sqrtd(2)) prods(1) = sqrtd(end); prods(2) = prods(1)^2; prods(3) = prods(2)^2 finalNum = prods(3)^2
sqrtd = 2.64575131106459 1.62657656169779 1.27537310685845 prods = 1.27537310685845 1.62657656169779 2.64575131106459 finalNum = 7.00000000000001
Find the difference between finalNum and num and compare to the proper relative floating-point accuracy:
diffNum = finalNum-num acc = eps(num)*6 diffNum-acc
diffNum =
5.329070518200751e-015
acc =
5.329070518200751e-015
ans =
0
Why did I use 6 when I calculated acc? Because I performed 6 floating point operations to get the my final answer, 3 sqrt applications, followed by squaring the answers 3 times.
Typical MATLAB Pitfall
The most common MATLAB pitfall I run across is when users check equality of values generated from the : operator. The code starts innocently enough. I've actually turned the logic around here a bit, so we will print out values in the loop where the loop counter may not be exactly what the user expects.
format short for ind = 0:.1:1; if ind ~= fix(10*ind)/10 disp(ind - fix(10*ind)/10) end end
5.5511e-017 1.1102e-016 1.1102e-016
And then the question is, why do any of the values print out? It's because computers can't represent all numbers exactly given a fixed storage size such as double precision.
Conclusions
This was a very simplified explanation about floating point, but one, I hope, that is easily understood. Let me know.
References
Here are some pointers to more resources.
Published with MATLAB® 7.2
Hi,
I have just found this floating point precision “bug” in an extremely simple operation:
e.g. if i type:
>> 0.3 == 0.1*3
ans =
0
It would be extremely useful if Matlab would give out a warning for such cases, because this can really throw you off track.
Ilya,
This is not a bug, but due to limited accuracy. What kind of warning would you expect, one in the editor perhaps?
–Loren
I have similar floating-point precision problem in my program. For instance;
a=0.1+0.2-0.3
a =
5.551115123125783e-017
I want “a” to be “0″ not that precise. How am I going to fix this problem in my code? If anyone could help me, I really appreciate.
Bashar-
You are asking for the impossible from floating point. Please read the links in the post to the FAQ, etc.
–Loren
Bashar, you could define a tolerance then round your computations accordingly. You might check out File Exchange submission 4261.
http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=4261&objectType=file
(test HTML link: <a href=”http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=4261&objectType=file”round2)
matlab must be doing some pretty strange stuff internally for the following: 0.1+0.2-0.3 = 5.551115123125783e-017
My question is why have the matlab dev team not sorted this out in 2008! Mathematica, Java, C# and C++ do not suffer such FP errors for such simple operations (with more complex operations then a certain degree of error is perhaps acceptable?) - Matlab seems to be alone in still have these really basic errors.
Matlab dev team care to comment?
Hugh-
This is not an error. If MATLAB were built to get that answer to be 0, we’d have to do a bunch of operations to sort the data in the “right” order. How would MATLAB actually know that IS the way the values should be accumulated from the real-world? You can have MATLAB perform calculations in whatever order you wish by using parentheses. Or introduce a function that sorts the data appropriately and then does the sums and differences.
As long as there is no way to perfectly represent some values, floating point calculations are subject to round-off error. As for the other languages, it depends on how the code is written so they get exactly 0. If the constants are built into the program rather than user-supplied at run-time, it’s possible that a C++ compiler will reorder the operations on constants. And Mathematica might be doing the calculations using extended precision (but I don’t know that for certain).
–Loren
Hugh—The C program below prints “Not equal” when compiled with gcc and run on Linux. An equivalent Java program also prints “Not equal”.
#include <stdio.h> int main(int nargin, char *argv[]) { double a = 0.1; double b = 0.3; if ( (a+a+a) != b ) { fprintf(stdout, "Not equal\n"); } else { fprintf(stdout, "Equal\n"); } return(0); }And the same for C++. This program prints “Not equal”:
#include <iostream> int main() { double a = 0.1; double b = 0.3; if ( (a+a+a) != b ) { std::cout << "Not equal\n"; } else { std::cout << "Equal\n"; } return(0); }