Example number
30
Example slug
example_30_condition_number_predicts_sensitivity
Background

Numerical linear algebra developed stability theory to make computation reliable under finite precision; ML inherits these issues at scale. The concept of condition number emerged in the 1960s when numerical analysts quantified the relationship between input perturbations and output errors for ill-posed problems. John G. F. Francis and Vera N. Kublanovskaya developed the QR algorithm, revealing how matrix properties (singular value distribution, eigenvalue spacing) control computational difficulty. Gene Golub and colleagues formalized backward error analysis, showing that even stable algorithms cannot correct for inherent problem ill-conditioning. The fundamental insight is that conditioning is a property of the problem, not the solver: a perfectly implemented algorithm cannot overcome poor conditioning; it merely achieves backward stability (computed solution satisfies a nearby problem). In machine learning, ill-conditioning appears in multiple contexts: (1) Feature scaling: unscaled features with vastly different ranges create ill-conditioned Hessians in neural networks, slowing training. (2) Near-collinearity: correlated features in linear regression inflate the normal equation condition number. (3) Penalty methods: interior-point optimization solvers for constrained problems generate increasingly ill-conditioned systems as the barrier parameter shrinks. Modern practice addresses conditioning through preprocessing (normalization, whitening), regularization (ridge, early stopping, dropout), and algorithm design (gradient descent with preconditioning, natural gradient exploiting Riemannian geometry).

Purpose

Show how ‘hard to train’ often means ‘ill-conditioned’ and how factorization choices change outcomes. The condition number $\kappa(A)$ is a single number that quantifies solution sensitivity: it predicts how much a small perturbation in the input $b$ will be amplified into the output $x$ when solving $Ax = b$. In machine learning, ill-conditioning (large $\kappa$) manifests as optimization difficulties—erratic gradient updates, slow convergence, training instability—even when the algorithm is implemented correctly. Recognizing this distinction between algorithm bugs and problem conditioning is crucial for effective debugging and intervention. This example demonstrates the perturbation bound $\|\delta x\| / \|x\| \lesssim \kappa(A) \|\delta b\| / \|b\|$, showing empirically that condition number accurately predicts sensitivity. Well-conditioned problems ($\kappa \approx 1$) are numerically benign; perturbations pass through unchanged. Ill-conditioned problems ($\kappa \gg 1$) amplify noise catastrophically. Understanding conditioning equips you to diagnose training failures, design regularization schemes (ridge regression, early stopping), and choose appropriate solvers (iterative methods with preconditioning for ill-conditioned systems). Condition number is not a solver-specific metric—it’s a property of the problem itself. No algorithm can overcome poor conditioning without additional information (regularization, constraints, prior knowledge).

Problem

Solve a well-conditioned and an ill-conditioned system and compare sensitivity to perturbations.

Solution (Math)

Large $\kappa(A)$ implies small input perturbations can cause large output changes. Compare systems with different singular value spreads.

We use:

  • Data matrix $X\in\mathbb{R}^{n\times d}$ (rows are examples).
  • Vectors are column vectors by default.
  • $\|x\|_2$ is Euclidean norm; $\langle x,y\rangle=x^Ty$.
Solution (Python)

Runnable diagnostic comparing condition number to empirical sensitivity. Prints cond(A) and relative solution change for well/ill-conditioned cases.


import numpy as np

def trial(A):
    b = np.array([1.,1.])
    db = np.array([0.,1e-6])
    x = np.linalg.solve(A,b)
    x2 = np.linalg.solve(A,b+db)
    return np.linalg.cond(A), np.linalg.norm(x2-x)/np.linalg.norm(x)

A_good = np.array([[1.,0.],[0.,0.9]])
A_bad  = np.array([[1.,0.],[0.,1e-6]])

for name,A in [("good",A_good),("bad",A_bad)]:
    condA, rel_dx = trial(A)
    print(name, "cond:", condA, "rel_dx:", rel_dx)
Code Introduction

This code demonstrates how matrix conditioning directly controls solution sensitivity: small perturbations in the right-hand side cause large solution changes when the condition number is large. The workflow consists of three steps: (1) define a trial function that solves a system twice—once with the original right-hand side $b$ and once with a tiny perturbation $b + \delta b$—then measure the fractional solution change relative to the condition number, (2) construct two diagonal test matrices with vastly different condition numbers (well-conditioned with $\kappa \approx 1.11$ and ill-conditioned with $\kappa = 10^6$), and (3) run the trial on both matrices and print the results side-by-side for empirical comparison. The key insight is the perturbation bound: relative solution error $\lesssim \kappa(A) \times$ relative input error. For well-conditioned matrices, small perturbations produce small solution changes—stable, predictable behavior. For ill-conditioned matrices, the same perturbation magnitude amplifies into a large solution change—the system is inherently sensitive to input noise. The trial function encapsulates this experiment in a reusable way: calling trial(A) returns both the condition number and the empirical relative solution change, enabling direct comparison of theory (perturbation bound) and practice (measured sensitivity). This pattern—diagnose conditioning, compare to theoretical predictions—is essential for debugging numerical algorithms and understanding why some systems are harder to solve than others. Condition number is a property of the problem, not the solver: even the most stable algorithm cannot correct for poor conditioning, though it can compute the solution accurately relative to the problem’s intrinsic sensitivity. Understanding conditioning diagnostics is crucial for identifying when a training difficulty reflects poor problem conditioning (requiring regularization or preprocessing) versus algorithm failure (requiring algorithmic correction) or implementation bugs (requiring code debugging). This code is the diagnostic pattern used throughout numerical computing: compute condition number, compare actual sensitivity to theoretical predictions, and adjust algorithms (regularization, preconditioning, iterative solving) accordingly to maintain numerical stability and computational efficiency.

Numerical Implementation Details

Numerical and Shape Notes

  • Shapes first: Declare shapes (e.g., $X \in \mathbb{R}^{n imes d}$, $w \in \mathbb{R}^{d}$, $b \in \mathbb{R}^{n}$). Vectors are column by convention; keep row/column usage consistent.
  • Axis discipline: Be explicit with axis in reductions and normalizations. For attention-like ops, softmax over keys (row-wise) so rows sum to ≈1.
  • Broadcasting: Check that broadcasts are intended (e.g., (n,1) with (n,d)). Prefer reshape/expand-dims to make semantics clear.
  • Stability eps: Add $arepsilon$ for divisions/logs and $arepsilon I$ (jitter) for SPD solves; use log-sum-exp for softmax.
  • Masking preserves shape: Masks should broadcast to the score/activation tensor; verify masked outputs keep the same shape and zero out excluded entries.
  • Dtype choices: Use float64 for clarity in scripts; with mixed precision, keep reductions/factorizations in float32/float64 to avoid under/overflow.
  • Sanity checks: Print shapes and residuals (e.g., ||Ax-b||, reconstruction error, row-sum ≈ 1). Assert finiteness and expected monotonicity where applicable.

Numerical and Implementation Notes

  • Dtype & precision: Prefer float64 for clarity; if using mixed precision, keep reductions (norms, softmax sums, factorizations) in float32/float64. Avoid explicit inverses; use solve, lstsq, Cholesky/QR/SVD.
  • Shapes & broadcasting: Annotate shapes (e.g., $X \in \mathbb{R}^{n imes d}$); vectors are column by default. Verify axes for reductions (axis) and ensure broadcasts are intended.
  • Stability: Use log-sum-exp for softmax; add small diagonal $arepsilon I$ (jitter) for SPD solves; prefer QR/SVD for ill-conditioned least squares.
  • Conditioning: Inspect np.linalg.cond(A) when solutions look unstable; regularize (ridge) or rescale features to improve conditioning.
  • Reproducibility: Set NumPy seed for random data; print shapes and residuals (e.g., ||Ax-b||, reconstruction errors) and assert finiteness.
  • Complexity & memory: Matmul ~ $O(n^3)$ for factorizations, $O(n^2)$ for triangular solves/products. Prefer vectorization over Python loops; avoid materializing large intermediates.
  • Masking & indexing: Use boolean masks that broadcast to target shapes; for attention-like ops, add $-\infty$ before softmax or zero-out after, then verify rows sum to ~1.
  • Sanity checks: Compare against references (e.g., lstsq vs. solve), check orthogonality (U.T @ U ≈ I), PSD (x.T @ A @ x > 0), and residual norms within tolerance (~1e-12 for float64).
  1. Trial function design: The trial(A) function encapsulates a controlled perturbation experiment: given a matrix $A$, it solves two systems with slightly different right-hand sides, then measures solution sensitivity. Wrapping this in a function enables systematic comparison across multiple matrices without code duplication.
  2. Baseline solution: b = np.array([1., 1.]) and x = np.linalg.solve(A, b) compute the reference solution. The choice of $b = [1, 1]^\top$ is arbitrary but standardizes the comparison—using $\|b\|_2 = \sqrt{2}$ as the baseline norm.
  3. Perturbation magnitude: db = np.array([0., 1e-6]) adds a tiny perturbation ($\sim 10^{-6}$) to only the second component. This is not a measurement error but a controlled input change to test system response. The magnitude $10^{-6}$ is small enough that perturbation theory applies (linear approximation $\delta x \approx A^{-1} \delta b$ is valid).
  4. Relative perturbation size: The relative perturbation in $b$ is $\|\delta b\|_2 / \|b\|_2 = 10^{-6} / \sqrt{2} \approx 0.7 \times 10^{-6}$. This quantifies input noise as a fraction of signal.
  5. Perturbed solution: x2 = np.linalg.solve(A, b + db) computes the solution with the perturbed right-hand side. Computing both solutions (baseline and perturbed) measures the actual solution change directly, avoiding accumulated errors from finite-difference approximations.
  6. Relative solution change: np.linalg.norm(x2 - x) / np.linalg.norm(x) computes the fractional solution change as a percentage. For well-conditioned problems, this is approximately equal to the relative input perturbation; for ill-conditioned problems, it’s much larger.
  7. Condition number extraction: np.linalg.cond(A) computes the spectral condition number $\kappa(A) = \sigma_{\max}(A) / \sigma_{\min}(A)$ using SVD internally. For diagonal matrices, this is simply the ratio of largest to smallest absolute diagonal entry.
  8. Empirical validation: Printing both condA and rel_dx for well-conditioned and ill-conditioned matrices demonstrates the predicted scaling $\text{rel\_dx} \approx \kappa(A) \times \text{rel\_db}$—validating the theory through numerical evidence.
What This Example Demonstrates

Pedagogical Significance

  • Learning goals: Build intuition for when and why this tool is used in ML, not just how to compute it.
  • ML-first framing: Tie the concept to a concrete task pattern (fit / project / decompose / solve / measure) to anchor understanding.
  • Shape discipline: Habitually annotating dimensions prevents silent bugs and reinforces linear map thinking.
  • Numerical habits: Prefer stable factorizations over inverses; check residuals and condition numbers to separate bugs from ill-conditioning.
  • Transfer: Reuse the same pattern across models (e.g., projection in PCA, orthogonalization in regressions, attention as weighted sums).
  • Assessment ideas: Quick checks: predict sensitivity from $\kappa(A)$, verify projection properties, or compare solver outputs within tolerance.

ML Examples and Patterns

  • Fit: Linear/logistic regression via least squares or softmax; regularization (ridge) improves conditioning and generalization.
  • Project: PCA/SVD for dimensionality reduction; orthogonal projections to subspaces for denoising and feature extraction.
  • Decompose: Eigen/SVD factorizations to expose structure (low rank, PSD) used in recommender systems, LSA, and spectral clustering.
  • Solve: Stable solves without inversion (Cholesky/QR/SVD; CG for SPD) for optimization steps and kernel methods.
  • Measure: Norms, angles, and condition number $\kappa(A)$ to diagnose sensitivity, stability, and training difficulty.
  • Condition number as amplification factor: The condition number $\kappa(A) = \sigma_{\max}(A) / \sigma_{\min}(A)$ (ratio of largest to smallest singular value) quantifies the worst-case amplification of input perturbations in the solution of $Ax = b$.
  • Perturbation bound: For small input perturbation $\delta b$, the relative solution change satisfies $\|\delta x\| / \|x\| \lesssim \kappa(A) \|\delta b\| / \|b\|$—the condition number acts as a multiplier on relative input error.
  • Diagonal structure isolation: For diagonal matrices $A = \text{diag}(\lambda_1, \ldots, \lambda_n)$, the condition number is $\kappa(A) = |\lambda_{\max}| / |\lambda_{\min}|$, making it easy to visualize: large disparities in eigenvalues cause large condition numbers.
  • Empirical validation: Running the trial function on well-conditioned (eigenvalues 1 and 0.9) vs. ill-conditioned (eigenvalues 1 and $10^{-6}$) matrices shows that relative solution change is approximately $\kappa \times$ relative input perturbation magnitude.
  • Problem vs. solver distinction: Condition number is intrinsic to the problem, not the algorithm. Even the most numerically stable solver (e.g., SVD-based least squares) cannot overcome ill-conditioning—it merely computes the solution to near machine precision, which may still be far from the truth due to input perturbations.
  • ML implications: In training, ill-conditioning manifests as optimization difficulty (erratic loss curves, slow convergence) or sensitivity to data noise. Regularization (ridge, dropout) improves conditioning by modifying the problem, not by better numerics.
Notes

Part 1: The Trial Function and Perturbation Protocol
The trial function measures how a tiny perturbation in the right-hand side $b$ propagates to the solution $x$ of the system $Ax = b$. Starting with nominal $b = [1., 1.]^\top$, a perturbation $\delta b = [0., 10^{-6}]^\top$ is added (perturbing only the second component by one millionth). The function solves both systems: x = solve(A, b) for the baseline and x2 = solve(A, b + db) for the perturbed case. The return value is the relative change in solution: $\frac{\|x_2 - x\|_2}{\|x\|_2}$, alongside the condition number $\kappa(A) = \|A\|_2 \|A^{-1}\|_2$, which is np.linalg.cond(A) (the ratio of largest to smallest singular value for general matrices, or largest to smallest eigenvalue magnitude for symmetric matrices).

Part 2: Well-Conditioned vs. Ill-Conditioned Matrices
Two test matrices illustrate the extremes: Well-conditioned case: $A_{\text{good}} = \begin{bmatrix} 1 & 0 \\ 0 & 0.9 \end{bmatrix}$ is a diagonal matrix with eigenvalues 1 and 0.9. The condition number is $\kappa(A_{\text{good}}) = 1.0 / 0.9 \approx 1.11$—nearly 1, which is the best possible (perfectly conditioned). When run with this matrix, the relative solution change is approximately $10^{-6}$, matching the relative perturbation magnitude in $b$. This is stable: the solution responds predictably to input changes. Ill-conditioned case: $A_{\text{bad}} = \begin{bmatrix} 1 & 0 \\ 0 & 10^{-6} \end{bmatrix}$ has eigenvalues 1 and $10^{-6}$. The condition number is $\kappa(A_{\text{bad}}) = 1.0 / 10^{-6} = 10^6$—huge, indicating severe ill-conditioning. When run with this matrix, the relative solution change is approximately $10^{-6} \times 10^6 = 1.0$ or larger—the solution can change by 100% from a one-millionth perturbation. This is unstable: small input noise gets amplified into large solution noise.

Part 3: The Fundamental Perturbation Bound
The phenomenon is explained by the condition number perturbation bound. For $Ax = b$ with small perturbation $\delta b$, the relative solution change satisfies $\frac{\|\delta x\|_2}{\|x\|_2} \lesssim \kappa(A) \frac{\|\delta b\|_2}{\|b\|_2}$. The condition number acts as an amplification factor: relative solution error is at most $\kappa(A)$ times the relative input error. When $\kappa(A) \approx 1$ (well-conditioned), perturbations pass through nearly unchanged. When $\kappa(A) \gg 1$ (ill-conditioned), perturbations are amplified dramatically. For $A_{\text{bad}}$ with $\kappa = 10^6$ and $\|\delta b\|_2 / \|b\|_2 \approx 10^{-6}$, the predicted worst-case relative solution change is $10^6 \times 10^{-6} = 1$, matching the observed behavior.

Why This Matters for ML and Numerics
Conditioning is a hidden killer in machine learning and scientific computing: (1) Training instability: Neural networks with poorly scaled features or loss landscapes (high Hessian condition number) can exhibit erratic gradient updates and slow convergence. Preprocessing—normalization, standardization, dimensionality reduction—reduces conditioning. (2) Regularization as conditioning control: Ridge regression $(X^\top X + \lambda I)^{-1}$ adds $\lambda$ to small singular values, improving conditioning at the cost of bias. As $\lambda$ increases, solutions become more stable but more biased (bias-variance trade-off). (3) Least squares sensitivity: For overdetermined systems $Ax = b$ (linear regression), the normal equations $X^\top X w = X^\top y$ square the condition number: $\kappa(X^\top X) = \kappa(X)^2$. This is why SVD-based solvers like np.linalg.lstsq are preferred—they work with $X$ directly, avoiding the squaring penalty. (4) Solver choice: For well-conditioned SPD systems, Cholesky decomposition is fast and stable. For ill-conditioned systems, iterative solvers (conjugate gradient with preconditioning) or regularized approaches (Tikhonov, ridge) are necessary.

Numerical and Shape Notes
Computing $\kappa(A)$ via np.linalg.cond(A) uses SVD internally (cost $O(n^3)$), so don’t call it repeatedly in loops. A rule of thumb: if $\kappa(A) > 10^{10}$ on a system with double precision (16 significant digits), expect to lose about 6 digits of accuracy due to rounding. Beyond $\kappa > 10^{12}$, the problem is effectively singular for practical purposes. Shape discipline: $A$ can be rectangular (SVD-based condition number works for all shapes); for square matrices, np.linalg.cond with default p=2 gives the spectral condition number ($\sigma_{\max}/\sigma_{\min}$). Printing both the condition number and the relative solution change together creates an empirical check: large condition number + large relative error = expected and not a solver bug. Gotchas: (1) Ill-conditioning is a property of the problem, not the solver—no algorithm can overcome it without regularization. (2) Attempted inversion via x = np.linalg.inv(A) @ b is even less stable than solving; always use solve. (3) For iterative solvers, preconditioning reduces effective condition number by exploiting problem structure (e.g., incomplete LU preconditioning for sparse systems).

History and Applications

Condition number and perturbation theory emerged as central concepts in numerical analysis during the 20th century, driven by the computational revolution. In the 1950s-1960s, as electronic computers became mainstream, numerical analysts recognized that two different algorithms could produce wildly different results for the same problem, not due to correctness but due to finite-precision arithmetic. Gene Golub and Wilkinson developed the theory of backward error analysis, showing that a computed solution $\tilde{x}$ is the exact solution to a nearby problem: $A \tilde{x} = b + \delta b$ where $\|\delta b\|$ is small. The condition number quantifies the worst-case scenario: for a given algorithm’s backward error, condition number determines the forward error in the solution. John G. F. Francis and Vera N. Kublanovskaya’s QR algorithm (independently developed in the late 1950s) was revolutionary because it had superb backward stability properties—computed eigenvalues are exact for a matrix with tiny perturbations. This motivated a paradigm shift: prefer structured algorithms with proven stability over mathematically equivalent formulas that are numerically unstable (e.g., normal equations vs. SVD for least squares). In machine learning, condition number became essential for understanding training dynamics. The Hessian matrix’s condition number controls gradient descent convergence: large $\kappa(H)$ produces slow convergence (many small eigenvalues drag the process along, while one large eigenvalue rushes ahead, creating oscillation and plateaus). Momentum, adaptive learning rates (Adam, RMSprop), and second-order methods (BFGS, Newton) all attempt to improve conditioning implicitly—momentum accelerates along slow directions; adaptive rates shrink large gradients, flattening the loss landscape; second-order methods approximate the inverse Hessian, effectively preconditioning. In modern deep learning, conditioning appears throughout: batch normalization reduces internal covariate shift (improves Jacobian conditioning in intermediate layers), layer normalization stabilizes training in Transformers, and natural gradient descent (using the Fisher information matrix, a curvature proxy) exploits problem structure to achieve better conditioning. Ill-conditioning is not a solver bug—it reflects fundamental properties of the problem. Addressing it requires domain insight: feature preprocessing (normalization, whitening, PCA), regularization (L2 penalty, dropout, label smoothing), or problem reformulation (change of variables, constraints). The historical arc—from Golub’s backward error analysis to modern ML’s conditioning-aware training techniques—shows that understanding numerical sensitivity is not a niche mathematical concern but central to reliable, efficient computation at scale.

Connection to Broader Examples

Conditioning synthesizes concepts from multiple chapters:

  • Chapter 10 (SVD): The condition number is defined using singular values: $\kappa(A) = \sigma_1 / \sigma_n$ where $\sigma_i$ are the singular values in decreasing order. Large gaps between singular values (steep decay) indicate ill-conditioning.
  • Chapter 12 (Least Squares): In linear regression, the normal equations $X^\top X w = X^\top y$ have condition number $\kappa(X^\top X) = \kappa(X)^2$, squaring the sensitivity. This is why SVD-based solvers like np.linalg.lstsq are preferred—they work with $X$ directly, avoiding the squaring penalty.
  • Chapter 13 (Solving Systems): Different factorizations have different stability properties. Cholesky (for SPD) is unconditionally stable; LU (general) requires pivoting. Condition number predicts how much these stability properties matter: well-conditioned $\Rightarrow$ even unstable methods work fine; ill-conditioned $\Rightarrow$ must use stable methods.
  • Chapter 14 (Conditioning & Stability): This example demonstrates the core of Chapter 14—the relationship between condition number and perturbation sensitivity.
  • Optimization in ML: The Hessian matrix $H = \nabla^2 f(w)$ of the loss function controls gradient descent convergence. Large condition number $\kappa(H)$ means slow convergence (eigenvalues spread over many orders of magnitude, creating plateaus in the loss landscape).
  • Regularization as conditioning control: Ridge regression modifies the system $(X^\top X + \lambda I) w = X^\top y$, improving conditioning by pushing small singular values up toward $\lambda$.

Comments