Unexpected numerical errors in matrix/vector multiplication

Hi,

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)
er =
   6.3080e-05

(By the way, norm(y1) = 1.2665e+06). Componentwise, the error is

>> y1(1)-y2(1)
ans = 
   1.5918
>> y1(2)-y2(2)
ans =
   -0.0898

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?

Thanks,

Henrique

[José Goulart]

It seems that this behavior somehow holds also for random vectors of "sufficiently large" size. For instance, if I generate random vectors t, x1 and x2 of dimension 8, no error arises, while for vectors of dimension 16 there is a small relative error (which is close to machine precision, though). Thus, it seems that the contents of my data vector t are somehow contributing for an aggravation of this behavior, since the relative error I get is much larger than eps.

[dpb]

I can't reproduce this here w/ R2012b...

>> t=rand(1000,1);
>> x=rand(2,1000);
>> z1=x*t;
>> z2=zeros(2,1);
>> z2(1)=x(1,:)*t;
>> z2(2)=x(2,:)*t;
>> z1-z2
ans =
   0
   0
>> all(z1==z2)
ans =
   1
>>

[José Goulart]

Thank you for your comment, dpb. I have run the same example code you have written, and here's my output:

>> t=rand(1000,1);
>> x=rand(2,1000);
>> z1=x*t;
>> z2=zeros(2,1);
>> z2(1)=x(1,:)*t;
>> z2(2)=x(2,:)*t;
>> z1-z2
ans =
   1.0e-12 *
   -0.2274
    0.1137
>> all(z1==z2)
ans =
     0

Thus, this seems to depend on the platform version. Mine is 8.1.0.604 (R2013a).

[Adam]

I also get a difference in R2014a. Indeed if I use rather larger numbers e.g.

t=rand(1000,1) * 1e20;
x=rand(2,1000) * 1e20;

then I end up with:

z1 - z2

ans =

        2.78536508839211e+27
        4.95176015714152e+27

[John D'Errico]

Adam - you do realize that difference IS machine precision?

Take the product of numbers on the order of 1e20. Add them in different order (because this is what the blas does, it controls the sequence of adds in big multiplies for efficiency.)

The error should be on the order of...

eps(1e40*1000)
ans =
   1.2379e+27

So completely expected.

[arich82]

I get something really, really weird.

If I start a new session of Matlab, the first time I run the code, I get no difference, but the second (and subsequent) time I run the code, I get non-trivial errors:

>> t=rand(1000,1);
x=rand(2,1000);
z1=x*t;
z2=zeros(2,1);
z2(1)=x(1,:)*t;
z2(2)=x(2,:)*t;
z1-z2

ans =

       0
       0

>> t=rand(1000,1);
x=rand(2,1000);
z1=x*t;
z2=zeros(2,1);
z2(1)=x(1,:)*t;
z2(2)=x(2,:)*t;
z1-z2

ans =

     1.0e-12 *

     -0.1990
      0.0568

If I restart Matlab, same thing happens. This is in 8.3.0.532 (R2014a) running in Ubuntu 14.04, with dual Intel E5-2690 v2.

Answer 1

[Oleg Komarov]

The dot product implementation might differ according to the input, i.e. dot product between two vectors as opposed between two matrices. Also, a quick search on the web exposes a few implementations of the dot product in the BLAS routines themselves and Matlab might be picking one or another according to different needs.

The example suggests zdotu.f

edit([matlabroot '/extern/examples/refbook/dotProductComplex.c']);

However, this doesn't tell us much about the dot product between matrices.

The difference in the results after startup calculation seems to be just chance. In fact, if you vary the number of random draws, the difference is not 0 anymore. To reset the state of the random generator without restarting, simply use

rng('default');

To address such a big discrepancy, we need to see what are the values of the dot products.

[José Goulart]

Indeed, that must be the cause.

Also, since the elements of the vectors span several orders of magnitude, adding them in different orders should produce fairly different results.

Thanks, everyone.