Norms, singular values and conditioning

Feed every unit vector to A. The outputs trace an ellipse (a line segment if A is singular). The longest output is the size of A, written |A|: the most that A can stretch anything. It is the right yardstick for a map because it bounds every single output: no vector is stretched by more than |A|.

The inequality. For z = 0 it is trivial. Otherwise z = |z|u with u a unit vector, so |Az| = |z||Au| ≤ |z||A|. Products. |CAx| ≤ |C||Ax| ≤ |C||A||x|; take the maximum over unit x. Transpose. (Ax, v) = (x, A^Tv), and |A| is the maximum of (Ax, v) over unit x and v (Lax p.90, from Theorem 2 of Chapter 7); the right-hand side gives |A^T| in the same way.

Theorem 15. Put C = A − B, so B = A − C = A(I − A^-1C), and |A^-1C| ≤ |A^-1||C| < 1. If Bx = 0 then x = A^-1Cx, so |x| ≤ |A^-1C||x|, which forces x = 0. So B has no kernel, and a square matrix with no kernel is invertible. The argument is Lax’s.

The eigenvalue formula. |Ax|² = (A^TAx, x). By the spectral theorem of Chapter 8, A^TA is symmetric, so this is at most λ|x|², where λ is its largest eigenvalue, with equality at a corresponding eigenvector. For 2×2, A^TA has trace F and determinant det(A^T)det(A) = D², so λ is the larger root of t² − Ft + D².

Lean: opNorm_le and opNorm_comp_le are Mathlib’s operator-norm lemmas, stated for real normed spaces; compiled clean (axioms: propext, Classical.choice, Quot.sound only). The eigenvalue formula, Theorem 15 and |A^T| = |A| are not stated in Lean. The script checks, exactly: for all 625 integer 2×2 matrices with entries −2..2 and 80 integer vectors (50,000 cases), |Ax|² ≤ λ|x|², comparing without square roots; where λ is an integer (273 matrices) it is attained at an eigenvector of A^TA; A A^T and A^TA have the same F and the same D². It also checks, with the 1-norm (the largest column sum, a different norm, for which the product rule still holds because the proof used nothing special about the Euclidean one), |CA| ≤ |C||A| for all 390,625 pairs, and Theorem 15 in 40,176 cases. The bound in Theorem 15 cannot be improved: A = I and B = I − E₁₁ has |A − B| = 1 = 1/|A^-1| and B is singular; the script finds 568 such equality cases among its inputs.

A real number is a sign times a size. A matrix is a stretch times a turn: first rotate or reflect (an orthogonal map U, which preserves lengths), then stretch along perpendicular axes (a symmetric map R with non-negative eigenvalues). Choosing the axes, once in the input space and once in the output space, makes A diagonal: A = W D V, with W and V orthogonal and D diagonal. The diagonal entries of D are the singular values, the semi-axes of the ellipse of section 1.

Polar decomposition, A invertible (real case). (AA^Tx, x) = |A^Tx|² > 0 for x ≠ 0, so AA^T is a positive symmetric matrix and has a unique positive square root R, R² = AA^T (Lax cites his Theorem 1 for this on p.169). Put U = R^-1A. Then UU^T = R^-1AA^TR^-1 = R^-1R²R^-1 = I, so U is orthogonal, and A = RU. For singular A Lax builds the isometry on the range of R and extends it (p.170); that argument is not reproduced here.

Singular values. By the spectral theorem R = WDW^T with W orthogonal and D diagonal. Then A = WDW^TU = WDV with V = W^TU orthogonal. The eigenvalues of R² = AA^T are the squares of the entries of D.

Exercise 14. If A^TAx = λx with λ ≠ 0, then Ax ≠ 0 (else A^TAx = 0) and AA^T(Ax) = λAx. So A carries the λ-eigenspace of A^TA injectively into that of AA^T; by the same argument A^T carries it back, so the two eigenspaces have equal dimension.

Lean: gram_det states det(A^TA) = det(A)², the product of the squared singular values. gram_same_nonzero_spectrum states t^m·charpoly(A^TA) = tⁿ·charpoly(AA^T) for any real m×n matrix, from Mathlib’s lemma for products in both orders: it says the two have the same non-zero roots with multiplicity; compiled clean, with the same axioms. The existence of the decomposition is not stated in Lean. The script builds A = W D V^T from the rotations with cosine and sine 3/5 and 4/5, 5/13 and 12/13, 8/17 and 15/17, 7/25 and 24/25 (and the identity), and D = diag(s₁, s₂) with 0 ≤ s₂ ≤ s₁ ≤ 5, so every entry is a fraction. For all 525 matrices it checks exactly A^TA V = V D², AA^T W = W D², A = RU with R = W D W^T symmetric and U = W V^T orthogonal. That checks the theorem on a family where it can be written in fractions; it is not a proof. Exercise 14 is checked on all 729 3×2 matrices with entries −1, 0, 1: the characteristic polynomial of AA^T is t times that of A^TA. And det(A^TA) = det(A)² on all 625 integer 2×2 matrices.

Solve Ax = b where b carries a small error. If A squeezes the circle into a thin ellipse, then an error in b along the short axis comes from a large change in x, because A made that direction small. The worst case is a solution x that A stretches the most, together with an error in b along the direction A stretches the least. The magnification is then the longest axis divided by the shortest: K = |A||A^-1|. It is at least 1, and equal to 1 only for a map that treats all directions alike.

Sensitivity. Subtracting the two equations gives Aδx = δb, so δx = A^-1δb and |δx| ≤ |A^-1||δb|. Also |b| = |Ax| ≤ |A||x|. Multiply: |δx||b| ≤ |A||A^-1||δb||x|, which is the claim with the denominators cleared. Attained: take x with |Ax| = |A||x| and δb with |A^-1δb| = |A^-1||δb|; both inequalities become equalities.

(i) I = AA^-1 and |I| = 1, so 1 ≤ |A||A^-1|. (ii) If Ax = μx then |μ||x| = |Ax| ≤ |A||x|, so every eigenvalue has |μ| ≤ |A|; in particular β ≤ |A|. The eigenvalues of A^-1 are the 1/μ, so the same argument gives 1/α ≤ |A^-1|. Multiply. (iii) From A = WDV, A^-1 = V^TD^-1W^T has singular values 1/σ_i, so |A| = σ_max and |A^-1| = 1/σ_min.

(iv) Write λ₁ ≥ λ₂ for the eigenvalues of A^TA, so λ₁ + λ₂ = F and λ₁λ₂ = D², and K² = λ₁/λ₂. Then (K + 1/K)² = K² + 2 + K^-2 = (λ₁ + λ₂)²/(λ₁λ₂) = F²/D². Since K + 1/K ≥ 2 with equality only at K = 1, K = 1 exactly when F = 2|D|, and by the two sum-of-squares identities that happens exactly when (a − d)² + (b + c)² = 0 (so D > 0) or (a + d)² + (b − c)² = 0 (so D < 0), with A not zero. (v) A^TA = V^TD²V has singular values σ_i², so its ratio is (σ_max/σ_min)².

Lean: condition_ge_one states 1 ≤ |A||B| whenever B is a left inverse of A (B·A = id) on a non-trivial real normed space. sensitivity states |δx||Ax| ≤ (|A||B|)(|Aδx||x|) with b = Ax and δb = Aδx, so it carries no division and holds when x or b is zero. frobenius_det_bounds states the two sum-of-squares identities and 2|D| ≤ F. trace_sq_gram states trace((A^TA)²) = F² − 2D² written out in the four entries, which is what turns (iv) into (v) for 2×2. All compiled clean. The script: K from the definition |A||A^-1| equals σ_max/σ_min for all 2112 invertible integer matrices with entries −3..3 (2401 matrices in all: 289 singular and 2112 invertible), in 60-digit decimal because it needs square roots; K + 1/K = F/|D| for all of them; K ≥ 1 for all, and 96 have K = 1 exactly where the integer test says; trace((A^TA)²) = F² − 2D² exactly and K(A^TA) = K(A)²; (ii) is checked exactly with the 1-norm on all 2112. Sensitivity is checked exactly (1-norm, cross-multiplied) on 285,696 cases, and the bound is attained for A = diag(1, 1/m), x = e₁, δb along e₂, with ratio exactly m = K, for m = 1 to 20. One caution from Lax p.248: the bound is a worst case. With k-digit arithmetic the relative error in x can be as large as 10^-kK(A); it is often smaller.

A tall matrix has no inverse, but it has a best substitute: the map that sends b to the least-squares answer of Chapter 8. That map is the pseudoinverse A⁺. In the singular value decomposition it is easy to write: undo each stretch that is not zero, and leave the zero directions alone. The classical way to compute the answer, solving A^TAc = A^Tb, squares every singular value, and so squares the condition number: half the digits are gone before the arithmetic starts.

Independent columns. A^TA is invertible (Chapter 8), so (A^TA)^-1A^TA = I. Then AA⁺ = A(A^TA)^-1A^T is the orthogonal projection onto the range of A (Chapter 8, section 3): it is symmetric and fixes every vector of the range, so AA⁺A = A. Since A⁺A = I, which is symmetric, A⁺AA⁺ = A⁺. Rank one. AA⁺A = uv^Tvu^Tuv^T/(|u|²|v|²) = u(v^Tv)(u^Tu)v^T/(|u|²|v|²) = uv^T. Condition. Section 3 (v): the singular values of A^TA are the squares of those of A. The SVD formula for A⁺ is stated here as the standard one; its proof is not given on this page.

Lean: left_inverse_normal states that (A^TA)^-1A^T is a left inverse of A when det(A^TA) is a unit; compiled clean. The Penrose conditions and the SVD formula are not stated in Lean. The script checks, exactly: the four Penrose conditions for 200 random integer rank-one 3×2 matrices with the rank-one formula; for 200 random integer 4×2 matrices of full column rank, that (A^TA)^-1A^T is a left inverse and satisfies all four conditions; and, on the 525 matrices A = W D V^T of section 2 (including the singular and zero ones), that V D⁺ W^T satisfies all four. The chapter 8 four-point fit gives A⁺b = (9/10, 9/10), the fitted line y = 0.9 + 0.9t. For that matrix A^TA = [[4,6],[6,14]] has F = 284 and D = 20, so K(A^TA) = 14.13 and K(A) = 3.76: here the squaring costs about one digit, which is harmless. It costs half of all digits when K(A) is large.

The Hilbert matrix H_n has entries 1/(i + j − 1), the integrals of x^i+j-2 over [0, 1]. It looks harmless and has an inverse made of whole numbers, but its columns are nearly parallel, so its ellipse is almost flat and its condition number grows very fast.

By direct computation in exact fractions: H₃H₃^-1 = I, the column sums of H₃^-1 in absolute value are 75, 408 and 390, and those of H₃ are 11/6, 13/12 and 47/60; 11/6 × 408 = 748. Nothing deeper is claimed; the growth for general n is observed here, not proved.

Lean: control_hilbert3_inverse states that the integer matrix above times H₃ is the identity over the rationals; a direct computation; compiled clean on the third try: the first proof left vector entries unevaluated, a first replacement left !![1,0,0;0,1,0;0,0,1] = 1 unsolved, and the second replacement closed it. The script verifies all of the above exactly for n = 1 to 8. It also solves H_nx = H_n(1, …, 1) in ordinary floating point (numpy, double precision, not exact) for n = 3, 5, 8, where the 2-norm condition numbers are about 5×10², 5×10⁵ and 1.5×10¹⁰. The relative errors of the computed solution are of order 10^-14, 10^-12 and 10^-7, and the script asserts only that each is below 2.3 × 10^-15 times K, a stand-in for Lax’s 10^-kK with double precision (about 16 digits) and a margin of 10. These floating-point digits depend on the machine and the library, so they are illustrations, not part of the results.