4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Federico Poloni Jun 15 '18 at 14:01
  • $\begingroup$ 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$. $\endgroup$ – Lukk Jun 15 '18 at 14:28
  • 2
    $\begingroup$ 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? $\endgroup$ – Federico Poloni Jun 15 '18 at 14:29
3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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