Error Estimates

An error estimate is a reasoned statement of how far a numerical answer may be from the intended mathematical or physical quantity. It may come from a theorem, a residual, a refinement study, a statistical error bar, an exactly solvable benchmark, or a comparison between independent methods.

In quantum mechanics, error estimates are not optional polish. They decide whether an energy splitting is resolved, whether a wave packet has propagated accurately, whether a Monte Carlo signal is larger than noise, and whether a small symmetry breaking is physical or numerical.

Different checks control different errors. Convergence Tests gives the practical refinement workflow. Detailed mechanisms live in Floating-Point Arithmetic, Conditioning and Stability, Discretization, and Monte Carlo Basics.

Helpful background. Review conditioning and stability to interpret a residual separately from the error in a solution.

From a physical model to a computed answer

A computation usually passes from a physical model to a continuum mathematical problem, then to a finite numerical problem, and finally to a computed answer.

An error estimate should say which gap it addresses. A small eigensolver residual does not estimate model error. A Monte Carlo standard error does not estimate finite-size bias. A grid-refinement trend does not estimate a coding mistake.

For a reported quantity $Q_{\mathrm{comp}}$, a schematic decomposition is

$$Q_{\mathrm{comp}}-Q_{\mathrm{phys}} = e_{\mathrm{model}} + e_{\mathrm{domain}} + e_{\mathrm{disc}} + e_{\mathrm{alg}} + e_{\mathrm{round}} + e_{\mathrm{stat}} + e_{\mathrm{impl}}.$$

Not every term appears in every calculation, and the terms need not be independent. The decomposition is a checklist: it prevents one small diagnostic from being mistaken for a full uncertainty analysis.

Truncation and Discretization Error

Truncation error comes from replacing an infinite or continuum object by a finite approximation. Examples include:

• grid spacing $h$ in finite differences;
• time step $\Delta t$ in time evolution;
• basis cutoff $N$ in a spectral or variational method;
• finite box size $L$ for a problem on the line;
• quadrature order or number of nodes;
• angular momentum, occupation-number, or energy cutoffs.

If a method has leading error $Ch^p$, then

$$Q(h) = Q + Ch^p + O(h^{p+1})$$

or sometimes

$$Q(h) = Q + Ch^p + O(h^{p+2}),$$

depending on symmetry and the method. The order statement is meaningful only when the solution is smooth enough and the boundary implementation has the same accuracy as the interior scheme.

For a second-order finite-difference eigenvalue calculation, one expects low resolved levels to change by about a factor of $4$ when $h$ is halved, after box-size error and algebraic solver error are under control.

Richardson Error Estimate

Suppose

$$Q(h)=Q+Ch^p+O(h^{p+1})$$

and the same calculation is repeated at $h/2$. Then

$$Q(h/2)-Q(h) = C\left(2^{-p}-1\right)h^p + O(h^{p+1}).$$

The finer-grid error is approximately

$$Q(h/2)-Q \approx \frac{Q(h)-Q(h/2)}{2^p-1}.$$

Thus an error estimate is

$$\epsilon_{h/2} \approx \frac{ \lvert Q(h/2)-Q(h)\rvert }{2^p-1}.$$

This estimate is useful only in the asymptotic refinement regime. If the observed ratios are inconsistent, another error source is probably dominating or the assumed order has not been reached.

Domain Truncation Error

A finite box can approximate an infinite-domain bound state, but the box position is an error source. If a localized wavefunction is not negligible at the boundary, refining the grid inside the same box will not fix the result.

Domain truncation checks include:

• increase the box size while holding resolution fixed;
• inspect wavefunction amplitude or probability near boundaries;
• compare boundary conditions when they should be irrelevant;
• estimate exponential tails in forbidden regions;
• avoid interpreting finite-box continuum levels as continuum energies.

For scattering and continuum states, box-size dependence is often the signal rather than a nuisance. The finite box changes the spectrum, and extracting physical observables requires additional analysis.

Algebraic Solver Error

After discretization, a finite problem still has to be solved. Eigensolvers, linear solvers, nonlinear root finders, and optimization algorithms introduce algebraic error.

For a linear system $Ax=b$, the residual is

$$r = b-A\tilde x.$$

The forward error satisfies

$$x-\tilde x = A^{-1}r$$

when $A$ is invertible. Therefore a small residual is reassuring only after conditioning is considered.

For an eigenpair of a Hermitian matrix, a common residual is

$$r = H\tilde\psi - \tilde E\tilde\psi.$$

The residual norm measures how well the finite matrix equation is solved. It does not include grid error, basis truncation error, or physical model error.

Roundoff Error

Roundoff error comes from finite-precision arithmetic. A standard local model is

$$\operatorname{fl}(x\circ y) = (x\circ y)(1+\delta), \qquad \lvert\delta\rvert\le u,$$

where $u$ is the unit roundoff and $\circ$ is a correctly rounded basic operation, assuming no overflow or underflow and a nonzero exact result.

Roundoff becomes visible when:

• many operations accumulate error;
• cancellation removes leading digits;
• derivative formulas divide by small powers of $h$;
• matrices are ill conditioned;
• nearly degenerate subspaces are compared vector by vector;
• tiny probabilities, splittings, or tunneling amplitudes are inferred from large intermediate quantities.

For a centered finite-difference second derivative, a schematic balance is

$$\text{error} \sim C_1h^2 + C_2\frac{u}{h^2}.$$

The first term decreases with refinement; the second can grow. The best grid spacing is not always the smallest grid spacing.

Statistical Error

Statistical error appears when an answer is estimated from random samples: measurement shots, Monte Carlo integration, stochastic trajectories, or randomized numerical algorithms.

For independent samples $Y_1,\dots,Y_N$ with sample mean

$$\bar Y_N = \frac1N \sum_{k=1}^N Y_k,$$

the estimated standard error is

$$\widehat{\operatorname{SE}} = \frac{s_N}{\sqrt N},$$

where $s_N$ is the sample standard deviation.

If samples are correlated, replace $N$ by an effective sample size. For Markov chains, a common approximation is

$$N_{\mathrm{eff}} \approx \frac{N}{2\tau_{\mathrm{int}}},$$

where $\tau_{\mathrm{int}}=\tfrac12+\sum_{k\ge1}\rho(k)$ for the sampled observable, assuming a stationary process with summable autocorrelations. Independent samples have $\tau_{\mathrm{int}}=1/2$.

Statistical error bars do not automatically include bias from equilibration, time-step error, finite-size effects, trial-wavefunction choices, or sign and phase problems.

Bias Versus Variance

Variance is random scatter. Bias is systematic displacement.

More samples reduce the variance of an unbiased estimator:

$$\operatorname{SE} \propto N^{-1/2}.$$

More samples do not automatically reduce bias. Examples of bias include:

• finite time step in an imaginary-time path integral;
• finite population or finite walker bias;
• incomplete equilibration in a Markov chain;
• variational bias from a restricted ansatz;
• finite box or finite basis truncation;
• regularization or cutoff choices.

A result with a tiny statistical error bar can still be wrong if the systematic error is larger.

Combining Error Estimates

Independent statistical errors are often combined in quadrature:

$$\epsilon_{\mathrm{tot}} \approx \left( \epsilon_1^2+\epsilon_2^2+\cdots \right)^{1/2}.$$

Deterministic systematic errors are less friendly. If their signs and correlations are unknown, a conservative bound adds magnitudes:

$$\epsilon_{\mathrm{sys}} \le \epsilon_1+\epsilon_2+\cdots.$$

In practice, quote the dominant known contributions separately when possible:

$$E = -0.5000003 \pm 2\times10^{-7}_{\mathrm{grid}} \pm 1\times10^{-8}_{\mathrm{solver}}.$$

The point is not to decorate the answer. The point is to make clear which uncertainty is controlled and which one remains the limiting source.

Error in Observables

Suppose a normalized exact state $\psi$ is approximated by another normalized state $\tilde\psi$ with

$$\lVert \tilde\psi-\psi\rVert \le \epsilon.$$

For a bounded observable $A$, a simple bound is

$$\left\lvert \langle\tilde\psi,A\tilde\psi\rangle - \langle\psi,A\psi\rangle \right\rvert \le 2\lVert A\rVert\epsilon.$$

The bound follows by adding and subtracting $\langle\tilde\psi,A\psi\rangle$ and applying Cauchy–Schwarz to the two differences. It is exact for normalized states; many quantum observables are unbounded in the continuum. Still, it gives the right warning: a small state-vector error must be interpreted in the norm and operator scale relevant to the observable.

For eigenstates, energies can converge faster than wavefunctions in some variational settings, while local observables may converge more slowly. Estimate the error in the quantity you actually report.

Small Splittings and Resolved Digits

If two computed energies are

$$E_1\pm\epsilon_1, \qquad E_2\pm\epsilon_2,$$

then the splitting

$$\Delta E=E_2-E_1$$

has uncertainty that depends on correlations between the two estimates. Common errors can cancel, while errors of opposite signs can add. If the errors are independent statistical errors, one may estimate

$$\epsilon_{\Delta} \approx \left( \epsilon_1^2+\epsilon_2^2 \right)^{1/2}.$$

For systematic discretization errors, cancellation is possible but must be demonstrated. A small printed $\Delta E$ is not resolved unless it is larger than the relevant uncertainty and stable under refinement.

Reporting Numerical Results

To make a numerical result interpretable, state:

• the quantity being approximated;
• the numerical method and discretization parameters;
• the refinement or residual checks performed;
• the dominant error estimate;
• whether the error is statistical, deterministic, or heuristic;
• the number of significant digits justified by the estimate;
• known uncontrolled errors.

Avoid reporting

$$E=-0.499999999731$$

if the grid-refinement uncertainty is $10^{-5}$. A better report is

$$E=-0.50000\pm 0.00001.$$

Fewer honest digits are more informative than many decorative digits.

Independent checks

Useful checks include:

• exact solutions such as the harmonic oscillator, infinite well, and two-level system;
• residuals for finite equations;
• conservation of norm, energy, trace, positivity, or symmetry labels;
• grid, basis, box-size, and time-step refinement;
• comparison between finite-difference, spectral, and variational methods;
• independent random seeds and autocorrelation analysis;
• higher precision or alternative summation for cancellation-prone calculations.

No single check is universal. The best error estimate triangulates from several checks whose failure modes are different.

For standard exact and controlled test cases, see Benchmark Problems.

Common Mistakes

• Reporting Monte Carlo error bars while ignoring systematic bias.
• Treating an eigensolver residual as a discretization error estimate.
• Quoting the formal order of a method without showing observed refinement.
• Refining $h$ while leaving the domain size or time step fixed and unconverged.
• Assuming roundoff is negligible because double precision was used.
• Reporting more digits than the uncertainty supports.
• Combining systematic errors in quadrature without justification.
• Using a conserved norm as the only accuracy test for time evolution.
• Claiming a tiny splitting without showing it survives all relevant error estimates.

References

• N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002.
• R. J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations, SIAM, 2007.
• C. P. Robert and G. Casella, Monte Carlo Statistical Methods, 2nd ed., Springer, 2004.
• J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, 3rd ed., Springer, 2002.
• J. M. Thijssen, Computational Physics, 2nd ed., Cambridge University Press, 2007.
• L. N. Trefethen and D. Bau, Numerical Linear Algebra, SIAM, 1997.

Exercises

1. Richardson estimate for a second-order method.

Suppose a quantity is computed as $Q(h)=1.024$ and $Q(h/2)=1.006$ with leading error $Ch^2$. Estimate the error in the finer value.

Solution

For $p=2$,

$$\epsilon_{h/2} \approx \frac{\lvert Q(h/2)-Q(h)\rvert}{2^2-1} = \frac{\lvert1.006-1.024\rvert}{3} = 0.006.$$

So the finer value is estimated as $1.006\pm0.006$ from this two-level refinement model.

2. Optimal grid spacing in a schematic error balance.

Minimize

$$E(h)=C_1h^2+C_2\frac{u}{h^2}$$

over $h\gt0$.

Solution

Differentiate:

$$\frac{dE}{dh} = 2C_1h - 2C_2u h^{-3}.$$

Set this to zero:

$$C_1h = C_2u h^{-3}, \qquad h^4 = \frac{C_2u}{C_1}.$$

Thus

$$h_{\mathrm{opt}} = \left( \frac{C_2u}{C_1} \right)^{1/4}.$$

The formula is schematic, but it shows why making $h$ arbitrarily small can increase roundoff-dominated error.

3. Monte Carlo standard error.

A Monte Carlo estimate uses $N=40{,}000$ independent samples with sample standard deviation $s_N=3$. Estimate the standard error of the mean.

Solution

Use

$$\widehat{\operatorname{SE}} = \frac{s_N}{\sqrt N} = \frac{3}{200} = 0.015.$$

4. Splitting uncertainty.

Two independently estimated energies have uncertainties $\epsilon_1=2\times10^{-6}$ and $\epsilon_2=3\times10^{-6}$. Estimate the statistical uncertainty in $\Delta E=E_2-E_1$.

Solution

For independent statistical uncertainties, combine in quadrature:

$$\epsilon_\Delta \approx \left( (2\times10^{-6})^2 + (3\times10^{-6})^2 \right)^{1/2} = \sqrt{13}\times10^{-6} \approx 3.6\times10^{-6}.$$

5. Why is a small residual not a full error estimate?

In particular, a normalized computed eigenvector has residual norm $10^{-12}$ in the chosen energy units, but halving the grid spacing changes its eigenvalue by $10^{-3}$ in the same units. Is it reasonable to report twelve correct decimal places for the continuum energy?

Solution

A residual says how well the computed answer solves the finite equations. It does not by itself include conditioning, discretization error, domain truncation, model error, roundoff amplification, or statistical uncertainty. For an ill-conditioned linear system, even a small residual can correspond to a large forward error because the error is $A^{-1}r$.

The eigenvalue example therefore does not justify twelve correct decimal places. The finite matrix equation is solved much more accurately than the grid approximates the continuum problem. Continue refining the grid, test the domain size, and estimate the remaining error from the observed convergence. One difference between grids is a useful warning, not yet a reliable error bound.