I am trying to calculate the following equation efficiently:
$\rho(t) =exp(-iHt/\hbar)\rho(0)exp(+iHt/\hbar)$
where $H$ and $\rho$ are some matrices of varying size (from 2x2 up to 128x128). It isn't terribly slow to "brute force it" which would be to just explicitly calculate it (which is easily done via MATLAB), and it takes something on the order of 1.5 ms when H and p are of size 128 by 128. However I am computing this on the order of millions of times, thus the computation can take days or weeks. I understand the most obvious way to improve this would be to optimize my code, however I'm reasonably confident I've done this and have parallelized it accordingly, thus the only thing left is to make the bulk of the computation more efficient. I am not a mathematician, thus I was wondering if anyone had any insight into how to more efficiently perform this calculation it would be greatly appreciated. My intuition says that it may be possible to replace the two matrix multiplications with a single matrix multiplication, resulting in a two-fold speedup, but I'm not sure. Any help would be greatly, appreciated.