Solving Normal Equations

What the System Solves

Given $A \in \mathbb{R}^{m \times n}$ and $b \in \mathbb{R}^m$ , a linear least-squares problem asks for the parameter vector that makes the residual as small as possible:

\min_x \|Ax-b\|_2^2

Every minimizer satisfies the normal equation

A^T A\hat{x}=A^T b

This equation is $always consistent$ . The minimizer is unique exactly when $\operatorname{rank}(A)=n$ , which requires $m\ge n$ . In that full-column-rank case,

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

When $A$ is rank deficient, the fitted vector $A\hat{x}$ is still unique, but the parameter vector is not. If $\hat{x}$ is one solution, then every solution is in

\hat{x}+\ker(A)

The column space of $A$ is a finite-dimensional closed subspace of $\mathbb{R}^m$ . Therefore, $b$ has a $unique orthogonal projection$ onto it. Any $\hat{x}$ satisfying $A\hat{x}=\operatorname{proj}_{\operatorname{col}(A)}b$ is a least-squares solution.

The normal equation expresses the same geometry because

A^T(A\hat{x}-b)=0

This equation says that the residual $A\hat{x}-b$ is orthogonal to every column of $A$ .

Moore Penrose Pseudoinverse

The Moore-Penrose pseudoinverse is defined for every real matrix, including rectangular and rank-deficient matrices. Write its singular value decomposition as

A=U\Sigma V^T

The displayed factorization is a singular value decomposition. To form the pseudoinverse, replace every nonzero singular value $\sigma_i$ with $1/\sigma_i$ and leave every zero singular value at zero. This gives

\Sigma^{\dagger}_{ii}=\begin{cases}1/\sigma_i,&\sigma_i>0,\\0,&\sigma_i=0,\end{cases}

The vector $x^{\dagger}$ is $the unique least-squares solution with minimum Euclidean norm$ . When $A$ has full column rank, the general definition simplifies to

A^{\dagger}=(A^T A)^{-1}A^T

The pseudoinverse is uniquely characterized by the four Penrose identities:

AA^{\dagger}A=A

Direct Solution with QR Decomposition

For a full-column-rank matrix with $m\ge n$ , use the thin $\mathrm{QR}$ decomposition

A=Q_1R_1,\qquad Q_1^TQ_1=I

The least-squares objective separates into

\begin{aligned} \|Ax-b\|_2^2 &=\|R_1x-Q_1^Tb\|_2^2 \\ &\quad{}+\|(I-Q_1Q_1^T)b\|_2^2 \end{aligned}

The second term does not depend on $x$ . Therefore, the minimizer is obtained by solving the $upper-triangular system$

R_1\hat{x}=Q_1^Tb

This route solves the original least-squares problem directly. It does not form $A^TA$ and then solve the normal equation.

Worked Polynomial Fit Example

The numerical factors and coefficients below are $rounded to five decimal places$ . The exact matrix equalities hold before rounding.

Fit $y=at^2+bt+c$ to seven observations. The design matrix and observation vector are

A=\begin{pmatrix}9&-3&1\\4&-2&1\\1&-1&1\\0&0&1\\1&1&1\\4&2&1\\9&3&1\end{pmatrix}

Each row of $A$ is $[t_i^2,t_i,1]$ . One valid thin $\mathrm{QR}$ decomposition is

Q_1=\begin{pmatrix}-0.64286&-0.56695&0.16496\\-0.28571&-0.37796&-0.24744\\-0.07143&-0.18898&-0.49487\\0&0&-0.57735\\-0.07143&0.18898&-0.49487\\-0.28571&0.37796&-0.24744\\-0.64286&0.56695&0.16496\end{pmatrix}

Signs of matching columns of $Q_1$ and rows of $R_1$ may change between valid $\mathrm{QR}$ implementations. The product $Q_1R_1$ and final fit do not change.

First compute

Q_1^Tb=\begin{pmatrix}-9.71429\\16.02569\\3.48060\end{pmatrix}

Then use $back substitution$ :

\begin{aligned} &-1.73205x_3=3.48060 \\ &\Longrightarrow \\ &x_3=-2.00952 \end{aligned}

The calculation gives:

\hat{x}=\begin{pmatrix}0.98095\\3.02857\\-2.00952\end{pmatrix}

Least-squares fit

\hat{y}(t)=0.98095t^2+3.02857t-2.00952

The cyan curve is the fitted polynomial, while the orange points are the seven observations in

b

$x$: $[-3, 3]$
$y$: $[-4.35, 15.9]$

Choosing a Numerical Method

Method	Approximate dense cost for $m\ge n$	Main strength	Main limitation
Normal equations and Cholesky	form $A^TA$ in $O(mn^2)$ , then factor in about $\frac{1}{3}n^3$ flops	fast when the problem is well conditioned and full rank	$\kappa_2(A^TA)=\kappa_2(A)^2$
Householder $\mathrm{QR}$	about $2mn^2-\frac{2}{3}n^3$ flops	avoids forming $A^TA$ and is usually more stable	requires a rank-aware variant when columns are dependent
$\mathrm{SVD}$	more work than $\mathrm{QR}$	determines rank and gives the $minimum-norm solution$	highest cost of the three

Use $\mathrm{QR}$ as the usual choice for dense, full-rank problems. Use $\mathrm{SVD}$ when the matrix is rank deficient, severely ill conditioned, or requires a minimum-norm solution. Form the normal matrix only when the speed gain is worth the loss of numerical stability. Do not compute an inverse explicitly in production code.