5
$\begingroup$

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)?

Additional information:

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).

Improve this question
$\endgroup$
17
  • 1
    $\begingroup$ 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. $\endgroup$
    – Kirill
    Commented Oct 21, 2016 at 19:26
  • 1
    $\begingroup$ 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. $\endgroup$
    – Geoff Oxberry
    Commented Oct 21, 2016 at 21:33
  • $\begingroup$ Out of curiosity, how do you construct U from the underlying Hamiltonian? (diagonalize H...?) $\endgroup$
    – roygvib
    Commented Oct 21, 2016 at 22:02
  • $\begingroup$ @roygvib Matrix exponential of the Hamiltonian (with some factor). $\endgroup$
    – The Quantum Physicist
    Commented Oct 21, 2016 at 22:05
  • $\begingroup$ 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... $\endgroup$
    – roygvib
    Commented Oct 21, 2016 at 22:13

2 Answers 2

Reset to default
3
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – Ka Wa Yip
    Commented Mar 7, 2017 at 2:20
  • $\begingroup$ @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. $\endgroup$
    – Federico Poloni
    Commented Mar 7, 2017 at 21:58
0
$\begingroup$

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.

Improve this answer
$\endgroup$
17
  • $\begingroup$ I'm sorry. This is not a valid solution for me. I'm looking for a solution without CUDA or PHI. $\endgroup$
    – The Quantum Physicist
    Commented Oct 21, 2016 at 16:48
  • $\begingroup$ Ok. I ask here, the following thing, because I can not comment on your question. Does your BLAS come from Atlas ? $\endgroup$
    – Mauro Vanzetto
    Commented Oct 21, 2016 at 17:56
  • $\begingroup$ My BLAS is OpenBLAS. $\endgroup$
    – The Quantum Physicist
    Commented Oct 21, 2016 at 18:01
  • 2
    $\begingroup$ Are you ordering the two matrix multiplies in the more efficient way? en.wikipedia.org/wiki/Matrix_chain_multiplication $\endgroup$
    – Bill Greene
    Commented Oct 21, 2016 at 19:03
  • 1
    $\begingroup$ @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? $\endgroup$
    – Mauro Vanzetto
    Commented Oct 22, 2016 at 13:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.