Krylov Time Evolution

Krylov time evolution approximates the action of a large matrix exponential on one state by exponentiating a small projected matrix. It can preserve the norm of the projected trajectory very accurately while still missing dynamics outside the search space. An error or defect estimate is therefore essential; normalization is not an adequate stopping rule.

Helpful background. Lanczos provides the Hermitian recurrence, and matrix exponentials provides the small-matrix calculation.

Projected exponential action

For fixed Hermitian $H$ and normalized initial state $v$, build $\mathcal K_m(H,v)=\operatorname{span}\{v,Hv,\ldots,H^{m-1}v\}$. Let $Q_m$ be an orthonormal basis with first column $v$ and $T_m=Q_m^\dagger HQ_m$. The approximation is

$$\psi_m(t)=Q_m e^{-iT_mt/\hbar}e_1.$$

It reproduces the action of powers through the represented Krylov directions. Lanczos gives a tridiagonal $T_m$; Arnoldi gives an upper-Hessenberg projection for a general non-Hermitian generator.

Because $T_m$ is Hermitian, $\|\psi_m(t)\|=1$ in exact arithmetic. At $m=1$, however, the method only multiplies the initial vector by the phase of its Rayleigh quotient. This may be very inaccurate despite exact norm conservation.

The defect and a useful bound

Using the Lanczos relation $HQ_m=Q_mT_m+\beta_mq_{m+1}e_m^T$, define the differential defect

$$d_m(t)=\dot\psi_m(t)+\frac{i}{\hbar}H\psi_m(t).$$

The projected terms cancel, leaving

$$d_m(t)=\frac{i\beta_m}{\hbar}q_{m+1} e_m^T e^{-iT_mt/\hbar}e_1.$$

Variation of constants and unitarity of the exact evolution give

$$\|\psi(t)-\psi_m(t)\| \leq\int_0^t\|d_m(s)\|\,ds.$$

The integral involves only a small projected exponential and the final recurrence coefficient. A single endpoint defect can accidentally vanish and is not, by itself, this integral bound. Practical estimators approximate the integral or compare enlarged spaces; document which criterion is used.

Loss of orthogonality degrades both the projected norm argument and the residual formula. Reorthogonalize when needed and occasionally compare against a larger Krylov dimension or independent propagator.

Restarts and time-dependent Hamiltonians

For a long interval, divide time into steps and rebuild the Krylov space from the new state. Increase $m$ or shorten the step when the estimated local error is too large. For exact unitary step operators, global state error is bounded by the sum of local errors through a telescoping argument; allocating a total tolerance across steps is therefore more meaningful than fixing an unexplained local threshold.

A large spectral width generally demands more Krylov vectors for a fixed step. A scalar shift $H\mapsto H-cI$ leaves $\mathcal K_m$ unchanged in exact arithmetic and replaces $T_m$ by $T_m-cI$. It can help the small-matrix exponential calculation, but it does not reduce the intrinsic Krylov projection error: after restoring the global phase, the projected trajectory is the same.

For $H(t)$, a Krylov action of the midpoint-frozen Hamiltonian solves only that frozen problem. Tightening the Krylov tolerance does not remove the time-ordering error. Use time-step refinement or a Magnus construction outside the exponential-action solver.

A two-level exactness check

Take $H=\hbar\Omega\sigma_x/2$ and $v=|0\rangle$. A one-vector space has $T_1=0$, so $\psi_1(t)=|0\rangle$ for all times, while the exact transition probability is $\sin^2(\Omega t/2)$. Its defect norm is $|\Omega|/2$.

With two independent vectors, the Krylov space is the full state space, so the projected propagation is exact up to floating-point error. This tests construction, projection, and reconstruction with no ambiguity about representation error.

A vanishing endpoint defect with a large state error

Use dimensionless energy and time units with $\hbar=1$. In the three-level Lanczos example, $H=\operatorname{diag}(0,1,2)$ and $v=(1,1,1)^T/\sqrt3$. The two-vector projection is $T_2=I+b\sigma_x$, where $b=\sqrt{2/3}$, and $\beta_2=1/\sqrt3$. Thus

$$e_2^T e^{-iT_2t}e_1=-i e^{-it}\sin(bt),\qquad \|d_2(t)\|=\frac{|\sin(bt)|}{\sqrt3}.$$

At $t_*=\pi/b$, this defect vanishes exactly. Nevertheless,

$$\psi_2(t_*)=-e^{-it_*}v,\qquad \psi(t_*)=\frac1{\sqrt3} \begin{pmatrix}1\\e^{-it_*}\\e^{-2it_*}\end{pmatrix}.$$

Their state error is

$$\|\psi(t_*)-\psi_2(t_*)\| =\sqrt{\frac83+\frac43\cos t_*}\simeq1.28534.$$

Both states have unit norm. Their overlap is $-(1+2\cos t_*)/3\simeq0.17395$, so a global phase adjustment cannot fix the discrepancy. By contrast, the integrated defect bound is $\sqrt2$ and correctly rejects this as an accurate step. An endpoint estimate such as $t_*\|d_2(t_*)\|$ would report zero.

This deliberately oversized step illustrates why a local asymptotic estimator must not be treated as an unconditional bound. Increase the Krylov dimension, shorten the interval, or control the defect integral including its quadrature error. Jawecki analyzes the distinction between integral bounds and quadrature estimates when the defect oscillates. Here a third independent Krylov vector spans the full space and removes the projection error.

Exercises

Why can a vanishing Lanczos coefficient terminate propagation exactly for the current initial state even if the full Hilbert space is much larger?

Solution

Exact breakdown means the current Krylov space is invariant under $H$. Every power of $H$ acting on the initial vector remains in it, and so does the exponential series. The full Hilbert space can contain other invariant sectors that the initial state never reaches.

Evaluate the defect integral for the three-level example up to $t_*=\pi/b$. Why is sampling only the two endpoints insufficient even though the defect is smooth?

Solution

On this interval $\sin(bs)\geq0$, so

$$\int_0^{t_*}\|d_2(s)\|\,ds =\frac1{\sqrt3}\int_0^{\pi/b}\sin(bs)\,ds =\frac{2}{\sqrt3 b}=\sqrt2.$$

Both endpoint values are zero, while the defect is positive in the interior. A single endpoint or endpoint-only trapezoidal estimate misses its entire contribution. Smoothness alone does not certify a quadrature at an unresolved scale; additional samples and error control are required.

References

• Al-Mohy, A. H., and N. J. Higham, “Computing the action of the matrix exponential, with an application to exponential integrators,” SIAM Journal on Scientific Computing 33, 488–511 (2011), author-hosted record — an alternative action method and comparison framework.
• Jawecki, T., “A study of defect-based error estimates for the Krylov approximation of φ-functions,” Numerical Algorithms 90, 323–361 (2022), doi:10.1007/s11075-021-01190-x, Sections 3 and 4.1 — integral error bounds and limitations of quadrature estimates for oscillatory defects.
• Saad, Y., “Analysis of some Krylov subspace approximations to the matrix exponential operator,” SIAM Journal on Numerical Analysis 29, 209–228 (1992), doi:10.1137/0729014 — projection and error analysis.