Identical math operations in C++ and matlab differ by ~1% (10^11 in this case!)

Hello,

I have a single (VERY LONG) line of code written in both C++ and matlab that evaluates a mathematical expression. The expression is just under 31,000 characters long, so it's pretty darn complicated.

When I evaluate this expression in matlab vs. C++ I get very different numbers:

MATLAB: 31856235311612.7

C++: 31745938896594.422

However, as you'll probably notice, they're not different enough to think I made a mistake in the transfer. I basically just copied and pasted the expression from C++, so no issues there.

I have two similar expressions of lower order that I also ported; one on the order of 10^9 (differs by about 1%), and one of order 10^7 which does not suffer the same magnitude of error:

MATLAB=3087973.34551207

C++=3087973.3455120716

The question is, which is "correct", and from where are these differences arising? I'm happy to provide code, but the first entry code is a 31,000 character expression!

FYI, in C++ all variables are stored as doubles.

Answer 1

[Jan]

A much simpler example:

MATLAB uses FDLIBM for the calculation of ACOS. I've added this library to a C-mex file also and compared the results of:

a = [ 0.82924747834055368, ...
      0.55849526207902467, ...
     -0.020776474703722032];
b = [-0.82889260818746668, ...
     -0.55895647212964961, ...
      0.022465670623305872];
acos(dot(a, b))

I found differences between the values calculated by C and Matlab up to 32768 EPS, depending on the C-compiler and the input values.

Let's look on the sensitivities:

b_var = b + [0, eps, 0];
acos(dot(a, b)) - acos(dot(a, b_var))
>> 1.2434e-013

This means that tiny differences in the input are amplified by the factor 560. This seems to be small. But together with some cancellation effects, such differences can explode: [EDITED]

r = 1 / (3.13980602793225 - acos(dot(a, b)));
>>  r = 2.2518e+014

r_var = 1 / (3.13980602793225 - acos(dot(a, b_var)));
>>  r_var = 7.7648e+012

Now the analysis of the sensitivity would reveal, that the algorithm is such instable, that it cannot be used to determine a single siginificant digit.

Conclusion: Calculating values using a 31'000 characters term without an analysis of the sensitivity is not meaningful. For a high sensitivity the results have a large larger random part.

[M S]

Thanks for the response. Are you asking that I add some epsilon to each parameter and see how much the expression changes? I guess if the expressions are already different when their parameters are identical, I don't quite understand why further perturbing the inputs will inform on what's going wrong in the calculation.

The other thing is, I *expect* differences between any two languages in numerical computation. The error you find, Jan, in your acos is fairly small though compared to the order of the numbers inputted. For my numbers, the error is on the order of 1 magnitude below the order of magnitude of the numbers, which in my mind is insanely high and unacceptable.

[Jan]

@M S: See the EDITED part. "Small" is a very relative term in numerics: subtract two "small" numbers from each other and build the reciprocal - suddenly it looks really large.

The variation of the inputs allows to estimate, if the solution is chaotic or near to chaotic in the first order approximation. If a solution is very sensitive, I hesitate to call it "solution". Imagine I simulate a pencil standing on the tip at first. Would you trust my work, if I claim that my simulation shows without doubt, that the pencil will fall to the left?!

No. A numerical result without sensitivity analysis is not scientifically meaningful. If the programming language has such a big influence, I doubt that the solution is stable. I do not think that you can analyse a 31'000 character equation symbolically, a numerical sensitivity test is obligatory.

[M S]

OK, so for the big one. if I allow .01 absolute error to one parameter:

=== MATLAB ===

+0.01 = 31816475818623.8

-0.01 = 31675519324512.4

=== C++ ===

+0.01 = 31816475818623.801

-0.01 = 31675519324512.441

So I think the answer is, they are moderately unstable. But you'll notice the answers are identical. Yup, I swapped 2 very similar parameter inputs, and when I fixed the code I now get identical solutions between C++ and matlab. I am so sorry to have wasted your time!! That said, I would not have discovered this bug had I not done your sensitivity analysis, so thank you.

Answer 2

[James Tursa]

Differences can arise from different order of operations, temporary results done in 80-bit vs 64-bit arithmetic, different rounding modes, etc, etc. As for which one is "correct", maybe neither is. Unless you examine the expression in detail to look for cancellation error etc who knows? How the heck do you expect to debug a 31000 character expression? I doubt anyone on this site will be willing to do so for you.

[M S]

As I indicated, the expression is now just copied and pasted from C++ to matlab, thus there are no order of operation mistakes or the like, and no need to debug the expression. The two expressions are identical. And besides, I agree with you, it would be ridiculous to ask anyone to debug this expression if I thought I had made a mistake in the transfer! I would never ask for that.

[James Tursa]

Just because you copied the same expression from C++ to MATLAB does *not* guarantee that there are no order of operation differences. That was the point of my comment ... the compiler might be reordering things around without you being aware of it.