Theory

Time Evolution using Krylov Methods

In this section, we will discuss some of the mathematics behind time-evolution using Krylov methods.

First, a refresher on a Krylov subspace. Given an initial state $v$, and matrix $A$, the $n^{\text{th}}$ order Krylov subspace is given by: $$\begin{array}{c}{\mathcal{K}}_n\equiv \text{span}\left\{v,Av,A^{2}v,...,A^{n-1}v\right\}\end{array}$$

where span refers to the space spanned by the vectors. In practice, this list of vectors is converted into an orthonormal basis by (modified-)Gram-Schmidt methods, Givens rotations or Householder reflections.

For the evaluation of the time-propagator, we seek approximations of the form:

$$\begin{array}{c}{e}^{A}v\approx p_{n-1}(A)v\end{array}$$

that is, that we can expand the effect of our evolution operator as a polynomial of order $n-1$ (in $A$). To do so, first we define our orthonormal basis using the Arnoldi method:

compute $v_0:=v/||v|{|}_2$
for $k=1,2,\dots ,n-1$:
    compute $v_k:=A{v}_{k-1}$
    for $j$ in $0,\dots ,k-1$:
        $h_{j,k-1}:=v_j^{†}{v}_k$
        $v_k:=v_k-h_{j,k-1}$
    $h_{k,k-1}:=||v_k|{|}_2$
    $v_k:=v_k/h_{k,k-1}$

then our basis $V_n$ is comprised of the vectors $v_k$ and the subscript $n$ refers to the fact that this is the basis formed from the $n^{\text{th}}$ order Krylov subspace.

The Arnoldi method has the useful property, that upon applying the basis vectors to the original matrix, the resulting smaller matrix $H_n\in {\mathbb{C}}^{n\times n}$, $$\begin{array}{c}{V}_n^{†}A{V}_n=H_n,\end{array}$$ is tridiagonal when $A$ is Hermitian, and is upper-Hessenberg otherwise. Both are more easily exponentiated, but the tridiagonal case especially so.

Now that we have our basis, we seek to find an approximate solution $e^{A}v\approx V_n{y}_0\in {\mathcal{K}}_n$. More precisely, we are looking for a solution of the minimization problem $$\begin{array}{c}{y}_0=\underset{y}{\mathrm{argmin}}{\Vert e^{A}v-V_{n}y\Vert }_2.\end{array}$$

Since the columns of $V_n$ are orthonormal, we can split the identity into the projection onto the range of $V$ and its orthogonal complement via $I=V_n{V}_n^{†}+(I-V_n{V}_n^{†})$. Thus,

$$\begin{array}{c}{\Vert e^{A}v-V_{n}y\Vert }_2={\Vert (V_n{V}_n^{†}+(I-V_n{V}_n^{†}))(e^{A}v-V_{n}y)\Vert }_2\end{array}$$

Expanding, using the pythagorean theorem, and simplifying gives

$$\begin{array}{c}{\Vert e^{A}v-V_{n}y\Vert }_2={\Vert V_n{V}_n^{†}{e}^{A}v-V_{n}y{\Vert }_2+\Vert e^{A}v-V_n{V}_n^{†}{e}^{A}v\Vert }_2.\end{array}$$

Note, that the second norm on the right does not include the variable $y$, while choosing

$$\begin{array}{c}y=V_n^{†}{e}^{A}v,\end{array}$$

makes the first norm on the right equal to zero. Since norms are non-negative, this expression is, therefore, the desired minimizing argument $y_0$.

This equation, however, still involves the original matrix exponential. By noting in the orthogonalization process of the Krylov subspace, that $$\begin{array}{c}v\equiv V_n{e}_0\end{array}$$ ($e_0$ being the unit vector at index $0$), we rewrite our solution as $$y_0=V_n^{†}{e}^A{V}_n{e}_0.$$

Hence

$$\begin{array}{c}{y}_0=V_n^{†}{e}^A{V}_n{e}_0=V_n^{†}\sum _{k=0}^{\infty }\frac{A^k}{k!}{V}_n{e}_0=\sum _{k=0}^{\infty }\frac{V_n^{†}{A}^k{V}_n}{k!}{e}_0.\end{array}$$

We would like to move the matrice $V_n$ into the exponentiation, because then we would be able to compute $y_0$ by

$$\begin{array}{c}\sum _{k=0}^{\infty }\frac{(V_n^{†}A{V}_n{)}^k}{k!}{e}_0=\sum _{k=0}^{\infty }\frac{H_n^k}{k!}{e}_0=e^{H_n}{e}_0,\end{array}$$

which only involves the small matrix $H_n$. However,

$$\begin{array}{c}(V_n^{†}A{V}_n{)}^k=V_n^{†}A{V}_n{V}_n^{†}A{V}_n\cdots V_n^{†}A{V}_n\ne V_n^{†}(A{)}^k{V}_n,\end{array}$$ since $V_n{V}_n^{†}\ne I$. Let us assume for a moment that the range of $V_n$ is invariant under the operation of $A$, i.e., $Av\in \mathcal{R}(V_n)$ for $v\in \mathcal{R}(V_n)$. Then, $V_n{V}_n^{†}A{V}_n=A{V}_n$, since $V_n{V}_n^{†}$ is the orthogonal projector onto $\mathcal{R}(V_n)$. Thus, under this assumption the desired equation does hold, and if we assume that the range of $V_n$ is almost invariant under the action of $A$ then the equation is still a good approximation.

$$\begin{array}{c}{e}^{A}v\approx V_n{y}_0\approx V_n{e}^H{e}_0.\end{array}$$

As a final point, if we want to compute the time evolution of a quantum system defined by a Hamiltonian $A$, we want to compute $e^{-\mathrm{i}A\tau }v$ for some value of $\tau \in \mathbb{R}$. We note that as the decomposition is invariant with respect to (non-zero) scalar multiples, we may modify $A\to -\mathrm{i}A\tau$ to give $$\begin{array}{c}{e}^{-\mathrm{i}A\tau }v\approx V_n{e}^{-\mathrm{i}{H}_n\tau }{e}_0.\end{array}$$

For more details on related methods see the following literature:

E. Gallopoulos, Y. Saad "On the parallel solution of parabolic equations." Proceedings of the 3rd international conference on Supercomputing, 1989 17-28. doi: 10.1145/318789.318793

E. Gallopoulos, Y. Saad "Efficient solution of parabolic equations by Krylov approximation methods." SIAM journal on scientific and statistical computing 13(5) 1992: 1236-1264. doi: 10.1137/0913071

Y. Saad, "Analysis of some Krylov subspace approximations to the matrix exponential operator." SIAM Journal on Numerical Analysis 29(1) 1992: 209-228. doi: 10.1137/0729014