## Loren on the Art of MATLABTurn ideas into MATLAB

Posted by Loren Shure,

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

### Contents

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


### One More Combination

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

### Note

Will replied on : 1 of 13

The error I always get in trouble with is switching between single and double precision values. While comparing results between different normalization methods from microarray utilities I spent days trying to track down the differences. Turns out there was some accumulation of error resulting from single-precision vs double-precision values.

Oscar replied on : 2 of 13

“All of the answers are correct! Welcome to the delights of floating point arithmetic.”

They are not all correct though, this is the equivalent of saying 1+1=POTATO.

I feel it would be much more helpful if the above post were a segue in the topic of machine precision and subtractive cancellation.

The machine epsilon is the smallest number we can add to 1 be able to distinguish it from 1. Matlab has a built in function that will tell you what your particular machine’s epsilon is. Try this:

>> eps

ans =

2.2204e-016


also,

a=1;
b=1+eps;
a<b

ans =

1


so our computer can still see a difference here. Now consider

a=1;
b=1+2.2204e-017;
a<b

ans =

0


BUMMER!

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.

Petter replied on : 3 of 13

matt fig replied on : 4 of 13

Loren,

I think you may have pointed to the wrong thread on the newsgroup.

Bruno Luong replied on : 5 of 13

This topic interests me greatly for many reasons. There is a gap between knowing a sum computed depend on the order and knowing *HOW* the order could affect the result.
As numerical developer, sometime it is crucial for me to know
(1) in which order the sum is carried out;
(2) The result should be stable if the sums are carried out on the same array at two different instants by the same program on the same computer.

The first point is important in order to perform a good *absolute* accuracy of the sum, taking into account the relative magnitudes and oscillation of the elements. The second point is even more important when the sum results are later mutually compared – for example when compute numerical Jacobian, gradient or Hessian.

In v.2009a, (2) does not meet http://www.mathworks.com/support/bugreports/532399, fortunately this “drawback” has been corrected from 2009B.

Still for a developer it remains for him/her to know which order the sum is carried out. Before multithreading is introduced, the sum is carried out sequentially; it is quite simple and straightforward and avoid some headache. When multithreading kicks in, things get much more complicate: I made various tests and on my computer with 2009A it seems the splitting occur somewhere in the range of 90000th element. On 2009B, it seems to be the mono-threading sum engine split artificial the array in order to imitate provide the same result as multi-threading sum. This is all good, but still very fuzzy. All I can say is it not as simple as before. I give up trying to predict on my own guessing what the sum-order is when carried out by Matlab engine. Let alone when building a reliable backward compatible code that can run consistently from version to version. Lacking any description of the sum routine, it is very frustrating. I ended up with my own MEX sum routine where I know exactly what is going on.

I believe an official document file that describes clearly the behavior of sum and various other multithreading calculations would be a big help for us.

Loren replied on : 6 of 13

Sorry for the bad link. I think it’s fixed. I have written a number of posts on numerical accuray. Here’s one of them. And here is a link to posts in the numerical accuracy category.

Bruno- In some sense, why does it matter to you which order the calculations are done? If you care, you can control that by writing the appropriate code. Absent that, the answers (not the earlier bug in sum) are all “correct”. If your overall calculation is that sensitive to order, you should, by all means, control that.

Oscar- I couldn’t agree more with you about the details of eps but I do still stand on my statement that all the answers are correct, especially to within eps(1e30) which is what you need to compare to in the example. Floating point precision is not something to ignore, and there are ways to minimize effects of cancellation between large numbers, for example – but it behooves the writer of an algorithm to understand when terms might come close to cancelling and potentially rewrite the algorithm to be less sensitive to that. In MATLAB we could have decided to sort all the elements in a vector and accumulate the positive values and negative values separately from smallest magnitude to largest. But that wouldn’t be fast, nor would it probably be what people were expecting. We just sum values as they come, in a single-threaded environment. This leads to possibly different results from a multithreaded answer, but neither answer is incorrect.

–Loren

Bruno Luong replied on : 7 of 13

Loren – yes I want to control the order, not always, but sometime. Before (read monothread Matlab) I can do so by sorting my array in appropriate order, then call SUM. Now (read multithread) how can I do? A for-loop?

Bobby Cheng replied on : 8 of 13

Oscar — I think you understand the issue. The important point here, I believe, is to understand that floating point computation has its limits. When the limits are hit, it is up to us to work around it. The problem

a=1;
b=1+2.2204e-017;
a<b  %returns true



is 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.

OysterEngineer replied on : 9 of 13

Loren,

I think it is great for you to bring up this topic regularly. Even though you’ve written about it many times, I think you need to continue to return to it, maybe once a year or so.

And, I really love the perspective you bring to this when you say that all the answers are correct. I certainly understand what you mean in that given the details of how MatLab operates, all the answers are correct.

However, for folks like me, MatLab is a tool that we use to represent reality. And, in that reality, only the 2nd answer is correct.

Merry Christmas!

Loren replied on : 10 of 13

OysterEngineer-

And Merry Christmas and happy holidays to all.

Just wondering if you have computations that are numerically very sensitive. What is usually required is what Bobby (and I) said in comments – if so, you may need to rewrite the algorithm to avoid cancellation between 2 large quantities and instead find a different way to just calculate the differences.

–Loren

Bobby Cheng replied on : 11 of 13

More importantly, if there is a need for this kind of special algorithms let us know or post here. If it looks like a good building block, we would consider adding it to MATLAB or help us prioritize work. :)

Functions like LOG1P for log(1+x), and EXPM1 for exp(x)-1 when x is very small are added for MATLAB this way.

Bruno Luong replied on : 12 of 13

Personally I don’t need a stock function that can perform an accurate sum; I already have developed several for my own uses, more or less slow, more or less accurate. I do however need a basic sum function (or the ability to set it in an compatible mode) where I do *know* the order of the additions that are carried out.

Most mathematicians spend their favorite times trying to bound the errors of the calculation. With an unspecified order of the sum, a crucial information is removed, making their works much more harder.

Kotya replied on : 13 of 13

> 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)?

At least I always try.

> 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?

I sometimes work from home, accessing my office desktop remotely. At home I have 32 bit Matlab under Mac OS, at work—64 bit Matlab under Win 7. The code is in the Dropbox folder, so I have the same code on both computers all the time.

Once I was debugging a code calculating some integrals resulting in tiny numbers. I was not interested in the numbers itself, I only wanted to have some checksum value to make sure the answer doesn’t change when I change something in the code. Due to this, I managed not to notice that I have encountered such small numbers. I was running the code on Mac, obtaining some checksum, then running the same code on Windows to make sure everything’s all right and was getting a wrong check-sum. It took quite some time before I realized that two computers give different tiny numbers.