High dimension matrix summation: a + b != b + a

I have two matrices X1 and X2. Both are 10x15x10x10 dimensions. There are three scalars: ,and .

Let i={i1,i2,i3,i4} denote the index of an element. For elements with the same index in X1 and X2, say X1(i) and X2(i), I want to compute a simple weighted sum:

+ *X1(i) + *X2(i)

Since the γ's are all scalars, I think I can do this multiplication and summation at the matrix level and get the result as a 10x15x10x10 matrix. However, when I change the order of summation, the results are also different. That is:

A = + *X1 + *X2;

B = + *X2 + *X1;

The differences in all entries of A and B sum up to 2.3856e-12. However, since this is just element-wise summation, A and B should be exactly the same.

Can someone explain what is going wrong?

I do notice that if I take out , I get the same answer... Why does adding a constant affect the final result?

Please see attached for the sample data and codes.

Thank you for your help!!

Answer 1

[Wan Ji]

It is a problem of machine error. If you have learned numerical analysis lessons, you will notice that big number will 'eat' small number when they are added.

a = 1;
b = eps;
for i = 1:10000
a = a+b;
end
c = a+10000*eps
c-a

Well the result is

ans =

     2.220446049250313e-12
     

When you follow this code, error will not occur!

A =  gamma0 + (gamma1*X1 + gamma2*X2); 
B =  gamma0 + (gamma2*X2 + gamma1*X1); 
A = A(:);
B = B(:);
sum(sum(B-A))

ans = 

0

This is beacuse the add operator order now is now the same for two expressions

[Tian]

This is super helpful!! Thank you so much for the explanation!!

[Wan Ji]

So, when you do numerical analysis, avoid adding a large number with a small number or dividing a large number by small number. The error will explode if you do so. Also, you cannot minus a number far too close to the minuend.

An example will also show why:

calculate where x=1000001, y = 1000000; (use single precision)

We use

a = sqrt(single(1000001)) - sqrt(single(1000000))

The answer is

  single

   4.8828125e-04

While with the following transformation

a = single(1)/(sqrt(single(1000001)) + sqrt(single(1000000)))

The answer is

a = single(1)/(sqrt(single(1000001)) + sqrt(single(1000000)))

a =

  single

   4.9999985e-04
   

The answer with double precision is

a =

     4.999998750000625e-04
     

You can see what I mean.

[Tian]

Thank you for the comments. I am very new to numerical analysis and this is great to know! My problem is to simulate results for a set of given coefficients. Some coefficients are large. To prevent errors from exploding:

(1) Are there any rule of thumbs to follow?

(2) How do we know whether the numbers are right or "exploded"?

I think I can make sure small numbers come first when I sum them...