Orthogonality and least squares

Fifty measurements never lie exactly on one line, so the system that would fit a line to them has no solution, and the useful question becomes which line comes closest.

The previous lesson established that $A x = b$ is solvable exactly when $b$ lies in the column space of $A$ . When there are more equations than unknowns the column space is a small subspace of a large space, and a $b$ built from measurements will miss it. This lesson replaces $b$ by the nearest point of the column space, and the entire method comes from one geometric fact: the shortest route from a point to a subspace meets it at a right angle.

Projection onto a line, again

The lesson on vectors projected one vector onto another by requiring the leftover to be perpendicular. Given $b$ and a direction $a$ , the projection is $c a$ where the condition $a \cdot (b - c a) = 0$ forces

c = \frac{a \cdot b}{a \cdot a}

That single equation is the seed of everything below, and it is worth naming what it accomplishes: it finds the point of the line ${c a}$ closest to $b$ . Closest, because any other point $c^{'} a$ on the line has $‖ b - c^{'} a ‖^{2} = ‖ b - c a ‖^{2} + ‖ (c - c^{'}) a ‖^{2}$ by Pythagoras, the two pieces being perpendicular, and the second piece is positive unless $c^{'} = c$ .

Example. Project $b = (4, 0, 0)$ onto the line through $a = (1, 1, 1)$ .

$a \cdot b = 4$ and $a \cdot a = 3$ , so $c = 4 / 3$ and the projection is $(1.3333, 1.3333, 1.3333)$ . The error is $(2.6667, - 1.3333, - 1.3333)$ , whose dot product with $(1, 1, 1)$ is $2.6667 - 1.3333 - 1.3333 = 0$ , confirming the right angle.

Now you. Project $b = (0, 0, 6)$ onto the line through $a = (1, 2, 2)$ .

Answer

$a \cdot b = 12$ and $a \cdot a = 9$ , so $c = 4 / 3 = 1.3333$ and the projection is $(1.3333, 2.6667, 2.6667)$ . The error $(- 1.3333, - 2.6667, 3.3333)$ dotted with $(1, 2, 2)$ gives $- 1.3333 - 5.3333 + 6.6667 = 0$ .

Projection onto a subspace, and the normal equations

Now replace the line by the column space of a matrix $A$ . The problem is: find the vector $p$ in $C (A)$ closest to $b$ . Since every element of the column space is $A \hat{x}$ for some $\hat{x}$ , the unknown is $\hat{x}$ , and the error vector is $e = b - A \hat{x}$ .

The condition is the same as before, generalised: the error must be perpendicular to the whole subspace, which means perpendicular to every column of $A$ . Perpendicularity to each column is a dot product being zero, and the dot products of all the columns with $e$ at once are the components of $A^{T} e$ . So the condition is

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

which rearranges into the normal equations:

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

This is a square system, of size $n$ by $n$ whatever the shape of $A$ , and it is always consistent. When the columns of $A$ are independent, $A^{T} A$ is invertible and the solution is unique. The projection is then $p = A \hat{x}$ , and the quantity being minimised is $‖ b - A x ‖^{2}$ , the sum of the squared errors, which is where the name least squares comes from.

Note what has happened. An unsolvable problem became a solvable one, of smaller size, by a change of question: not "which $x$ satisfies the equations" but "which $x$ misses by the least". The method was published by Legendre in 1805 and used by Gauss to recover the orbit of Ceres from a handful of observations in 1801, after the asteroid was lost behind the sun.

Fitting a straight line

Here is the whole machine on real numbers. A ball rolls along a track and its position is measured at four times: at $t = 0, 1, 2, 3$ seconds the readings are $y = 1, 3, 4, 6$ metres. Fit $y = c + m t$ .

Demanding exactness gives four equations in two unknowns: $c = 1$ , $c + m = 3$ , $c + 2 m = 4$ , $c + 3 m = 6$ . In matrix form $A x = b$ with $A = [(1, 1, 1, 1), (0, 1, 2, 3)]$ , $x = (c, m)$ and $b = (1, 3, 4, 6)$ . The four points are not collinear, so no solution exists: $b$ is a vector in $R^{4}$ and the column space is only a plane.

Form the normal equations. $A^{T} A$ has first row $((1, 1, 1, 1) \cdot (1, 1, 1, 1), (1, 1, 1, 1) \cdot (0, 1, 2, 3)) = (4, 6)$ and second row $(6, 14)$ , since $0 + 1 + 4 + 9 = 14$ . And $A^{T} b = ((1, 1, 1, 1) \cdot b, (0, 1, 2, 3) \cdot b) = (14, 29)$ , since $0 + 3 + 8 + 18 = 29$ . So

4 c + 6 m = 14 6 c + 14 m = 29

The determinant is $56 - 36 = 20$ , so a solution exists and is unique: $m = (4 \times 29 - 6 \times 14) / 20 = 32 / 20 = 1.6$ and $c = (14 - 1.6 \times 6) / 4 = 4.4 / 4 = 1.1$ . The fitted line is $y = 1.1 + 1.6 t$ , so the best estimate of the ball's speed is $1.6$ m s⁻¹.

The fit can be audited. The predicted positions are $1.1$ , $2.7$ , $4.3$ and $5.9$ , so the residuals are $- 0.1$ , $+ 0.3$ , $- 0.3$ and $+ 0.1$ . Their sum is zero, which is the first normal equation, that the error is perpendicular to the column of ones. Their dot product with $(0, 1, 2, 3)$ is $0 + 0.3 - 0.6 + 0.3 = 0$ , which is the second. The sum of squares is $0.01 + 0.09 + 0.09 + 0.01 = 0.20$ , and no other line achieves a smaller value: that is what the perpendicularity guaranteed.

Example. Fit $y = c + m t$ to the three points $(1, 2)$ , $(2, 3)$ and $(3, 5)$ .

The sums are $\sum t = 6$ , $\sum y = 10$ , $\sum t^{2} = 14$ and $\sum t y = 2 + 6 + 15 = 23$ . The normal equations are $3 c + 6 m = 10$ and $6 c + 14 m = 23$ . Their determinant is $42 - 36 = 6$ , so $m = (3 \times 23 - 6 \times 10) / 6 = 9 / 6 = 1.5$ and $c = (10 - 9) / 3 = 0.3333$ . The fitted values are $1.8333$ , $3.3333$ and $4.8333$ , with residuals $+ 0.1667$ , $- 0.3333$ and $+ 0.1667$ , summing to zero as required, and a squared error of $0.1667$ .

Now you. Fit a horizontal line $y = c$ to the four readings $1, 3, 4, 6$ . What is $c$ , and what is it in ordinary language?

Answer

Now $A$ is the single column $(1, 1, 1, 1)$ , so the normal equation is $4 c = 14$ , giving $c = 3.5$ . Fitting a constant by least squares gives the arithmetic mean, which is what the mean is for: the number whose squared deviations from the data are smallest.

The projection matrix

Substituting $\hat{x} = (A^{T} A)^{- 1} A^{T} b$ into $p = A \hat{x}$ gives the projection as a matrix acting on $b$ :

P = A (A^{T} A)^{- 1} A^{T}

Two properties identify it, and both can be read off the geometry rather than computed. $P^{2} = P$ , because projecting a point that is already in the subspace leaves it alone. And $P^{T} = P$ : projection matrices are symmetric, which follows from the formula since $(A^{T} A)^{- 1}$ is symmetric and the outer factors transpose into each other. Symmetry is not obvious geometrically, and the final lesson explains why projections belong to the best-behaved family of matrices.

Two special cases check the formula. If $A$ is square and invertible, its column space is everything, nothing needs projecting, and indeed $P = A (A^{- 1} (A^{T})^{- 1}) A^{T} = I$ . If $A$ is a single column $a$ , then $A^{T} A$ is the number $a \cdot a$ and $P b = a (a \cdot b) / (a \cdot a)$ , the line formula from the start of the lesson.

Orthonormal bases

Projection got complicated because the columns of $A$ were not perpendicular. When they are, everything simplifies at once.

A set of vectors is orthonormal when each has length one and any two distinct ones are perpendicular. If $q_{1}, \dots, q_{k}$ are orthonormal then the coordinates of any vector in their span are simply dot products: writing $v = c_{1} q_{1} + \dots + c_{k} q_{k}$ and dotting both sides with $q_{i}$ kills every term but one, leaving $c_{i} = v \cdot q_{i}$ . No system needs solving.

Take $q_{1} = (0.6, 0.8)$ and $q_{2} = (- 0.8, 0.6)$ , which are unit vectors, since $0.36 + 0.64 = 1$ , and perpendicular, since $- 0.48 + 0.48 = 0$ . The coordinates of $(2, 3)$ in this basis are $(2, 3) \cdot (0.6, 0.8) = 1.2 + 2.4 = 3.6$ and $(2, 3) \cdot (- 0.8, 0.6) = - 1.6 + 1.8 = 0.2$ . Reconstructing: $3.6 (0.6, 0.8) + 0.2 (- 0.8, 0.6) = (2.16 - 0.16, 2.88 + 0.12) = (2, 3)$ .

A square matrix $Q$ whose columns are orthonormal is called an orthogonal matrix, and it satisfies $Q^{T} Q = I$ , since the entry in row $i$ and column $j$ of $Q^{T} Q$ is $q_{i} \cdot q_{j}$ . So its inverse is its transpose, which costs nothing to compute. Such matrices preserve every length and every angle, because $(Q u) \cdot (Q v) = u^{T} Q^{T} Q v = u \cdot v$ , and in the plane they are exactly the rotations, of determinant $+ 1$ , and the reflections, of determinant $- 1$ . Their determinant is always $\pm 1$ , since $\det (Q)^{2} = \det (Q^{T} Q) = 1$ .

Any independent set can be turned into an orthonormal one by Gram-Schmidt: keep the first vector, normalised; from the second, subtract its projection onto the first and normalise what remains; from the third, subtract its projections onto both and normalise; and so on. Starting from $(1, 1, 0)$ and $(1, 0, 1)$ : the first normalises to $q_{1} = (0.7071, 0.7071, 0)$ . The second has $a_{2} \cdot q_{1} = 0.7071$ , so subtracting $0.7071 q_{1} = (0.5, 0.5, 0)$ leaves $(0.5, - 0.5, 1)$ , of length $\sqrt{1.5} = 1.2247$ , giving $q_{2} = (0.4082, - 0.4082, 0.8165)$ . Their dot product is $0.2887 - 0.2887 + 0 = 0$ .

Example. Verify that $Q = [(0.6, 0.8), (- 0.8, 0.6)]$ is orthogonal, and find the coordinates of $(5, 0)$ in its columns.

Each column has length one and the two are perpendicular, as computed above, so $Q^{T} Q = I$ . The coordinates are dot products: $(5, 0) \cdot (0.6, 0.8) = 3.0$ and $(5, 0) \cdot (- 0.8, 0.6) = - 4.0$ . Checking: $3 (0.6, 0.8) - 4 (- 0.8, 0.6) = (1.8 + 3.2, 2.4 - 2.4) = (5, 0)$ . Note that $3^{2} + (- 4)^{2} = 25 = 5^{2}$ , so the length is the same computed in either basis, which is what preserving length means.

Now you. What is the inverse of $Q = [(0.6, 0.8), (- 0.8, 0.6)]$ , and what transformation is it?

Answer

The inverse is the transpose, $Q^{T} = [(0.6, - 0.8), (0.8, 0.6)]$ . Since $\det Q = 0.36 + 0.64 = 1$ , $Q$ is a rotation, by the angle whose cosine is $0.6$ and sine is $0.8$ , about $53.13 °$ . Its inverse is the rotation by $- 53.13 °$ , which is what the transpose describes.

Honest limits

The normal equations are the right derivation and the wrong algorithm. Forming $A^{T} A$ squares the condition number of $A$ , the ratio of its largest stretching factor to its smallest, and squaring a condition number of $10^{6}$ gives $10^{12}$ , which exceeds what double precision can absorb. Numerical libraries fit by QR factorisation, in which $A$ is written as an orthogonal matrix times a triangular one, or by the singular value decomposition, and they never form $A^{T} A$ . The geometry above is exactly right; the arithmetic route to it is not.

Least squares is also a choice, not a law. Squaring the residuals means one badly wrong measurement counts as much as many slightly wrong ones, so a single outlier can drag the whole line. It is the correct choice when errors are independent and normally distributed with equal variance, which Gauss proved and which is sometimes even true. When it is not, other measures of misfit are used, at the price of losing the linear algebra: minimising the sum of absolute errors has no formula like the normal equations.

Finally, a fitted line is a description of the data, not a licence to extrapolate. The ball's fitted position at $t = 0$ is $1.1$ metres while the measurement was $1.0$ , and nothing in the method says which is closer to the truth.

Coordinates in an orthonormal basis turned out to be free, which raises a question left open two lessons ago. If a basis can be chosen, what does a change of basis do to the matrix of a transformation, and how much can a good choice simplify it? That is the next lesson, and it leads directly to the last third of the course.