There continue to be a flurry of queries suggesting that MATLAB has bugs when it comes to certain operations like addition with more than values. Sometimes what prompts this is a user noticing that he or she gets different answers when using MATLAB on multiple cores or processors and comparing those results to running on a single core. Here is one example exchange in the MATLAB newsgroup.
In actuality, you can always get different results doing floating point sums in MATLAB, even sticking to a single thread. Most languages do not guarantee an order of execution of instructions for a given line of code, and what you get might be compiler dependent (this is true for C). Having multiple threads simply compounds the issue even more since the compiler has more options for how it will control the statement execution.
In MATLAB, as you add more elements to a variable, you may see different results depending on how many threads are being used. The reason is numerical roundoff and is not a bug. Fortunately, it is easy to simulate having more than one core or thread for some simple computations by forcing the computational order using parentheses. Let's try out a small example. We'll compute the following expression several ways and compare the results.
1e-16 + 1e-16 + 1e30 - 1e30
One a single, the sum would get evaluate from left to right. Using parentheses, we can force this to happen.
sequentialSum = ((1e-16 + 1e-16) + 1e30) - 1e30
sequentialSum = 0
The most logical way a vector would get split for between two threads is the first half on one and the second on the other. We simulate that here.
twoThreadSum = (1e-16 + 1e-16) + (1e30 - 1e30)
twoThreadSum = 2e-016
Let's try taking every other element instead. That would essentially be
alternateSum = (1e-16 + 1e30) + (1e-16 - 1e30)
alternateSum = 0
All of the answers are correct! Welcome to the delights of floating point arithmetic. Have you encountered a situation similar to one of these that puzzled you? How did you figure out that the issue was really a non-issue? Do you program defensively when it matters (e.g., don't check for equality for floating point answers, but test that the values lie within some limited range)? Let me know here.
Get the MATLAB code
Published with MATLAB® 7.9
Comments are closed.
13 CommentsOldest to Newest
>> eps ans = 2.2204e-016also,
a=1; b=1+eps; a<b ans = 1so our computer can still see a difference here. Now consider
a=1; b=1+2.2204e-017; a<b ans = 0BUMMER! You'll come across this problem in things as simple as the quadratic formula and Taylor Expansions. You get around it by being aware of the numerical behavior of your algorithm.
a=1; b=1+2.2204e-017; a<b %returns trueis interesting because it does not cause any issue 99% of the time and yet it points out a potential flaw in the algorithm. Moreover, one cannot reimplement "+" differently. The calculation of b essentially asking the computer to store a number that is impossible to represent exactly in IEEE double precision. The computer tries its best but for some this still represents a failure. That leaves only one thing. People must intervene, aka no free lunch here. First it is important to note that it is too late by the time we do a < b. Since the root cause is the computation of b, we must be careful not to throw away information. Adding the two numbers with huge differences in order is throwing away information. It is not uncommon to see code rewritten to cater for extreme cases
a=1; b=1, c=2.2204e-017; (a < b && c >= 0) || (a == b && c > 0)Of course, this code does not work for all values of a, b, c. But I think one can see the point. More importantly, this is still very narrow in scope of consideration. I assume c is computed correctly in the first place. Should we make MATLAB smarter to extend the precision when needed? Maybe. But in a lot of cases, it will merely be masking the fundamental issues in the algorithm and delay the problem. To me, the only long term solution is to write carefully crafted numerical algorithms.