The condition number, and how many digits a problem throws away

Some problems amplify error

Give a computation an input that is wrong in its seventh digit. Some computations return an answer wrong in the seventh digit; others return an answer wrong in the first. The second kind are ill-conditioned, and the condition number is how much they amplify the error. Nothing about the algorithm or the hardware changes this. It is a property of the problem.

For a matrix, module 2's singular values give the number directly:

κ = σ_max / σ_min

the ratio of the largest stretch the matrix applies to the smallest. A matrix that stretches one direction by 1,000 and another by 0.001 has κ = 10^6. The identity has κ = 1. A singular matrix, which flattens some direction entirely, has κ = ∞.

The rule of thumb

Solving a linear system, fitting a least-squares model, or inverting a matrix with condition number κ loses about log10(κ) digits of accuracy. Float32 has seven; float64 has sixteen.

κ = 10^2   float32 keeps 5 digits   float64 keeps 14
κ = 10^6   float32 keeps 1 digit    float64 keeps 10
κ = 10^8   float32 keeps nothing    float64 keeps 8
κ = 10^16  float64 keeps nothing

The Hilbert matrix, whose entries are 1/(i + j − 1), is the classic case: the 5 × 5 version has κ ≈ 5 × 10^5, the 10 × 10 version κ ≈ 10^13. Solving a 10 × 10 Hilbert system in float64 gives three correct digits; in float32, garbage with no warning. Module 2's advice to solve rather than invert stands, but it does not rescue a problem this badly conditioned. Nothing does, except changing the problem.

Where it comes from in practice: scale

The commonest ill-conditioned matrix in machine learning is the one you build from unscaled features. Take a regression with one feature in the range 0 to 1 and another in the range 0 to 1,000,000. The matrix XᵀX from module 2's normal equations has one direction stretched by about 10^12 relative to the other: κ ≈ 10^12. Solving it in float64 leaves four digits. Solving it in float32 leaves none, and the coefficient on the small feature comes out as noise.

Standardise both features to mean 0 and spread 1, and κ drops to the ratio set by their correlation, usually under 10. Same data, same model, twelve orders of magnitude of conditioning recovered by subtracting a mean and dividing by a spread. This is the numerical reason behind the feature-scaling advice in machine-learning-foundations; the statistical reasons are there, and this is the arithmetic one.

Polynomial features are worse. Fitting x, x², …, x^10 on x between 0 and 1,000 gives columns ranging from 10^3 to 10^30, a condition number beyond float64's reach, and coefficients that are pure rounding error. Rescale x to [−1, 1] first, or use an orthogonal polynomial basis, which is the same fix by another name.

Gradient descent feels it too

Module 3's curvature lesson showed that plain gradient descent zigzags when the loss is much more curved in one direction than another, and that the ratio of curvatures sets the number of steps. That ratio is the condition number of the Hessian. A quadratic loss with κ = 10^4 needs on the order of 10^4 steps of plain gradient descent to converge to a fixed accuracy; with κ = 10 it needs about ten. Feature scaling, normalisation layers, and adaptive optimisers such as Adam are all, from this angle, ways of reducing the condition number the optimiser sees. The numerical and the optimisation problems are one problem.

Correlated features

Two columns that are nearly copies of each other, a price in dollars and the same price in cents, make XᵀX nearly singular, because the matrix cannot tell them apart. The condition number rises without limit as the correlation approaches 1, and the fitted coefficients become huge, opposite in sign and meaningless, while the predictions stay fine. The L2 penalty of module 6 adds λ to every singular value, which bounds κ by σ_max / λ and is the numerical reason ridge regression was invented, before anyone described it as a prior.

Computing it

import numpy as np
X = np.column_stack([np.random.rand(1000), 1e6 * np.random.rand(1000)])
print(np.linalg.cond(X))                 # about 1e6 for X, 1e12 for X.T @ X
Xs = (X - X.mean(0)) / X.std(0)
print(np.linalg.cond(Xs))                # a few

np.linalg.cond computes the singular-value ratio. Print it for any matrix you are about to solve with, and read log10 of the result as the digits you are about to lose. Above 10^8 in float32 or 10^15 in float64, the answer the solver returns is not an approximation of the truth; it is an artefact of rounding, and the solver will not say so.