Floating Point and Conditioning

Floating Point Representation

Real numbers form an infinite continuum; a computer's memory is finite. The IEEE 754 double-precision standard (the default in NumPy, MATLAB, Julia, and most scientific software) represents every nonzero number in the form:

IEEE 754 Double Precision

x = \pm 2^e \cdot (1 + m), \quad m \in [0,1), \quad e \in [-1022,\, 1023]

A 64-bit double stores: 1 sign bit, 11 exponent bits (biased by 1023), and 52 mantissa bits. The mantissa bits represent the fractional part m ∈ [0,1), so the significand 1.m lies in [1,2). Together this gives about 15–17 significant decimal digits. Numbers outside [≈5×10⁻³²⁴, ≈1.8×10³⁰⁸] overflow or underflow to ±∞ or ±0.

Because the mantissa has a finite number of bits, most real numbers cannot be represented exactly — they must be rounded to the nearest representable float. The key quantity measuring this granularity is machine epsilon ε_mach: the smallest positive number such that 1 + ε_mach ≠ 1 in floating point arithmetic. For double precision, ε_mach = 2⁻⁵² ≈ 2.2 × 10⁻¹⁶.

Machine epsilon characterizes the relative rounding error in a single floating-point operation. When you compute fl(x), the floating-point representation of a real x, the relative error satisfies |fl(x) − x| / |x| ≤ ε_mach. This bound is tight: numbers of the form 1 + δ where δ ≈ ε_mach/2 round to exactly 1.0. The set of representable doubles has gaps that grow with magnitude — the spacing between adjacent doubles near 10¹⁵ is approximately 10¹⁵ · ε_mach ≈ 0.2, so you cannot distinguish 10¹⁵ from 10¹⁵ + 0.1 in double precision.

Catastrophic Cancellation

The most dangerous source of floating-point error is catastrophic cancellation: subtracting two nearly equal numbers. If a = 1.000000000000001 and b = 1.000000000000000 (both stored to 16 significant digits), then a − b = 1 × 10⁻¹⁵, which has only one significant digit — fifteen digits of precision were annihilated in a single subtraction. The relative error in the result can be enormous even though neither input had any error. The classic example is evaluating the quadratic formula: for a nearly zero discriminant b² − 4ac, direct subtraction is numerically disastrous; the numerically stable alternative uses the relation x₁x₂ = c/a to compute the smaller root after finding the larger one accurately.

A practical rule: if you compute a − b and |a − b| ≪ |a|, you have lost approximately log₁₀(|a|/|a−b|) decimal digits of precision. Restructuring the computation to avoid such subtractions — using the formula log(1+x) instead of log(1+x) computed naively for small x, or using Kahan summation for long sums — is standard practice in numerical software.

The Condition Number

Consider the linear system Ax = b. Suppose the right-hand side b has a small perturbation δb (perhaps from measurement noise or rounding). The perturbed system A(x + δx) = b + δb gives δx = A⁻¹δb. The relative change in the solution compared to the relative change in the data is bounded by:

Condition Number Bound

\frac{\|\delta x\|}{\|x\|} \leq \kappa(A)\,\frac{\|\delta b\|}{\|b\|}, \quad \kappa(A) = \|A\|\cdot\|A^{-1}\|

The condition number κ(A) = ‖A‖ · ‖A⁻¹‖ is the worst-case amplification factor for relative errors. Using the 2-norm, κ₂(A) = σ_max / σ_min — the ratio of the largest to smallest singular value. A condition number of 10⁶ means a relative perturbation of 10⁻¹⁶ in b (machine epsilon) can cause a relative error of up to 10⁻¹⁰ in x — losing 6 decimal digits of accuracy.

The condition number κ(A) does not depend on b; it is a property of the matrix A alone. Intuitively, it measures how "squished" the matrix is: a well-conditioned matrix maps a sphere to a (nearly) sphere; an ill-conditioned matrix maps a sphere to a very elongated ellipsoid, with the ratio of longest to shortest axis equal to κ(A). When σ_min ≈ 0, the matrix is nearly singular, and κ(A) → ∞ — the matrix almost collapses a whole direction to zero, and any noise in that direction in b gets massively amplified in x.

The same analysis applies to perturbations in A itself. If we solve (A + δA)x̃ = b, the relative error in x̃ compared to the true solution is bounded by approximately κ(A) · ‖δA‖/‖A‖. In floating-point arithmetic, the matrix A is itself rounded when stored, so δA is roughly ε_mach · ‖A‖. The backward error (the size of the perturbation needed to make the computed answer exact) and forward error (the actual error in the answer) are related by the condition number: forward error ≲ κ(A) × backward error.

Computing the Condition Number

In MATLAB/NumPy, cond(A) computes κ₂(A) = σ_max/σ_min via SVD. This costs O(mn min(m,n)) for an m×n matrix — often too expensive just to diagnose a system. A cheaper estimate is the reciprocal condition number estimate rcond(A), which uses an O(n²) heuristic based on LU decomposition. As a rule of thumb: if κ(A) ≈ 10^k, you lose approximately k decimal digits of accuracy relative to what the arithmetic could provide. For double precision with 16 digits, a condition number of 10¹² leaves only ~4 significant digits in the solution — often unacceptably poor.

Condition Number via SVD

\kappa_2(A) = \frac{\sigma_{\max}}{\sigma_{\min}}

For a square matrix A with SVD A = UΣVᵀ, the 2-norm condition number is the ratio of the largest to smallest singular value. For a symmetric positive definite matrix (eigenvalues λ₁ ≥ ... ≥ λₙ > 0), κ₂(A) = λ_max/λ_min. The Frobenius norm condition number uses ‖A‖_F instead of ‖A‖₂ — more expensive but sometimes more informative for structured matrices.

Well-Conditioned vs. Ill-Conditioned Problems

A well-conditioned problem is one where small changes in the input lead to small changes in the output: κ(A) is small (say, κ < 100 for most purposes). An ill-conditioned problem has κ(A) ≫ 1 — a tiny perturbation in the data can cause a huge change in the solution. Ill-conditioning is a property of the mathematical problem, not of the algorithm used to solve it. A stable algorithm applied to an ill-conditioned problem will still return a poor answer — it just won't make things worse than the conditioning dictates.

The canonical example of an ill-conditioned matrix family is the Hilbert matrix: Hᵢⱼ = 1/(i+j−1). The 10×10 Hilbert matrix has condition number κ₂(H₁₀) ≈ 1.6 × 10¹³ — solving H₁₀x = b with double-precision arithmetic yields at most 16 − 13 = 3 reliable decimal digits in x, no matter which algorithm you use. The Hilbert matrix arises naturally in polynomial regression (fitting a polynomial of degree n−1 to n equally-spaced points using the monomial basis), which is why polynomial regression with high degree is numerically treacherous.

Sources of Ill-Conditioning in Practice

Ill-conditioned linear systems appear throughout applied mathematics and engineering:

Multicollinear features: In regression, if two columns of the design matrix X are nearly parallel (highly correlated features), XᵀX has a tiny eigenvalue, and its condition number is large. Ridge regularization (adding λI) lifts all eigenvalues by λ, capping the condition number at (σ_max² + λ)/λ.
Polynomial fitting in the monomial basis: The Vandermonde matrix V with Vᵢⱼ = xᵢʲ⁻¹ is extremely ill-conditioned for high degrees and closely-spaced nodes. Switching to an orthogonal polynomial basis (Legendre, Chebyshev) reduces the condition number dramatically.
Discretized differential equations: Finite-difference or finite-element discretizations of PDEs produce matrices whose condition number grows as O(1/h²) where h is the mesh spacing. Preconditioning — multiplying by an approximate inverse — is essential to make iterative solvers converge.
Near-singular systems: If a physical system is near resonance, near bifurcation, or nearly degenerate (e.g., two circuits with nearly equal frequencies), the matrix governing the system is nearly singular and numerically sensitive.

Numerical Stability of Algorithms

Two algorithms that solve the same problem mathematically can have very different numerical behavior. A numerically stable algorithm returns a result that is the exact answer to a slightly perturbed problem (small backward error). An unstable algorithm may introduce errors much larger than the condition number requires — it amplifies rounding errors unnecessarily.

Gaussian elimination without pivoting is unstable: if a small pivot element is encountered, dividing by it amplifies rounding errors explosively. Partial pivoting — always using the largest element in the current column as the pivot — prevents this and makes LU decomposition backward stable in practice. Similarly, computing A⁻¹ explicitly and then forming A⁻¹b to solve Ax = b is both wasteful (O(n³) for the inverse vs. O(n³/3) for LU) and less stable than directly solving via LU or QR decomposition.

Backward Error Analysis

(A + \delta A)\tilde{x} = b + \delta b, \quad \frac{\|\delta A\|}{\|A\|},\, \frac{\|\delta b\|}{\|b\|} = O(\varepsilon_{\text{mach}})

A backward-stable algorithm for solving Ax = b finds x̃ that exactly satisfies (A + δA)x̃ = b + δb, where ‖δA‖/‖A‖ and ‖δb‖/‖b‖ are of order ε_mach. The forward error — how much x̃ differs from the true x — is then bounded by κ(A) × ε_mach. Backward stability means the algorithm does not introduce errors beyond what the problem's conditioning dictates.

The QR decomposition approach to least-squares (solve Rw = Qᵀy instead of forming XᵀX) is a prime example of stability improvement. The normal equations XᵀXw = Xᵀy square the condition number: κ(XᵀX) = κ(X)². If κ(X) = 10⁸ (a moderately ill-conditioned matrix), κ(XᵀX) = 10¹⁶ — essentially singular in double precision. The QR approach works directly with X and has condition number κ(X) = 10⁸, allowing 8 accurate digits where the normal equations would give 0. This is why all serious software (scipy.linalg.lstsq, NumPy's lstsq, R's lm()) uses QR or SVD internally, never the normal equations directly.

When to Worry About Numerical Accuracy

A practical checklist for numerical reliability:

Check the condition number first: if κ(A) > 10^(16/p) where p is the desired significant digits, you cannot reliably achieve p-digit accuracy in double precision.
Avoid forming AᵀA explicitly for least-squares; use QR or SVD directly.
Never invert a matrix to solve a linear system; use a factored form (LU, QR, Cholesky) and solve by substitution.
Normalize your data in regression: standardizing columns of X brings its singular values closer together and reduces κ(XᵀX).
Use regularization (ridge, Tikhonov) when the problem is inherently ill-conditioned; the regularization parameter λ caps the effective condition number.
Monitor residuals: after computing x̃, evaluate ‖Ax̃ − b‖/‖b‖ (the relative residual). A small residual means the backward error is small; a large forward error despite small residual means the problem is ill-conditioned.

Engineering Implication

In NumPy: np.finfo(float).eps gives ε_mach ≈ 2.2e-16. np.linalg.cond(A) computes κ₂(A). np.linalg.solve(A,b) uses LU with partial pivoting (backward stable). np.linalg.lstsq(A,b) uses SVD and sets a threshold on small singular values. MATLAB's backslash operator automatically chooses LU or QR based on matrix structure. When your system seems to give garbage results despite a correct algorithm, compute the condition number first — the answer often has more to do with the problem's inherent sensitivity than with any bug in your code.

Key Takeaways

IEEE 754 double precision has machine epsilon ε_mach ≈ 2.2×10⁻¹⁶ — the relative rounding error in a single float operation. Catastrophic cancellation (subtracting nearly equal numbers) can destroy many digits at once. The condition number κ(A) = σ_max/σ_min bounds how much relative error in b is amplified into error in the solution of Ax = b: ‖δx‖/‖x‖ ≤ κ(A) ‖δb‖/‖b‖. A condition number of 10^k costs k decimal digits of accuracy. Ill-conditioning is a property of the problem, not the algorithm; a stable algorithm (like LU with pivoting or QR) cannot do better than the condition number allows, but an unstable algorithm can do far worse. Always solve Ax = b via factorization, never by computing A⁻¹; for least squares, use QR or SVD to avoid squaring the condition number.