Matrix exponential by eigenvectors - implementation issues

Question

I posted a similar question yesterday but I deleted it since I think that I had to reformulate it after some insights.

I want to calculate $$ \exp(-i\Delta t\,\mathcal{H}) = V\,\mathrm{diag}(\{\exp(-i\Delta t\,\epsilon_k\})\,V^{-1} ~~~~~~~~~~~~~(*) $$ where $\mathcal{H}$ is a hermitian matrix (the hamilton matrix) and $V$, $\epsilon_k$ are its eigenvectors and -values.

---

Since I'm implementing this in a given DFT package, I make use of the existing eigensolver routine. All numerics are done by ScaLAPACK. It first computes the Cholesky factorization of the overlap matrix $S_{ij}=\langle\phi_i|\phi_j\rangle$: $$ S = U^TU $$ where $U$ is the upper triangular matrix. In the next step, $U$ is inverted for further usage. It then solves the eigenvalue problem for the transformed hamilton matrix $$ \mathcal{H}' = U^{-\dagger}\mathcal{H}\,U^{-1} $$ which yields the eigenvectors $V'$. This is done since we have a generalized eigenvalue problem because a non-orthogonal basis set is used. The eigenvectors are afterwards backtransformed via $$ V = U^{-1} V'. $$

I tried to calculate $(*)$ by inverting $V$ which needs a LU decomposition since $V$ is a non-hermitian matrix. This turned out to be unstable after some iterations although it initially produced correct results.

I found a different way to circumvent inverting $V$ which uses the not backtransformed eigenvectors $V'$. Here, the exponential is explicitly calculated via $$ (U^{-1}V')\,\mathrm{diag}(\{\exp(-i\Delta t\,\epsilon_k\})(U^{-1}V')^{-1}=U^{-1}V'\,\mathrm{diag}(\{\exp(-i\Delta t\,\epsilon_k\})V'^\dagger U $$ where $V^{-1}=V'^{\dagger}$ can be applied since $V'$ is unitary. Unfortunately, this yields oscillating solutions around the right solution (e.g. total energy).

I'm not sure why this isn't working and I'm quite new to the field. I tried $S^{-1}$ and $S$ in above formula instead of $U^{-1}$ and $U$ but it didn't help, also I tried a lot more modifications but I just can't figure out where the problem could be.

I would really appreciate some help. Thanks and best regards!

Answer 1

(Heavily edited because I made a big mistake in the original post)

We write the Time Dependent Kohn-Sham equations as:

$Hc = i\hbar S\frac{\delta c}{\delta t} $

where $c$ are the coefficients of your wavefunction. This is also the same as:

$S^{-1}Hc = i\hbar\frac{\delta c}{\delta t} $

and so we can update the wavefunctions as:

$c(t+\Delta t) = exp(-i\hbar S^{-1}H\Delta t)c(t)$.

So now we want to compute $exp(-i\hbar S^{-1}H\Delta t)$, through the eigendecomposition as you said.

$S^{-1}H = X^{-1}DX$

Let's write that in terms of a familiar problem:

$S^{-1}HX = DX$

$HX = DSX$

Using the cholesky decomposition, as you suggested, we can turn this into the standard eigenvalue problem.

$S = LL^{T}$

$(L^{-1}HL^{-T})(L^{T}X) = D(L^{T}X)$

$L^T X = X'$

$L^{-1}HL^{-T} = H'$

and then just solve

$H'X' = DX'$.

Now let's work our way back.

$X = L^{-T}X'$

$X^{-1} = X'^{-1}L^{T} = L^{T}X'^{T}$

(since $X'$ is an orthogonal matrix).

So now we know D, X, and $X^{-1}$, and thus can compute:

$exp(-i\hbar S^{-1}H\Delta t) = exp(-i\hbar X^{-1}DX\Delta t) = X^{-1}exp(-i\hbar D\Delta t)X$.

Let's make sure I'm correct with some source code:

import scipy.linalg
import numpy.random
import numpy

# Setup
N = 10
H = 0.01*numpy.random.rand(N,N)
H = H+H.T # symmetric
S = numpy.random.rand(N,N)
S = S+S.T # symmetric
S = S.dot(S) # positive definite
Sinv = scipy.linalg.inv(S)
I = numpy.identity(N)

# Reference answer
ref = scipy.linalg.expm(Sinv.dot(H))

# Our solution
## first our cholesky factors
L = scipy.linalg.cholesky(S,lower=True)
Linv = scipy.linalg.solve_triangular(L, I, lower=True)

## now solve the normal eigenvalue problem
Hprime = Linv.dot(H).dot(Linv.T)
D, Xprime = scipy.linalg.eigh(Hprime)
D = numpy.diag(D)

## back subsite for X and X^{-1}
X = Linv.T.dot(Xprime)
Xinv = Xprime.T.dot(L.T)

our_answer = X.dot(scipy.linalg.expm(D)).dot(Xinv)

print(scipy.linalg.norm(ref - our_answer))