Reputation: 11
I encountered some numerical questions when running simulation on MatLab. Here please find the questions:
I found that A*A'
(a matrix times its transpose) is not guaranteed to be symmetric in MatLab. Can I know what is the reason? And because I will have A*C*A'
, where C
is a symmetric matrix, and I would like to keep A*C*A'
as symmetric. Is there any method to fix the numerical difference created by the transpose operation?
I implemented a for loop in Matlab to compute a set of matrices. Small numerical difference (around 10^(-10)
) in each round accumulates to the next run, and it finally diverges after around 30 rounds. Is there any method to fix small error in each run and do not affect the result at the same time.
Thank you for reading my questions!
Upvotes: 0
Views: 74
Reputation: 2636
"I found that A*A' (a matrix times its transpose) is not guaranteed to be symmetric in MatLab."
I would dispute that statement as written. The MATLAB parser is smart enough to recognize that the operands of A*A' are the same and call a symmetric BLAS routine in the background to do the work, and then manually copy one triangle into the other resulting in an exactly symmetric result. Where one usually gets into trouble is by coding something that the parser cannot recognize. E.g.,
A = whatever;
B = whatever;
X = A + B;
(A+B) * (A+B)' <-- MATLAB parser will call generic BLAS routine
X * X' <-- MATLAB parser will call symmetric BLAS routine
In the first matrix multiply above, the MATLAB parser may not be not smart enough to recognize the symmetry so a generic matrix multiply BLAS routine (e.g., dgemm) could be called to do the work and the result is not guaranteed to be exactly symmetric. But in the second matrix multiply above the MATLAB parser does recognize the symmetry and calls a symmetric BLAS matrix multiply routine.
For the ACA' case, I don't know of any method to force MATLAB to generate an exact symmetric result. You could manually copy one resulting triangle into the other after the fact. I suppose you could also factor C into two parts X*X' and then regroup but that seems like too much work for what you are trying to do.
Upvotes: 1