error in big matrix multiplication

I have two matrices as Q a 4-by-4 matrix and x as a 4-by-10000 matrix . When I compute the multiplication as

Ans1 = Q * x;

the answer is different than having a for loop like

for i = 1 :10000
  Ans2(:,i) = Q * x(:,i);
end

The order of error (Ans1 - Ans2) is 1e-16 for my special case of Q and x. I know it is small but I was wondering what is the source of this error and how I can eliminate it.

Answer 1

[Walter Roberson]

The 4*10000 case would be delegated to the highly optimized multi-threaded libraries such as LINPACK or BLAS. The additions could end up being done in a different order, leading to a slightly different round-off result. Also, potentially the highly optimized libraries could use an instruction such as Fused Multiply And Add, which does one fewer rounding operations, potentially leading to a more accurate -- but different -- result.

[TOSA2016]

Is using GPU makes this worse due to data structure it is using?

[Walter Roberson]

What I described above does not use GPU.

If you are using gpuarray() or calling into gpu kernels you have built, then the result can certainly differ compared to serial operations.

Serial operations are not necessarily any more accurate than parallel operations. Even just for addition, you need special computation routines to achieve the effect of the rounding coming out as if the precision of the operands was indefinitely long. Knuth has a routine in The Art Of Computer Programming; unfortunately I never did understand how that routine worked.