Numerical Error, Conditioning, and Stability

A stable algorithm controls the numerical errors it introduces; a well-conditioned problem does not amplify small input changes excessively. These are different properties. A solver can accurately solve a nearby matrix problem while the requested physical quantity remains sensitive to discretization, uncertain parameters, or a nearly closed spectral gap.

Required background. Precision diagnostics explains roundoff diagnostics; basis choice defines the overlap metric used by a finite representation.

Forward error and backward error

For an exact task $y=f(x)$ and computed result $\widehat y$, forward error measures $\widehat y-y$. Backward error asks how much the input must change so that $\widehat y=f(x+\delta x)$. A backward-stable algorithm produces a small $\delta x$ relative to the input scale. The forward error also depends on the sensitivity of $f$.

For a nonsingular linear system $Ax=b$, define $r=b-A\widehat x$. Since $A(x-\widehat x)=r$,

$$\frac{\|x-\widehat x\|}{\|x\|} \le \kappa(A)\frac{\|r\|}{\|b\|}, \qquad \kappa(A)=\|A\|\|A^{-1}\|.$$

The inequality uses $\|b\|\le\|A\|\|x\|$ and assumes $b\ne0$. A small relative residual guarantees a small forward error only with a useful condition bound. The scaled quantity $\|r\|/(\|A\|\|\widehat x\|+\|b\|)$ is a common normwise backward-error diagnostic, with the exact interpretation depending on allowed perturbations.

For $A=\operatorname{diag}(1,10^{-8})$ and $b=(1,10^{-8})^T$, the exact solution is $(1,1)^T$. The approximation $(1,0)^T$ has residual norm $10^{-8}$ but an order-one error in the second component. Tightening the stopping tolerance can reduce the residual, but cannot change the conditioning of this linear system.

Hermitian eigenpair residuals

Let $H=H^\dagger$, $\|v\|_2=1$, and $r=Hv-\theta v$. Expanding in an orthonormal eigenbasis gives

$$\|r\|_2^2=\sum_j |c_j|^2|E_j-\theta|^2, \qquad \sum_j|c_j|^2=1.$$

Consequently, at least one eigenvalue lies within $\|r\|_2$ of $\theta$. This statement does not identify which eigenvalue, nor does it prove $v$ approximates one particular eigenvector. If the target eigenspace is separated from every other eigenvalue by a distance at least $\delta$ from $\theta$, then

$$\|(I-P)v\|_2\le\frac{\|r\|_2}{\delta},$$

where $P$ projects onto that eigenspace. Near degeneracies, compare subspaces and observables rather than individual eigenvectors. A tiny residual may coexist with a large change of eigenvector direction when the gap is tiny.

These are finite-matrix statements. They do not control errors between $H$ and the intended continuum Hamiltonian. A nonorthogonal basis also changes which residual norm belongs in the bound.

Residuals in a nonorthogonal trial space

Let $H_{ij}=\langle\phi_i|\hat H|\phi_j\rangle$ and $S_{ij}=\langle\phi_i|\phi_j\rangle$, with $H=H^\dagger$ and $S$ positive definite. For nonzero coefficients $c$, define $r=Hc-\theta Sc$. These residual entries are inner products with the trial functions, not expansion coefficients of a state. Their ordinary Euclidean norm therefore depends on how the basis functions are scaled.

Use the Cholesky reduction, $S=LL^\dagger$, $\widetilde H=L^{-1}HL^{-\dagger}$, and $y=L^\dagger c$. It gives

$$\widetilde H y-\theta y=L^{-1}r, \qquad \eta=\frac{\|L^{-1}r\|_2}{\sqrt{c^\dagger Sc}} =\sqrt{\frac{r^\dagger S^{-1}r}{c^\dagger Sc}}.$$

Thus $\operatorname{dist}(\theta,\sigma(\widetilde H))\le\eta$. The inverse metric appears because $r$ contains overlaps; the corresponding coefficient residual is $S^{-1}r$, whose $S$-norm is $\sqrt{r^\dagger S^{-1}r}$. Evaluate $L^{-1}r$ by a triangular solve, without forming $S^{-1}$. Uncertainty in $H$, $S$, or the factorization remains a separate numerical error.

For an explicit scaling check, take the physical two-level Hamiltonian $\hat H=\operatorname{diag}(0,2)$ in fixed energy units and the normalized state $\psi=(\sqrt3/2,1/2)^T$. Its Rayleigh quotient is $\theta=1/2$, and $\|(\hat H-\theta)\psi\|=\sqrt3/2$. Now use trial functions $\phi_i=\epsilon e_i$, with $\epsilon>0$. Then

$$S=\epsilon^2I,\qquad H=\epsilon^2\hat H,\qquad c=\psi/\epsilon,\qquad r=\epsilon(\hat H-\theta)\psi.$$

At $\epsilon=10^{-8}$, the raw norm is only $8.66\times10^{-9}$, although the state and its spectral error have not improved. The invariant residual is still $\eta=\sqrt3/2$; indeed the nearest eigenvalue remains $1/2$ away from $\theta$. Here $\kappa_2(S)=1$: the misleading raw residual is caused by basis scaling, even without poor overlap conditioning.

When the trial space is truncated, $\eta$ measures only its projected residual. If $P$ is its orthogonal Hilbert-space projector and $\psi=\sum_i c_i\phi_i$ is normalized and in the operator domain, then

$$\|(\hat H-\theta)\psi\|^2 =\eta^2+\|(I-P)\hat H\psi\|^2.$$

The omitted component is orthogonal to the trial space. An exactly solved Ritz problem has $\eta=0$ but can still have a nonzero full residual. Neither this missing component nor continuum accuracy is certified by the finite-basis solve; use basis refinement and appropriate independent bounds. For trial functions only in the quadratic-form domain, the full operator residual need not even exist.

Non-Hermitian sensitivity

For a simple eigenvalue of a diagonalizable non-Hermitian matrix, let $Hx=\lambda x$ and $y^\dagger H=\lambda y^\dagger$. First-order perturbation theory gives

$$\delta\lambda=\frac{y^\dagger(\delta H)x}{y^\dagger x} +O(\|\delta H\|^2).$$

For unit-norm left and right vectors, a small $|y^\dagger x|$ signals large sensitivity. The Hermitian residual-distance bound does not transfer unchanged to such a problem. Complex absorbing potentials, effective Hamiltonians, and Liouvillians can be nonnormal; inspect left/right conditioning or pseudospectral behavior when apparently stable eigenvalues are suspect.

At a defective eigenvalue, the simple-eigenvalue expansion itself fails. More digits and a smaller residual do not justify interpreting every decimal place of a near-exceptional-point spectrum as robust physics.

Stability of propagation

For a time-step map $U_h$ repeatedly applied to a state, a local error may accumulate or be amplified. In units with $\hbar=1$, explicit Euler for $i\dot c=Ec$ gives $c_{n+1}=(1-iEh)c_n$, with amplification modulus $\sqrt{1+E^2h^2}>1$. Norm growth follows even though the differential equation is unitary.

Conversely, an exactly unitary step can have the wrong phase. Stability and invariant preservation are valuable constraints, but accuracy still requires step refinement or comparison with a known solution. Report both the conservation defect and the error in the target observable.

Exercises

For $H=\operatorname{diag}(0,\epsilon)$ with $\epsilon>0$ and $v=(1,1)^T/\sqrt2$, evaluate the Rayleigh quotient and residual. Does a small residual for small $\epsilon$ make $v$ close to the ground-state eigenvector?

Solution

$\theta=\epsilon/2$ and $r=(-\epsilon/2,\epsilon/2)^T/\sqrt2$, so $\|r\|=|\epsilon|/2$. Nevertheless the overlap with $(1,0)^T$ is always $1/\sqrt2$. The gap shrinks with the residual; the two-dimensional low-energy subspace is the stable object.

A three-level Hamiltonian is $\hat H=\begin{pmatrix}0&0&g\\0&2&0\\g&0&4\end{pmatrix}$, with real $g\ne0$. Restrict the trial space to the first two coordinate vectors and take $\psi=e_1$, $\theta=0$. Compare the Ritz residual, full residual, and exact ground energy. Does solving the trial problem more accurately remove the discrepancy?

Solution

The trial Hamiltonian is $\operatorname{diag}(0,2)$, so its residual is exactly zero. In the full space, $(\hat H-\theta)\psi=g e_3$ has norm $|g|$. The coupled first/third block has eigenvalues $2\pm\sqrt{4+g^2}$, hence the ground energy is $2-\sqrt{4+g^2}=-g^2/4+O(g^4)$ at small $g$. The error comes from the omitted direction; tightening the trial eigensolver cannot remove it.

References

• Anderson, E., et al., LAPACK Users’ Guide, 3rd ed., SIAM, 1999, Accuracy and Stability and Generalized Symmetric Definite Eigenproblems — spectral bounds and the Cholesky reduction used for the overlap-metric residual.
• Higham, N. J., Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, publisher edition — conditioning and stability.
• Trefethen, L. N., and M. Embree, Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators, Princeton University Press, 2005 — nonnormal sensitivity and transient behavior.