# Smart way to multiply 3 matrices

I have a quantum mechanics simulation where I need to multiply three matrices that look like this:

$$\rho(t_1)=U^\dagger \rho(t_0) \, U$$

where $U^\dagger$ is the hermitian conjugate of $U$. This evolves the density matrix $\rho$ from one time point to another time point. $U$ is called the evolution operator.

My question is: Is there a smart way to reduce the time required to do this matrix multiplication in C++? Is there anything better than using BLAS's zgemm 2 times (or zhemm 2 times, since the matrices are Hermitian)?

What I mean with smart is not only using a different library, but also finding a way to mathematically or computationally reduce the number of multiplications that have to be done in that operation.

The sizes of the matrices I deal with range from $2^7$ to $2^{15}$ in side-length (All these matrices are square matrices).

• I know this doesn't directly answer you, but it's often possible to avoid storing unitary matrices directly. Unitary matrices (e.g., netlib.org/lapack/lug/node128.html) are more often stored internally as products of elementary unitary matrices, which can help. – Kirill Oct 21 '16 at 19:26
• Another comment that doesn't directly answer you, but sometimes only the action of matrices are needed on some known vectors, rather than the matrices themselves. If you can't get away with DGEMVs, you might be stuck with DGEMMs, or you could consider implementing a BLAS-like operation that does this sort of thing in one step and see if you can get any speedup. – Geoff Oxberry Oct 21 '16 at 21:33
• Out of curiosity, how do you construct U from the underlying Hamiltonian? (diagonalize H...?) – roygvib Oct 21 '16 at 22:02
• @roygvib Matrix exponential of the Hamiltonian (with some factor). – The Quantum Physicist Oct 21 '16 at 22:05
• I initially thought of split-operator things, but it may be not effective for general discrete models... Another thought is that if rho(t=0) can be represented somehow as \sum_k |k><k| with a relatively few states {|k>}, then the calculation may be reduced to matrix-vector things... – roygvib Oct 21 '16 at 22:13

Have you considered working with the Cholesky (or low-rank) factor of $\rho(t_0)$ rather than with the matrix itself? This might reduce the number of products that you need to make, and it has the additional benefit of preserving positive semidefiniteness across your computations.

It is already cheaper to use Cholesky factors for this computation alone, if I am not making mistakes with the costs: let $n$ be the side of each matrix appearing here; computing the upper triangular Cholesky factor $R$ such that $\rho(t_0)=R^\dagger R$ costs $\frac13 n^3$ (counting with the traditional model addition=multiplication=1); computing $RU$ costs $n^3$ (because $R$ is triangular), and computing $(RU)^\dagger (RU)$ costs $n^3$ (because you only need to compute half of the entries). OTOH, computing $U^\dagger \rho(t_0)$ costs $2n^3$, and then computing half of the entries of $(U^\dagger \rho(t_0))U$ costs $n^3$.

If $\rho(t_0)$ has low rank (which in the examples I saw when I studied quantum mechanics happened quite often for initial states), then one can start from a rectangular $R$ coming from the QR of its low-rank factor rather than a Cholesky decomposition, and this method is even cheaper. If you have to make more computations with your matrices after this triple product, then there may be additional savings: for instance, to compute products and solve linear system you can work with the low-rank factors directly rather than forming the last product.

• Thanks. I wonder if MATLAB (another language known for fast matrix operations) somehow does this in its own matrix multiplication scheme, or do I have to implement this myself? – Ka-Wa Yip Mar 7 '17 at 2:20
• @kyle No, as far as I know Matlab does not do this. The only similar-looking optimization that it does is computing $A^*A$ by computing only its lower or upper triangular -- but slightly more complex tasks are not similarly optimized. – Federico Poloni Mar 7 '17 at 21:58

Considering only a computational point of view you could think of using cuBLAS the CUDA version of BLAS.

At the link is an example of use. In the last part there is a example in C++, with thrust.

• I'm sorry. This is not a valid solution for me. I'm looking for a solution without CUDA or PHI. – The Quantum Physicist Oct 21 '16 at 16:48
• Ok. I ask here, the following thing, because I can not comment on your question. Does your BLAS come from Atlas ? – Mauro Vanzetto Oct 21 '16 at 17:56
• My BLAS is OpenBLAS. – The Quantum Physicist Oct 21 '16 at 18:01
• Are you ordering the two matrix multiplies in the more efficient way? en.wikipedia.org/wiki/Matrix_chain_multiplication – Bill Greene Oct 21 '16 at 19:03
• @TheQuantumPhysicist In openBLAS site there is a page with extension. In particular there is a function gemm3m, from Intel description: The ?gemm3m routines perform a matrix-matrix operation with general complex matrices. These routines are similar to the ?gemm routines, but they use fewer matrix multiplication operations In the section Note application they write resulting in significant savings in compute time for large matrices. Have you seen this function? – Mauro Vanzetto Oct 22 '16 at 13:15