SVD and Least Squares

In practice, $A\mathbf{x}=\mathbf{b}$ is often inconsistent — more equations than unknowns, and $\mathbf{b}\notin\mathcal{C}(A)$. We can't solve exactly, but we can do the next best thing: find the $\hat{\mathbf{x}}$ that minimizes the residual's squared length:

$\hat{\mathbf{x}} = \arg\min_{\mathbf{x}} \|A\mathbf{x} - \mathbf{b}\|^2.$

This is least squares. Geometrically: $A\hat{\mathbf{x}}$ is the orthogonal projection of $\mathbf{b}$ onto the column space $\mathcal{C}(A)$. The residual $\mathbf{b}-A\hat{\mathbf{x}}$ is orthogonal to every column of $A$:

$A^\top(\mathbf{b} - A\hat{\mathbf{x}}) = \mathbf{0}.$

The closest point in a subspace is always the orthogonal foot of the perpendicular — this geometric fact is the entire engine of least squares. The projection $\hat{\mathbf{b}}=A\hat{\mathbf{x}}$ is the "best approximation to $\mathbf{b}$ that $A$ can produce."

Rearranging $A^\top(\mathbf{b}-A\hat{\mathbf{x}})=\mathbf{0}$ gives the normal equations:

$A^\top A\,\hat{\mathbf{x}} = A^\top\mathbf{b}.$

If $A$ has full column rank, $A^\top A$ is invertible and the unique least-squares solution is:

$\hat{\mathbf{x}} = (A^\top A)^{-1}A^\top\mathbf{b}.$

The matrix $A^+ = (A^\top A)^{-1}A^\top$ is called the left pseudoinverse — it satisfies $A^+A=I$ when $A$ has full column rank.

Never solve the normal equations by forming $A^\top A$ and inverting it. $A^\top A$ squares the condition number of $A$: if $\kappa(A)=10^7$ then $\kappa(A^\top A)=10^{14}$, causing catastrophic loss of precision. In practice, use QR decomposition (np.linalg.lstsq) or SVD — both bypass $A^\top A$ entirely and are numerically stable.

The orthogonal projection of any vector $\mathbf{b}$ onto $\mathcal{C}(A)$ is $\hat{\mathbf{b}}=A\hat{\mathbf{x}}=A(A^\top A)^{-1}A^\top\mathbf{b}$. Define the projection matrix:

$P = A(A^\top A)^{-1}A^\top.$

$P$ has two defining properties:

• Idempotent: $P^2=P$ (projecting twice = projecting once).
• Symmetric: $P^\top=P$.

The complementary projector $I-P$ projects onto the orthogonal complement $\mathcal{N}(A^\top)$ (the left null space). For any $\mathbf{b}$:

$\mathbf{b} = \underbrace{P\mathbf{b}}_{\in\,\mathcal{C}(A)} + \underbrace{(I-P)\mathbf{b}}_{\in\,\mathcal{N}(A^\top)}.$

This is the orthogonal direct-sum decomposition of $\mathbb{R}^m$ — the four fundamental subspaces (Lesson 2) made concrete and computational. If $A=Q$ is already orthonormal, $P=QQ^\top$ and computing the projection is a single matrix multiply.

Every $m\times n$ matrix $A$ (no assumptions on shape or rank) factors as:

$A = U\Sigma V^\top$

where:

• $U\in\mathbb{R}^{m\times m}$ is orthogonal (left singular vectors — orthonormal basis of $\mathbb{R}^m$).
• $\Sigma\in\mathbb{R}^{m\times n}$ is diagonal, non-negative entries $\sigma_1\ge\sigma_2\ge\cdots\ge\sigma_r>0$ (singular values).
• $V\in\mathbb{R}^{n\times n}$ is orthogonal (right singular vectors — orthonormal basis of $\mathbb{R}^n$).

Geometric reading: $A$ performs three operations in sequence — $V^\top$ rotates in input space, $\Sigma$ stretches along $r$ axes (and kills the rest), $U$ rotates in output space. Every matrix is just stretch-and-rotate, once you find the right axes. The singular values $\sigma_i$ are the stretching factors; they are always real and non-negative even when $A$ is not symmetric.

The SVD delivers all four fundamental subspaces cleanly:

The SVD makes the abstract theory of Lesson 2 numerically concrete: just run np.linalg.svd(A) and read off the four subspaces from the columns of $U$ and $V$. The singular values are the "importance" of each direction — a near-zero singular value signals near-rank-deficiency, which is exactly what the determinant can't tell you reliably.

The best rank-$k$ approximation to $A$ in both the Frobenius norm and the spectral norm is:

$A_k = \sum_{i=1}^k \sigma_i \mathbf{u}_i \mathbf{v}_i^\top = U_k \Sigma_k V_k^\top$

where $U_k, V_k$ take the first $k$ columns and $\Sigma_k$ the top-left $k\times k$ block. The approximation error is:

$\|A - A_k\|_F = \sqrt{\sigma_{k+1}^2 + \cdots + \sigma_r^2}.$

This is the Eckart-Young theorem (1936) — no rank-$k$ matrix can approximate $A$ better.

Image compression. A $1000\times 1000$ grayscale image needs $10^6$ numbers. Taking $k=50$ singular values: store $U_{1000\times 50}$, $\Sigma_{50\times 50}$, $V_{50\times 1000}$ — just $105{,}000$ numbers, a $9.5\times$ compression. Visual quality is often indistinguishable from the original if the singular values decay rapidly.

The Moore-Penrose pseudoinverse $A^+$ exists for any matrix and satisfies four conditions (Moore-Penrose axioms):

$AA^+A=A, \quad A^+AA^+=A^+, \quad (AA^+)^\top=AA^+, \quad (A^+A)^\top=A^+A.$

Via SVD it's clean: flip the roles of $U$ and $V$, invert the nonzero singular values:

$A^+ = V\Sigma^+U^\top, \quad \Sigma^+ = \text{diag}(1/\sigma_1, \dots, 1/\sigma_r, 0, \dots, 0).$

Unified behavior:

• $A$ square invertible: $A^+ = A^{-1}$.
• $A$ full column rank: $A^+ = (A^\top A)^{-1}A^\top$ (left inverse).
• $A$ full row rank: $A^+ = A^\top(AA^\top)^{-1}$ (right inverse).
• $A$ rank-deficient: $A^+\mathbf{b}$ gives the minimum-norm least-squares solution — the unique solution with smallest $\|\hat{\mathbf{x}}\|$.

np.linalg.pinv(A) computes this via SVD with automatic tolerance for small singular values.

import numpy as np

# Overdetermined system: 5 equations, 2 unknowns
rng = np.random.default_rng(42)
m, n = 5, 2
A = rng.standard_normal((m, n))
true_x = np.array([3.0, -1.0])
b = A @ true_x + 0.1 * rng.standard_normal(m)  # noisy observations

# Stable least squares (QR internally)
x_hat, residuals, rank, sv = np.linalg.lstsq(A, b, rcond=None)
print(f"Estimated x: {x_hat}")     # close to [3, -1]
print(f"Singular values: {sv}")

# Full SVD
U, S, Vt = np.linalg.svd(A, full_matrices=True)
print(f"Singular values from SVD: {S}")

# Low-rank approximation of a matrix
B = rng.standard_normal((100, 80))
Ub, Sb, Vtb = np.linalg.svd(B, full_matrices=False)  # economy SVD
k = 10
Bk = (Ub[:, :k] * Sb[:k]) @ Vtb[:k, :]             # rank-k approx
error = np.linalg.norm(B - Bk, 'fro')
print(f"Frobenius error (rank-{k} approx): {error:.4f}")
print(f"Theoretical min: {np.sqrt(np.sum(Sb[k:]**2)):.4f}")

Use full_matrices=False (economy SVD) when $m\gg n$ — it skips the $m-n$ columns of $U$ that are orthogonal to $\mathcal{C}(A)$ and never used.

Two landmark uses of SVD in engineering:

Image / data compression. The singular values of natural images decay rapidly (smooth images have most energy in few singular vectors). Truncated SVD at rank $k$ gives near-perfect perceptual quality at a fraction of the storage. JPEG-style compression uses a discrete cosine transform (a fixed unitary basis), but SVD is optimal — it finds the data-specific best basis.

Recommender systems. Netflix / Spotify: build a users$\times$items rating matrix $R$ (mostly missing). Approximate $R\approx U_k\Sigma_k V_k^\top$ — the $k$ latent factors represent taste dimensions. User row $i$ of $U_k$ is a "taste fingerprint"; item row $j$ of $V_k$ is a "genre fingerprint". Predicted rating $\hat{R}_{ij}=\mathbf{u}_i^\top\Sigma_k\mathbf{v}_j$. The Netflix Prize was won with SVD-based matrix factorization.

The unifying theme: singular values rank-order importance. Chop below a threshold, keep what matters, discard the noise. This is information compression at its most principled.

Across this four-lesson arc you've built the complete engineer's linear algebra toolkit:

The SVD is the capstone: it reveals the eigenstructure of any matrix (not just square symmetric ones), delivers all four fundamental subspaces numerically, provides the best low-rank approximation provably, and unifies the pseudoinverse and least squares in a single computation. Every serious numerical algorithm in machine learning and scientific computing — PCA, linear regression, image compression, recommender systems — traces back to these ideas.