I have observed an apparently very strange behavior when running a Matlab script. I have a very long data vector t (with 4560000 elements) and several other vectors x1,...,xM, all with the same dimension as t. Now, I basically want to compute the scalar product between each xm and t. I have done that using two procedures: the first one consists in computing each scalar product between a vector xm and t inside a for loop, while the other simply consists in forming a matrix X whose columns are the vectors xm and then computing the product of X transposed and t. Surprisingly, the results are very different!
For simplicity, let us suppose M=2. If I compute:
>> y1 = [x1 x2].'*t;
>> y2 = zeros(2,1);
>> y2(1) = x1.'*t;
>> y2(2) = x2.'*t;
Then I get a large relative normalized error
>> er = norm(y2-y1)/norm(y1)
(By the way, norm(y1) = 1.2665e+06). Componentwise, the error is
In principle, these alternatives should always yield the same results, regardless of the contents of those vectors, since they should be translated into machine code corresponding to equivalent operations. So, it seems odd to me that some kind of numerical issue arises.
Has anyone seen a similar behavior before and knows what might be causing that? Am I missing something with respect to the internal implementation of such operations?