Why indexing vs not indexing the 2D matrix will lead to different result?
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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!
2 commentaires
the cyclist
le 4 Oct 2021
Modifié(e) : the cyclist
le 4 Oct 2021
I couldn't duplicate the behavior here (R2021b) ...
load s.mat
s(:,:,1)*s(:,:,1)' - s*s'
but I did on my machine (R2020a) ...
and I didn't even need to index into the 3rd dimension to do it.
Réponse acceptée
Steven Lord
le 4 Oct 2021
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.
Plus de réponses (1)
James Tursa
le 6 Oct 2021
Modifié(e) : James Tursa
le 8 Oct 2021
Some links related to this known behavior:
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:
0 commentaires
Voir également
Catégories
En savoir plus sur Matrix Indexing dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!