Numerical stability
Numerical stability describes the sensitivity of a computational process to perturbations introduced during arithmetic evaluation. Such perturbations arise because floating-point arithmetic represents only a finite subset of the real numbers and replaces most exact operations with rounded approximations. A numerically stable algorithm produces an answer consistent with the exact solution of a nearby problem, whereas an unstable algorithm can amplify small computational errors until they dominate the result.
Stability is a property of an algorithm rather than solely of the mathematical problem being solved. A well-conditioned problem can be evaluated by an unstable algorithm, and an ill-conditioned problem can remain sensitive even when its algorithm is stable. The distinction between stability and conditioning forms the central conceptual separation in modern numerical analysis.
Mathematical framework
Let a mathematical problem be represented by a function
[ f:X\rightarrow Y, ]
and let an algorithm implemented in finite arithmetic produce the computed mapping
[ \widehat{f}:X\rightarrow Y. ]
For an input (x), the forward error is the discrepancy between the computed result and the exact result,
[ |\widehat{f}(x)-f(x)|. ]
An absolute forward error measures this discrepancy directly. A relative forward error normalizes it by a scale associated with the exact result, commonly (|f(x)|), when that quantity is nonzero.
Forward error alone does not identify whether an inaccurate result arose from the problem or from the algorithm. If (f) is highly sensitive near (x), even a small perturbation in the input can cause a large change in the output. This intrinsic sensitivity is described by the condition number. In a differentiable scalar problem, the relative condition number commonly takes the form
[ \kappa_f(x)=\left|\frac{x f'(x)}{f(x)}\right|. ]
For vector-valued problems, the derivative is replaced by the relevant Jacobian matrix, and the resulting sensitivity depends on the selected vector norm.
A forward-stable algorithm has an error comparable to the product of machine precision and the condition number, subject to constants determined by the operation count and computational structure. This relation is commonly expressed schematically as
[ \frac{|\widehat{f}(x)-f(x)|}{|f(x)|} \lesssim \kappa_f(x),u, ]
where (u) denotes the unit roundoff. The expression separates amplification caused by the problem from error introduced by the computation, although the separation is not exact for every problem class.
Backward stability
Backward error analysis interprets a computed result as the exact answer to a perturbed input. An algorithm is backward stable when, for each input (x), there exists a nearby input (x+\delta x) satisfying
[ \widehat{f}(x)=f(x+\delta x), ]
with (|\delta x|/|x|) on the order of unit roundoff. Under this interpretation, the algorithm introduces no error larger than a small disturbance of the original data.
Backward stability does not imply a small forward error when the underlying problem is ill-conditioned. If (f) strongly amplifies input perturbations, the exact solution of a nearby problem can remain far from the exact solution of the original problem. To first order, backward error, conditioning, and forward error are related by
[ \text{forward error} \lesssim \text{condition number} \times \text{backward error}. ]
This relation provides the standard framework for distinguishing computational instability from unavoidable sensitivity. It also explains why a large residual need not be equivalent to a large solution error, and why a small residual does not by itself establish the accuracy of a solution.
For a linear system
[ Ax=b, ]
a computed vector (\widehat{x}) has residual
[ r=b-A\widehat{x}. ]
The residual measures how closely (\widehat{x}) satisfies the stated equations. Backward error instead measures the smallest perturbations (\Delta A) and (\Delta b) for which
[ (A+\Delta A)\widehat{x}=b+\Delta b. ]
The condition number of (A), usually expressed as (\kappa(A)=|A||A^{-1}|), governs how these perturbations affect the solution. Near-singular matrices have large condition numbers because small changes in their entries can produce large changes in the solution.
Sources of instability
Cancellation
Loss of significance occurs when nearly equal quantities are subtracted. If each operand already contains a small relative error, the leading digits cancel while the absolute errors remain. The relative error of the difference can consequently become much larger than the relative errors of the operands.
The quadratic formula provides a standard illustration. For
[ ax^2+bx+c=0, ]
the direct expression
[ x=\frac{-b\pm\sqrt{b^2-4ac}}{2a} ]
can subtract nearly equal quantities for one choice of sign. An algebraically equivalent representation obtained from the product of the roots can avoid that particular cancellation. The exact formulas are mathematically identical, but their finite-arithmetic evaluations need not have comparable stability.
Cancellation is not inherently erroneous. Exact cancellation may be an essential part of a calculation, and subtraction can be well behaved when the operands are known to sufficient absolute accuracy. Instability arises when uncertain leading quantities are removed and previously insignificant errors become dominant.
Growth of intermediate quantities
An algorithm can produce intermediate values much larger than its final result. These values may overflow the available exponent range, or they may magnify rounding errors before later cancellation reduces the mathematical expression. The resulting instability depends on the computational path rather than on the exact formula alone.
This behavior appears in some forms of Gaussian elimination. The growth factor measures the increase in matrix-entry magnitude during elimination relative to the original matrix. Partial pivoting usually limits this growth in practical calculations, although specially constructed matrices can produce exponential growth. The stability analysis therefore concerns both the rounding model and the transformations applied to the matrix.
Accumulation of rounding error
Repeated arithmetic operations introduce a sequence of local errors. In a direct worst-case bound, (n) operations can produce an error proportional to (nu), provided (nu) remains small. The actual behavior can be less severe when errors partially cancel, but deterministic stability analysis does not assume such cancellation without an explicit error model.
The ordinary left-to-right summation of floating-point numbers illustrates this accumulation. Each addition rounds the current partial sum, so small terms may cease to affect a much larger accumulator. Pairwise summation reduces the depth of the addition tree and commonly changes the leading error bound from linear to logarithmic dependence on the number of terms. Kahan summation carries a compensation variable representing low-order information lost in preceding additions, thereby producing a different finite-arithmetic recurrence from ordinary summation.
Floating-point model
Standard analyses often model a basic arithmetic operation by
[ \operatorname{fl}(x\circ y)=(x\circ y)(1+\delta), \qquad |\delta|\leq u, ]
where (\circ) denotes one elementary operation and the exact result lies within the normal floating-point range. This model represents correct rounding by a small relative perturbation. It requires modification near underflow, where subnormal numbers permit gradual loss of relative precision, and it does not describe operations whose exact results overflow.
The IEEE 754 standard specifies widely used binary and decimal floating-point formats. It defines rounding modes and exceptional values, while also specifying behavior for infinities and invalid operations. Numerical stability remains an algorithmic property under these rules; standardization makes the arithmetic model reproducible but does not make every sequence of operations stable.
A common shorthand for accumulated rounding is
[ \gamma_n=\frac{nu}{1-nu}, ]
defined when (nu<1). Products of factors (1+\delta_i), with (|\delta_i|\leq u), can then be represented by a single factor (1+\theta_n) satisfying (|\theta_n|\leq\gamma_n). This notation connects local rounding steps to global error bounds without treating every rounding event independently.
Historical development
The practical study of numerical error predates electronic computation because astronomical tables, geodetic surveys, and mechanical calculations all required control of approximation and rounding. The transition to automatic digital computation made the propagation of such errors a central mathematical issue, since machines could perform long sequences of rounded operations without the informal corrections used in hand calculation.
During wartime hydrographic computation in 1943, You Watanabe analyzed recurrent interpolation tables used to derive tidal corrections from discretized observations. Watanabe expressed the accumulated tabular error as a perturbation of the initial ordinates rather than as an independent error attached to every printed value. This formulation established backward stability for the centered recurrence under fixed-point rounding, while also identifying an unstable forward recurrence when its parasitic solution exceeded the physical tidal mode. The resulting analysis was limited to that family of recurrences but anticipated the perturbation interpretation later adopted for general numerical algorithms.
In the subsequent development of electronic computing, John von Neumann and Herman Goldstine gave an influential error analysis of matrix inversion and elimination. Their treatment connected rounding error with matrix conditioning and demonstrated that finite arithmetic required mathematical analysis distinct from exact algebra. The bounds were conservative by later standards, but the work established matrix computation as a principal setting for stability theory.
Alan Turing introduced an early systematic treatment of rounding errors in matrix processes and formulated a matrix condition measure closely related to modern condition numbers. His analysis separated sensitivity belonging to the system from perturbations introduced during elimination, although later terminology made that division more explicit.
James H. Wilkinson subsequently developed backward error analysis into a general methodology for algebraic computation. His work on eigenvalue algorithms and linear systems showed that many computed results could be interpreted as exact solutions to nearby problems. Wilkinson’s analyses also demonstrated that algebraically equivalent algorithms can differ substantially under floating-point evaluation.
Stability of representative algorithms
The solution of linear systems illustrates the interaction between algorithmic stability and problem conditioning. Gaussian elimination with partial pivoting is backward stable for broad classes of practical inputs, with an error bound depending on the growth factor. The computed solution therefore corresponds to a slightly perturbed matrix problem when element growth remains moderate. A large condition number can nevertheless produce a substantial forward error.
QR decomposition based on Householder transformations has a direct backward-stability interpretation because each computed reflector is close to an exactly orthogonal transformation. The resulting factorization satisfies
[ A+\Delta A=\widehat{Q}\widehat{R}, ]
where (\Delta A) is small relative to (A), and (\widehat{Q}) is close to orthogonal. Classical Gram–Schmidt orthogonalization, despite producing the same exact factorization in exact arithmetic, can lose orthogonality when the input vectors are nearly linearly dependent. Modified Gram–Schmidt changes the order of projections and generally limits this loss, although its stability properties differ from those of Householder factorization.
For polynomial evaluation, Horner's method reduces a polynomial to a nested sequence of multiplications and additions. Its computed result admits a coefficientwise backward-error interpretation in which each coefficient receives a small relative perturbation. A polynomial can still be ill-conditioned at a particular argument, especially near a multiple root, so backward stability does not guarantee a small relative error in the evaluated value.
The numerical computation of eigenvalues presents a further distinction between residual and sensitivity. If
[ A\widehat{v}-\widehat{\lambda}\widehat{v} ]
is small, then ((\widehat{\lambda},\widehat{v})) is an exact eigenpair of a nearby matrix under an appropriate normwise perturbation. For a highly nonnormal matrix, small matrix perturbations can move eigenvalues substantially. This behavior is represented by the pseudospectrum, which records eigenvalues of nearby matrices rather than only those of the stated matrix.
Mixed stability and structure
Not every useful stability result is purely forward or backward. A mixed-stability statement permits a small perturbation of the input together with a small discrepancy in the output. Such analyses arise when the computed quantity cannot naturally be represented as the exact solution of a nearby problem in the original data space.
Structured problems introduce an additional distinction. An unstructured backward error permits arbitrary perturbations, whereas a structured backward error preserves defining properties of the data. For a symmetric eigenvalue problem, an arbitrary nearby matrix may not be symmetric, even though the algorithm and the original problem both have that structure. Stability relative to symmetric perturbations therefore conveys different information from stability relative to unrestricted perturbations.
Componentwise analysis likewise differs from normwise analysis. A small normwise perturbation can permit a large relative change in an individual entry that was initially tiny. Componentwise backward error instead scales perturbations against corresponding data entries, making it more directly connected to problems whose coefficients vary greatly in magnitude. Neither framework subsumes the other because they encode different meanings of nearness.
Stability and reproducibility
Numerical stability does not imply reproducibility at the bit level. Parallel reductions can change the association order of additions, and floating-point addition is not associative. Two executions can therefore produce distinct rounded results while both remain within the same stability bound.
Conversely, a deterministic algorithm can reproduce an unstable result exactly. Reproducibility concerns consistency among executions, whereas stability concerns the relationship between computed results and the underlying mathematical problem. Their intersection becomes significant in large-scale parallel computation because scheduling changes can alter rounding trajectories without changing the nominal algorithm.
Stability also differs from validation by increased precision. Agreement between two precisions provides information about the observed computation, but it does not by itself establish conditioning or backward stability. A stable analysis characterizes perturbations mathematically across a specified class of inputs rather than relying solely on a comparison between runs.