Singular Value Decomposition and the Pseudoinverse

Purpose

The SVD is the decomposition that works for every matrix: no squareness, no symmetry, no full-rank requirement. This note derives it from the eigendecomposition of $A^{T}A$, defines the pseudoinverse through it, and states Eckart-Young, which is the license behind every truncated-SVD compression and PCA plot. It leans on Eigenvalues and Diagonalization and finishes the least-squares story from Orthogonality, Projections, and Least Squares.

Deriving the SVD from $A^{T}A$

Let $A\in {\mathbb{R}}^{m\times n}$. The matrix $A^{T}A\in {\mathbb{R}}^{n\times n}$ is symmetric and positive semidefinite, so by the spectral theorem it has an orthonormal eigenbasis $v_1,\dots ,v_n$ with real eigenvalues ${\lambda }_1\ge \cdots \ge {\lambda }_n\ge 0$. Define the singular values ${\sigma }_i=\sqrt{{\lambda }_i}$, and for each ${\sigma }_i>0$ define

These $u_i$ are orthonormal: $u_i^T{u}_j=\frac{v_i^T{A}^{T}A{v}_j}{{\sigma }_i{\sigma }_j}=\frac{{\lambda }_j{v}_i^T{v}_j}{{\sigma }_i{\sigma }_j}$, which is 1 when $i=j$ and 0 otherwise. Completing $\{u_i\}$ to an orthonormal basis of ${\mathbb{R}}^m$ and collecting everything in matrices gives

with $U\in {\mathbb{R}}^{m\times m}$ and $V\in {\mathbb{R}}^{n\times n}$ orthogonal and $\Sigma$ diagonal with ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge 0$. This is Strang’s construction in 18.06 lecture 29. The rank $r$ of $A$ is the number of nonzero singular values, and the compact form $A={\sum }_{i=1}^r{\sigma }_i{u}_i{v}_i^T$ writes $A$ as a sum of $r$ rank-one pieces, largest first.

Sanity check both eigenproblems: $A^{T}A=V{\Sigma }^T\Sigma V^T$ and $A{A}^T=U\Sigma {\Sigma }^T{U}^T$, so the right singular vectors diagonalize $A^{T}A$ and the left ones diagonalize $A{A}^T$, sharing the nonzero spectrum ${\sigma }_i^2$.

Geometry

The SVD says every linear map is rotate, stretch, rotate. $V^T$ rotates the input so the axes line up with $v_1,\dots ,v_n$, $\Sigma$ stretches axis $i$ by ${\sigma }_i$ (killing axes past rank $r$), and $U$ rotates into output space. The image of the unit sphere under any $A$ is a hyperellipse whose semi-axes have lengths ${\sigma }_1,\dots ,{\sigma }_r$ pointing along $u_1,\dots ,u_r$ (Trefethen and Bau open the whole subject with this picture in lecture 4). ${\sigma }_1=\Vert A{\Vert }_2$ is the largest stretch factor any unit vector experiences, and ${\sigma }_r$ the smallest nonzero one; their ratio ${\sigma }_1/{\sigma }_r$ is the condition number.

flowchart LR
    X["input x"] --> VT["V transpose: rotate input axes onto right singular vectors"]
    VT --> SIG["Sigma: stretch axis i by sigma i, drop axes past rank r"]
    SIG --> U["U: rotate onto left singular vectors in output space"]
    U --> Y["output Ax"]
    style VT fill:#e3f2fd
    style SIG fill:#e8f5e9
    style U fill:#e3f2fd

The Moore-Penrose pseudoinverse

Invert what is invertible and leave the rest alone:

$A^{+}$ is the unique matrix satisfying the four Penrose conditions: $A{A}^{+}A=A$, $A^{+}A{A}^{+}=A^{+}$, and both $A{A}^{+}$ and $A^{+}A$ symmetric (Golub and Van Loan, ch. 5). When $A$ is square and invertible, $A^{+}=A^{-1}$. When $A$ has full column rank, $A^{+}=(A^{T}A{)}^{-1}{A}^T$, the least-squares operator from the normal equations.

The general statement: $\hat{x}=A^{+}b$ is always a least-squares solution of $Ax=b$, and among all least-squares solutions it is the one with minimum norm. Rank deficiency makes the minimizer non-unique, adding any null-space component leaves the residual unchanged, and $A^{+}b$ is the choice with zero null-space component. This is exactly what np.linalg.lstsq returns, since its SVD-based driver applies a truncated pseudoinverse (numpy.linalg.lstsq docs).

SVD, pseudoinverse, and least squares are one package

Undo the three-step geometry in reverse: rotate with $U^T$, divide by the nonzero ${\sigma }_i$ (zero out the rest), rotate back with $V$. That is $A^{+}=V{\Sigma }^{+}{U}^T$, and it solves least squares in every rank case at once:

• Square invertible $A$: $A^{+}=A^{-1}$, the exact solution.
• Full column rank: $A^{+}=(A^{T}A{)}^{-1}{A}^T$, the normal-equations solution.
• Rank deficient: $A^{+}b$ is the minimum-norm least-squares solution, the unique minimizer with no null-space component.

The normal equations only cover the middle case; the SVD route degrades gracefully through all three.

Eckart-Young: best low-rank approximation

Truncate the rank-one sum after $k$ terms: $A_k={\sum }_{i=1}^k{\sigma }_i{u}_i{v}_i^T$. The Eckart-Young-Mirsky theorem says no rank-$k$ matrix does better:

with $A_k$ attaining both (Mirsky extended the result to every unitarily invariant norm). The singular value spectrum is therefore a complete answer to “how compressible is this matrix”: if the ${\sigma }_i$ decay fast, a few rank-one terms capture almost everything, and ${\sigma }_{k+1}$ tells you exactly the error of stopping at $k$.

PCA is this theorem applied to data. Center the data matrix $X\in {\mathbb{R}}^{n\times p}$ (subtract column means), take its SVD, and the right singular vectors $v_i$ are the principal components; the variance along component $i$ is ${\sigma }_i^2/(n-1)$, and the fraction of variance explained by the top $k$ components is ${\sum }_{i\le k}{\sigma }_i^2/{\sum }_i{\sigma }_i^2$ (Stanford STATS 305C). The repo has a small image-compression demo of exactly this truncation in content/math/linear-algebra/pca_image_compression.py.

NumPy

import numpy as np
 
rng = np.random.default_rng(0)
A = rng.standard_normal((6, 4)) @ np.diag([10.0, 3.0, 0.5, 0.01])
 
U, s, Vh = np.linalg.svd(A, full_matrices=False)
print(s)                                  # descending: [25.57  5.64  0.95  0.01]
print(np.allclose(U @ np.diag(s) @ Vh, A))  # True
 
k = 2                                     # rank-2 truncation
Ak = U[:, :k] @ np.diag(s[:k]) @ Vh[:k]
print(np.linalg.norm(A - Ak, 2), s[k])    # equal, per Eckart-Young
 
Apinv = Vh.T @ np.diag(1/s) @ U.T         # pseudoinverse by hand
print(np.allclose(Apinv, np.linalg.pinv(A)))  # True

Conventions from the numpy.linalg.svd docs that trip people up: the function returns Vh $=V^T$, not $V$; singular values come back sorted descending as a 1-D array, so reconstruction needs np.diag; and full_matrices=True (the default) returns square $U$ and $V^T$ padded with basis vectors for the null spaces, while full_matrices=False gives the compact economy shapes. np.linalg.pinv applies an rcond cutoff before inverting singular values, because dividing by a tiny ${\sigma }_i$ that is really numerical noise amplifies error by $1/{\sigma }_i$; treating near-zero singular values as exactly zero is the stable choice for ill-conditioned systems.

Sources

• MIT 18.06 Linear Algebra, Spring 2010 (Strang), lecture 29
• Interactive Linear Algebra (Margalit and Rabinoff)
• numpy.linalg.svd documentation
• Stanford STATS 305C, PCA lecture
• Trefethen and Bau, Numerical Linear Algebra (SIAM, 1997), lectures 4-5
• Golub and Van Loan, Matrix Computations (4th ed.), chapter 5

Related notes

• Eigenvalues, Eigenvectors, and Diagonalization
• Orthogonality, Projections, and Least Squares
• Matrix Theory