Least Squares and the Pseudoinverse

⚡ Key Takeaways

Least squares finds the vector $\hat{x}$ minimizing $\|Ax-b\|^2$ when $Ax=b$ has no exact solution. The unique minimum-norm solution is $\hat{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).

💡 The Core Insight

You cannot always solve $Ax=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.

DefinitionThe Least Squares Problem

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

1

\quad \hat{x}=\arg\min_{x\in\mathbb{R}^n}\;\|Ax-b\|_2^{\,2}

Derivation hint: This is the definition; $\|v\|_2^2=\sum_i v_i^2=v^\top v$ is the squared length of the residual vector $v=Ax-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=Ax-b$ is the residual; $\|\cdot\|_2$ is the Euclidean norm.

DefinitionThe Normal Equations

Any minimizer $\hat{x}$ of $\|Ax-b\|^2$ satisfies the normal equations:

2

\quad A^\top A\,\hat{x}=A^\top b

Derivation hint: Set the gradient $\nabla_x\|Ax-b\|^2=2A^\top(Ax-b)$ equal to zero. Here $A^\top A\in\mathbb{R}^{n\times n}$ is symmetric and positive semidefinite. If $A$ has full column rank ( $\operatorname{rank}A=n$ ), then $A^\top A$ is invertible and the solution is unique:

3

\quad \hat{x}=(A^\top A)^{-1}A^\top b

Derivation hint: Left-multiply [2] by $(A^\top A)^{-1}$ , which exists exactly when $A$ 's columns are linearly independent.

DefinitionThe Moore–Penrose Pseudoinverse

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

4

\quad AA^{+}A=A,\quad A^{+}AA^{+}=A^{+},\quad (AA^{+})^\top=AA^{+},\quad (A^{+}A)^\top=A^{+}A

Derivation hint: These four equations uniquely determine $A^+$ ; the first two are "weak inverse" axioms, the last two force the products $AA^+$ and $A^+A$ to be orthogonal projectors (symmetric). The general least-squares / minimum-norm solution is then

5

\quad \hat{x}=A^{+}b

Derivation hint: $A^+$ projects $b$ onto $\operatorname{col}(A)$ and inverts on the row space, automatically selecting the smallest-norm minimizer.

TheoremExistence, Uniqueness, and Minimum Norm

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

TheoremPseudoinverse via the SVD

If $A=U\Sigma V^\top$ is a singular value decomposition with nonzero singular values $\sigma_1,\dots,\sigma_r$ , then

6

\quad A^{+}=V\,\Sigma^{+}U^\top,\qquad \Sigma^{+}=\operatorname{diag}\!\left(\tfrac{1}{\sigma_1},\dots,\tfrac{1}{\sigma_r},0,\dots,0\right)

Derivation hint: Substitute $A=U\Sigma V^\top$ into the four Penrose conditions [4] and verify, using $U^\top U=I$ , $V^\top V=I$ , and that $\Sigma^+$ inverts the nonzero block while zeroing the rest.

ProofThe Normal Equations from Orthogonality

Step 1. Define $f(x)=\|Ax-b\|^2=(Ax-b)^\top(Ax-b)=x^\top A^\top A x-2b^\top A x+b^\top b$ .

Step 2. Differentiate with respect to $x$ :

\nabla f(x)=2A^\top A x-2A^\top b.

Since $f$ is convex (its Hessian $2A^\top A\succeq 0$ ), any stationary point is a global minimum.

Step 3. Set $\nabla f=0$ : $A^\top A x=A^\top b$ , which is equation [2].

Step 4 (geometric meaning). Rewrite as $A^\top(b-A\hat{x})=0$ , i.e. $A^\top r=0$ . This says the residual $r$ is orthogonal to every column of $A$ , hence orthogonal to $\operatorname{col}(A)$ . Thus $A\hat{x}$ is exactly the orthogonal projection of $b$ onto $\operatorname{col}(A)$ — the closest point in the column space to $b$ , confirming it minimizes distance.

Step 5 (full rank). If $\operatorname{rank}A=n$ , the columns are independent, so $A^\top 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(\hat{x}+z)=A\hat{x}$ , so $\hat{x}+z$ is also a minimizer. Decompose any solution as $\hat{x}=\hat{x}_{\text{row}}+\hat{x}_{\text{null}}$ with $\hat{x}_{\text{row}}\in\operatorname{row}(A)=\ker(A)^\perp$ . Since these components are orthogonal, $\|\hat{x}\|^2=\|\hat{x}_{\text{row}}\|^2+\|\hat{x}_{\text{null}}\|^2$ , minimized by taking $\hat{x}_{\text{null}}=0$ . That unique minimum-norm choice lies in the row space and equals $A^{+}b$ . $\blacksquare$

∎

CF

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.

🧠 Intuition

Think of $\operatorname{col}(A)$ as a flat plane (or hyperplane) sitting inside $\mathbb{R}^m$ . The vector $b$ is a point in space that, generically, floats off the plane — that's why $Ax=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\hat{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^\top 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).

Figure

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

ExampleOverdetermined: Fitting a Line Through Three Points

Given: Fit $y=c_0+c_1t$ to data $(t,y)=(0,1),(1,3),(2,4)$ . Three equations, two unknowns.

A=\begin{pmatrix}1&0\\1&1\\1&2\end{pmatrix},\quad b=\begin{pmatrix}1\\3\\4\end{pmatrix},\quad x=\begin{pmatrix}c_0\\c_1\end{pmatrix}.

Step 1 — form $A^\top A$ and $A^\top b$ :

A^\top A=\begin{pmatrix}3&3\\3&5\end{pmatrix},\qquad A^\top b=\begin{pmatrix}8\\11\end{pmatrix}.

Step 2 — solve the normal equations $\begin{pmatrix}3&3\\3&5\end{pmatrix}x=\begin{pmatrix}8\\11\end{pmatrix}$ . Determinant $=15-9=6$ .

\hat{x}=\frac{1}{6}\begin{pmatrix}5&-3\\-3&3\end{pmatrix}\begin{pmatrix}8\\11\end{pmatrix}=\frac{1}{6}\begin{pmatrix}40-33\\-24+33\end{pmatrix}=\frac{1}{6}\begin{pmatrix}7\\9\end{pmatrix}.

Step 3 — read off: $c_0=7/6\approx1.167$ , $c_1=9/6=1.5$ .

Residual check: fits are $\hat 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 $\operatorname{col}A$ ).

\boxed{\;y=1.167+1.5\,t\quad(\text{best-fit line})\;}

ExampleUnderdetermined: Minimum-Norm Solution via Pseudoinverse

Given: One equation, two unknowns: $x_1+2x_2=5$ , so $A=(1\;\;2)$ , $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^\top(AA^\top)^{-1}$ .

AA^\top=1^2+2^2=5,\qquad (AA^\top)^{-1}=\tfrac15.

Step 2 — compute:

A^{+}=\begin{pmatrix}1\\2\end{pmatrix}\cdot\tfrac15=\begin{pmatrix}0.2\\0.4\end{pmatrix},\qquad \hat{x}=A^{+}b=\begin{pmatrix}0.2\\0.4\end{pmatrix}\cdot5=\begin{pmatrix}1\\2\end{pmatrix}.

Step 3 — verify minimum norm. Check $\hat x$ satisfies the constraint: $1+2(2)=5$ ✓. Its norm is $\sqrt{1+4}=\sqrt5\approx2.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).

\boxed{\;\hat{x}=(1,\,2),\quad \|\hat{x}\|=\sqrt5\;(\text{minimum-norm solution})\;}

⚠️ Common mistake

Do not invert $A^\top A$ in practice. Forming and solving the normal equations [3] squares the condition number ( $\kappa(A^\top A)=\kappa(A)^2$ ), amplifying rounding errors badly when $A$ is ill-conditioned. Professional solvers use QR decomposition ( $A=QR\Rightarrow R\hat x=Q^\top b$ ) or the SVD, which work directly with $\kappa(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 $AA^{+}\neq I$ ; instead these are orthogonal projectors onto $\operatorname{row}(A)$ and $\operatorname{col}(A)$ respectively.

Counterexample

Least squares is not invariant to rescaling variables or weighting equations. Consider doubling the second equation in a system: replacing row $a_2^\top x=b_2$ with $2a_2^\top x=2b_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 $\hat c=10/3\approx3.33$ ; if we instead believe the third point is twice as reliable and weight it doubly, $\hat 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.)

connection

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

Orthogonal Projection — $A\hat x=P\,b$ where $P=AA^{+}$ is the projector onto $\operatorname{col}(A)$ ; least squares is projection.
<a class="wiki-link" href="/wiki/cs-ml-singular-value-decomposition">Singular Value Decomposition</a> — provides the universal formula [6] for $A^{+}$ and reveals how singular values control conditioning and stability.
<a class="wiki-link" href="/wiki/qr-decomposition">QR Decomposition</a> — the numerically preferred route to $\hat x$ via Gram–Schmidt or Householder reflections.
<a class="wiki-link" href="/wiki/linear-regression">Linear Regression</a> — 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^\top A+\lambda I)^{-1}A^\top$ to tame ill-posed problems.
Four Fundamental Subspaces — the pseudoinverse maps $\operatorname{col}(A)\to\operatorname{row}(A)$ bijectively and kills the orthogonal complements.

🔬 Where this appears

⏳ History

- 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.

open-problem

The pseudoinverse is discontinuous at points where the rank of $A$ changes: a tiny perturbation that lifts a zero singular value to $\varepsilon$ makes $A^{+}$ jump by $\sim 1/\varepsilon$ . 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 $\lambda$ 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 ( $\ell_1$ , Huber loss), sparse recovery ( $\ell_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 $\ell_2$ question definitively; for other norms there is no closed-form analogue.

Formal Statement

Why It Holds / Derivation

Building Intuition

Worked Examples

Common Misconceptions

Connections

Real-World Applications

Historical Development

Subtleties & Limits

See Also