Least Squares and the Pseudoinverse

Least squares finds the vector

\overset{x}{^}

minimizing

∥ A x - b ∥^{2}

when

A x = b

has no exact solution. The unique minimum-norm solution is

\overset{x}{^} = A^{+} b

, where

A^{+}

is the Moore–Penrose pseudoinverse. Applies to overdetermined/underdetermined linear systems, data fitting, regression, and signal processing. Originated with Gauss and Legendre (c. 1795–1805); pseudoinverse formalized by Moore (1920) and Penrose (1955).

Given a (usually rectangular or rank-deficient) matrix $A \in R^{m \times n}$ and a vector $b \in R^{m}$ , the system $A x = b$ generally has no solution when $m > n$ (too many constraints) and infinitely many when $m < n$ (too few). The least squares method replaces the demand for an exact solution with the demand for the best approximate one: the $x$ that makes the residual $A x - b$ as small as possible in Euclidean length. The Moore–Penrose pseudoinverse $A^{+}$ is the single matrix that, applied to $b$ , always returns the correct least-squares answer — and among all least-squares answers, the one of smallest norm. The surprising fact is that this one object simultaneously inverts square invertible matrices, solves overdetermined systems by projection, and solves underdetermined systems by minimum-norm selection — a unified generalization of matrix inversion that always exists and is unique.

The Core Insight

You cannot always solve $A x = b$ , but you can always project $b$ onto the column space of $A$ and solve that. The residual of the best fit is perpendicular to everything $A$ can produce — this orthogonality is the whole theory in one geometric picture. The pseudoinverse is just the bookkeeping that makes "project, then solve, then pick the smallest" into a single matrix.

Formal Statement

The Least Squares Problem

Given $A \in R^{m \times n}$ and $b \in R^{m}$ , the least squares problem is to find $\overset{x}{^} \in R^{n}$ minimizing the squared Euclidean residual:

\overset{x}{^} = ar g x \in R^{n} min ∥ A x - b ∥_{2}^{2}

⟳This is the definition;

∥ v ∥_{2}^{2} = \sum_{i} v_{i}^{2} = v^{⊤} v

is the squared length of the residual vector

v = A x - b

Symbols: $A$ is the design/coefficient matrix ( $m$ rows = equations, $n$ columns = unknowns); $b$ is the observation vector; $x$ is the parameter vector; $r = A x - b$ is the residual; $∥ \cdot ∥_{2}$ is the Euclidean norm.

The Normal Equations

Any minimizer $\overset{x}{^}$ of $∥ A x - b ∥^{2}$ satisfies the normal equations:

A^{⊤} A \overset{x}{^} = A^{⊤} b

⟳Set the gradient

\nabla_{x} ∥ A x - b ∥^{2} = 2 A^{⊤} (A x - b)

equal to zero.

Here $A^{⊤} A \in R^{n \times n}$ is symmetric and positive semidefinite. If $A$ has full column rank ( $rank A = n$ ), then $A^{⊤} A$ is invertible and the solution is unique:

\overset{x}{^} = (A^{⊤} A)^{- 1} A^{⊤} b

⟳Left-multiply [2] by

(A^{⊤} A)^{- 1}

, which exists exactly when

A

's columns are linearly independent.

The Moore–Penrose Pseudoinverse

The Moore–Penrose pseudoinverse $A^{+} \in R^{n \times m}$ of any $A \in R^{m \times n}$ is the unique matrix satisfying the four Penrose conditions:

A A^{+} A = A, A^{+} A A^{+} = A^{+}, (A A^{+})^{⊤} = A A^{+}, (A^{+} A)^{⊤} = A^{+} A

⟳These four equations uniquely determine

A^{+}

; the first two are "weak inverse" axioms, the last two force the products

A A^{+}

and

A^{+} A

to be orthogonal projectors (symmetric).

The general least-squares / minimum-norm solution is then

\overset{x}{^} = A^{+} b

⟳

A^{+}

projects

b

onto

col (A)

and inverts on the row space, automatically selecting the smallest-norm minimizer.

Existence, Uniqueness, and Minimum Norm

For every $A$ and $b$ :

A least-squares minimizer of $∥ A x - b ∥^{2}$ always exists.
The set of minimizers is ${\overset{x}{^} + z : z \in ker A}$ — a single point iff $A$ has full column rank.
Among all minimizers, $A^{+} b$ is the unique one of smallest Euclidean norm $∥ x ∥$ .

Pseudoinverse via the SVD

If $A = U Σ V^{⊤}$ is a singular value decomposition with nonzero singular values $σ_{1}, \dots, σ_{r}$ , then

A^{+} = V Σ^{+} U^{⊤}, Σ^{+} = diag (\frac{1}{σ _{1}}, \dots, \frac{1}{σ _{r}}, 0, \dots, 0)

⟳Substitute

A = U Σ V^{⊤}

into the four Penrose conditions [4] and verify, using

U^{⊤} U = I

V^{⊤} V = I

, and that

Σ^{+}

inverts the nonzero block while zeroing the rest.

Why It Holds / Derivation

The Normal Equations from Orthogonality

Step 1. Define

f (x) = ∥ A x - b ∥^{2} = (A x - b)^{⊤} (A x - b) = x^{⊤} A^{⊤} A x - 2 b^{⊤} A x + b^{⊤} b

. Step 2. Differentiate with respect to

x

\nabla f (x) = 2 A^{⊤} A x - 2 A^{⊤} b .

Since

f

is convex (its Hessian

2 A^{⊤} A ⪰ 0

), any stationary point is a global minimum. Step 3. Set

\nabla f = 0

A^{⊤} A x = A^{⊤} b

, which is equation [2]. Step 4 (geometric meaning). Rewrite as

A^{⊤} (b - A \overset{x}{^}) = 0

, i.e.

A^{⊤} r = 0

. This says the residual

r

is orthogonal to every column of $A$ , hence orthogonal to

col (A)

. Thus

A \overset{x}{^}

is exactly the orthogonal projection of

b

onto

col (A)

— the closest point in the column space to

b

, confirming it minimizes distance. Step 5 (full rank). If

rank A = n

, the columns are independent, so

A^{⊤} A

is symmetric positive definite and invertible, giving the closed form [3]. Step 6 (rank-deficient / minimum norm). If

ker A \neq = {0}

, any

z \in ker A

gives

A (\overset{x}{^} + z) = A \overset{x}{^}

, so

\overset{x}{^} + z

is also a minimizer. Decompose any solution as

\overset{x}{^} = \overset{x}{^}_{row} + \overset{x}{^}_{null}

with

\overset{x}{^}_{row} \in row (A) = ker (A)^{⊥}

. Since these components are orthogonal,

∥ \overset{x}{^} ∥^{2} = ∥ \overset{x}{^}_{row} ∥^{2} + ∥ \overset{x}{^}_{null} ∥^{2}

, minimized by taking

\overset{x}{^}_{null} = 0

. That unique minimum-norm choice lies in the row space and equals

A^{+} b

■

Carl Friedrich Gauss | 1777–1855

Gauss claimed to have used least squares as early as 1795 to fit astronomical and geodetic data, famously predicting the reappearance of the asteroid Ceres in 1801 after it was lost behind the Sun — by fitting an orbit to a handful of noisy observations. He published the method in Theoria Motus Corporum Coelestium (1809) and, crucially, gave it a probabilistic justification: least squares yields the most probable parameters when measurement errors follow what we now call a Gaussian (normal) distribution. This connection between minimizing squared error and maximizing likelihood under normal noise remains the deep reason squared error dominates statistics.

Building Intuition

Think of $col (A)$ as a flat plane (or hyperplane) sitting inside $R^{m}$ . The vector $b$ is a point in space that, generically, floats off the plane — that's why $A x = b$ has no solution. The best you can do is drop a perpendicular from $b$ straight down onto the plane; the foot of that perpendicular is $A \overset{x}{^}$ , and the dropped segment is the residual $r$ , which by construction makes a right angle with the plane. "Least squares" literally means "shortest perpendicular," and the right-angle condition $A^{⊤} r = 0$ is the normal equations.

The pseudoinverse adds a second idea for the underdetermined case. If the plane is too big (many $x$ map to the same point), then once you've landed on the plane there are infinitely many "addresses" $x$ for the same destination — they differ by a null-space direction in which $A$ is blind. The pseudoinverse chooses the address closest to the origin, refusing to add any invisible null-space component. So $A^{+}$ does two clean geometric acts: project (kill the part of $b$ outside the column space) and trim (kill the part of $x$ inside the null space).

Least squares as orthogonal projection: b is projected onto the column space of A, leaving a perpendicular residual.

Worked Examples

Overdetermined: Fitting a Line Through Three Points

Given: Fit

y = c_{0} + c_{1} t

to data

(t, y) = (0, 1), (1, 3), (2, 4)

. Three equations, two unknowns.

A = 111012, b = 134, x = (c_{0} c_{1}) .

Step 1 — form $A^{⊤} A$ and $A^{⊤} b$ :

A^{⊤} A = (3335), A^{⊤} b = (811) .

Step 2 — solve the normal equations

(3335) x = (811)

. Determinant

= 15 - 9 = 6

\overset{x}{^} = \frac{1}{6} (5 - 3 - 3 3) (811) = \frac{1}{6} (40 - 33 - 24 + 33) = \frac{1}{6} (79) .

Step 3 — read off:

c_{0} = 7/6 \approx 1.167

c_{1} = 9/6 = 1.5

. Residual check: fits are

\overset{y}{^} = (1.167, 2.667, 4.167)

; residuals

(- 0.167, 0.333, - 0.167)

sum to zero (as they must, since the column of ones is in

col A

y = 1.167 + 1.5 t (best-fit line)

Underdetermined: Minimum-Norm Solution via Pseudoinverse

Given: One equation, two unknowns:

x_{1} + 2 x_{2} = 5

, so

A = (12)

b = (5)

. Infinitely many solutions lie on a line; we want the smallest-norm point. Step 1 — for a full-row-rank wide matrix,

A^{+} = A^{⊤} (A A^{⊤})^{- 1}

A A^{⊤} = 1^{2} + 2^{2} = 5, (A A^{⊤})^{- 1} = \frac{1}{5} .

Step 2 — compute:

A^{+} = (12) \cdot \frac{1}{5} = (0.2 0.4), \overset{x}{^} = A^{+} b = (0.2 0.4) \cdot 5 = (12) .

Step 3 — verify minimum norm. Check

\overset{x}{^}

satisfies the constraint:

1 + 2 (2) = 5

✓. Its norm is

1 + 4 = 5 \approx 2.236

. Any other solution, e.g.

(5, 0)

, has norm

5

— larger. The solution

(1, 2)

is parallel to the row

(1, 2)

, confirming it lies in the row space (zero null-space component).

\overset{x}{^} = (1, 2), ∥ \overset{x}{^} ∥ = 5 (minimum-norm solution)

Common Misconceptions

Do not invert $A^{⊤} A$ in practice. Forming and solving the normal equations [3] squares the condition number (

κ (A^{⊤} A) = κ (A)^{2}

), amplifying rounding errors badly when

A

is ill-conditioned. Professional solvers use QR decomposition (

A = QR \Rightarrow R \overset{x}{^} = Q^{⊤} b

) or the SVD, which work directly with

κ (A)

. The normal equations are pedagogically clean but numerically fragile. $A^{+}$ is not generally a true inverse. Only when

A

is square and invertible does

A^{+} = A^{- 1}

. In general

A^{+} A \neq = I

and

A A^{+} \neq = I

; instead these are orthogonal projectors onto

row (A)

and

col (A)

respectively. Least squares is not invariant to rescaling variables or weighting equations. Consider doubling the second equation in a system: replacing row

a_{2}^{⊤} x = b_{2}

with

2 a_{2}^{⊤} x = 2 b_{2}

leaves the exact solution set unchanged, but it changes the least-squares answer, because the squared residual of that equation now counts four times as heavily. Concretely, fitting

(0, 0), (1, 0), (1, 10)

to a constant

y = c

gives

\overset{c}{^} = 10/3 \approx 3.33

; if we instead believe the third point is twice as reliable and weight it doubly,

\overset{c}{^}

shifts toward

5

. Least squares optimizes the given metric — choosing the wrong weights silently distorts the fit. (This is why weighted least squares and generalized least squares exist.)

Connections

Least squares sits at a crossroads of linear algebra, geometry, and statistics:

Orthogonal Projection — $A \overset{x}{^} = P b$ where $P = A A^{+}$ is the projector onto $col (A)$ ; least squares is projection.
Singular Value Decomposition — provides the universal formula [6] for $A^{+}$ and reveals how singular values control conditioning and stability.
QR Decomposition — the numerically preferred route to $\overset{x}{^}$ via Gram–Schmidt or Householder reflections.
Linear Regression — ordinary least squares is the statistical estimator; under Gaussian noise it is the maximum-likelihood and (Gauss–Markov) best linear unbiased estimator.
Tikhonov Regularization (ridge regression) — replaces $A^{+}$ with $(A^{⊤} A + λ I)^{- 1} A^{⊤}$ to tame ill-posed problems.
Four Fundamental Subspaces — the pseudoinverse maps $col (A) \to row (A)$ bijectively and kills the orthogonal complements.

Real-World Applications

Orbit determination & GPS — Gauss's original use; modern GPS receivers solve an overdetermined nonlinear least-squares problem (more satellites than the 4 unknowns: position + clock bias) to pin location to meters.
Curve and surface fitting / CAD — fitting polynomials, splines, or planes to measured points in metrology and reverse engineering.
Econometrics & statistics — ordinary least squares regression underlies essentially all empirical economics and the social sciences.
Computer vision — camera calibration, homography estimation, and bundle adjustment minimize reprojection error, often a (nonlinear) least-squares problem.
Machine learning — linear regression, the closed-form solution of ridge regression, and the normal equations seed gradient-based training; pseudoinverse appears in minimum-norm interpolation theory.
Geodesy & seismic tomography — inverting noisy travel-time data for subsurface structure via regularized least squares, where $A^{+}$ tames rank-deficiency.
Signal processing — Wiener filtering, system identification, and deconvolution all reduce to least-squares estimates of model parameters.

Historical Development

1795 — Gauss (age 18) reportedly uses least squares privately for geodetic computations, but does not publish.
1801 — Gauss recovers the lost asteroid Ceres by fitting an orbit to sparse observations, dramatically demonstrating the method's power.
1805 — Adrien-Marie Legendre publishes least squares first, in an appendix to Nouvelles méthodes pour la détermination des orbites des comètes, coining the term méthode des moindres carrés.
1809 — Gauss publishes Theoria Motus, claiming priority and adding the probabilistic (Gaussian-error) justification, igniting a priority dispute with Legendre.
1821–1823 — Gauss proves the Gauss–Markov theorem: least squares gives the best linear unbiased estimator regardless of the error distribution (given finite, equal variances).
1900–1909 — Erik Ivar Fredholm and others develop integral-equation analogues; foundations for generalized inverses laid.
1920 — E. H. Moore defines a "general reciprocal" (the pseudoinverse) for arbitrary matrices, in obscure notation that goes largely unnoticed.
1951 — Arne Bjerhammar independently rediscovers the generalized inverse in a geodesy context.
1955 — Roger Penrose gives the clean four-condition characterization [4] in "A generalized inverse for matrices," fusing the concept with modern matrix theory; the object becomes the Moore–Penrose pseudoinverse.
1965 — Golub & Kahan develop stable SVD algorithms, making pseudoinverse computation routine and numerically sound.

Subtleties & Limits

The pseudoinverse is discontinuous at points where the rank of $A$ changes: a tiny perturbation that lifts a zero singular value to $ε$ makes $A^{+}$ jump by $\sim 1/ ε$ . This makes naive rank determination numerically delicate — how small is "zero"? — and motivates truncated SVD and regularization, which trade a little bias for vastly better stability. Choosing the truncation rank or the regularization parameter $λ$ optimally (cross-validation, L-curve, discrepancy principle) is an active, problem-dependent art with no universal answer.

A deeper subtlety: squared error is optimal only under Gaussian noise. Real data with outliers or heavy tails make least squares non-robust — a single bad point can dominate. This drives ongoing research into robust regression ( $ℓ_{1}$ , Huber loss), sparse recovery ( $ℓ_{1}$ /LASSO, compressed sensing), and the broader question of which loss function the geometry of a given data-generating process actually warrants. The pseudoinverse answers the $ℓ_{2}$ question definitively; for other norms there is no closed-form analogue.