Matrix Exponentials Numerically

Matrix exponentials appear whenever a finite-dimensional linear system is solved exactly over a time interval. In quantum mechanics, the central example is time-independent Schrödinger evolution:

$$\psi(t+\Delta t) = \exp\left( - \frac{i}{\hbar} H\Delta t \right) \psi(t).$$

The mathematical definition of matrix functions is reviewed in Matrix Functions and Exponentials. This page focuses on numerical choices: diagonalization, scaling and squaring, and Krylov methods for applying an exponential to a vector.

Full Exponential Versus Action

There are two different computational tasks:

1. form the full matrix $e^A$;
2. compute the action $e^A v$ on one or more vectors.

For large quantum systems, the second task is usually the right one. Even if $A$ is sparse, $e^A$ is usually dense. Forming the full exponential can destroy the memory advantage of Sparse Matrices.

Time evolution of a state usually needs

$$v_{\mathrm{new}} = e^A v, \qquad A = - \frac{i}{\hbar} H\Delta t,$$

not the full matrix $e^A$.

Exact Diagonalization Approach

If a Hermitian Hamiltonian is small enough to diagonalize densely, write

$$H = V\Lambda V^\dagger, \qquad \Lambda = \operatorname{diag}(E_1,\dots,E_n).$$

Then

$$\exp\left( - \frac{i}{\hbar} Ht \right) = V \operatorname{diag} \left( e^{-iE_1t/\hbar}, \dots, e^{-iE_nt/\hbar} \right) V^\dagger.$$

This method is accurate and transparent when the full spectrum is already available. It is also a useful benchmark for smaller versions of a problem.

The cost is the dense eigensolve and the storage of eigenvectors. For large sparse Hamiltonians, diagonalizing the whole matrix just to propagate one state is usually wasteful.

General Diagonalizable Matrices

For a diagonalizable non-Hermitian matrix,

$$A = X\Lambda X^{-1}, \qquad e^A = Xe^\Lambda X^{-1}.$$

This formula is mathematically simple but can be numerically dangerous when $X$ is ill conditioned. Non-normal matrices can have sensitive eigenvectors and large transient behavior even when their eigenvalues look harmless.

For closed quantum systems with Hermitian $H$, the unitary diagonalization is much better conditioned than a generic eigenvector decomposition. For non-Hermitian effective Hamiltonians, absorbing layers, and Liouvillians, eigenvector conditioning must be treated seriously.

Scaling and Squaring

Scaling and squaring is a standard dense-matrix method for computing $e^A$. Choose an integer $s$ so that $A/2^s$ has a manageable norm, approximate $e^{A/2^s}$ accurately, and then square repeatedly:

$$e^A = \left( e^{A/2^s} \right)^{2^s}.$$

The small exponential is often approximated by a Padé approximant. This is the basis of robust general-purpose matrix exponential routines.

Scaling and squaring is excellent for moderate dense matrices when the full matrix exponential is required. It is less attractive when $A$ is huge and sparse but only $e^A v$ is needed, because repeated squaring tends to create dense matrices.

Krylov Exponential Action

For large sparse quantum Hamiltonians, a common goal is to approximate $e^A v$ without forming $e^A$. Build a Krylov subspace

$$\mathcal K_m(A,v) = \operatorname{span} \{v,Av,A^2v,\dots,A^{m-1}v\}.$$

Let $Q_m$ be an orthonormal basis for this subspace, with first vector

$$q_1 = \frac{v}{\lVert v\rVert}.$$

Project $A$ to the smaller matrix

$$A_m = Q_m^\dagger A Q_m.$$

Then approximate

$$e^A v \approx \lVert v\rVert Q_m e^{A_m}e_1,$$

where $e_1$ is the first coordinate vector. The small exponential $e^{A_m}$ can be computed densely, because $m$ is much smaller than the full dimension.

This is the exponential-action analogue of the Krylov ideas in Sparse Eigensolvers. The difference is that here the goal is a time-evolved vector, not selected eigenpairs.

Lanczos for Hermitian Hamiltonians

For quantum real-time evolution with Hermitian $H$,

$$A = - \frac{i}{\hbar} H\Delta t$$

is anti-Hermitian. A Lanczos basis built from $H$ gives a small Hermitian projected Hamiltonian $T_m$. The Krylov approximation becomes

$$e^{-iH\Delta t/\hbar}v \approx \lVert v\rVert Q_m \exp\left( - \frac{i}{\hbar} T_m\Delta t \right) e_1.$$

Because $T_m$ is Hermitian, the small exponential is unitary inside the Krylov subspace. The remaining error is truncation error from using a finite Krylov dimension and any loss of orthogonality in finite precision.

Time Evolution and Unitarity

If $H$ is Hermitian, the exact propagator

$$U(\Delta t) = \exp\left( - \frac{i}{\hbar} H\Delta t \right)$$

is unitary. Numerically, one should monitor

$$\lVert \psi_{n+1}\rVert - \lVert \psi_n\rVert$$

and physically relevant observables. Norm preservation alone is not enough, but norm drift is an immediate warning for closed-system dynamics.

The general time-step context is Time-Stepping Methods. Matrix exponential actions are one way to build a time step with better structural behavior than generic explicit methods.

Time-Dependent Hamiltonians

When $H(t)$ changes during the step, the exact evolution is time ordered. A common approximation freezes or samples the Hamiltonian over a short interval:

$$\psi_{n+1} \approx \exp\left( - \frac{i}{\hbar} H(t_n+\Delta t/2)\Delta t \right) \psi_n.$$

This midpoint exponential can be useful, but it is not the full time-ordered solution unless commutator effects are negligible. If

$$[H(t_1),H(t_2)] \ne 0,$$

then decreasing $\Delta t$ and comparing with a time-ordering-aware method is essential. See Time Ordering.

Conditioning and Error

The matrix exponential can be sensitive to perturbations, especially for non-normal matrices. Two matrices with nearby entries can have exponentials that differ substantially over long times or in strongly non-normal problems.

For anti-Hermitian quantum generators, the exact exponential has norm $1$ in the Hilbert-space norm. This helps, but it does not remove discretization error, Krylov truncation error, roundoff error, or time-dependent approximation error.

Useful diagnostics include:

• norm drift for closed-system states;
• residual or error estimates from the exponential-action routine;
• convergence under Krylov dimension or tolerance changes;
• convergence under time-step refinement;
• comparison with dense diagonalization on a smaller problem;
• conservation of energy for time-independent $H$;
• sensitivity to basis, grid, and boundary choices.

Choosing a Method

Task	Common method	Main caution
small Hermitian $H$, many times	dense diagonalization	full eigenvector storage
moderate dense general $A$	scaling and squaring	non-normal sensitivity
large sparse $H$, one state	Krylov exponential action	truncation and orthogonality
many right-hand sides	block Krylov or reuse structure	memory growth
time-dependent $H(t)$	short-step exponential or Magnus-type method	time ordering

No method removes the need for physical convergence checks. The numerical exponential solves the finite matrix problem; it does not prove that the finite matrix is the right approximation to the continuum system.

Common Mistakes

• Forming $e^A$ when only $e^A v$ is needed.
• Assuming $e^A$ is sparse because $A$ is sparse.
• Using a non-Hermitian eigenvector decomposition without checking conditioning.
• Treating norm conservation as a complete accuracy test.
• Ignoring time ordering for a time-dependent Hamiltonian.
• Comparing long-time wavefunctions without separating global phase, relative phase, and shape errors.
• Forgetting that dense diagonalization solves the finite matrix exactly only after discretization errors have already entered.
• Using a Krylov dimension that is too small for the time step or spectral width.

Cross-Links

• Matrix Functions and Exponentials
• Time-Stepping Methods
• Matrix Diagonalization
• Sparse Matrices
• Sparse Eigensolvers
• Conditioning and Stability
• Floating-Point Arithmetic
• Unitary Operators
• Time-Evolution Operator
• Time Ordering
• Unitary Time Evolution

References

• N. J. Higham, Functions of Matrices: Theory and Computation, SIAM, 2008.
• A. H. Al-Mohy 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.
• R. B. Sidje, “Expokit: A software package for computing matrix exponentials”, ACM Transactions on Mathematical Software 24, 130-156, 1998.
• Y. Saad, Numerical Methods for Large Eigenvalue Problems, revised ed., SIAM, 2011.
• G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, 2013.
• D. J. Tannor, Introduction to Quantum Mechanics: A Time-Dependent Perspective, University Science Books, 2007.

Exercises

1. Let $H=V\Lambda V^\dagger$ be Hermitian diagonalization. Derive the formula for $e^{-iHt/\hbar}$.

Solution

Using the power series and $V^\dagger V=I$,

$$H^k = V\Lambda^kV^\dagger.$$

Therefore

$$e^{-iHt/\hbar} = V e^{-i\Lambda t/\hbar} V^\dagger.$$

Since $\Lambda$ is diagonal, $e^{-i\Lambda t/\hbar}$ is diagonal with entries $e^{-iE_jt/\hbar}$.

2. Show that if $A^\dagger=-A$, then $e^A$ is unitary.

Solution

For an anti-Hermitian $A$,

$$(e^A)^\dagger = e^{A^\dagger} = e^{-A}.$$

Because $A$ commutes with $-A$,

$$(e^A)^\dagger e^A = e^{-A}e^A = I.$$

3. Why is $e^A v$ often a better target than $e^A$ for a large sparse Hamiltonian?

Solution

The state update needs only the action of the exponential on the current state. The full matrix exponential is usually dense even when $A$ is sparse, so forming it can be far more expensive in memory and time. Krylov action methods exploit sparse matrix-vector products and avoid constructing the dense propagator.

4. In the Krylov approximation $e^A v\approx\lVert v\rVert Q_m e^{A_m}e_1$, what is the role of the small matrix $A_m$?

Solution

$A_m=Q_m^\dagger A Q_m$ is the projection of the large matrix onto the Krylov subspace. The exponential is computed on this small matrix and then lifted back with $Q_m$. The approximation is accurate when the Krylov subspace captures the part of the dynamics generated by repeatedly applying $A$ to $v$.