Matrix exponential of hermitian matrix with eigenvectors from generalized eigenvalue problem

Question

I want to calculate the following expression $$ \exp(-i\Delta t\mathbf{H}) $$ where $\mathbf{H}\in\mathbb{C}^{n\times n}$ is a hermitian matrix. Since I have a highly optimized eigensolver in the code I'm working with (ScaLAPACK & Fortran 90) I was thinking about utilizing it by calculating the matrix exponential with the eigenvector ansatz $$ \exp(-i\Delta t\mathbf{H}) = \mathbf{C}\,\mathrm{diag}\Big[\exp(-i\Delta t \epsilon)\Big]\,\mathbf{C}^{-1} $$ where $\mathbf{C}\in\mathbb{C}^{n\times n}$ contains the eigenvectors of $\mathbf{H}$ and $\epsilon$ indicate its eigenvalues. In this case, solving the generalized eigenvalue problem $$ \mathbf{H}\mathbf{C} = \epsilon\mathbf{S}\mathbf{C} $$ yields the eigenvectors and -values. The generalized EVP is reduced to an ordinary EVP via a Cholesky transform of $\mathbf{S}=\mathbf{U}^T\mathbf{U}$, \begin{align} &~\mathbf{H}\mathbf{C}=\epsilon\mathbf{U}^T\mathbf{U}\mathbf{C}\\ \Leftrightarrow & ~\mathbf{U}^{-T}\mathbf{H}\mathbf{C} = \epsilon\mathbf{U}\mathbf{C}\\ \Leftrightarrow & ~\mathbf{U}^{-T}\mathbf{H}\mathbf{U}^{-1}\mathbf{U}\mathbf{C} = \epsilon\mathbf{U}\mathbf{C} \end{align} and by setting $\mathbf{H}'\equiv \mathbf{U}^{-T}\mathbf{H}\mathbf{U}^{-1}$ and $\mathbf{C}'\equiv \mathbf{U}\mathbf{C}$: $$ \mathbf{H}'\mathbf{C}'=\epsilon\mathbf{C}'. $$ The eigenvectors are backtransformed via $\mathbf{C}=\mathbf{U}^{-1}\mathbf{C}'$ when the eigenproblem is solved. Not sure about the general case, but here I have non-orthogonal eigenvectors so I can not make use of $$ \mathbf{C}^{-1} = \mathbf{C}^{T}. $$ Considering the paper of Moler and Van Loan (19 dubious ways ...), the eigenvector approach should always work out for a hermitian matrix if I understood it right. But I am not sure how to handle this situation. My next step would be the straightforward numerical calculation of $\mathbf{C}^{-1}$ but I am a bit afraid that I may miss a point here or that there is a better (i.e. cheaper) solution.

Answer 1

I guess your problem is the following. Your working in a non-orthogonal basis $\phi_i$, and solved the time-independent Schrödinger equation

$$\mathbf H \mathbf C = E \mathbf S \mathbf C$$

Here $\mathbf H$ is the Hamiltonian, and $\mathbf S$ is the (positive-definite) overlap matrix with $S_{ij} = \langle{\phi_i},\phi_j\rangle$.

Now you want to propagate the thus obtained initial bound state via the time-dependent Schrödinger equation, which is given by

$$ \mathbf S \, i\partial_t \mathbf C = \mathbf H \mathbf C$$

Note that the overlap matrix $\mathbf S$ is also present in the time-dependent equation. By employing the same Cholesky decomposition as in your OP, the previous equation becomes

$$ i\partial_t \mathbf C^\prime = \underbrace{(U^{-T}\mathbf H U^{-1})}_{=:\mathbf H^\prime}\mathbf C^\prime$$

The solution to this equation is

$$\mathbf C^\prime(t) = \exp(-i\mathbf H^\prime t) \, \mathbf C^\prime(0)$$

By solving the time-independent equation, you already have the eigenvalues and eigenvectors of $\mathbf H^\prime$. Use these to calculate the exponential (--I guess you know how), and finally use another backsubstitution step to obtain $\mathbf C(t)$ from $\mathbf C^\prime(t)$ as $\mathbf C(t) = \mathbf U^{-1} \mathbf C^\prime(t)$.

So, as in the time-independent problem, the clue is again not to work in the original, non-orthogonal basis, but in the Cholesky-transformed orthogonal basis.