Why indexing vs not indexing the 2D matrix will lead to different result?

I encounted an indexing problem: I have a 2x2 Cholesky factor and I want to recover the original SPD matrix.

load s.mat

s(:,:,1)*s(:,:,1)' - s*s'

and the s.mat file is attached in this question.

My question is, why s(:,:,1)*s(:,:,1)' and s*s' are different? I am working on the Bayesian optimization and the small different will be enlarged during the calculation.

Thanks in advance for the help!

[the cyclist]

I couldn't duplicate the behavior here (R2021b) ...

load s.mat
s(:,:,1)*s(:,:,1)' - s*s'

ans = 2×2

0 0 0 0

but I did on my machine (R2020a) ...

and I didn't even need to index into the 3rd dimension to do it.

[Rupeng Li]

Hi the cyclist, thanks for your quick reply.

Yes I found this too... not sure what is the reason...

s(:,:) == s

returns all ture, but

s(:,:)*s(:,:) == s*s

does not.

So I use s=rand([2,2]) and confirmed that the problem does not exist, which means the problem is in the matrix I uploaded. But what makes this matrix special?

By the way, the entries of the matrix are output from the L-BFGS-B wrapper in Matlab central [1], written in C or Fortran? So maybe there are rounding issues across platforms. Given such reason I still don't know s(:,:)*s(:,:) does no equals s .

[1] LBFGSB (L-BFGS-B) mex wrapper - File Exchange - MATLAB Central (mathworks.com)

Answer 1

[Steven Lord]

For a long time [Bobby Cheng and Peter Perkins call out "a few release" and "several releases" in a thread from 2007 on the old MATLAB newsgroup (Thread title "Outer prod" started by Marcus M. Edvall on Sep 12, 2007, 2:52:12 AM if that link doesn't work)] MATLAB has been able to detect operations of the form A*A' or A'*A and compute them in such a way that they are exactly symmetric.

But this optimization only works if MATLAB can detect that you're multiplying a matrix and the transpose of itself. In the case s*s' it can be certain that s and s are the same variable and so it can apply the optimization. In the case s(:, :)*s(:, :) it cannot be certain that s(:, :) and s(:, :) are the same variable (those are each the results of an indexing operation on a variable) so it cannot be certain that optimization is valid. Therefore it doesn't apply that optimization.

Sometimes the optimized and unoptimized answers will be the same. Other times they may differ slightly.

There are other older threads (Thread title "An interesting Matlab[sic] bug..." started by Jordan Rosenthal on June 14, 2006, 1:10:13 PM) and this other one (Thread title "Some unstable phenomenon of matrix operation, what is the reason?" started by Sui Huang on Aug 29, 2005, 1:47:58 PM) also discuss this phenomenon. But I don't remember offhand when this optimization was introduced.

[Rupeng Li]

Hi Steven, thanks for your reply. This is really a good answer. I looked through the email correspondence now I have an idea how to keep results consistent with the best effort. Indeed, I tested on 4000 samples of s and found that less than 20 have such issues, which added the difficulty of dubugging.

Answer 2

[James Tursa]

Some links related to this known behavior:

https://www.mathworks.com/matlabcentral/answers/599818-computational-complexity-of-matrix-multiplication?s_tid=srchtitle

https://www.mathworks.com/matlabcentral/answers/723553-matrix-multiplication-accuracy-problem?s_tid=srchtitle

https://www.mathworks.com/matlabcentral/answers/319587-is-it-possible-to-declare-the-matrix-is-symmetric-in-advance-before-multiplication?s_tid=srchtitle

https://www.mathworks.com/matlabcentral/answers/86390-matrix-transpose-multiplication-precision?s_tid=srchtitle

Basically, if you use the same variable in a matrix multiply operation (or at least variables with the same data pointers), MATLAB can recognize that they are the same and call symmetric BLAS library routines in the background. Otherwise generic BLAS library routines are called. Symmetric routines only do about 1/2 the work and guarantee symmetric results. Generic routines do not since they do the calculations in a different order and don't do the calculations for the lower vs upper triangle in the same order.

In the case of an indexed expression such as X(:,:,1) * X(:,:,1)', realize that each individual X(:,:,1) expression generates a separate deep copy of the data in a new memory location. I.e., the first X(:,:,1) data pointer is different from the second X(:,:,1) data pointer. When MATLAB then does the matrix multiply operation it sees the two different data pointers and can't recognize the symmetry. E.g.,

X(:,:,1) * X(:,:,1)' % causes two separate data copies, generic BLAS routine called, result not guaranteed symmetric

X1 = X(:,:,1); X1 * X1' % only one data copy involved, symmetric BLAS routine called, result guaranteed symmetric

All of the above discussion is based on MATLAB historical behavior. It is possible that in the future MATLAB may implement pointer-based sub-array capabilities where contiguous indexed expressions do not generate deep data copies, but I have yet to see any examples demonstrating that this behavior occurs for current versions. A mex routine that mimics this behavior can be found here, but this routine has not been updated for the latest MATLAB versions yet:

https://www.mathworks.com/matlabcentral/fileexchange/65842-sharedchild-creates-a-shared-data-copy-of-contiguous-subset