0

I am using matlab. I have a function f(x) and I want to apply f(x) to a set of values. So I wrote 2 codes. The first code is a simple for loop. At one point x0, inside this for loop, I find that f(x0)=1.0000 and then I use f(x0)-1=-4.7684e-07.

My second code is using arrayfun on f(x). And at the same input value x0, I found that the results is 1.0000 but arrayfun(f,x0)-1=4.7684e-07!

This error 4.7684e-07 looks tiny. But the for loop gives me a number below 1 and the arrayfun gives me a number above 1. This is really a big difference in my work, as my subsequent computations largely hinges on whether this number is below 1 or above 1, as this number is supposed to be a probability.

Now my question is: why arrayfun has such problem? There is no random numbers in my code, why arrayfun generates a different result than for loop? which one should I trust? Is there a way to avoid this kind of precision problem? Note that in this code, all of my variables are in single type. Is this causing the problem?

KevinKim
  • 1,382
  • 3
  • 18
  • 34
  • this seems indeed odd as `arrayfun` works like a loop under the hood. How about a [MVE](https://stackoverflow.com/help/minimal-reproducible-example)? – max Apr 04 '20 at 05:50
  • 1
    @max let me try to see if I can modify my original code. – KevinKim Apr 04 '20 at 06:01
  • Maybe the difference is already earlier and one of the x0 isn't exactly the same? When debugging your code, make sure to look at 17 significant digits. This is more than the 15 `format long` will show. Code you could use to print: https://stackoverflow.com/a/35626253/2732801 – Daniel Apr 04 '20 at 08:25
  • 1
    is there a reason why you work with `single`? This is indeed not much precision. – Robert Seifert Apr 04 '20 at 10:00
  • 1
    Please read [mre]. We cannot guess at what the problem is if we don’t know what you are doing. Your description is rich, but not complete. Show what you did, don’t tell what you did! – Cris Luengo Apr 04 '20 at 17:13

1 Answers1

0

addition and multiplication are not commutative under limited precision, and results depends on the order. For example, adding (a+b+c+d) is expected to give you different answers, up to the floating point round-off error, compared to (a+c+d+b). This is especially an issue in parallel computing where the orders of the accumulation from multiple parallel thread is not ensured.

You can either force the operation order or use higher precision, like double, to reduce this error.

FangQ
  • 1,444
  • 10
  • 18