52 views (last 30 days)

I made a function that calculates autocorrelation using partial matrices. When I tried to verify it by comparing with autocorrelation using whole matrices, I thought I had made a mistake, but then I found out Matlab gives different results depending on how the matrix is separated into parts. However, there should be no difference between the calculations!

Now, do not give me mere floating point arithmetic as the explanation: It does NOT explain this difference.

Here is a very simple test that shows the problem:

a1=randn(100,1); a2=randn(100,1);

A=[a1 a2];

R0=A'*A;

R1=[a1 a2]'*[a1 a2];

R2=[a1'*a1 a1'*a2; a2'*a1 a2'*a2];

[R0-R1 R2-R1]

ans =

1.0e-13 *

0.1421 -0.0089 0.2842 -0.0089

-0.0089 -0.2842 -0.0089 0

If you know anything about matrix algebra, you should know that R0 is internally calculated EXACTLY in the same way as R1 and R2 above. So, floating point arithmetic as such does not explain the difference.

So, obviously Matlab must be rounding the intermediate results differently in these cases, but I cannot understand WHY they are rounded differently? There is no sensible reason in these cases to do it differently. Mathworks, why you are doing this? It just makes all verification more difficult when the results are not what they are supposed to be.

James Tursa
on 21 Jan 2021

Edited: James Tursa
on 27 Jan 2021

This difference is this:

R0=A'*A;

In the above expression, you are using the exact same matrix A on both sides of the operation. The MATLAB parser can see this and thus call a symmetric BLAS matrix multiply routine, which fills in about 1/2 the matrix result and then MATLAB copies the results into the other half. This is faster than a generic matrix multiply and R0 is guaranteed to be exactly Hermitian/symmetric in this case.

R1=[a1 a2]'*[a1 a2];

In the above expression, you are building the two [a1 a2] matrices on the fly. The MATLAB parser is not smart enough to recognize that the two operands are in fact built the exact same way (it just sees two different operands with two different data pointers), so it calls a generic BLAS matrix multiply routine to get the R1 result. The operations are done in a different order than the symmetric routine, so in general you will get a slightly different result. And the result is not guaranteed to be exactly Hermitian/symmetric.

So here again, this is just the normal differences that would be expected when doing floating point arithmetic in a different order. Nothing else is going on, and MATLAB isn't trying to do anything to make your life more difficult. In fact the parser is trying to help you by giving you a way to speed up your calculations and guarantee an exact Hermitian/symmetric result. There are ways to force the calculation to use either the symmetric or generic routine in the background. My advice is simply to recognize what is going on and then write code accordingly. Assigning the operands to a single variable and then using that variable in a way that the parser will recognize the symmetry would probably be the best approach.

James Tursa
on 21 Jan 2021

John D'Errico
on 21 Jan 2021

The Mathworks is NOT rounding numbers just to make your life difficult. They are not rounding numbers at all, at least not beyond that which happens in normal floating point arithmetic. And apparently, you are telling us that we cannot tell you the answer, since you already "know" what it is.

What you do not seem to appreciate is that MATLAB uses the blas to do linear algebra. This makes your computations more efficient. At the same time, the blas can change the sequence of the computations involved. And as soon as you play with the sequence of the adds in a computation, you incur the possibility that you will see a problem.

So I'm sorry, but what you "know" about matrix agebra does not always apply. Mathematical proofs do not apply down in the least significant bits of floating point numbers.

John D'Errico
on 22 Jan 2021

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 1 Comment

## Direct link to this comment

https://es.mathworks.com/matlabcentral/answers/723553-matrix-multiplication-accuracy-problem#comment_1278553

⋮## Direct link to this comment

https://es.mathworks.com/matlabcentral/answers/723553-matrix-multiplication-accuracy-problem#comment_1278553

Sign in to comment.