# Matrix exponential of hermitian matrix with eigenvectors from generalized eigenvalue problem

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.

• I lost you at the third equation. Why do you need to solve a generalized eigenvalue problem to compute the eigenvectors and eigenvalues of a Hermitian matrix? What's wrong with a normal eigensolver? Where did that S pop out of? Jun 15 '18 at 14:01
• The reason why a generalized EVP has to be solved is intrinsic to this problem, I should have said that. It has a physical background: basically, the Schrödinger eqn. is solved for the case that the ansatz wave function is represented via a non-orthogonal basis set $\{\phi_i\}$. This yields the overlap matrix $S_{ij}=\langle\psi_i|\psi_j\rangle$.
– Lukk
Jun 15 '18 at 14:28
• So what is the function that you really have to compute? It's not the one in the first equation, is it? Where does $S$ appear in it? Jun 15 '18 at 14:29

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.