Numerical methods for machine learning: conditioning, factorizations, and stability

Why this subject matters now

For most of the last decade the dominant narrative in machine learning was that numerical care did not matter. Throw float32 at a large enough model, let Adam absorb the conditioning problems, and the loss goes down. That narrative broke in two places at once. Training moved to float16 and bfloat16, where the exponent range is small enough that a plain sum of gradients can overflow or vanish, and inference moved to int8 and 4-bit formats, where every rounding decision is visible in the output. A practitioner today is expected to know why a matrix multiply in bfloat16 loses three decimal digits relative to float64, why loss scaling exists, why the normal equations are never formed when a QR factorization is available, and why a covariance matrix that is mathematically positive definite can fail a Cholesky factorization. These are not exotic questions. They are the difference between a training run that diverges at step ten thousand and one that does not.

The subject itself is old and settled. The rounding model was standardized in IEEE 754 in 1985 and revised in 2008 and 2019. The analysis of Gaussian elimination is due to Wilkinson in the 1960s, the conjugate gradient method dates to Hestenes and Stiefel in 1952, and the singular value decomposition and its optimality were understood by Eckart and Young in 1936. What is new is the audience. The results that numerical analysts proved for scientific computing now govern whether a hundred-billion-parameter model trains, because a transformer forward pass is a long chain of matrix multiplications, softmaxes, and layer norms, each of which is a place where finite precision either behaves or does not. This page develops the classical theory and then connects each piece to where it bites in modern training. The optimization side of the story, gradient descent and its convergence, momentum, Newton and L-BFGS as algorithms, is developed in the combinatorial optimization page. The statistical reading of least squares and shrinkage lives in statistical learning. Here the single concern is the arithmetic, what the machine actually computes and how far it is from what the mathematics asked for.

Core theory

Floating point and the rounding model

A binary floating-point number in IEEE 754 double precision is stored as a sign bit \( s \), an 11-bit exponent \( e \), and a 52-bit fraction \( f \), representing

The leading 1 is implicit, so the significand carries 53 bits of precision. The set of representable numbers is not uniform. Between consecutive powers of two \( [2^k, 2^{k+1}) \) the spacing is constant at \( 2^{k-52} \), and it doubles at each power of two. The relative spacing is therefore roughly constant everywhere, and that constant is the quantity that governs everything downstream. Define the unit roundoff \( u = 2^{-53} \approx 1.11 \times 10^{-16} \), and machine epsilon \( \varepsilon_{\text{mach}} = 2^{-52} \approx 2.22 \times 10^{-16} \) as the gap between \( 1 \) and the next representable double. Machine epsilon has an operational definition that does not mention the format at all. It is the smallest \( \varepsilon \) such that \( \text{fl}(1 + \varepsilon) \neq 1 \). Searching for it by repeated halving,

terminates at exactly \( 2.220446049250313 \times 10^{-16} \), which agrees to the last digit with the value NumPy reports from the format description. The same loop in single precision returns \( 1.1920929 \times 10^{-7} = 2^{-23} \), the 24-bit significand of float32. The rounding rule that produces \( u \) is round-to-nearest-even. A real number is rounded to the nearest representable value, with ties going to the one with an even last bit, which removes the statistical bias that round-half-up would introduce over a long computation.

The single fact that makes error analysis possible is the standard model of floating-point arithmetic. For any two representable numbers \( x, y \) and any of the four basic operations \( \circ \in \{+, -, \times, \div\} \), the computed result satisfies

In words, the machine computes the exact result and then rounds it once, and rounding perturbs the answer by at most a relative \( u \). IEEE 754 guarantees this for \( +, -, \times, \div \) and for the square root, because it mandates that these operations be correctly rounded. The returned value is the exactly-rounded true result. The model is the lever for every bound that follows. It says nothing about the absolute error, only the relative error, which is why the trouble is never a single rounding but the way relative errors in intermediate quantities get amplified into the final answer. That amplification has a name, and it is the whole subject.

Catastrophic cancellation

Rounding a single number loses at most one unit in the last place. Cancellation is different. It does not create new error, it exposes error that was already present. When two nearly equal numbers are subtracted, the leading digits agree and cancel, and the result is built from the low-order digits where the previous rounding errors live. The relative error of the difference can be very large even though each operand was accurate to full precision. Concretely, if \( \hat{a} = a(1 + \delta_a) \) and \( \hat{b} = b(1 + \delta_b) \) with \( |\delta_a|, |\delta_b| \le u \), then

and the amplification factor \( |a| / |a - b| \) blows up precisely when \( a \approx b \). The canonical example is \( (1 - \cos x)/x^2 \), which tends to \( \tfrac12 \) as \( x \to 0 \). For small \( x \), \( \cos x \) is close to \( 1 \), so \( 1 - \cos x \) subtracts two nearly equal numbers. Running it in double precision shows the digits leave one by one.

At \( x = 10^{-8} \) the naive expression returns exactly zero. The value \( \cos(10^{-8}) \) rounds to \( 1 \) in double precision because the true value \( 1 - 5 \times 10^{-17} \) is closer to \( 1 \) than the next double below it, so \( 1 - \cos x \) is \( 0 \) and the answer has no correct digits. The stable form uses the identity \( 1 - \cos x = 2\sin^2(x/2) \), which never subtracts nearly-equal quantities and stays accurate to the last digit. The cure is always the same. Rewrite the expression so the subtraction of near-equal terms disappears. The quadratic formula is the other textbook case. The roots of \( x^2 - (10^8 + 10^{-8})x + 1 = 0 \) are \( 10^8 \) and \( 10^{-8} \). Computing the small root as \( (-b - \sqrt{b^2 - 4ac})/2a \) subtracts two numbers that agree to sixteen digits and returns \( 1.4901 \times 10^{-8} \), a relative error of \( 0.49 \). Computing it instead as \( c/(a \, r_1) \), using the fact that the product of the roots is \( c/a \), returns \( 10^{-8} \) exactly. Cancellation also destroys the textbook one-pass variance formula \( \operatorname{Var} = \tfrac{1}{n-1}(\sum x_i^2 - (\sum x_i)^2/n) \). For the five numbers \( 10^8 + 1, \dots, 10^8 + 5 \) the true variance is \( 2.5 \), the numerically stable two-pass formula returns \( 2.5 \), and the one-pass formula returns \( 2.0 \), because \( \sum x_i^2 \) and \( (\sum x_i)^2/n \) are each near \( 5 \times 10^{16} \) and their difference lives entirely in the digits that rounding already corrupted.

Summation and Kahan compensation

Adding \( n \) numbers naively accumulates rounding error that grows with \( n \). Let \( s_k \) be the exact partial sums and \( \hat{s}_k \) the computed ones. Each addition contributes a relative perturbation, \( \hat{s}_k = (\hat{s}_{k-1} + x_k)(1 + \delta_k) \), and unrolling the recursion gives a computed total

Using \( \prod (1 + \delta_j) = 1 + \theta \) with \( |\theta| \le (n-1)u/(1 - (n-1)u) \approx (n-1)u \) for \( (n-1)u \ll 1 \), the forward error is bounded by

The bound is proportional to \( n \) and to the sum of magnitudes, not to the magnitude of the answer, so summing many terms of the same sign in single precision is where it hurts. Summing \( 10^6 \) copies of \( 0.1 \) sequentially in float32 returns \( 100958.34 \) against a true value of \( 100000 \), an absolute error of \( 958 \), because \( 0.1 \) is not representable exactly and each of a million additions rounds again against a partial sum that has grown far larger than the term being added, so the term's low bits fall off the end.

Kahan's compensated summation recovers those lost bits by carrying a running correction. After forming \( t = s + y \), the quantity \( (t - s) - y \) is, in exact arithmetic, zero, but in floating point it equals the part of \( y \) that did not fit into \( t \), with the sign flipped. Feeding it back removes the error at the next step.

The analysis replaces the \( (n-1)u \) factor with a bound of the form \( 2u + O(nu^2) \), so the error no longer grows with \( n \) to first order. On the same million-element float32 sum, Kahan returns exactly \( 100000.0 \), an absolute error of zero to the precision printed. This is why every serious reduction, including the ones inside cuBLAS and the accumulation of a bfloat16 matmul into an fp32 register, does the accumulation in a wider type or with compensation rather than naively in the storage type. NumPy sidesteps the problem differently, with pairwise summation, which sums in a balanced binary tree so the error grows like \( u \log n \) rather than \( u n \). Its plain .sum() over ten million float32 copies of \( 0.1 \) returns \( 999989.44 \) against \( 1000000 \), a relative error of \( 10^{-5} \), far better than the sequential loop but not as good as Kahan.

Conditioning, or how sensitive the problem is

Stability is a property of an algorithm. Conditioning is a property of the problem, independent of any algorithm. The condition number measures how much a relative perturbation of the input can be amplified into a relative perturbation of the output. For a differentiable scalar function \( f \), perturb \( x \) to \( x(1 + \delta) \) and expand.

For \( f(x) = \log x \) this is \( 1/|\log x| \), which is large near \( x = 1 \), matching the cancellation seen above. For \( f(x) = \sqrt{x} \) it is \( \tfrac12 \), well conditioned everywhere. For subtraction \( f(a,b) = a - b \) with \( a \approx b \) the relative condition number is \( (|a| + |b|)/|a - b| \), which is exactly the cancellation amplification factor. For a matrix-vector product \( y = Ax \), the relative condition number with respect to perturbations in \( x \) is

where the last inequality uses \( \|x\| = \|A^{-1} A x\| \le \|A^{-1}\| \|Ax\| \). The quantity \( \kappa(A) = \|A\| \|A^{-1}\| \) is the condition number of the matrix. In the 2-norm it equals the ratio of the largest to the smallest singular value, \( \sigma_{\max}/\sigma_{\min} \), so a matrix is ill-conditioned exactly when its smallest singular value is small relative to its largest, that is, when it is close to singular. For the linear system \( Ax = b \), perturbing \( b \) to \( b + \Delta b \) changes the solution by \( \Delta x = A^{-1} \Delta b \), and

with an analogous bound for perturbations in \( A \). The condition number is the amplifier. A system with \( \kappa(A) = 10^{10} \) can turn a relative input error of \( 10^{-16} \) into a relative output error of \( 10^{-6} \), losing ten of the sixteen available digits before any algorithm is even chosen.

Forward error, backward error, and stability

An algorithm computes \( \hat{y} \) where the true answer is \( y = f(x) \). The forward error is \( \|\hat{y} - y\| / \|y\| \), the thing one actually wants small. It is usually impossible to bound directly. The backward error asks a different question. For what perturbed input would the computed answer be exactly correct? It is the smallest \( \Delta x \) such that \( \hat{y} = f(x + \Delta x) \), measured as \( \|\Delta x\|/\|x\| \). An algorithm is backward stable if it always produces an answer with backward error on the order of the unit roundoff. The computed answer is then the exact answer to a problem negligibly different from the one posed. This is the strongest guarantee one can ask of finite-precision arithmetic, because rounding the input already perturbs it by \( u \). An algorithm cannot do better than solve a nearby problem exactly. The two errors are tied together by the condition number in one line.

This is the master inequality of the field. A backward-stable algorithm has backward error \( \approx u \), so its forward error is at most \( \kappa u \). If that product is large, the fault lies with the problem, not the algorithm. No method solving an ill-conditioned problem in \( u \)-precision arithmetic can promise a small forward error. The practical program is therefore split cleanly. Choose a backward-stable algorithm, so the backward error is guaranteed to be \( \approx u \). The forward error is then whatever the conditioning allows, and nothing can improve it short of higher precision or a better-conditioned reformulation of the problem.

The Hilbert matrix \( H_{ij} = 1/(i + j - 1) \) is the standard demonstration because its condition number grows rapidly with size while every entry is a simple rational. Solving \( H x = b \) with \( x = (1, \dots, 1) \), using a backward-stable LU factorization, gives the results below.

Read the table across. The backward error stays pinned at machine precision for every size, which is the signature of a backward-stable solver. The computed \( \hat{x} \) exactly solves a system \( (H + \Delta H)\hat{x} = b \) with \( \|\Delta H\|/\|H\| \approx 10^{-16} \), and the residual \( \|H\hat{x} - b\| \) confirms it. The forward error, by contrast, tracks \( \kappa u \) faithfully. At \( n = 12 \), \( \kappa \approx 10^{16} \), and the solution has no correct digits at all, precisely because \( \kappa u \approx 1 \). The solver did nothing wrong. The problem is the difficulty. These residuals were cross-checked by explicitly forming \( H\hat{x} - b \) and confirming it is at machine precision, the sanity check that this host's history of a broken CPU BLAS makes mandatory.

The squared condition number, or why the normal equations are avoided

The least-squares problem \( \min_x \|Ax - b\|_2 \) for a tall \( A \in \R^{m \times n} \) has the exact solution characterized by the normal equations \( A^\top A \, x = A^\top b \). Forming \( A^\top A \) and solving it looks efficient and is a trap, because it squares the condition number. The singular values of \( A^\top A \) are the squares of the singular values of \( A \), so

The forward error of the computed least-squares solution therefore scales like \( \kappa(A)^2 u \) when the problem is squared, against \( \kappa(A) u \) for a method that works with \( A \) directly. Take a matrix with \( \kappa(A) = 10^8 \) and a consistent right-hand side so the exact answer is known. Solving the normal equations by Cholesky gives a relative error of \( 0.12 \), close to the \( \kappa^2 u = 2.2 \) that the theory predicts as an upper bound and consistent with having lost roughly sixteen of sixteen digits. Solving the same problem by a Householder QR factorization, which never forms \( A^\top A \), gives a relative error of \( 4.3 \times 10^{-10} \), close to \( \kappa u = 2.2 \times 10^{-8} \). The difference is eight orders of magnitude, and it comes entirely from refusing to square the condition number. This is why LAPACK's least-squares driver and every well-written regression solver use QR or the SVD, not the normal equations. The statistical reading of the same fact, that \( X^\top X \) should not be inverted, is in the statistical learning page.

LU factorization and pivoting

Gaussian elimination factors a square \( A \) into a unit lower-triangular \( L \) and an upper-triangular \( U \). Without care it is unstable, because a small pivot forces a large multiplier \( \ell_{ik} = a_{ik}/a_{kk} \) that magnifies rounding errors from the rows it combines. Partial pivoting fixes this by swapping the row with the largest-magnitude entry in the current column to the pivot position, which forces every multiplier to satisfy \( |\ell_{ik}| \le 1 \). The factorization becomes \( PA = LU \) with a permutation \( P \), and the backward error is governed by the growth factor

the ratio of the largest entry appearing anywhere during elimination to the largest entry of the original matrix. Wilkinson's backward-error bound for LU with partial pivoting is \( \|\Delta A\|_\infty \le c\, n\, \rho\, u \, \|A\|_\infty \) for a modest constant \( c \), so the method is backward stable precisely when \( \rho \) is not large. With partial pivoting \( \rho \le 2^{n-1} \), and this bound is attained by a specific adversarial matrix. The one with \( -1 \) below the diagonal, \( 1 \) on it, and \( 1 \) in the last column doubles the last column at every step. For \( n = 20 \) that construction produces a growth factor of exactly \( 524288 = 2^{19} \), confirming the bound is tight. The reassuring fact, established over decades of practice and partly explained by average case analysis, is that this exponential growth essentially never happens on matrices that arise in applications. The observed growth factor is almost always \( O(\sqrt{n}) \) or smaller, which is why Gaussian elimination with partial pivoting remains the default dense solver despite the pathological worst case.

Cholesky factorization

When \( A \) is symmetric positive definite, the factorization specializes to \( A = L L^\top \) with \( L \) lower triangular and positive diagonal. Equating entries of \( A \) and \( L L^\top \) gives the recurrence directly. The \( (j, j) \) entry yields the diagonal, and the \( (i, j) \) entry for \( i > j \) yields the off-diagonal.

Cholesky needs no pivoting for stability, which is a genuinely special property. Symmetric positive definiteness guarantees that every quantity under the square root is positive and that the entries of \( L \) cannot grow, so the growth factor is bounded by \( 1 \) and the method is unconditionally backward stable. It also costs half of LU, \( \tfrac13 n^3 \) flops against \( \tfrac23 n^3 \), because it exploits symmetry and computes only one triangle. One further property gets used constantly. Because the algorithm fails, by hitting a negative number under the square root, exactly when the matrix is not positive definite, attempting a Cholesky factorization is the standard cheap test for positive definiteness. Checking eigenvalues costs \( O(n^3) \) with a large constant, while a Cholesky attempt costs \( \tfrac13 n^3 \) and either succeeds or throws. The matrix \( \left[\begin{smallmatrix} 1 & 2 \\ 2 & 1 \end{smallmatrix}\right] \) has eigenvalues \( 3 \) and \( -1 \), and a Cholesky attempt on it raises immediately at the second diagonal, where \( 1 - 2^2 = -3 < 0 \) has no real square root. This is how well-written code checks that a covariance estimate is usable before trusting it. In machine learning the failure is common. A sample covariance from fewer samples than dimensions is positive semidefinite at best and numerically indefinite in practice, which is exactly why ridge and its cousins add \( \lambda I \) to push the smallest eigenvalue safely positive.

Householder QR

The QR factorization writes \( A = QR \) with \( Q \) orthogonal and \( R \) upper triangular. It is the numerically preferred route to least squares because orthogonal transformations preserve the 2-norm exactly, \( \|Qx\| = \|x\| \), so they cannot amplify error. A product of Householder reflections has a backward error of \( O(u) \) regardless of the conditioning of \( A \). The workhorse is the Householder reflector, the matrix that reflects a vector across a hyperplane so as to zero out all but its first component. Given a vector \( x \), the goal is an orthogonal \( H \) with \( Hx = \pm\|x\| e_1 \). Reflection across the hyperplane orthogonal to a unit vector \( v \) is

Here \( H \) is symmetric and orthogonal, since \( H^\top H = (I - 2vv^\top)^2 = I - 4vv^\top + 4v(v^\top v)v^\top = I \) when \( v^\top v = 1 \). To see it maps \( x \) onto \( \pm\|x\| e_1 \), choose \( u = x - \alpha e_1 \) with \( \alpha = \pm\|x\| \). Then \( \|u\|^2 = \|x\|^2 - 2\alpha x_1 + \alpha^2 = 2(\|x\|^2 - \alpha x_1) \) and \( u^\top x = \|x\|^2 - \alpha x_1 = \tfrac12\|u\|^2 \), so

The sign of \( \alpha \) is chosen opposite to \( x_1 \), that is \( \alpha = -\text{sign}(x_1)\|x\| \), so that \( u = x - \alpha e_1 \) never subtracts near-equal quantities. Choosing the same sign would reintroduce catastrophic cancellation in forming \( u_1 = x_1 - \alpha \). Applying \( n \) such reflectors, each zeroing the subdiagonal of one column, reduces \( A \) to upper-triangular \( R \), and the accumulated product of reflectors is \( Q \). Reflectors are never formed as explicit matrices. Each is applied as \( Hy = y - 2v(v^\top y) \), a matrix-vector product plus a rank-one update, so the whole factorization costs \( 2n^2(m - n/3) \) flops for an \( m \times n \) matrix and stores only the vectors \( v \). Gram-Schmidt produces the same factorization in exact arithmetic but is not backward stable in its classical form, because the computed columns of \( Q \) lose orthogonality as the algorithm proceeds. Modified Gram-Schmidt repairs this partially, but Householder is the method LAPACK uses.

The singular value decomposition

Every matrix \( A \in \R^{m \times n} \) factors as \( A = U \Sigma V^\top \), with \( U \in \R^{m \times m} \) and \( V \in \R^{n \times n} \) orthogonal and \( \Sigma \) diagonal with nonnegative entries \( \sigma_1 \ge \sigma_2 \ge \cdots \ge 0 \). Existence follows from the spectral theorem applied to \( A^\top A \), which is symmetric positive semidefinite and therefore has an orthonormal eigenbasis \( v_1, \dots, v_n \) with eigenvalues \( \lambda_i \ge 0 \). Set \( \sigma_i = \sqrt{\lambda_i} \) and, for \( \sigma_i > 0 \), \( u_i = A v_i / \sigma_i \). These \( u_i \) are orthonormal because \( u_i^\top u_j = v_i^\top A^\top A v_j / (\sigma_i \sigma_j) = \lambda_j v_i^\top v_j / (\sigma_i\sigma_j) = \delta_{ij} \), and \( A v_i = \sigma_i u_i \) rearranges to \( AV = U\Sigma \), hence \( A = U\Sigma V^\top \). The decomposition lays bare the four fundamental subspaces. The columns of \( U \) with \( \sigma_i > 0 \) span the column space, the remaining columns of \( U \) span the left null space, the columns of \( V \) with \( \sigma_i > 0 \) span the row space, and the remaining columns of \( V \) span the null space. The rank is the number of nonzero singular values, and the condition number in the 2-norm is \( \sigma_1/\sigma_r \). Numerically the SVD is the gold standard because the singular values reveal exactly how far the matrix is from each lower rank.

That last point is the Eckart-Young-Mirsky theorem, the reason the SVD underlies principal component analysis, low-rank adaptation, and every spectral compression scheme. Truncating the SVD to its top \( k \) terms, \( A_k = \sum_{i=1}^k \sigma_i u_i v_i^\top \), gives the best rank-\( k \) approximation in both the spectral and Frobenius norms, and the errors are the discarded singular values.

Proof for the spectral norm. First, \( \|A - A_k\|_2 = \|\sum_{i>k} \sigma_i u_i v_i^\top\|_2 = \sigma_{k+1} \), since the tail is itself an SVD with largest singular value \( \sigma_{k+1} \). Now let \( B \) be any matrix of rank at most \( k \). Its null space has dimension at least \( n - k \). The subspace spanned by \( v_1, \dots, v_{k+1} \) has dimension \( k + 1 \). Two subspaces of \( \R^n \) whose dimensions sum to more than \( n \) must intersect nontrivially, so there is a unit vector \( w \) in both, with \( Bw = 0 \) and \( w = \sum_{i=1}^{k+1} c_i v_i \) where \( \sum c_i^2 = 1 \). Then

using orthonormality of the \( u_i \) and \( \sigma_i \ge \sigma_{k+1} \) for \( i \le k+1 \). Hence no rank-\( k \) matrix beats \( A_k \), and \( A_k \) attains the bound, proving optimality. The Frobenius case follows the same argument with Mirsky's extension to all unitarily invariant norms. Running this on an \( 80 \times 60 \) matrix with a geometric spectrum \( \sigma_i = 2^{-i/4} \) confirms the theorem to the digit. At \( k = 5 \) the measured \( \|A - A_k\|_2 = 0.4204482 \) equals \( \sigma_6 = 0.4204482 \) exactly, and the Frobenius error \( 0.7768870 \) equals \( (\sum_{i>5} \sigma_i^2)^{1/2} \) exactly. Reconstruction of the full matrix from its SVD had error \( 8 \times 10^{-13} \), the machine-precision check that the factorization was computed correctly on this host.

The power method

Iterative eigenvalue methods matter because the largest matrices are never factored densely. The power method finds the dominant eigenvector by repeated multiplication and normalization. Start from a random \( v_0 \) and set \( v_{k+1} = A v_k / \|A v_k\| \). Expand \( v_0 \) in the eigenbasis, \( v_0 = \sum_i c_i x_i \), assuming \( |\lambda_1| > |\lambda_2| \ge \cdots \). Then

Every term after the first decays like \( |\lambda_i/\lambda_1|^k \), so \( v_k \) converges to \( x_1 \) at the linear rate \( r = |\lambda_2/\lambda_1| \), the ratio of the two largest eigenvalues in magnitude. The closer the top two eigenvalues, the slower the convergence. On a symmetric matrix with eigenvalues \( 10, 6, 3, 1 \), the ratio is \( r = 0.6 \), and the eigenvector error indeed falls by a factor near \( 0.6 \) per step. The eigenvalue estimate, taken as the Rayleigh quotient \( v_k^\top A v_k \), converges faster, at rate \( r^2 = 0.36 \), because for a symmetric matrix the Rayleigh quotient is stationary at an eigenvector and so its error is quadratic in the eigenvector error. The measured error ratios settle to exactly \( 0.36 \) within ten iterations. The power method is the conceptual seed of PageRank, of the largest-eigenvalue estimates used to set learning rates, and of the randomized SVD, which runs a few steps of a block power method on \( A \) times a random test matrix.

Conjugate gradients

For a symmetric positive definite \( A \), solving \( Ax = b \) is equivalent to minimizing the convex quadratic \( \phi(x) = \tfrac12 x^\top A x - b^\top x \), whose gradient is the residual \( \nabla\phi = Ax - b = -r \). Steepest descent takes the residual as the search direction and converges at the rate set by the condition number \( \kappa = \lambda_{\max}/\lambda_{\min} \), with error contracting by a factor \( (\kappa - 1)/(\kappa + 1) \) per step. For large \( \kappa \) that is \( 1 - 2/\kappa \), so the number of iterations grows linearly in \( \kappa \). Conjugate gradients does far better by choosing search directions that are \( A \)-orthogonal, \( p_i^\top A p_j = 0 \) for \( i \neq j \), rather than merely orthogonal. The key structural fact is that after \( k \) steps CG has minimized \( \phi \) exactly over the \( k \)-dimensional Krylov subspace

and it does so with a three-term recurrence that stores only a handful of vectors and needs one matrix-vector product per iteration. Because the error after \( k \) steps is the minimum over all degree-\( k \) polynomials \( q \) with \( q(0) = 1 \) of \( \|q(A) e_0\|_A \), and such polynomials can be made small on the spectrum of \( A \) using Chebyshev polynomials, the \( A \)-norm error obeys

The decisive feature is \( \sqrt{\kappa} \) where steepest descent has \( \kappa \). The iteration count to reach a fixed tolerance grows like \( \sqrt{\kappa} \) rather than \( \kappa \), a quadratic improvement. On an SPD system of size \( 200 \) with \( \kappa = 1000 \), CG reaches a relative residual of \( 10^{-8} \) in \( 173 \) iterations while steepest descent needs \( 5119 \), a ratio of about thirty, close to \( \sqrt{\kappa} = 31.6 \). CG is the reason large sparse and implicit systems, the normal equations of a large regression, the linear solves inside a Gauss-Newton step, the Hessian-vector products of a Newton-CG optimizer, are solved without ever forming a matrix. All CG needs is a routine that multiplies by \( A \).

The QR algorithm for eigenvalues

The dense eigenvalue workhorse is the QR algorithm, not to be confused with the QR factorization it uses as a subroutine. The idea is simple to state. Factor \( A_k = Q_k R_k \), then form \( A_{k+1} = R_k Q_k \) by multiplying the factors back in the opposite order. Each step is an orthogonal similarity transformation, \( A_{k+1} = Q_k^\top A_k Q_k \), so the eigenvalues are preserved, and under mild conditions the sequence \( A_k \) converges to upper-triangular (Schur) form with the eigenvalues on the diagonal, ordered by magnitude. The unshifted iteration is exactly a disguised power method on a whole basis at once, which explains both why it works and why it would be slow. Practical implementations first reduce \( A \) to upper Hessenberg form by Householder reflections, which makes each QR step cost \( O(n^2) \) instead of \( O(n^3) \), and then apply spectral shifts, replacing \( A_k \) by \( A_k - \mu_k I \) with a carefully chosen \( \mu_k \) that accelerates convergence to cubic near an eigenvalue. The shifted, Hessenberg-reduced, deflating QR algorithm is what LAPACK's dgeev and dsyev run, and it computes all eigenvalues of a dense matrix in \( O(n^3) \) with the backward stability that orthogonal transformations guarantee.

Forward-mode automatic differentiation

Training reduces to computing gradients of a scalar loss with respect to millions of parameters, and automatic differentiation is what makes that exact and cheap. It is neither symbolic differentiation, which produces unwieldy expression trees, nor numerical differentiation by finite differences, which trades a step-size bias against exactly the cancellation error analyzed above. Automatic differentiation applies the chain rule to the elementary operations of a program, and it comes in two modes. Forward mode is cleanest through dual numbers. Augment every real \( a \) to \( a + b\varepsilon \) with \( \varepsilon^2 = 0 \), where \( b \) carries the derivative. Arithmetic on duals follows from the nilpotency rule.

The multiplication rule is the product rule and the function rule is the chain rule, both falling out of a Taylor expansion truncated by \( \varepsilon^2 = 0 \). Evaluate a program on \( x + 1\cdot \varepsilon \) and the \( \varepsilon \)-component of the output is \( f'(x) \) exactly, to machine precision, with no step size and no cancellation. Differentiating \( g(x) = \sin(x^2) + e^x \) at \( x = 1.2 \) by seeding a dual returns \( 3.6331338237 \), matching the analytic \( 2x\cos(x^2) + e^x \) to every printed digit. The cost of forward mode is one extra number carried per input direction, so it computes a directional derivative, a Jacobian-vector product, in one pass. That makes it efficient when the number of inputs is small and the number of outputs large, the opposite of the regime that dominates deep learning.

Reverse-mode automatic differentiation

A neural network has one scalar output, the loss, and millions of inputs, the parameters. Forward mode would need one pass per parameter. Reverse mode gets all the partials in a single backward pass, which is why it, backpropagation, is the algorithm of the field. Reverse mode records the computation as a directed graph of elementary operations, then propagates derivatives from the output back to the inputs. For each intermediate \( v \), define the adjoint \( \bar{v} = \partial L / \partial v \), the sensitivity of the final output to that node. The chain rule says a node's adjoint is the sum over its consumers of the consumer's adjoint times the local partial derivative of the consumer with respect to the node.

Seed the output with \( \bar{L} = 1 \), visit the nodes in reverse topological order so every consumer is processed before its inputs, and accumulate. Each elementary operation contributes a local rule. An addition passes its adjoint unchanged to both inputs, a multiplication \( w = uv \) sends \( \bar{u} \mathrel{+}= \bar{w} v \) and \( \bar{v} \mathrel{+}= \bar{w} u \), a \( \sin \) sends \( \bar{u} \mathrel{+}= \bar{w}\cos u \), and so on. The total cost of the backward pass is a small constant times the cost of the forward pass, independent of the number of inputs, which is the property that makes gradient descent on billions of parameters affordable. The implementation section below builds this from scratch in fewer than fifty lines and checks it against torch.autograd and jax.grad on a nontrivial function. The three agree to \( 1.1 \times 10^{-16} \), one unit in the last place.

Reverse mode's cost is memory. Every intermediate value produced in the forward pass must be kept alive until the backward pass consumes it, because the local partials depend on the forward values. For a deep network the activation memory scales with depth times the layer width, and it is the reason training a large model needs far more memory than inference. Gradient checkpointing trades that memory for recomputation. It stores only a subset of activations, the checkpoints, and during the backward pass recomputes the intermediate values between checkpoints on demand from the nearest stored one. Checkpointing every \( \sqrt{n} \)-th layer of an \( n \)-layer network reduces activation memory from \( O(n) \) to \( O(\sqrt{n}) \) at the cost of one extra forward pass, an asymptotic sweet spot analyzed by Chen and colleagues. It is standard in every large-model training stack.

Newton and quasi-Newton methods

Newton's method for solving \( g(x) = 0 \) linearizes. Write \( g(x + \Delta) \approx g(x) + g'(x)\Delta \), set to zero to get \( \Delta = -g(x)/g'(x) \), and iterate \( x_{k+1} = x_k - g(x_k)/g'(x_k) \). For minimizing \( f \) this is applied to \( g = \nabla f \), giving the step \( x_{k+1} = x_k - [\nabla^2 f(x_k)]^{-1} \nabla f(x_k) \). Its signature property is quadratic convergence near a root. Let \( e_k = x_k - x^\star \) and Taylor-expand \( g \) about the root, using \( g(x^\star) = 0 \).

Dividing by \( g'(x_k) \) and substituting the Newton update \( x_{k+1} = x_k - g(x_k)/g'(x_k) \) gives

so the error is squared at each step. The number of correct digits roughly doubles per iteration once the iterate is close enough, provided \( g'(x^\star) \neq 0 \). The catch is the cost and conditioning of the Hessian. Forming and factoring \( \nabla^2 f \) is \( O(n^3) \) per step and impossible for a large model, and a near-singular Hessian makes the step unstable. Quasi-Newton methods, of which BFGS is the archetype, build an approximation to the inverse Hessian from successive gradient differences using the secant condition \( \nabla^2 f \, (x_{k+1} - x_k) \approx \nabla f(x_{k+1}) - \nabla f(x_k) \), and achieve superlinear convergence without ever forming the Hessian. L-BFGS keeps only the last few gradient-difference pairs and so needs \( O(n) \) memory. The convergence theory and the BFGS update itself are derived in the combinatorial optimization page. What matters numerically is that quasi-Newton buys most of Newton's fast local convergence while sidestepping the ill-conditioned, expensive Hessian solve.

Mixed precision arithmetic

Modern accelerators are fast in low precision and slow in high precision, and the gap is large. The H100 measured in this repository reaches \( 51.4 \) fp32 TFLOP/s on a \( 8192^3 \) matmul but \( 728.7 \) bf16 TFLOP/s and \( 700.8 \) fp16 TFLOP/s on the same shape, a factor of roughly fourteen, and \( 409.7 \) TFLOP/s in TF32, NVIDIA's 19-bit tensor-core format. Exploiting that speed means training in 16-bit formats, and the two available ones make opposite tradeoffs, which the standard model \( \text{fl}(x) = x(1 + \delta) \) quantifies through the size of \( \delta \).

float16 keeps ten significand bits but only five exponent bits, so its dynamic range stops at \( 6.55 \times 10^4 \) and its smallest normal number is \( 6.1 \times 10^{-5} \). bfloat16 sacrifices three significand bits, keeping only seven, to preserve the full eight-bit exponent of float32, so it spans the same \( 10^{38} \) range with coarser precision. The consequence for training is that fp16 gradients underflow. A gradient of magnitude \( 10^{-8} \), entirely plausible deep in a network, rounds to exactly zero in fp16, since it is below even the subnormal range. Loss scaling is the fix. Multiply the loss by a large constant \( S \) before the backward pass, which scales every gradient by \( S \) and lifts small ones back into the representable range, then divide the gradients by \( S \) before the optimizer step. Scaling by \( S = 2^{16} \) turns the \( 10^{-8} \) gradient into \( 6.55 \times 10^{-4} \), well inside fp16's range, and unscaling recovers \( 9.997 \times 10^{-9} \), the original value to four digits. bfloat16 needs no loss scaling because its exponent range already covers the gradient, though at the cost of precision. The same \( 10^{-8} \) gradient survives in bf16 as \( 3.004 \times 10^{-8} \), coarse but nonzero. On a \( 1024^3 \) matmul against an fp64 reference, the measured relative errors are \( 5.7 \times 10^{-7} \) in fp32, \( 3.6 \times 10^{-4} \) in fp16, and \( 2.9 \times 10^{-3} \) in bf16, exactly the ordering the significand bit counts predict. The standard recipe, from Micikevicius and colleagues, is to keep a master copy of the weights in fp32, compute the forward and backward passes in 16-bit, accumulate the matmuls into fp32 registers, and apply loss scaling for fp16. The accumulate-in-fp32 step is the same Kahan-style insight that a reduction must be done in wider precision than its inputs.

Worked problems

Implementation

The first block runs the CPU experiments in NumPy. It covers machine epsilon by bisection, catastrophic cancellation, Kahan summation, the Hilbert-matrix conditioning table, and the normal-equations-versus-QR comparison. Every printed number in the comments was produced by executing this code with a NumPy build whose bundled OpenBLAS is numerically reliable, and each factorization was cross-checked by a residual or a reconstruction, because this host once shipped a broken reference BLAS that returned plausible-looking wrong results for large factorizations.

import numpy as np
from scipy.linalg import qr, cholesky, hilbert, solve_triangular

# ---- machine epsilon by bisection ----------------------------------
eps = 1.0
while 1.0 + eps / 2.0 != 1.0:            # stop when the halving stops mattering
    eps /= 2.0
print(eps)                               # 2.220446049250313e-16 == 2**-52

# ---- catastrophic cancellation: (1 - cos x) / x^2 -> 1/2 -----------
for x in [1e-2, 1e-4, 1e-6, 1e-8]:
    naive  = (1.0 - np.cos(x)) / x**2            # subtracts near-equal terms
    stable = 2.0 * np.sin(x / 2.0)**2 / x**2     # 1 - cos x = 2 sin^2(x/2)
    print(x, naive, stable)                      # naive -> 0.0 at 1e-8, all digits gone

# ---- Kahan compensated summation in float32 ------------------------
def kahan(arr):
    s = np.float32(0.0); c = np.float32(0.0)
    for v in arr:
        y = v - c                 # bring in the running correction
        t = s + y                 # add; low bits of y may be lost
        c = (t - s) - y           # recover exactly the lost part
        s = t
    return s
small = np.full(1_000_000, 0.1, dtype=np.float32)
naive = np.float32(0.0)
for v in small:
    naive = np.float32(naive + v)
print(float(naive), float(kahan(small)))   # 100958.34 vs 100000.0 (true 100000)

# ---- Hilbert conditioning: forward error tracks kappa * eps --------
for n in [6, 10, 12]:
    H = hilbert(n)
    kappa = np.linalg.cond(H, 2)
    xtrue = np.ones(n)
    b = H @ xtrue
    xhat = np.linalg.solve(H, b)
    fwd = np.linalg.norm(xhat - xtrue) / np.linalg.norm(xtrue)
    resid = np.linalg.norm(H @ xhat - b) / (np.linalg.norm(H) * np.linalg.norm(xhat))
    print(n, f"{kappa:.2e}", f"{fwd:.2e}", f"{resid:.2e}")
    # 12 1.64e+16 1.29e-01 9.32e-17  -- forward error ~ kappa*eps, backward ~ eps

# ---- normal equations vs QR: cond(A)^2 vs cond(A) ------------------
m, n = 50, 12
rng = np.random.default_rng(3)
U0, _ = qr(rng.standard_normal((m, n)), mode="economic")
V0, _ = qr(rng.standard_normal((n, n)))
svals = np.logspace(0, 8, n)                    # condition number 1e8
A = (U0 * svals) @ V0.T
xtrue = rng.standard_normal(n)
b = A @ xtrue
# normal equations via Cholesky of A^T A -- squares the condition number
L = cholesky(A.T @ A, lower=True)
x_ne = solve_triangular(L.T, solve_triangular(L, A.T @ b, lower=True))
# QR -- works with A directly, condition number NOT squared
Q, R = qr(A, mode="economic")
x_qr = solve_triangular(R, Q.T @ b)
print(np.linalg.norm(x_ne - xtrue) / np.linalg.norm(xtrue))   # 1.23e-01
print(np.linalg.norm(x_qr - xtrue) / np.linalg.norm(xtrue))   # 4.26e-10

# ---- SVD low-rank approximation and Eckart-Young -------------------
U, s, Vt = np.linalg.svd(A, full_matrices=False)
for k in [2, 5, 10]:
    Ak = (U[:, :k] * s[:k]) @ Vt[:k, :]
    print(k, np.linalg.norm(A - Ak, 2), s[k])   # spectral error == sigma_{k+1}

The second block is the automatic differentiation comparison. A reverse-mode engine is built from scratch as a small tape of operations, each recording its parents and local derivatives, and a backward pass accumulates adjoints in reverse topological order. It is checked against PyTorch and JAX on the same nontrivial scalar function \( f(x, y, z) = \sin(xy) + e^{x - z}/\log(1 + y^2) - x/z \). All three agree to a unit in the last place. PyTorch and JAX are shown side by side because they represent the two dominant designs, an eager tape and a traced-and-compiled functional transform.

import math

# A minimal reverse-mode autodiff: a Wengert tape of Var nodes.
class Var:
    def __init__(self, value, parents=(), local_grads=()):
        self.value = float(value)
        self.parents = parents          # the inputs this node was built from
        self.local_grads = local_grads  # d(self)/d(parent) for each parent
        self.grad = 0.0                 # accumulated adjoint dL/d(self)
    def __add__(s, o):
        o = o if isinstance(o, Var) else Var(o)
        return Var(s.value + o.value, (s, o), (1.0, 1.0))
    def __sub__(s, o):
        o = o if isinstance(o, Var) else Var(o)
        return Var(s.value - o.value, (s, o), (1.0, -1.0))
    def __mul__(s, o):
        o = o if isinstance(o, Var) else Var(o)
        return Var(s.value * o.value, (s, o), (o.value, s.value))   # product rule
    def __truediv__(s, o):
        o = o if isinstance(o, Var) else Var(o)
        return Var(s.value / o.value, (s, o), (1.0 / o.value, -s.value / o.value**2))
    __radd__ = __add__

def vsin(x): return Var(math.sin(x.value), (x,), (math.cos(x.value),))
def vexp(x): return Var(math.exp(x.value), (x,), (math.exp(x.value),))
def vlog(x): return Var(math.log(x.value), (x,), (1.0 / x.value,))

def backward(node):
    topo, seen = [], set()
    def build(v):                       # reverse topological order
        if id(v) not in seen:
            seen.add(id(v))
            for p in v.parents:
                build(p)
            topo.append(v)
    build(node)
    node.grad = 1.0                     # seed dL/dL = 1
    for v in reversed(topo):
        for p, g in zip(v.parents, v.local_grads):
            p.grad += v.grad * g        # chain rule: adjoint flows to parents

x, y, z = Var(0.7), Var(1.3), Var(2.1)
out = vsin(x * y) + vexp(x - z) / vlog(Var(1.0) + y * y) - x / z
backward(out)
print(out.value)                        # 0.7053737397126023
print(x.grad, y.grad, z.grad)           # 0.5708823315  0.1862105559  -0.0904731746

import torch
torch.set_default_dtype(torch.float64)

x = torch.tensor(0.7, requires_grad=True)
y = torch.tensor(1.3, requires_grad=True)
z = torch.tensor(2.1, requires_grad=True)

out = torch.sin(x * y) + torch.exp(x - z) / torch.log(1 + y * y) - x / z
out.backward()                          # reverse-mode pass over the autograd graph

print(out.item())                       # 0.7053737397126023
print(x.grad.item(), y.grad.item(), z.grad.item())
# 0.5708823315009642 0.18621055585527035 -0.0904731746258266
# max abs difference from the from-scratch engine: 1.11e-16 (one ulp)

import jax
jax.config.update("jax_enable_x64", True)   # match float64 precision
import jax.numpy as jnp
from jax import grad

def f(v):
    x, y, z = v
    return jnp.sin(x * y) + jnp.exp(x - z) / jnp.log(1 + y * y) - x / z

g = grad(f)(jnp.array([0.7, 1.3, 2.1]))     # functional transform, not a tape
print(g)   # [ 0.57088233  0.18621056 -0.09047317 ]
# with x64 enabled, matches torch to 1.11e-16; the default float32 gives ~1e-7

The third block demonstrates mixed precision on the H100. It shows the fp16 underflow of a small gradient, its rescue by loss scaling, and the accuracy of matmul in each format measured against an fp64 reference. It is written for both frameworks. The reference is computed in float64 on the GPU, which goes through cuSOLVER and cuBLAS and is reliable, unlike this host's historical CPU BLAS.

import torch
dev = "cuda"

# format properties: significand precision vs dynamic range
for dt in (torch.float32, torch.bfloat16, torch.float16):
    fi = torch.finfo(dt)
    print(dt, fi.eps, fi.max, fi.tiny)
# fp16 eps 9.77e-04  max 6.55e+04  tiny 6.10e-05   (5 exponent bits)
# bf16 eps 7.81e-03  max 3.39e+38  tiny 1.18e-38   (8 exponent bits, coarse)

# a small gradient underflows fp16 to exactly zero
g = 1e-8
print(torch.tensor(g, dtype=torch.float16, device=dev).item())    # 0.0

# loss scaling lifts it back into range, then unscale
S = 2.0**16
scaled = torch.tensor(g * S, dtype=torch.float16, device=dev)
print(scaled.item() / S)                                          # 9.997e-09 recovered

# bf16 keeps the range, so no underflow, but coarse precision
print(torch.tensor(g, dtype=torch.bfloat16, device=dev).item())   # 3.004e-08

# matmul accuracy vs an fp64 reference (H100, modest 1024^3, cross-checked)
torch.manual_seed(0)
A = torch.randn(1024, 1024, device=dev, dtype=torch.float64)
B = torch.randn(1024, 1024, device=dev, dtype=torch.float64)
ref = A @ B
for dt in (torch.float32, torch.float16, torch.bfloat16):
    C = (A.to(dt) @ B.to(dt)).double()
    print(dt, (torch.norm(C - ref) / torch.norm(ref)).item())
# fp32 5.7e-07   fp16 3.6e-04   bf16 2.9e-03   (ordered by significand bits)

import jax, jax.numpy as jnp
import numpy as np

# format properties via ml_dtypes / numpy finfo
for dt in (jnp.float32, jnp.bfloat16, jnp.float16):
    fi = jnp.finfo(dt)
    print(dt, fi.eps, fi.max, fi.tiny)

# fp16 underflow of a small gradient and the loss-scaling rescue
g = 1e-8
print(jnp.asarray(g, dtype=jnp.float16))                 # 0.0
S = 2.0**16
print(jnp.asarray(g * S, dtype=jnp.float16) / S)         # ~9.997e-09 recovered
print(jnp.asarray(g, dtype=jnp.bfloat16))                # ~3.004e-08

# matmul accuracy vs fp64 reference on the accelerator
key = jax.random.PRNGKey(0)
A = jax.random.normal(key, (1024, 1024), dtype=jnp.float64)
B = jax.random.normal(jax.random.PRNGKey(1), (1024, 1024), dtype=jnp.float64)
ref = A @ B
for dt in (jnp.float32, jnp.float16, jnp.bfloat16):
    C = (A.astype(dt) @ B.astype(dt)).astype(jnp.float64)
    print(dt, float(jnp.linalg.norm(C - ref) / jnp.linalg.norm(ref)))
# same ordering: fp32 < fp16 < bf16 in relative error

How it is done in practice

The gap between these derivations and a deployed system is filled by libraries that have absorbed decades of numerical care so their users need not repeat it. Dense linear algebra runs through LAPACK and the BLAS. LAPACK's dgesv for general solves calls a blocked LU with partial pivoting, dpotrf for Cholesky, dgeqrf for Householder QR, and dgesdd for a divide-and-conquer SVD, all built on the level-3 BLAS dgemm that vendors tune to the last cache line. The blocking matters as much as the arithmetic. Casting a factorization as a sequence of matrix-matrix products lets it run at the machine's peak throughput, which on the H100 in this repository is \( 51.4 \) fp32 TFLOP/s and \( 728.7 \) bf16 TFLOP/s on large matmuls, against a small fraction of that for an unblocked, memory-bound implementation. This is also where numerical reliability can silently fail. This host once shipped a reference LAPACK that returned wrong answers for factorizations of size fifty and up while matrix multiply stayed correct, which is the worst failure mode because the wrong answers look plausible. The habit that catches it is to verify every factorization by a residual or a reconstruction, exactly the checks embedded in the code above.

At training scale, the numerical decisions are made once and baked into the framework. Automatic mixed precision in PyTorch and the equivalent in JAX keep an fp32 master copy of the weights, run the forward and backward passes in bf16 or fp16, and accumulate every matmul into fp32 registers on the tensor cores, which is why bf16 training reaches near-fp16 speed without the accuracy collapse that a naive all-bf16 computation would suffer. fp16 training adds a dynamic loss scaler that raises the scale when no overflow is detected and backs off when an inf or nan appears in the gradients. bf16, with its full exponent range, has largely displaced fp16 for large-model pretraining precisely because it removes the loss-scaling failure mode, trading precision that the optimizer's own noise floor makes irrelevant. Optimizer state is another quiet numerical choice. Adam's moment estimates are kept in fp32 even when the weights are bf16, because accumulating a running average in bf16 would suffer the summation error analyzed above, the same reason Kahan summation exists. Gradient checkpointing is turned on for the deepest models to hold activation memory to \( O(\sqrt{n}) \), and the conjugate-gradient and Lanczos routines inside second-order and natural-gradient methods run matrix-free, never forming the curvature matrix they implicitly invert.

The current research frontier

The active questions cluster around ever-lower precision and the numerics of scale. The 8-bit floating-point formats E4M3 and E5M2, standardized in a joint proposal from NVIDIA, Arm, and Intel and now native on Hopper and Blackwell tensor cores, push training below sixteen bits. The open problem is which tensors, weights, activations, gradients, tolerate four or five significand bits and which need per-tensor or per-block scaling factors to survive, a line of work carried by NVIDIA's Transformer Engine and by the microscaling (MX) formats from a cross-industry consortium including Microsoft and AMD. Below that, 4-bit and even ternary weight formats for inference, from the GPTQ and AWQ quantization line and Microsoft's BitNet, turn the whole model into a quantization problem where the Hessian of the loss with respect to the weights decides which bits can be discarded, a direct application of the conditioning theory in this page to the weight space rather than the input space.

On the factorization side, randomized numerical linear algebra, developed by Halko, Martinsson, and Tropp and extended by groups at Oxford, Michigan, and Berkeley, replaces deterministic \( O(n^3) \) factorizations with sketch-and-solve methods that touch the matrix through a few products with random test matrices, exactly the block power method idea, and come with probabilistic error bounds tight enough for production. Communication-avoiding and mixed-precision iterative refinement, associated with Higham and Demmel and their collaborators, solve a system in low precision and correct it with a few residual computations in higher precision, recovering full accuracy at a fraction of the cost, a technique now shipping in LAPACK's mixed-precision routines and directly relevant to the fp8-with-fp32-correction pattern in training. The autodiff frontier has moved to higher-order and structured derivatives, such as forward over reverse for Hessian-vector products, the source-to-source transforms in JAX and in Enzyme, which differentiates LLVM IR and so handles code the tracing frameworks cannot, and the differentiation of implicit functions and optimization layers, where the derivative of an argmin is obtained from the implicit function theorem on the optimality conditions rather than by unrolling the solver. The common thread is that the classical results, conditioning, stability, the SVD, the chain rule, are being re-derived in each new regime, which is why they repay being understood from first principles.

Open source to read

• numpy/numpy: the reference array library. Read numpy/linalg/linalg.py for the LAPACK bindings and numpy/core/src/umath for how reductions like pairwise summation are actually implemented.
• scipy/scipy: scipy/linalg exposes the full LAPACK surface, including the QR, Cholesky, and SVD drivers with control over pivoting and economy mode that NumPy hides. The place to go when the default solve is not stable enough.
• Reference-LAPACK/lapack: the canonical Fortran implementations. Read SRC/dgeqrf.f for blocked Householder QR and SRC/dpotrf.f for Cholesky. The header comments are a numerical-analysis course in themselves.
• google/jax: the cleanest modern autodiff. Start with jax/_src/interpreters/ad.py for the forward- and reverse-mode transforms and jax/_src/custom_derivatives.py for how custom vector-Jacobian products are registered.
• pytorch/pytorch: the tape-based autograd engine lives in torch/csrc/autograd. torch/autograd/function.py shows the Python surface, and tools/autograd/derivatives.yaml lists the local derivative rule for every operation.
• HIPS/autograd: the original NumPy-transparent reverse-mode library, small enough to read end to end. autograd/core.py is the entire engine and is the clearest tutorial on how a tape and vector-Jacobian products fit together.
• NVIDIA/TransformerEngine: production fp8 and mixed-precision training. Read the recipe and scaling-factor logic to see loss scaling and per-tensor scaling as they are actually deployed on Hopper and Blackwell.

Common misconceptions

"Double precision is accurate to sixteen digits, so rounding error never matters." The sixteen digits are the precision of a single stored number, and they say nothing about a computation. A single catastrophic cancellation can leave zero correct digits, as \( (1 - \cos x)/x^2 \) at \( x = 10^{-8} \) returned exactly zero against a true value of one-half. Accuracy is a property of the algorithm and the conditioning, not of the storage format.

"A small residual means the solution is accurate." The residual measures backward error, not forward error. The Hilbert system at \( n = 12 \) had a residual at machine precision, \( 10^{-16} \), while the solution itself had no correct digits, because the forward error is the residual amplified by \( \kappa \approx 10^{16} \). A small residual only certifies that the computed answer solves a nearby problem. Whether that is the right answer depends on the condition number.

"Solving the normal equations is the efficient way to do least squares." Forming \( A^\top A \) squares the condition number and can double the number of lost digits. On a problem with \( \kappa(A) = 10^8 \), the normal equations gave a relative error of \( 0.12 \) while QR on the same data gave \( 4 \times 10^{-10} \). The flops saved are never worth the digits lost.

"bfloat16 is just a less accurate float16." They differ in what they sacrifice. bfloat16 keeps float32's eight exponent bits and drops significand bits, so it has the full dynamic range and coarse precision. float16 keeps more significand bits but only five exponent bits, so it is more precise but overflows at \( 65504 \) and underflows small gradients to zero. That is why fp16 needs loss scaling and bf16 does not.

"Automatic differentiation is just the finite-difference approximation done automatically." Finite differences approximate the derivative with a step size and suffer both truncation bias and cancellation error, losing about half the available digits at the optimal step. Automatic differentiation applies the chain rule to exact local derivatives and is accurate to machine precision. The from-scratch engine above matched PyTorch to one unit in the last place, which no finite difference can do.

"A positive-definite covariance matrix will always factor with Cholesky." In exact arithmetic yes, but a matrix that is positive definite on paper can be numerically indefinite when it is ill-conditioned or estimated from too few samples, and the Cholesky attempt then fails at a negative square root. That failure is a feature. It is the cheapest reliable test that a matrix is usable, and it is exactly why regularizers add \( \lambda I \) to lift the smallest eigenvalue.

"Conjugate gradients is just a smarter step size for gradient descent." It is a different algorithm with a different convergence class. Its iteration count scales with \( \sqrt{\kappa} \) rather than \( \kappa \) because its search directions are \( A \)-orthogonal and it minimizes over an expanding Krylov subspace, giving the quadratic improvement demonstrated by \( 173 \) CG iterations against \( 5119 \) for steepest descent on the same \( \kappa = 1000 \) system.

Self-check

References

1. Trefethen, L. N., & Bau, D. (1997). Numerical Linear Algebra. SIAM. The clearest modern treatment of conditioning, stability, and the factorizations, and the source for the backward-error framing used throughout this page.
2. Golub, G. H., & Van Loan, C. F. (2013). Matrix Computations (4th ed.). Johns Hopkins University Press. The comprehensive reference for LU, Cholesky, QR, the SVD, and iterative methods.
3. Higham, N. J. (2002). Accuracy and Stability of Numerical Algorithms (2nd ed.). SIAM. The definitive account of rounding-error analysis, including the summation and cancellation bounds derived here. doi:10.1137/1.9780898718027
4. Nocedal, J., & Wright, S. J. (2006). Numerical Optimization (2nd ed.). Springer. Newton and quasi-Newton methods, the conjugate gradient derivation, and the convergence theory. doi:10.1007/978-0-387-40065-5
5. Griewank, A., & Walther, A. (2008). Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation (2nd ed.). SIAM. The reference for forward and reverse mode and the checkpointing memory/recompute tradeoff. doi:10.1137/1.9780898717761
6. Demmel, J. W. (1997). Applied Numerical Linear Algebra. SIAM. Complementary treatment of the same material with a strong emphasis on the LAPACK algorithms and perturbation theory.
7. Wilkinson, J. H. (1965). The Algebraic Eigenvalue Problem. Oxford University Press. The origin of modern backward-error analysis and the growth-factor bounds for Gaussian elimination.
8. Saad, Y. (2003). Iterative Methods for Sparse Linear Systems (2nd ed.). SIAM. Krylov subspaces, conjugate gradients, and the Chebyshev-polynomial convergence bound. doi:10.1137/1.9780898718003
9. Goldberg, D. (1991). What Every Computer Scientist Should Know About Floating-Point Arithmetic. ACM Computing Surveys, 23(1), 5-48. The standard introduction to IEEE 754 and its rounding model. doi:10.1145/103162.103163
10. IEEE. (2019). IEEE Standard for Floating-Point Arithmetic (IEEE Std 754-2019). The normative definition of the formats, rounding modes, and correctly-rounded operations. doi:10.1109/IEEESTD.2019.8766229
11. Kahan, W. (1965). Further remarks on reducing truncation errors. Communications of the ACM, 8(1), 40. The compensated-summation algorithm derived here. doi:10.1145/363707.363723
12. Householder, A. S. (1958). Unitary Triangularization of a Nonsymmetric Matrix. Journal of the ACM, 5(4), 339-342. The reflector and the stable QR factorization. doi:10.1145/320941.320947
13. Eckart, C., & Young, G. (1936). The approximation of one matrix by another of lower rank. Psychometrika, 1(3), 211-218. The optimality of the truncated SVD. doi:10.1007/BF02288367
14. Hestenes, M. R., & Stiefel, E. (1952). Methods of Conjugate Gradients for Solving Linear Systems. Journal of Research of the National Bureau of Standards, 49(6), 409-436. The original conjugate gradient paper. doi:10.6028/jres.049.044
15. Francis, J. G. F. (1961-1962). The QR Transformation, Parts I and II. The Computer Journal, 4(3-4). The dense eigenvalue algorithm sketched here. doi:10.1093/comjnl/4.3.265
16. Golub, G., & Kahan, W. (1965). Calculating the Singular Values and Pseudo-Inverse of a Matrix. SIAM Journal on Numerical Analysis, 2(2), 205-224. The bidiagonalization route to the SVD. doi:10.1137/0702016
17. Baydin, A. G., Pearlmutter, B. A., Radul, A. A., & Siskind, J. M. (2018). Automatic Differentiation in Machine Learning: a Survey. Journal of Machine Learning Research, 18(153), 1-43. arXiv:1502.05767
18. Micikevicius, P., Narang, S., Alben, J., et al. (2018). Mixed Precision Training. International Conference on Learning Representations (ICLR). Loss scaling and the fp32 master-weight recipe. arXiv:1710.03740
19. Chen, T., Xu, B., Zhang, C., & Guestrin, C. (2016). Training Deep Nets with Sublinear Memory Cost. The square-root-checkpointing analysis for reverse-mode memory. arXiv:1604.06174
20. Halko, N., Martinsson, P. G., & Tropp, J. A. (2011). Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions. SIAM Review, 53(2), 217-288. The foundation of randomized numerical linear algebra. doi:10.1137/090771806
21. Anderson, E., Bai, Z., Bischof, C., et al. (1999). LAPACK Users' Guide (3rd ed.). SIAM. The reference implementations behind NumPy, SciPy, and every dense solver discussed here.
22. Maclaurin, D., Duvenaud, D., & Adams, R. P. (2015). Autograd: Effortless Gradients in NumPy. ICML AutoML Workshop. The reverse-mode library whose core is small enough to read end to end, and the design ancestor of JAX. github.com/HIPS/autograd