10. Numerical Methods and Computational Stability

Introduction and Context

Machine-learning systems are often described through exact mathematical expressions: a gradient is computed, a linear system is solved, a Hessian is factorized, or a probability is evaluated. On a computer, however, these operations are performed with finite precision, limited memory, and finite computational time. The implemented algorithm therefore produces an approximation rather than an exact mathematical result. Numerical methods provide the principles and algorithms needed to make that approximation accurate, efficient, and reliable.

This issue is central to machine learning. Least-squares regression may require solving a large linear system; neural-network training repeatedly applies gradient-based updates; probabilistic models evaluate exponentials and logarithms; principal component analysis depends on eigenvalue computation; and scientific or industrial ML systems may process matrices with millions or billions of entries. A numerically unstable implementation can produce overflow, inaccurate gradients, slow convergence, or apparently inexplicable model failure even when the underlying mathematical formulation is correct.

Numerical reliability has three related dimensions. Approximation error arises because a continuous or exact problem is replaced by a finite computational procedure. Conditioning describes how sensitive the mathematical problem is to perturbations in its input. Stability describes how sensitively an algorithm behaves in the presence of computational errors. A well-designed ML implementation considers all three.

Learning Objectives

By the end of this chapter, a student should be able to:

1. distinguish truncation, discretization, round-off, and propagation errors;

2. explain floating-point representation, cancellation, overflow, underflow, and machine epsilon;

3. distinguish the conditioning of a problem from the stability of an algorithm;

4. analyse convergence using error measures, residuals, rates, and stopping criteria;

5. apply direct and iterative methods to linear systems arising in ML;

6. select appropriate approaches for numerical differentiation and gradient computation;

7. explain the role of scaling, preconditioning, sparse computation, and low-rank methods;

8. diagnose numerical failure in practical machine-learning pipelines.

Approximation, Error, and the Meaning of a Numerical Answer

Suppose the exact value of a quantity is x and its computed approximation is x-hat. The absolute error is

|x - x-hat|,

while the relative error, when x is nonzero, is

|x - x-hat| / |x|.

Relative error is often more informative because an absolute error of 0.001 is substantial when the true value is 0.002 but negligible when the true value is 10,000.

A numerical algorithm may introduce several kinds of error. Truncation error results when an infinite process is stopped after finitely many terms. For example, replacing an infinite Taylor series with its first few terms creates truncation error. Discretization error occurs when a continuous domain is represented by a finite grid or a finite set of samples. Round-off error is caused by representing numbers with finite precision. Model or data error arises when the mathematical problem itself is only an approximation of the real phenomenon.

Errors can be classified as forward or backward. Forward error measures the distance between the computed answer and the exact answer. Backward error asks a different question: for which slightly perturbed input would the computed answer be exact? An algorithm is often considered numerically stable if its output is the exact solution to a nearby problem. This viewpoint is powerful because small backward error, combined with a well-conditioned problem, generally implies small forward error.

Error propagation explains why a small local error may become significant. If a computation repeatedly transforms a value, errors can be amplified by derivatives, matrix norms, or recursive feedback. In an iterative ML algorithm, errors may also influence later parameters, changing the trajectory of optimization rather than merely perturbing one final result.

Floating-Point Arithmetic: The Computational Number System

Computers commonly represent a real number in floating-point form:

x = (-1)^s × m × b^e,

where s is a sign bit, m is a significand or mantissa, b is the base, and e is an exponent. In a finite format, only a limited number of significant digits and exponent values are available. Most real numbers therefore cannot be represented exactly.

When an operation produces a value that is rounded to the available format, the result may be written as fl(x). For a standard rounded floating-point system, the relative representation error is typically bounded by a small quantity related to machine epsilon, often denoted by εmach. Machine epsilon is the smallest number such that, in the working precision, 1 + εmach is distinguishable from 1. Its numerical value depends on the format; it is not a universal constant.

Several phenomena are especially important in ML implementations:

• Rounding: the exact result lies between representable numbers and is replaced by the nearest available value.
• Overflow: a magnitude exceeds the largest representable value, producing infinity or an exception.
• Underflow: a nonzero magnitude becomes too small to represent normally and may be rounded toward zero.
• Catastrophic cancellation: subtracting nearly equal numbers destroys leading significant digits.
• Absorption: adding a very small number to a much larger number may leave the larger number unchanged.

The expression sqrt(x² + y²) can overflow even when the final mathematical result is representable, because x² and y² may overflow separately. A scaled formulation, such as |x|sqrt(1 + (y/x)²) when |x| is dominant, is safer. Similarly, direct computation of exp(z) can overflow for large z, while softmax should be evaluated by subtracting the maximum input first:

softmax(z)i = exp(zi - zmax) / sumj exp(zj - zmax).

The subtraction does not change the mathematical result but substantially improves numerical safety. Log-sum-exp uses the same principle:

log sumj exp(zj) = zmax + log sumj exp(zj - zmax).

These transformations illustrate a general rule: algebraically equivalent expressions need not be computationally equivalent.

Summation is another important example. Adding many numbers in arbitrary order can accumulate round-off error, particularly when positive and negative terms differ greatly in magnitude. Pairwise summation, sorting by magnitude, or compensated methods such as Kahan summation can improve accuracy. In distributed ML, reduction order may vary across devices, so exact reproducibility may require deterministic reduction strategies.

Conditioning and Stability

The condition number measures the sensitivity of a problem to small changes in its input. For a scalar function y = f(x), a local relative condition measure is approximately

κ(x) = |x f'(x) / f(x)|,

when the quantities are meaningful. A large condition number indicates an ill-conditioned problem: even highly accurate input can lead to a substantially changed output. An ill-conditioned problem cannot be made reliable merely by selecting a different implementation; the underlying mathematical task itself is sensitive.

For a nonsingular matrix A, the usual condition number in a chosen norm is

κ(A) = ||A|| ||A⁻¹||.

If κ(A) is large, solving Ax = b is sensitive to perturbations in A or b. The condition number of a matrix depends on the norm and may be estimated through singular values. In the 2-norm, κ₂(A) = σmax/σmin. Nearly dependent columns make σmin small and the problem poorly conditioned.

Algorithmic stability is different. A stable algorithm controls the errors introduced during computation. The normal-equation approach to least squares, for example, solves

XᵀXβ = Xᵀy.

Although mathematically valid, forming XᵀX squares the condition number approximately: κ(XᵀX) ≈ κ(X)². QR factorization or singular value decomposition is usually more stable because it avoids this amplification. Thus, a difficult problem may be handled reasonably by a stable algorithm, while an avoidable unstable algorithm may fail even on a well-conditioned problem.

Scaling is often the first remedy. If one feature is measured in dollars and another in millimetres, their numerical ranges may differ by many orders of magnitude. Standardization, nondimensionalization, and appropriate rescaling can improve conditioning and optimization geometry. Scaling must be fitted using training data only when it is part of a predictive pipeline, otherwise information leakage may occur.

Convergence and Stopping Decisions

An iterative method produces a sequence x⁽⁰⁾, x⁽¹⁾, … intended to approach a solution x. The method converges if x⁽ᵏ⁾ approaches x as k tends to infinity. The error is e⁽ᵏ⁾ = x⁽ᵏ⁾ − x. In practice, x is unknown, so stopping decisions rely on observable quantities such as successive-step size, objective improvement, gradient norm, or residual.

A method has linear convergence if, asymptotically,

||e⁽ᵏ⁺¹⁾|| ≤ c||e⁽ᵏ⁾||, 0 < c < 1.

It has superlinear convergence if the ratio of successive errors tends to zero, and quadratic convergence if the next error is approximately proportional to the square of the current error. Newton-type methods can be quadratically convergent near a regular solution, but they may be expensive or unreliable far from it.

For a linear system Ax = b, the residual is r⁽ᵏ⁾ = b − Ax⁽ᵏ⁾. A small residual is necessary for a good solution but does not always guarantee a small solution error when A is ill-conditioned. A practical stopping rule may combine a relative residual condition,

||r⁽ᵏ⁾|| / ||b|| ≤ tolerance,

with a step-size condition and a maximum iteration limit. In optimization, stopping solely when the loss changes little can be misleading: the objective may be flat, poorly scaled, or affected by noise. Gradient norms, validation behaviour, and parameter changes should also be considered.

Finite precision creates a convergence floor. After sufficient iterations, round-off and accumulated error may prevent further improvement. An algorithm that continues indefinitely may waste resources or even become less accurate. Stopping criteria should reflect data scale, precision, conditioning, and the consequences of false convergence.

Solving Linear Systems in Machine Learning

Linear systems appear in regression, Gaussian models, kernel methods, constrained optimization, numerical PDEs, and second-order training methods. The preferred method depends on matrix structure, size, sparsity, and accuracy requirements.

Direct methods

Gaussian elimination with pivoting transforms a system into triangular form. Partial pivoting, which exchanges rows to use a sufficiently large pivot, reduces instability caused by division by small numbers. For symmetric positive-definite matrices, Cholesky factorization A = LLᵀ is efficient and usually stable. QR factorization A = QR is especially useful for least squares, while singular value decomposition A = UΣVᵀ provides robust treatment of rank deficiency and exposes ill-conditioning.

Direct methods can be computationally expensive for very large matrices and may destroy sparsity through fill-in. Nevertheless, factorization is often valuable when the same matrix must be solved against many right-hand sides.

Iterative methods

Iterative methods are attractive when matrices are large, sparse, or accessible only through matrix-vector products. The Jacobi and Gauss–Seidel methods update components using previously computed approximations. They are conceptually simple but may converge slowly or fail unless matrix conditions are favourable.

For symmetric positive-definite systems, conjugate gradient (CG) constructs mutually conjugate search directions and often converges much faster than basic stationary iterations. Its rate depends on the condition number and eigenvalue distribution. For general systems, GMRES minimizes the residual over a growing Krylov subspace, though memory costs may increase and restarting may be required. Stochastic gradient methods can also be viewed as iterative numerical solvers for optimization problems, although their gradients are noisy and their convergence analysis differs from deterministic linear algebra.

Preconditioning

A preconditioner M approximates A in a way that makes the transformed system easier to solve. Instead of solving Ax = b directly, one may solve M⁻¹Ax = M⁻¹b. A good preconditioner reduces the effective condition number without being so expensive that it eliminates the benefit of iteration. Examples include diagonal scaling, incomplete factorization, and problem-specific approximations. In ML, adaptive parameter scaling, feature normalization, and approximate curvature matrices play related roles.

Numerical Differentiation and Automatic Differentiation

Derivatives are fundamental to optimization, sensitivity analysis, and scientific ML. The forward finite-difference approximation is

f'(x) ≈ [f(x+h) − f(x)] / h.

Its truncation error is generally O(h), but choosing h too small increases cancellation and round-off error. Choosing h too large increases truncation error. A central difference,

f'(x) ≈ [f(x+h) − f(x−h)]/(2h),

has truncation error O(h²) under suitable smoothness assumptions, but requires two function evaluations.

Finite differences are useful for checking gradients, but they are often inefficient for high-dimensional models and can be unreliable when the function is noisy, discontinuous, or badly scaled. Automatic differentiation (AD) applies the chain rule to elementary computational operations. Forward mode is efficient when the number of input dimensions is small relative to outputs; reverse mode is efficient for scalar losses with many parameters and is the basis of backpropagation. AD is not the same as symbolic differentiation and does not eliminate floating-point error, but it usually avoids the step-size trade-off of finite differences.

A practical gradient check compares an AD gradient with a central-difference estimate on a small, well-scaled test problem. It should be performed away from nondifferentiable points and with a carefully selected step size. Agreement to machine precision is not expected; the comparison is diagnostic rather than a proof.

Scalable Matrix Computation

Modern ML frequently uses matrices too large for dense storage. Sparse matrices store only nonzero entries, reducing memory and computation when sparsity is genuine. Sparse matrix-vector multiplication is central to graph neural networks, recommender systems, and scientific computing. However, sparse formats introduce indexing overhead and may perform poorly when the matrix becomes moderately dense.

Low-rank approximation is another important strategy. Truncated SVD approximates a matrix by its leading singular components, reducing storage and accelerating downstream calculations. Randomized numerical linear algebra uses random projections to estimate these components efficiently for large data sets. The approximation rank must be selected according to error tolerance, spectral decay, and the task’s sensitivity to discarded components.

Blocked algorithms improve cache use and parallelism. GPUs and distributed systems accelerate matrix operations, but they also introduce communication costs, nondeterministic reduction order, and mixed-precision effects. A scalable implementation therefore considers not only asymptotic complexity but also memory bandwidth, data movement, synchronization, and numerical reproducibility.

Mixed-precision computation can accelerate training by using lower precision for matrix operations while retaining higher precision for sensitive accumulations, normalization statistics, or parameter updates. Loss scaling helps prevent underflow in low-precision gradients. Such methods require monitoring for overflow, gradient distortion, and unstable updates.

An Illustrative Case Study: An Unstable Logistic-Regression Pipeline

Background

A financial-analytics team builds a binary classifier for loan default using logistic regression. The data contain income measured in annual currency units, age in years, and a highly correlated set of transaction variables. The team implements the model using gradient descent and evaluates the sigmoid function directly.

Situation and problem

During training, the loss occasionally becomes NaN. In other runs, the loss decreases very slowly, and the estimated coefficients vary substantially with small changes in the data. The team initially suspects a software defect.

Analysis

Inspection reveals that the transaction variables have widely different scales and strong multicollinearity. The design matrix is ill-conditioned. The implementation computes sigmoid(z) = 1/(1 + exp(−z)); for large negative z, exp(−z) overflows. The loss also evaluates log(sigmoid(z)) directly, creating log(0) when probabilities underflow. Finally, the team solves the regularized least-squares subproblem through normal equations, amplifying the condition number.

Approach and solution

The features are standardized using training-set statistics, and redundant variables are examined. The stable softplus and log-sum-exp formulations are used for the loss. The optimization code clips or otherwise safely handles extreme intermediate values, while retaining mathematically meaningful gradients. For a Newton-style variant, a Cholesky or QR-based solve is used rather than explicitly forming an inverse or relying on normal equations. The team monitors gradient norms, relative loss changes, validation loss, and parameter norms. A small gradient check confirms the AD implementation on a controlled data subset.

Outcome

The NaN failures disappear, convergence becomes more consistent, and coefficient estimates become less sensitive to small perturbations. Regularization improves both statistical generalization and numerical behaviour, although it does not eliminate the need for stable computation.

Learning points

The case demonstrates that numerical failure can arise from the interaction of representation, conditioning, algorithm choice, and stopping rules. Stable elementary functions, feature scaling, appropriate factorizations, regularization, and diagnostic monitoring must be treated as parts of one engineering problem.

Common Failure Modes and Professional Practice

Several mistakes recur in numerical ML. Explicitly computing a matrix inverse is usually less stable and less efficient than solving a system. Forming XᵀX unnecessarily can square conditioning. Ignoring units and feature scales can make optimization appear defective. Comparing floating-point numbers with exact equality can cause incorrect branching. Using a fixed learning rate without considering curvature can produce divergence. Treating a small residual as proof of an accurate solution is unsafe for ill-conditioned systems.

Reliable practice includes the following principles:

• scale variables and inspect ranges before optimization;
• use stable mathematical identities for exponentials, logarithms, norms, and probabilities;
• prefer factorizations or iterative solvers to explicit inverses;
• estimate conditioning when sensitivity matters;
• use pivoting, regularization, or preconditioning appropriately;
• monitor residuals, gradients, objective values, and validation metrics;
• test algorithms on synthetic cases with known answers;
• use higher precision selectively for verification and sensitive stages;
• document tolerances, random seeds, precision modes, and solver settings;
• distinguish numerical instability from data quality, model misspecification, and statistical uncertainty.

Numerical reliability also has professional and governance implications. In high-stakes domains, an unstable computation may produce inconsistent decisions that are difficult to audit. Reproducibility requires recording software versions, hardware, precision, solver configuration, and preprocessing statistics. Approximation error should be included in uncertainty analysis when model outputs support financial, medical, legal, or safety-related decisions.

Current Developments

Research and practice increasingly combine numerical linear algebra with automatic differentiation, probabilistic programming, and hardware-aware computation. Differentiable solvers allow optimization through numerical algorithms, but their gradients may inherit solver inaccuracies. Neural operators and scientific ML depend on stable discretization and scalable iterative methods. Randomized algorithms reduce the cost of large matrix approximations, while communication-avoiding and asynchronous methods address distributed-scale computation.

The growth of specialized accelerators has also made mixed precision a standard design consideration. Future systems will increasingly select precision adaptively, balancing speed, energy, accuracy, and reproducibility. The central principle remains unchanged: computational scale must not be achieved by ignoring numerical behaviour.

Chapter Summary

Numerical methods turn mathematical ML formulations into reliable computational procedures. Approximation error includes truncation, discretization, and round-off effects, while floating-point arithmetic creates overflow, underflow, cancellation, and order-dependent summation. Conditioning measures the sensitivity of a problem; stability measures the behaviour of an algorithm. These concepts must be distinguished when diagnosing error.

Convergence analysis guides iterative optimization and linear-system solvers, but stopping rules must account for residuals, gradients, conditioning, and finite-precision floors. QR, Cholesky, SVD, pivoting, conjugate gradient, GMRES, and preconditioning provide different tools for different matrix structures. Stable formulas, automatic differentiation, feature scaling, sparse and low-rank computation, mixed precision, and careful monitoring are essential for scalable ML. Reliable numerical practice is simultaneously a technical, professional, and governance responsibility.

Key Terms

Machine epsilon: The smallest positive quantity that distinguishes 1 from a nearby floating-point number in a specified precision.

Round-off error: Error caused by representing or operating on numbers with finite precision.

Truncation error: Error introduced by terminating an otherwise infinite mathematical process.

Discretization error: Error caused by replacing a continuous problem with a finite representation.

Forward error: Difference between a computed result and the exact result.

Backward error: Perturbation of the input for which the computed result is exact.

Condition number: A measure of the sensitivity of a problem’s output to input perturbations.

Numerical stability: The ability of an algorithm to limit the effects of computational errors.

Catastrophic cancellation: Severe loss of significant digits when nearly equal quantities are subtracted.

Overflow: A result whose magnitude exceeds the representable range.

Underflow: A result whose magnitude is too small for normal representation.

Residual: The discrepancy b − Ax for an approximate solution of Ax = b.

Convergence: The approach of an iterative sequence toward a target solution.

Preconditioner: An operator that transforms a problem into one more favourable for iterative solution.

Cholesky factorization: Factorization of a symmetric positive-definite matrix as LLᵀ.

QR factorization: Representation of a matrix as an orthogonal matrix times an upper triangular matrix.

Singular value decomposition: Factorization A = UΣVᵀ revealing rank and conditioning.

Automatic differentiation: Systematic application of the chain rule through a computational graph.

Mixed precision: Use of more than one numerical precision within a computation.

Kahan summation: A compensated summation technique that reduces accumulated floating-point error.

Review Questions

Short Answer Questions

1. Distinguish forward error from backward error.

2. What is machine epsilon, and why does it depend on numerical format?

3. Explain the difference between conditioning and algorithmic stability.

4. Why can a small residual fail to imply a small solution error?

5. State one advantage of automatic differentiation over finite differences.

Long Answer Questions

1. Explain floating-point arithmetic and analyse overflow, underflow, cancellation, and absorption with ML-related examples.

2. Compare normal equations, QR factorization, and SVD for least-squares problems.

3. Discuss convergence rates and stopping criteria for iterative numerical methods.

4. Explain the role of scaling and preconditioning in improving computational reliability.

5. Describe how sparse, low-rank, randomized, and mixed-precision methods support scalable matrix computation.

Analytical / Application Questions

1. A softmax implementation produces NaN values when logits have large magnitude. Diagnose the numerical cause and propose a stable reformulation.

2. Two nearly identical features cause regression coefficients to change greatly across samples, although training loss is small. Analyse the likely conditioning issue and recommend computational and statistical remedies.

3. A conjugate-gradient solver reduces the residual rapidly at first but then stagnates. Develop a diagnostic plan involving conditioning, preconditioning, precision, stopping tolerances, and matrix properties.

22. Multiple Choice Questions with Answers

Question 1

Which statement best distinguishes conditioning from numerical stability?

A. Conditioning concerns algorithm speed, whereas stability concerns model accuracy.

B. Conditioning describes problem sensitivity to input perturbations, whereas stability describes an algorithm’s response to computational errors.

C. Conditioning applies only to nonlinear problems, whereas stability applies only to linear problems.

D. They are two names for the same property.

Correct Answer: B. Conditioning describes problem sensitivity to input perturbations, whereas stability describes an algorithm’s response to computational errors.

Question 2

Why is subtracting the maximum logit before computing softmax numerically useful?

A. It changes the predicted class deliberately.

B. It makes all logits positive.

C. It prevents unnecessary overflow or underflow without changing the mathematical softmax values.

D. It removes the need for normalization.

Correct Answer: C. It prevents unnecessary overflow or underflow without changing the mathematical softmax values.

Question 3

Which operation can approximately square the condition number of a least-squares problem?

A. Using QR factorization

B. Using singular value decomposition

C. Forming XᵀX before solving

D. Standardizing the columns of X

Correct Answer: C. Forming XᵀX before solving

Question 4

For a symmetric positive-definite linear system, which iterative method is especially appropriate?

A. Conjugate gradient

B. Newton interpolation

C. Forward finite difference

D. Kahan summation

Correct Answer: A. Conjugate gradient

Question 5

What is catastrophic cancellation?

A. Overflow during matrix multiplication

B. Loss of significant digits when nearly equal numbers are subtracted

C. Failure of a sparse matrix to remain sparse

D. The exact equality of two floating-point numbers

Correct Answer: B. Loss of significant digits when nearly equal numbers are subtracted

Question 6

A small residual in Ax = b does not necessarily imply a small solution error primarily when:

A. A is well-conditioned

B. b is stored in double precision

C. A is ill-conditioned

D. the system is triangular

Correct Answer: C. A is ill-conditioned

Question 7

Which method generally has second-order truncation error for a smooth scalar function?

A. Forward finite difference

B. Central finite difference

C. Direct symbolic substitution

D. Jacobi iteration

Correct Answer: B. Central finite difference

Question 8

What is the principal purpose of a preconditioner?

A. To increase the dimension of the system

B. To transform the system into one that is easier for an iterative solver

C. To replace all matrix-vector products with dense inverses

D. To guarantee exact arithmetic

Correct Answer: B. To transform the system into one that is easier for an iterative solver

Question 9

In reverse-mode automatic differentiation, computation is especially efficient when:

A. There are many outputs and very few inputs

B. There is a scalar loss and many model parameters

C. The function is discontinuous everywhere

D. No chain rule is required

Correct Answer: B. There is a scalar loss and many model parameters

Question 10

Which practice is most appropriate for a numerically reliable ML implementation?

A. Explicitly compute matrix inverses whenever possible

B. Use exact equality tests for floating-point stopping conditions

C. Monitor residuals or gradients, scale variables, and use stable factorizations or solvers

D. Continue iterations indefinitely to overcome finite-precision limits

Correct Answer: C. Monitor residuals or gradients, scale variables, and use stable factorizations or solvers