What are the efficient options to perform matrix multiplications of the form A*B*A and A*B*A' ? I am particularly intereted in the case when A and B are Symmetric or Hermitian and when these multiplications are performed repeatedly inside a loop, such as
sum(A*B*A for B in Bset),
where Bset is an array of matrices.
Some time ago i found a post on this topic, but i am not able to retrive it anymore.
If your matrices are positive definite as well you should have a look at https://github.com/JuliaStats/PDMats.jl
See also the discussion in
https://github.com/JuliaLang/julia/issues/16573
We’ve long wanted a specialized dot(x,M,y) function, but no one has gotten around to implementing it.
Wouldn’t that be the same as A*sum(Bset)*A? But I guess your real expression doesn’t allow such simplification…
No, in fact i rather need to compute expressions such as
sum(log(A*B*A) for B in Bset)
Unfortunately PDMats concerns real positive matrices only. I am interested in either real or complex positive matrices. As per the A*B*A' expression, PDmats decomposes the expression as A*L*L'A' using the Cholesky decomposition B=L*L' and uses low-level function rmul! in LinearAlgebra. I have tried and i don’t get a significant difference in computation time and exactly the same memory allocation.