Guides / Linear algebra

Least squares and the SVD: solve AᵀAx̂ = Aᵀb, check the residual, and read the rank off the singular values

When a system Ax = b has more equations than unknowns, there is usually no exact solution, and least squares finds the x that makes ‖Ax − b‖ as small as possible. It does that by solving the square system A^TAx̂ = A^Tb, and the leftover residual comes out perpendicular to every column of A. The singular value decomposition, A = UΣV^T, explains all of it and still works when the columns are dependent. The mistakes come from reading singular values as eigenvalues, either of A^TA or of A itself, and from inverting A^TA when it has no inverse.

Why it works

Every vector Ax lies in the column space of A. The closest such vector to b is the perpendicular projection of b onto that space, and what is left over, r = b − Ax̂, must be perpendicular to every column: A^T(b − Ax̂) = 0. Rearranged, that is exactly the normal equations. If the columns of A are independent, A^TA is invertible and x̂ is unique.

For the SVD, A^TA is symmetric with nonnegative eigenvalues, so it has orthonormal eigenvectors v_i. Then ‖Av_i‖² = v_i^TA^TAv_i = λ_i, so A stretches v_i by σ_i = √λ_i, and the images Av_i are perpendicular to each other too. Normalizing them gives the u_i. A maps one orthonormal basis to another, stretching each direction by its own σ_i.

Example 1: fitting a line by least squares

Fit y = c + mx to the points (0, 6), (1, 0), (2, 0).

1. A = 101112, b = 600One row per point: c · 1 + m · x = y. Three equations, two unknowns.
2. A^TA = 3335, A^Tb = 60Entries of A^TA are dot products of columns: N = 3, Σx = 3, Σx² = 5. A^Tb holds Σy = 6 and Σxy = 0.
3. 3c + 3m = 6, 3c + 5m = 0 ⇒ m = −3, c = 5Subtract the equations: 2m = −6.
4. Ax̂ = (5, 2, −1), r = b − Ax̂ = (1, −2, 1)The fitted values and the residual.
5. A^Tr = (1 − 2 + 1, 0 − 2 + 2) = (0, 0), ‖r‖² = 6The residual is perpendicular to both columns, which confirms x̂. No other line through these points has a sum of squared errors below 6.

Example 2: an SVD by hand

Find the SVD of A = 221−1.

1. A^TA = 5335Column dot products: 4 + 1 = 5, 4 − 1 = 3, 4 + 1 = 5.
2. λ₁ = 8, λ₂ = 2Trace 10, determinant 16: λ² − 10λ + 16 = (λ − 8)(λ − 2). The unit eigenvectors are v₁ = (1, 1)/√2 and v₂ = (1, −1)/√2.
3. σ₁ = √8 = 2√2, σ₂ = √2Square roots of the eigenvalues.
4. u₁ = Av₁/σ₁ = (1, 0)Av₁ = (4, 0)/√2, divided by 2√2.
5. u₂ = Av₂/σ₂ = (0, 1)Av₂ = (0, 2)/√2, divided by √2. The u's come out perpendicular, as they must.
6. A = 10012√200√2 · 1√2111−1Check the product: the first row is 2√2 · (1, 1)/√2 = (2, 2) and the second is √2 · (1, −1)/√2 = (1, −1).

Example 3: a rank-one matrix and the pseudoinverse

For A = 122436 and b = (1, 0, 1), find the rank, the null space and the shortest least-squares solution.

1. A^TA = 14282856, eigenvalues 70 and 0The second column is twice the first, so the determinant is 14 · 56 − 28² = 0.
2. σ₁ = √70, σ₂ = 0, rank 1One nonzero singular value. v₁ = (1, 2)/√5 and u₁ = (1, 2, 3)/√14.
3. null space: multiples of v₂ = (2, −1)/√5A(2, −1) = (0, 0, 0). Adding any multiple of it to x leaves Ax unchanged, so least squares has infinitely many solutions.
4. x⁺ = A⁺b = v₁(u₁ · b)/σ₁ = (2/35)(1, 2)Only the σ₁ term survives: u₁ · b = 4/√14, and (4/√14)/√70 · (1, 2)/√5 = 4(1, 2)/70.
5. Ax⁺ = (2/7)(1, 2, 3)The projection of b onto the line through (1, 2, 3), as least squares requires. x⁺ has no component along the null space, which makes it the shortest of all the solutions.

Example 4: the best rank-one approximation

Find the best rank-one approximation to the A of Example 2, and its error.

1. A₁ = σ₁u₁v₁^T = 2√2 · (1, 0) · (1, 1)/√2Keep the largest term of the SVD.
2. A₁ = 2200An outer product: column (1, 0) times row (2, 2).
3. A − A₁ = 001−1, ‖A − A₁‖₂ = √2 = σ₂The error is the part that was dropped, σ₂u₂v₂^T, and its size is σ₂. No rank-one matrix gets closer: that is the Eckart–Young theorem, the basis of PCA and low-rank compression.

Common mistakes, and the exact line they happen on

1. The eigenvalues of AᵀA taken as the singular values

The third step of Example 2.

1. A^TA has eigenvalues 8 and 2Right so far.
2. σ₁ = 8, σ₂ = 2This is the mistake. A^TA applies A twice, so its eigenvalues are the squared stretches. Check with v₁: ‖Av₁‖ = ‖(4, 0)/√2‖ = 2√2 ≈ 2.83, not 8. The singular values are √8 and √2.

2. The eigenvalues of A taken as the singular values

The same matrix, without forming AᵀA.

1. A = 221−1, det(A − λI) = λ² − λ − 4Right so far, as a calculation of eigenvalues.
2. σ = |(1 ± √17)/2| ≈ 2.56 and 1.56This is the mistake. For a matrix that isn't symmetric, the eigenvectors aren't perpendicular and the eigenvalues don't measure stretching. The largest stretch is σ₁ = 2√2 ≈ 2.83, more than 2.56. The two only agree, up to sign, for a symmetric matrix.

3. Inverting AᵀA when the columns are dependent

The least-squares problem of Example 3.

1. A^TA = 14282856, A^Tb = (4, 8)Right so far.
2. x̂ = (A^TA)⁻¹A^TbThis is the mistake. det(A^TA) = 14 · 56 − 28² = 0, so there is no inverse. The normal equations reduce to the single equation 14x₁ + 28x₂ = 4, a whole line of solutions. In floating point the "inverse" exists only through rounding error and gives huge, meaningless weights. Use the pseudoinverse, or np.linalg.lstsq, which gives the shortest solution (2/35)(1, 2).

Practice

1. Fit y = mx, a line through the origin, to (1, 2), (2, 3), (3, 7).

   Show answer

   A is the column (1, 2, 3), so m = (A · b)/(A · A) = (2 + 6 + 21)/14 = 29/14 ≈ 2.071.

2. Find the constant c that best fits the values 3, 5, 10 in the least-squares sense.

   Show answer

   A is a column of ones, so the normal equation is 3c = 18 and c = 6, the mean. The residual (−3, −1, 4) adds to 0, which is A^Tr = 0.

3. Find the singular values of 300−2.

   Show answer

   3 and 2. A^TA = diag(9, 4); singular values are never negative.

4. Find the singular values, rank and null space of 1111.

   Show answer

   A^TA has eigenvalues 4 and 0, so σ = 2 and 0, and the rank is 1. The null space is spanned by (1, −1).

5. Find the best rank-one approximation to 3001, and its error.

   Show answer

   Keep σ₁ = 3: 3000. The error is σ₂ = 1.