Numerical Linear Algebra: Conditioning, Stability, and Iteration

On paper, a matrix is invertible or it is not. On a computer working in finite precision, that binary is the wrong question. Every number carries a rounding error near $10^{-16}$, and the right-hand side of a physical solve is itself measured or approximated. The question becomes: how much does the answer move when the data moves a little? One number answers it.

The condition number

The singular value decomposition writes $A = U\Sigma V^{\top}$ with $\Sigma = \operatorname{diag}(\sigma_1 \ge \dots \ge \sigma_n \ge 0)$. Geometrically $A$ maps the unit sphere to an ellipsoid whose semi-axis lengths are the $\sigma_i$. The condition number is the axis ratio,

$$\kappa(A) = \frac{\sigma_{\max}}{\sigma_{\min}} \ge 1.$$

It bounds error amplification exactly. Let $A x = b$ and perturb the right-hand side to $b + \delta b$, giving $x + \delta x$. Then $A\,\delta x = \delta b$, so $\delta x = A^{-1}\delta b$ and $\lVert\delta x\rVert \le \sigma_{\min}^{-1}\lVert\delta b\rVert$, while $\lVert b\rVert = \lVert A x\rVert \le \sigma_{\max}\lVert x\rVert$. Multiplying the two,

$$\frac{\lVert\delta x\rVert}{\lVert x\rVert} \;\le\; \kappa(A)\,\frac{\lVert\delta b\rVert}{\lVert b\rVert}.$$

A relative error of $10^{-16}$ in the data can surface as a relative error of $\kappa \cdot 10^{-16}$ in the solution. When $\kappa \approx 10^{8}$, half your digits are gone before you start; when $\kappa$ is small the ellipse is nearly a circle and errors barely grow.¹

The condition number is an axis ratiointeractive

a = 1.00b = 1.00c = 1.00d = 1.03

κ(A) = σ_max/σ_min

135.4

σ_max

2.015

σ_min

0.015

δb_x = 0.00δb_y = 0.12

x = A⁻¹b

(1.00, 0.00)

δx

(-4.00, 4.00)

amplification

66.7 ≤ κ

κ > 100: the ellipse is a near-line and a tiny δb swings x hard.

A nearly singular system, by hand

Take two almost-parallel rows,

$$A = \begin{bmatrix} 1 & 1 \\ 1 & 1.0001 \end{bmatrix}, \qquad \det A = 0.0001.$$

The determinant is nonzero, so $A$ inverts, but its singular values are $\sigma_{\max} \approx 2.0000$ and $\sigma_{\min} \approx 5.0\times 10^{-5}$, giving

$$\kappa(A) = \frac{2.0000}{5.0\times 10^{-5}} \approx 4.0\times 10^{4}.$$

Solve $A x = b$ for $b = (2,\, 2.0001)^{\top}$. Row two minus row one gives $0.0001\,x_2 = 0.0001$, so $x_2 = 1$ and $x_1 = 1$:

$$x = \begin{bmatrix} 1 \\ 1 \end{bmatrix}.$$

Now nudge only the second entry of $b$ by $\delta b = (0,\, 0.0001)^{\top}$, a relative change of $\lVert\delta b\rVert / \lVert b\rVert \approx 3.5\times 10^{-5}$. The same elimination now reads $0.0001\,x_2 = 0.0002$, so

$$x' = \begin{bmatrix} 0 \\ 2 \end{bmatrix}, \qquad \delta x = \begin{bmatrix} -1 \\ 1 \end{bmatrix}.$$

The solution moved by $100\%$ ($\lVert\delta x\rVert / \lVert x\rVert = 1$). The amplification is $1 / (3.5\times 10^{-5}) \approx 2.8\times 10^{4}$, safely under $\kappa$, and the bound $1 \le \kappa \cdot 3.5\times 10^{-5} \approx 1.41$ is nearly tight because $\delta b$ points along the smallest singular direction. A stray bit in the fourth decimal of $b$ rewrites the answer.

Direct is too expensive: iterate

A dense $LU$ factorization solves $A x = b$ in $O(n^{3})$. The systems in a physics simulation are large and sparse, often $10^{5}$ unknowns with a handful of nonzeros per row, so $n^{3}$ is out of reach and factorization would destroy the sparsity anyway. Split $A = D + L + U$ into its diagonal, strictly lower, and strictly upper parts, and iterate.

Jacobi updates every component from the previous sweep; Gauss-Seidel reuses each freshly computed value at once:

$$x^{k+1} = D^{-1}\!\big(b - (L+U)x^{k}\big), \qquad x^{k+1} = (D+L)^{-1}\!\big(b - U x^{k}\big).$$

Each is $x^{k+1} = R\,x^{k} + c$ for an iteration matrix $R$ ($R_{\mathrm J} = -D^{-1}(L+U)$, $R_{\mathrm{GS}} = -(D+L)^{-1}U$). The error obeys $e^{k+1} = R\,e^{k}$, so the method converges for every starting guess if and only if the spectral radius $\rho(R) < 1$, and the smaller $\rho(R)$ the faster.² For a symmetric positive-definite $A$, conjugate gradient does far better: its error in the $A$-norm shrinks by $\big((\sqrt{\kappa}-1)/(\sqrt{\kappa}+1)\big)^{k}$ per step, so it reaches a fixed tolerance in about $\sqrt{\kappa}$ iterations rather than $\kappa$. Conditioning sets the iteration count directly.

The numerical backbone

This is the machinery under a contact solver. The velocity-level contact system $A\lambda = b$ has $A = J M^{-1} J^{\top}$, the Delassus operator, which is symmetric positive semidefinite: exactly the shape conjugate gradient and its projected variants want. Projected Gauss-Seidel, the sequential-impulse method of Box2D and Bullet, is the Gauss-Seidel split above with a $\lambda \ge 0$ clamp after each component. Its convergence rate is $\rho(R)$, and $\rho$ grows as the contacts couple more tightly, so a tall stack or a large mass ratio drives $\kappa$ up and the solver needs more sweeps. The condition number is not an abstraction here; it is the number of iterations per frame.

Footnotes

1. Trefethen & Bau, Numerical Linear Algebra (1997). ↩

2. Golub & Van Loan, Matrix Computations (4th ed., 2013). ↩