6.5 Least-squares problems

Big idea

An inconsistent system has no solution, and until now that ended the conversation. Real data makes inconsistent systems constantly: more measurements than unknowns, each with a little noise, and no vector satisfies all of them at once. Refusing to answer is not useful.

So change the question. Instead of demanding $A\mathbf{x} = \mathbf{b}$, ask for the $\mathbf{x}$ that makes $A\mathbf{x}$ as close to $\mathbf{b}$ as possible. Every $A\mathbf{x}$ lies in the column space of $A$, so the closest reachable vector is the orthogonal projection of $\mathbf{b}$ onto that column space, and the answer is any $\mathbf{x}$ that produces it.

That reduces a new problem to the previous section, but there is a better route than projecting first. The condition “the error is orthogonal to every column” can be written as one matrix equation in $\mathbf{x}$ alone, and solving it needs no orthogonal basis at all. That equation is the practical content of the section.

Decoder

A least-squares solution of an inconsistent system is a vector minimizing the Euclidean norm of the residual, and it need not be unique.

“Least squares” names the quantity being minimized: the sum of the squares of the entries of $\mathbf{b} - A\mathbf{x}$, which is the squared length of that vector. Minimizing the length and minimizing the sum of squares are the same problem, and squares are easier to differentiate and to write as a dot product.

The uniqueness warning is about $\mathbf{x}$, not about the projection. The closest point $\hat{\mathbf{b}}$ in the column space is always unique. If the columns of $A$ are dependent, several different $\mathbf{x}$ produce that same point, and all of them are least-squares solutions.

Definitions and results

Least-squares solution. For $A$ of size $m \times n$ and $\mathbf{b}$ in $\mathbb{R}^m$, a vector $\hat{\mathbf{x}}$ in $\mathbb{R}^n$ is a least-squares solution when $|\mathbf{b} - A\hat{\mathbf{x}}| \leq |\mathbf{b} - A\mathbf{x}|$ for every $\mathbf{x}$ in $\mathbb{R}^n$.

Reduction to projection. Let $\hat{\mathbf{b}}$ be the orthogonal projection of $\mathbf{b}$ onto $\text{Col}\,A$. Since $\hat{\mathbf{b}}$ is in the column space, the system $A\mathbf{x} = \hat{\mathbf{b}}$ is consistent, and its solutions are exactly the least-squares solutions of the original system.

The normal equations. $\hat{\mathbf{x}}$ is a least-squares solution exactly when

$$ A^TA\hat{\mathbf{x}} = A^T\mathbf{b} $$

This comes from the orthogonality condition directly. Requiring $\mathbf{b} - A\hat{\mathbf{x}}$ to be orthogonal to every column of $A$ says $A^T(\mathbf{b} - A\hat{\mathbf{x}}) = \mathbf{0}$, which rearranges into the display above. The normal equations are always consistent, whatever $A$ and $\mathbf{b}$ are, so a least-squares solution always exists.

When it is unique. The following conditions say the same thing: the columns of $A$ are linearly independent; $A^TA$ is invertible; the least-squares solution is unique. In that case

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

although in practice you solve the square system rather than form the inverse.

Residual. The vector $\mathbf{b} - A\hat{\mathbf{x}}$ is the residual, and its length is the least-squares error. It equals zero exactly when the original system was consistent to begin with. It is always orthogonal to every column of $A$, which is the cheapest available check on your arithmetic.

Projection matrix. When the columns are independent, $\hat{\mathbf{b}} = A(A^TA)^{-1}A^T\mathbf{b}$, so that product of three matrices projects onto the column space of $A$. It depends only on $A$.

The QR route. If $A = QR$ with independent columns, substituting into the normal equations and cancelling the invertible factors leaves the triangular system

$$ R\hat{\mathbf{x}} = Q^T\mathbf{b} $$

solved by back substitution. This is the numerically preferred method, because forming $A^TA$ squares the spread between the largest and smallest scales in the data and can destroy accuracy.

Worked examples

An inconsistent system, solved anyway

Take

$$ A = \begin{bmatrix} 1 & 0 \\ 0 & 1 \\ 1 & 1 \end{bmatrix}, \qquad \mathbf{b} = \begin{bmatrix} 4 \\ 3 \\ 1 \end{bmatrix} $$

The system is inconsistent: the first two rows force $x_1 = 4$ and $x_2 = 3$, but the third then demands $7 = 1$.

Form the normal equations. $A^TA$ is the $2 \times 2$ matrix of dot products of columns:

$$ A^TA = \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}, \qquad A^T\mathbf{b} = \begin{bmatrix} 4 + 1 \\ 3 + 1 \end{bmatrix} = \begin{bmatrix} 5 \\ 4 \end{bmatrix} $$

Solve $2x_1 + x_2 = 5$ and $x_1 + 2x_2 = 4$. Doubling the first and subtracting the second gives $3x_1 = 6$, so $x_1 = 2$ and then $x_2 = 1$. The columns of $A$ are independent, so this is the only least-squares solution.

The projection is $A\hat{\mathbf{x}} = (2, 1, 3)$ and the residual is

$$ \mathbf{b} - A\hat{\mathbf{x}} = \begin{bmatrix} 4 \\ 3 \\ 1 \end{bmatrix} - \begin{bmatrix} 2 \\ 1 \\ 3 \end{bmatrix} = \begin{bmatrix} 2 \\ 2 \\ -2 \end{bmatrix} $$

Check it against both columns: $(2,2,-2) \cdot (1,0,1) = 2 - 2 = 0$ and $(2,2,-2) \cdot (0,1,1) = 2 - 2 = 0$. The least-squares error is $\sqrt{4 + 4 + 4} = 2\sqrt{3}$, about $3.46$.

The same answer two ways

Take

$$ A = \begin{bmatrix} 1 & 2 \\ 1 & 0 \\ 0 & 1 \end{bmatrix}, \qquad \mathbf{b} = \begin{bmatrix} 4 \\ 0 \\ -1 \end{bmatrix} $$

Normal equations first:

$$ A^TA = \begin{bmatrix} 2 & 2 \\ 2 & 5 \end{bmatrix}, \qquad A^T\mathbf{b} = \begin{bmatrix} 4 + 0 \\ 8 + 0 - 1 \end{bmatrix} = \begin{bmatrix} 4 \\ 7 \end{bmatrix} $$

Subtracting the first equation from the second gives $3x_2 = 3$, so $x_2 = 1$, and then $2x_1 + 2 = 4$ gives $x_1 = 1$. So $\hat{\mathbf{x}} = (1, 1)$, $\hat{\mathbf{b}} = A\hat{\mathbf{x}} = (3, 1, 1)$, and the residual is $(1, -1, -2)$ with length $\sqrt{6}$. Check: $(1,-1,-2) \cdot (1,1,0) = 0$ and $(1,-1,-2) \cdot (2,0,1) = 2 - 2 = 0$.

Now the QR route, reusing the factorization built for this same matrix earlier in the chapter:

$$ Q = \begin{bmatrix} 1/\sqrt{2} & 1/\sqrt{3} \\ 1/\sqrt{2} & -1/\sqrt{3} \\ 0 & 1/\sqrt{3} \end{bmatrix}, \qquad R = \begin{bmatrix} \sqrt{2} & \sqrt{2} \\ 0 & \sqrt{3} \end{bmatrix} $$

Compute $Q^T\mathbf{b}$: the first entry is $(4 + 0)/\sqrt{2} = 2\sqrt{2}$ and the second is $(4 - 0 - 1)/\sqrt{3} = \sqrt{3}$. Back substitution in $R\hat{\mathbf{x}} = Q^T\mathbf{b}$ gives $\sqrt{3}\,x_2 = \sqrt{3}$, so $x_2 = 1$, and then $\sqrt{2}\,x_1 + \sqrt{2} = 2\sqrt{2}$, so $x_1 = 1$. Same answer, no $A^TA$ formed.

Reading the residual as a quality score

In the first example the error was $2\sqrt{3}$, about $3.46$, against a data vector of length $\sqrt{16 + 9 + 1} = \sqrt{26}$, about $5.10$. Most of the data vector lies outside the column space, so the model $A\mathbf{x}$ explains it poorly. In the second example the error was $\sqrt{6}$, about $2.45$, against $\|\mathbf{b}\| = \sqrt{17}$, about $4.12$; better, but the plane still misses.

Neither number is meaningful on its own. The residual is only interpretable next to the scale of the data, which is why the comparison above divides one by the other rather than quoting the error alone.

Practice

The projection half: find the point of a column space closest to a given vector, and the distance to it.

The algebra half: the normal equations form a square system with a unique solution whenever the columns are independent, so solving them is ordinary work.

Videos

Watch for the derivation of the normal equations from the orthogonality of the error, and for the remark on what goes wrong when the columns are dependent.

Quiz

Five items on setting up and solving normal equations and on reading the residual.