Orthogonality, least squares and SVD

The last chapter left us facing a cloud of points through which no line can be drawn. Let’s take a smaller cloud, of three points. Weights of $1$, $2$ and $3$ kg are hung on a spring, and the pointer on the scale stops at the marks $1$, $2$ and $2$. By Hooke’s law the reading grows linearly: $y = C + Dt$, where $t$ is the mass of the weight, $C$ is where the pointer would stand with no weight, and $D$ is how many marks each kilogram moves it. Three measurements give three equations in two unknowns:

$$C + D = 1, \qquad C + 2D = 2, \qquad C + 3D = 2.$$

The first two give $D = 1$ and $C = 0$, and then the third demands $3 = 2$. There’s no line through the three points. But the spring does have a stiffness; it’s just that a hand shook or the scale was marked unevenly. The problem changes: find the line that errs least of all. What “least of all” means will have to be agreed on, and an agreement two centuries old turned out to be more successful than anyone could have expected. The method of least squares will grow out of it, and out of that the decomposition of any matrix into a rotation, a stretch and another rotation.

The chapter is laid out as a lab notebook. It has six experiments with data, and between the experiments, the tools without which they couldn’t be carried out.

Experiment 1. Three measurements and not a single line

Write the system as a matrix, as in the chapter on Gaussian elimination:

$$\underbrace{\begin{pmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \end{pmatrix}}_{A}\underbrace{\begin{pmatrix} C \\ D \end{pmatrix}}_{\mathbf x} = \underbrace{\begin{pmatrix} 1 \\ 2 \\ 2 \end{pmatrix}}_{\mathbf b}.$$

Now look at it column by column. The product $A\mathbf x$ is $C\,\mathbf a_1 + D\,\mathbf a_2$, where $\mathbf a_1 = (1, 1, 1)$ and $\mathbf a_2 = (1, 2, 3)$ are the columns of $A$. All such combinations fill a plane in three-dimensional space passing through the origin: it’s the span of the columns, or the image of the matrix (chapter 37). Solving the system means reaching the point $\mathbf b = (1, 2, 2)$ without leaving this plane. And $\mathbf b$ lies outside it, so there’s no solution.

Since $\mathbf b$ can’t be reached, let’s reach the point of the plane nearest to it. Distance in space is the square root of the sum of the squares of the differences of coordinates, like the diagonal of a box from the chapter on Pythagoras’ theorem. And the coordinates of the difference $\mathbf b - A\mathbf x$ are exactly the misses of the line at the three measurements: $y_i - (C + Dt_i)$.

So “the nearest point of the plane” and “the line with the smallest sum of squared misses” are one and the same problem, stated in two languages. We know geometry better, so we’ll solve it in that language. True, there are usually hundreds of measurements rather than three, and the vector $\mathbf b$ lives in the space $\mathbb R^{100}$ or $\mathbb R^{1000}$. It can’t be drawn, so we need a ruler that doesn’t care how many dimensions the space has.

A ruler for any number of dimensions

In the chapter on vectors the dot product in the plane was computed from coordinates: $a_1b_1 + a_2b_2$. For vectors in $\mathbb R^n$ the formula is the same, only with $n$ terms:

$$\langle \mathbf u, \mathbf v\rangle = u_1v_1 + u_2v_2 + \dots + u_nv_n = \mathbf u^{\mathsf T}\mathbf v.$$

Through it are defined the length $|\mathbf u| = \sqrt{\langle \mathbf u, \mathbf u\rangle}$ (Pythagoras’ theorem applied $n - 1$ times) and the distance $|\mathbf u - \mathbf v|$. In the plane the dot product also equalled $|\mathbf u|\,|\mathbf v|\cos\varphi$. In $\mathbb R^n$ there’s no protractor, and the angle is defined precisely by this equation.

The angle between non-zero vectors, from $0^\circ$ to $180^\circ$. The dot product: the sum of pairwise products of the coordinates. The product of the lengths. Dividing by it removes the scale: double the length of either vector, and the angle doesn’t change. Example: for the columns $\mathbf a_1 = (1, 1, 1)$ and $\mathbf a_2 = (1, 2, 3)$ we get $\langle \mathbf a_1, \mathbf a_2\rangle = 6$, $|\mathbf a_1| = \sqrt3$, $|\mathbf a_2| = \sqrt{14}$ and $\cos\theta = \frac{6}{\sqrt{42}} \approx 0.926$, that is, $\theta \approx 22.2^\circ$. The columns of our matrix point almost the same way.

The definition is honest only if the fraction on the right never goes outside the interval $[-1, 1]$; otherwise there’s simply no such angle. That’s what the main inequality of linear algebra asserts.

The trick is that two vectors always lie in one plane, even in $\mathbb R^{1000}$. Drop a perpendicular from the tip of $\mathbf u$ to the line of the vector $\mathbf v$: the shadow can’t be longer than the vector itself. All the steps below are calculations with the dot product, so the picture suggests nothing beyond the calculations.

If $\mathbf v = \mathbf 0$, both sides are zero. Otherwise consider the vector $\mathbf p = \frac{\langle \mathbf u, \mathbf v\rangle}{\langle \mathbf v, \mathbf v\rangle}\,\mathbf v$ on the line of the vector $\mathbf v$. Its length is $\frac{|\langle \mathbf u, \mathbf v\rangle|}{|\mathbf v|^2}\,|\mathbf v| = \frac{|\langle \mathbf u, \mathbf v\rangle|}{|\mathbf v|}$. The remainder $\mathbf e = \mathbf u - \mathbf p$ is perpendicular to $\mathbf v$. The dot product is linear in each argument, so $\langle \mathbf e, \mathbf v\rangle = \langle \mathbf u, \mathbf v\rangle - \frac{\langle \mathbf u, \mathbf v\rangle}{\langle \mathbf v, \mathbf v\rangle}\langle \mathbf v, \mathbf v\rangle = 0$. The vector $\mathbf p$ is a multiple of $\mathbf v$, so $\langle \mathbf p, \mathbf e\rangle = 0$ as well. Pythagoras’ theorem in $\mathbb R^n$: expand $|\mathbf u|^2 = \langle \mathbf p + \mathbf e, \mathbf p + \mathbf e\rangle = |\mathbf p|^2 + 2\langle \mathbf p, \mathbf e\rangle + |\mathbf e|^2$. The middle term is zero by the previous step, leaving $|\mathbf u|^2 = |\mathbf p|^2 + |\mathbf e|^2$. The square $|\mathbf e|^2$ isn’t negative, so $|\mathbf p| \le |\mathbf u|$, that is, $\frac{|\langle \mathbf u, \mathbf v\rangle|}{|\mathbf v|} \le |\mathbf u|$. Multiplying by $|\mathbf v|$ gives the inequality. Equality means $\mathbf e = \mathbf 0$: the vector $\mathbf u$ coincides with its shadow $\mathbf p$ and is a multiple of $\mathbf v$.

For sums this inequality was published by Augustin-Louis Cauchy in 1821, and for integrals by Viktor Bunyakovsky in 1859 and independently by Hermann Schwarz in 1888; in English it’s usually called the Cauchy–Schwarz inequality, while some books call it the Cauchy–Bunyakovsky inequality. The integrals aren’t there by chance: in the chapter on Fourier series the dot product of functions was $\int f g\,dx$, and everything we prove about vectors holds for signals too.

The most useful angle is the right angle.

For which $a$ are the vectors $(1, 2, a, 1)$ and $(3, -1, 2, 1)$ in $\mathbb R^4$ orthogonal?

The dot product is $1 \cdot 3 + 2 \cdot (-1) + a \cdot 2 + 1 \cdot 1 = 2 + 2a$. It vanishes at $a = -1$.

A shadow on a plane

The proof of the inequality solved the problem of the nearest point on a line along the way. The shadow $\mathbf p = \frac{\langle \mathbf b, \mathbf a\rangle}{\langle \mathbf a, \mathbf a\rangle}\,\mathbf a$ of the vector $\mathbf b$ on the line of the vector $\mathbf a$ is a vector whose signed length is the projection from the chapter on vectors: the useful part of the barge haulers’ effort that pulled the barge along the river. Now we need a shadow not on a line but on a plane, and indeed on any subspace. The main property of the shadow doesn’t change.

The idea: join $\mathbf b$ to any other point $\mathbf y$ of the subspace, and you get a right triangle in which the perpendicular is shorter than the hypotenuse. The picture shows the plane containing $\mathbf b$, $\mathbf p$ and $\mathbf y$; in $\mathbb R^n$ there’s always such a plane.

Let $\mathbf p \in W$ with $\mathbf b - \mathbf p$ orthogonal to $W$. Take any other point $\mathbf y \in W$. The difference $\mathbf p - \mathbf y$ also lies in $W$: a subspace is closed under subtraction. So $\langle \mathbf b - \mathbf p, \mathbf p - \mathbf y\rangle = 0$, and the angle of the triangle at the vertex $\mathbf p$ is a right angle. Since $\mathbf b - \mathbf y = (\mathbf b - \mathbf p) + (\mathbf p - \mathbf y)$, expanding as in the proof of the Cauchy–Schwarz inequality gives Pythagoras’ theorem: $|\mathbf b - \mathbf y|^2 = |\mathbf b - \mathbf p|^2 + |\mathbf p - \mathbf y|^2$. The second term is positive for $\mathbf y \ne \mathbf p$, so $\mathbf y$ is further from $\mathbf b$ than $\mathbf p$ is. Uniqueness is proved at the same time. Conversely: suppose $\mathbf b - \mathbf q$ is not orthogonal to some vector $\mathbf w \in W$. Moving from $\mathbf q$ along $\mathbf w$, we go along a chord of the circle with centre $\mathbf b$ and radius $|\mathbf b - \mathbf q|$ that passes through $\mathbf q$. The chord isn’t a tangent, because the radius at $\mathbf q$ isn’t perpendicular to it, and inside the circle there are points of $W$ nearer to $\mathbf b$. In numbers: $|\mathbf b - \mathbf q - s\mathbf w|^2 = |\mathbf b - \mathbf q|^2 - 2s\langle \mathbf b - \mathbf q, \mathbf w\rangle + s^2|\mathbf w|^2$, and for a small $s$ of the right sign this is less than $|\mathbf b - \mathbf q|^2$. The nearest point of a subspace is the foot of the perpendicular dropped from $\mathbf b$, and only that.

Back to the spring. The subspace $W$ is the image of the matrix $A$, and its vectors have the form $A\mathbf x$. The difference $\mathbf b - A\hat{\mathbf x}$ is orthogonal to $W$ exactly when it’s orthogonal to each column: $\mathbf a_j^{\mathsf T}(\mathbf b - A\hat{\mathbf x}) = 0$. The rows $\mathbf a_j^{\mathsf T}$ are the rows of the transposed matrix, so all the conditions together are written as $A^{\mathsf T}(\mathbf b - A\hat{\mathbf x}) = \mathbf 0$.

A square symmetric matrix made of the dot products of the columns of $A$ with each other: cell $(i, j)$ holds $\langle \mathbf a_i, \mathbf a_j\rangle$. Its size is the number of unknowns, not the number of measurements. The best approximate solution: for it the sum of the squared residuals $|\mathbf b - A\mathbf x|^2$ is smallest. The dot products of the columns with the right-hand side. Example: for the spring $A^{\mathsf T}A = \left(\begin{smallmatrix} 3 & 6 \\ 6 & 14 \end{smallmatrix}\right)$ and $A^{\mathsf T}\mathbf b = (5, 11)$. The system $3C + 6D = 5$, $6C + 14D = 11$ gives $D = \frac12$, $C = \frac23$. The line $y = \frac23 + \frac t2$ predicts $\frac76$, $\frac53$, $\frac{13}{6}$, the residuals are $-\frac16$, $\frac13$, $-\frac16$, and the sum of their squares is $\frac16$. A check of perpendicularity: $-\frac16 + \frac13 - \frac16 = 0$ and $-\frac16 + \frac23 - \frac12 = 0$.

One question remains: do the normal equations always have a unique solution?

One and the same problem lives in two pictures at once. In the plane of the measurements we move the line and look at the misses. In the space where each axis is one measurement we drop a perpendicular from $\mathbf b$ to the plane of the columns. Compare them.

Experiment 2. The best line

Before deriving the general formula, try to beat the normal equations by hand. The squares in the picture are literally the squared residuals, and your job is to make their total area as small as possible.

For $n$ points $(t_i, y_i)$ the matrix $A$ has two columns: ones and the values $t_i$. The normal equations become a system of two equations that can be solved once and for all. Write averages with a bar: $\bar t = \frac1n\sum t_i$, $\bar y = \frac1n\sum y_i$.

The line $y = C + Dt$ is drawn through $50$ points by the method of least squares. What’s the sum of the residuals $e_1 + e_2 + \dots + e_{50}$?

Exactly zero, for any points The smallest possible number The sum of the squared residuals It depends on how the points are placed

The residual is orthogonal to every column of $A$, including the column of ones: $\langle (1, \dots, 1), \mathbf e\rangle = e_1 + \dots + e_n = 0$. The misses above and below balance each other. If the line is forced through the origin (without $C$), there’s no column of ones, and this property is lost.

Squares have a downside. A miss of $10$ marks weighs as much as a hundred misses of one mark, so a single gross error, a knocked sensor or a typo in a table, can swing a line built from a thousand good points. Try it in the widget. You can minimise the sum of the absolute values of the residuals instead: an outlier shifts that line much less, but it has no formula, and it’s sometimes not unique.

And one more limitation. We counted the misses vertically, because the mass of the weight is known exactly and it’s the scale that errs. If both coordinates are inaccurate, it’s fairer to measure the distance to the line along the perpendicular, and then the answer is different: experiment 5 will lead to it.

On the other hand, the word “linear” in the name of the method refers to the unknowns, not to the graph. A parabola $y = C + Dt + Et^2$ is fitted by the same normal equations, adding to $A$ a third column of the numbers $t_i^2$: the unknowns $C$, $D$, $E$ still enter each equation to the first power. An exponential $y = ke^{rt}$ is reduced to a line by taking logarithms: $\ln y = \ln k + rt$.

Four measurements: $(0, 1)$, $(1, 3)$, $(2, 4)$, $(3, 4)$. What value of $y$ at $t = 4$ does the least-squares line predict? Give the answer as a number.

The averages: $\bar t = 1.5$, $\bar y = 3$. The deviations of $t$: $-1.5$, $-0.5$, $0.5$, $1.5$; the deviations of $y$: $-2$, $0$, $1$, $1$. The numerator is $3 + 0 + 0.5 + 1.5 = 5$, the denominator $2.25 + 0.25 + 0.25 + 2.25 = 5$, so $D = 1$ and $C = 3 - 1.5 = 1.5$. At $t = 4$ the line $y = 1.5 + t$ gives $5.5$.

Projections, lines and normal equations can be practised.

A tool: right angles instead of equations

The normal equations have to be solved because the columns $\mathbf a_1$ and $\mathbf a_2$ aren’t perpendicular: there’s $22^\circ$ between them, and the shadow on one column “spills over” onto the other. If they were perpendicular and of unit length, the matrix $A^{\mathsf T}A$ would become the identity and there’d be nothing to solve: $\hat{\mathbf x} = A^{\mathsf T}\mathbf b$.

Coordinates in such a basis are computed without any systems. If $\mathbf v = c_1\mathbf q_1 + \dots + c_k\mathbf q_k$, take the dot product of both sides with $\mathbf q_j$: all the terms except the $j$th vanish, leaving $c_j = \langle \mathbf v, \mathbf q_j\rangle$. For the same reason the projection of any $\mathbf b$ onto the span equals $\langle \mathbf b, \mathbf q_1\rangle\mathbf q_1 + \dots + \langle \mathbf b, \mathbf q_k\rangle\mathbf q_k$: the difference between $\mathbf b$ and this sum is orthogonal to each $\mathbf q_j$. The Fourier coefficients from chapter 33 were computed in exactly the same way. It remains to learn how to straighten a basis.

We build the vectors one at a time: from each new one we subtract its shadows on those already built, and normalise the remainder. The picture shows the first two steps; after that everything repeats.

The vector $\mathbf a_1$ is non-zero, since otherwise the set would be dependent. Put $\mathbf q_1 = \mathbf a_1 / |\mathbf a_1|$, a unit vector with the same span. Subtract from $\mathbf a_2$ its shadow on $\mathbf q_1$: $\mathbf a_2' = \mathbf a_2 - \langle \mathbf a_2, \mathbf q_1\rangle\mathbf q_1$. The remainder is orthogonal to $\mathbf q_1$: $\langle \mathbf a_2', \mathbf q_1\rangle = \langle \mathbf a_2, \mathbf q_1\rangle - \langle \mathbf a_2, \mathbf q_1\rangle\langle \mathbf q_1, \mathbf q_1\rangle = 0$, since $\langle \mathbf q_1, \mathbf q_1\rangle = 1$. The remainder isn’t zero, since otherwise $\mathbf a_2$ would be a multiple of $\mathbf a_1$. Put $\mathbf q_2 = \mathbf a_2' / |\mathbf a_2'|$. The span hasn’t changed: $\mathbf q_1$, $\mathbf q_2$ are expressed through $\mathbf a_1$, $\mathbf a_2$ and vice versa, since $\mathbf a_2 = |\mathbf a_2'|\,\mathbf q_2 + \langle \mathbf a_2, \mathbf q_1\rangle\mathbf q_1$. At step $j$ we subtract the shadows on all the vectors built so far: $\mathbf a_j' = \mathbf a_j - \sum_{i < j}\langle \mathbf a_j, \mathbf q_i\rangle\mathbf q_i$. Taking the dot product with $\mathbf q_l$ for $l < j$, we get $\langle \mathbf a_j, \mathbf q_l\rangle - \langle \mathbf a_j, \mathbf q_l\rangle = 0$: only one term of the sum survives, the others are killed by the orthogonality of the $\mathbf q_i$. The remainder isn’t zero, because $\mathbf a_j$ doesn’t lie in the span of the preceding $\mathbf a$’s, which coincides with the span of the preceding $\mathbf q$’s. Normalise: $\mathbf q_j = \mathbf a_j' / |\mathbf a_j'|$. Induction on $j$ gives the whole set, and at each step the span is preserved for the same reasons as at the second.

Let’s straighten the columns of the spring problem. The first gives $\mathbf q_1 = \frac{1}{\sqrt3}(1, 1, 1)$. The shadow of the second on it is $\frac{6}{3}(1, 1, 1) = (2, 2, 2)$, and the remainder is $(1, 2, 3) - (2, 2, 2) = (-1, 0, 1)$. Subtracting the shadow on the column of ones means subtracting from the $t_i$ their average $\bar t = 2$: that’s where the deviations in the slope formula come from. Now $\mathbf q_2 = \frac{1}{\sqrt2}(-1, 0, 1)$, and the projection of $\mathbf b = (1, 2, 2)$ comes out without a single equation: $\frac53(1, 1, 1) + \frac12(-1, 0, 1) = \left(\frac76, \frac53, \frac{13}{6}\right)$, the same predictions as before.

The process bears the names of Jørgen Gram (1883) and Erhard Schmidt (1907). Written with matrices, it gives the decomposition $A = QR$: the columns of $Q$ are orthonormal, and $R$ is upper triangular, because $\mathbf a_j$ is expressed only through $\mathbf q_1, \dots, \mathbf q_j$.

In the plane there are very few orthogonal matrices. The first column is a unit vector $(\cos\varphi, \sin\varphi)$, the second is a unit vector perpendicular to it, and there are two of those: $(-\sin\varphi, \cos\varphi)$ and $(\sin\varphi, -\cos\varphi)$. In the first case we get a rotation through the angle $\varphi$ with determinant $1$, in the second a reflection with determinant $-1$. These are the familiar motions from the chapter on symmetry, the ones that keep the origin in place. For example, $\frac15\left(\begin{smallmatrix} 3 & -4 \\ 4 & 3 \end{smallmatrix}\right)$ rotates the plane through $\arctan\frac43 \approx 53.13^\circ$. Any permutation matrix, with a single one in each row and each column, is orthogonal too: it only rearranges the coordinates.

Experiment 3. A circle turns into an ellipse

A rectangular matrix has no eigenvectors: it takes vectors from one space to another, and $A\mathbf v$ can’t be compared with $\mathbf v$. Even a square matrix, a shear for example, may have few eigendirections. But let’s ask differently: what does a matrix turn the unit circle into?

Whatever matrix you choose, the circle becomes an ellipse (or a segment), and an ellipse has two perpendicular semi-axes. The most surprising thing is their preimages: the circle always has two perpendicular vectors $\mathbf v_1$, $\mathbf v_2$ that the matrix takes to these semi-axes, that is, to perpendicular vectors again. Even a shear, which has only one eigendirection, preserves one pair of right angles. This holds for any matrix.

It’s not $A$ that has eigenvectors but the symmetric matrix $A^{\mathsf T}A$, and for it the hunt from the last chapter always succeeds. Its eigenvectors will be the perpendiculars we need. On the left of the picture is the space where the $\mathbf v$’s live, on the right the space where the $A\mathbf v$’s live.

The $n \times n$ matrix $A^{\mathsf T}A$ is symmetric: $(A^{\mathsf T}A)^{\mathsf T} = A^{\mathsf T}(A^{\mathsf T})^{\mathsf T} = A^{\mathsf T}A$. By the spectral theorem it has an orthonormal basis of eigenvectors $\mathbf v_1, \dots, \mathbf v_n$ with eigenvalues $\lambda_1, \dots, \lambda_n$. Number them in decreasing order. All the $\lambda_i$ are non-negative: $\lambda_i = \lambda_i\langle \mathbf v_i, \mathbf v_i\rangle = \langle \mathbf v_i, A^{\mathsf T}A\mathbf v_i\rangle = \langle A\mathbf v_i, A\mathbf v_i\rangle = |A\mathbf v_i|^2$. Put $\sigma_i = \sqrt{\lambda_i} = |A\mathbf v_i|$: these are the lengths of the images. The images are perpendicular: for $i \ne j$ we have $\langle A\mathbf v_i, A\mathbf v_j\rangle = \langle \mathbf v_i, A^{\mathsf T}A\mathbf v_j\rangle = \lambda_j\langle \mathbf v_i, \mathbf v_j\rangle = 0$. Let $\sigma_1, \dots, \sigma_r$ be the non-zero ones among these numbers. The vectors $\mathbf u_i = A\mathbf v_i / \sigma_i$ for $i \le r$ are unit vectors and pairwise orthogonal; extend them to an orthonormal basis of $\mathbb R^m$ by the Gram–Schmidt process. For $i > r$ we have $|A\mathbf v_i| = 0$, that is, $A\mathbf v_i = \mathbf 0$. The vectors $A\mathbf v_1, \dots, A\mathbf v_n$ span the whole image of $A$, the non-zero ones among them are the first $r$, and non-zero pairwise orthogonal vectors are linearly independent. So $r$ equals the rank. The equations $A\mathbf v_i = \sigma_i\mathbf u_i$ are written all at once as $AV = U\Sigma$, where the columns of $V$ and $U$ are the vectors $\mathbf v_i$ and $\mathbf u_i$, and $\Sigma$ is a rectangular matrix with the numbers $\sigma_i$ on the diagonal and zeros elsewhere. The matrix $V$ is orthogonal, $V^{-1} = V^{\mathsf T}$, and $A = U\Sigma V^{\mathsf T}$. The point of the circle $\cos s\,\mathbf v_1 + \sin s\,\mathbf v_2$ goes to $\sigma_1\cos s\,\mathbf u_1 + \sigma_2\sin s\,\mathbf u_2$: that’s the ellipse with semi-axes $\sigma_1$ and $\sigma_2$. An orthogonal $n \times n$ matrix: a rotation, perhaps with a reflection. It takes $\mathbf v_i$ to $\mathbf e_i$: first the space is turned so that the directions we need lie along the axes. A rectangular diagonal $m \times n$ matrix: it stretches the $i$th axis by a factor of $\sigma_i$, and drops the extra coordinates or pads with zeros if $m \ne n$. An orthogonal $m \times m$ matrix: it turns the axes $\mathbf e_i$ into the directions $\mathbf u_i$. As always, read from right to left. Example: for $A = \left(\begin{smallmatrix} 3 & 0 \\ 4 & 5 \end{smallmatrix}\right)$ the matrix $A^{\mathsf T}A = \left(\begin{smallmatrix} 25 & 20 \\ 20 & 25 \end{smallmatrix}\right)$ has eigenvalues $45$ and $5$ with vectors $\mathbf v_1 = \frac{1}{\sqrt2}(1, 1)$, $\mathbf v_2 = \frac{1}{\sqrt2}(1, -1)$. So $\sigma_1 = 3\sqrt5 \approx 6.708$ and $\sigma_2 = \sqrt5 \approx 2.236$. The images: $A\mathbf v_1 = \frac{1}{\sqrt2}(3, 9)$ and $A\mathbf v_2 = \frac{1}{\sqrt2}(3, -1)$; they’re perpendicular, and $\mathbf u_1 = \frac{1}{\sqrt{10}}(1, 3)$, $\mathbf u_2 = \frac{1}{\sqrt{10}}(3, -1)$. Check: $\sigma_1\sigma_2 = 15 = |\det A|$.

The proof shows how to compute them: the singular values are the square roots of the eigenvalues of $A^{\mathsf T}A$. For a square matrix their product is $|\det A|$: rotations don’t change area, and $\Sigma$ stretches it by a factor of $\sigma_1\sigma_2\cdots\sigma_n$. For a symmetric matrix with non-negative eigenvalues the singular value decomposition coincides with the spectral one. And for the shear $\left(\begin{smallmatrix} 1 & 1 \\ 0 & 1 \end{smallmatrix}\right)$, from which the eigenvector hunter came away almost empty-handed, the singular values are the golden ratio $\varphi \approx 1.618$ and $1/\varphi \approx 0.618$: $A^{\mathsf T}A = \left(\begin{smallmatrix} 1 & 1 \\ 1 & 2 \end{smallmatrix}\right)$, and its characteristic equation $\lambda^2 - 3\lambda + 1 = 0$ has roots $\varphi^2$ and $\varphi^{-2}$.

Find the singular values of the matrix $\left(\begin{smallmatrix} 3 & -8 \\ 4 & 6 \end{smallmatrix}\right)$.

The columns $(3, 4)$ and $(-8, 6)$ are orthogonal: $-24 + 24 = 0$. So $A^{\mathsf T}A = \left(\begin{smallmatrix} 25 & 0 \\ 0 & 100 \end{smallmatrix}\right)$, the eigenvalues are $100$ and $25$, and the singular values are $10$ and $5$, the lengths of the columns. The right singular vectors here are simply $\mathbf e_2$ and $\mathbf e_1$ (the one stretched more comes first): the circle is stretched along the coordinate axes, and $U$ rotates the result. Check: $10 \cdot 5 = 50 = |3 \cdot 6 - (-8) \cdot 4|$.

Experiment 4. A photograph in layers

A black-and-white photograph is a table of brightnesses, that is, a matrix. Its singular value decomposition can be rewritten as a sum: the column $\mathbf u_i$ times the row $\mathbf v_i^{\mathsf T}$ gives a matrix of the same size as $A$, and

$$A = \sigma_1\mathbf u_1\mathbf v_1^{\mathsf T} + \sigma_2\mathbf u_2\mathbf v_2^{\mathsf T} + \dots + \sigma_r\mathbf u_r\mathbf v_r^{\mathsf T}.$$

Each term is a picture of rank $1$. In it all the rows are proportional to a single row $\mathbf v_i^{\mathsf T}$, and the brightness of a pixel is the product of a “row weight” and a “column weight”. That’s what tartan looks like: vertical stripes multiplied by horizontal ones. A photograph is a sum of tartans, and the terms come in decreasing order of loudness $\sigma_i$.

The loudest layer. For a photograph it’s usually a blurred “average” image: the overall lighting by rows and columns. The last layer kept. Everything after it is thrown away. Example: to store $A_k$ for an $m \times n$ photograph we need $k$ columns $\mathbf u_i$ of $m$ numbers, $k$ rows $\mathbf v_i^{\mathsf T}$ of $n$ numbers and $k$ numbers $\sigma_i$, $k(m + n + 1)$ in all. For a $1000 \times 1000$ picture and $k = 50$ that’s $100{,}050$ numbers instead of a million, about $10\,\%$.

Why should it be the quiet layers that are thrown away and not some others? To answer, we have to agree on how to measure the error of an approximation.

The trick is a count of dimensions: a matrix of low rank has a large kernel, and it’s bound to meet the directions that $A$ stretches most. The picture shows the case $2 \times 2$ and $k = 1$; in general the argument is the same.

First $A_k$ itself. The difference $A - A_k = \sigma_{k+1}\mathbf u_{k+1}\mathbf v_{k+1}^{\mathsf T} + \dots$ is already written as a singular value decomposition, and its largest stretch is $\sigma_{k+1}$. In the picture $A - A_1$ squashes the circle into a segment of length $2\sigma_2$ along $\mathbf u_2$. Now any $B$ of rank at most $k$. By the rank–nullity theorem (chapter 37) its kernel, the vectors that $B$ sends to zero, has dimension at least $n - k$. In the picture the kernel is a line; rotate it, and that’s the choice of $B$. The span of the vectors $\mathbf v_1, \dots, \mathbf v_{k+1}$ has dimension $k + 1$, and $(n - k) + (k + 1) > n$. Put together bases of the kernel and of the span: that’s at least $n + 1$ vectors in $\mathbb R^n$, and they’re linearly dependent. In the dependence, move the terms from the kernel’s basis to the left and the rest to the right: on the left is a vector of the kernel, on the right the same vector of the span. It isn’t zero, since otherwise both sides would be zero combinations of independent vectors and all the coefficients would be zero. Take this common vector with unit length and call it $\mathbf x$. In the plane $\mathbf v_1$ and $\mathbf v_2$ span everything, and $\mathbf x$ is simply a unit vector of the kernel. Since $B\mathbf x = \mathbf 0$, we get $(A - B)\mathbf x = A\mathbf x$. Expand $\mathbf x = c_1\mathbf v_1 + \dots + c_{k+1}\mathbf v_{k+1}$ with $c_1^2 + \dots + c_{k+1}^2 = 1$. Then $A\mathbf x = \sum c_i\sigma_i\mathbf u_i$ and $|A\mathbf x|^2 = \sum c_i^2\sigma_i^2 \ge \sigma_{k+1}^2$. In the picture $A\mathbf x$ lies on the ellipse, and the ellipse doesn’t go inside the circle of radius $\sigma_2$. So $\|A - B\| \ge |(A - B)\mathbf x| \ge \sigma_{k+1} = \|A - A_k\|$: no matrix of rank $k$ approximates $A$ better than the sum of the first $k$ layers.

And if the error is measured by the sum of squares over all the pixels

For a picture it’s more natural to add up the squared errors at all the pixels: $\|M\|_F^2 = \sum_{i,j} m_{ij}^2$ (the Frobenius norm). Multiplying by an orthogonal matrix on the left doesn’t change the lengths of the columns, and on the right the lengths of the rows, so $\|M\|_F^2 = \|\Sigma\|_F^2 = \sigma_1^2(M) + \sigma_2^2(M) + \dots$. For $A_k$ the error is $\sigma_{k+1}^2 + \dots + \sigma_r^2$, and this is also the best result.

We’ll need an estimate: for subspaces $U$ and $W$ of the space $\mathbb R^n$ the intersection has dimension at least $\dim U + \dim W - n$. It follows from the rank–nullity theorem applied to the map of pairs $(\mathbf y, \mathbf z) \mapsto \mathbf y - \mathbf z$ from $U \times W$ (dimension $\dim U + \dim W$) to $\mathbb R^n$: its kernel consists of the pairs $(\mathbf y, \mathbf y)$ with $\mathbf y \in U \cap W$, and its image is no bigger than $\mathbb R^n$.

Let’s prove that for any $B$ of rank at most $k$ and any $i \ge 1$, $\sigma_{k+i}(A) \le \sigma_i(A - B)$. Take three subspaces: the span of $\mathbf v_1, \dots, \mathbf v_{k+i}$ (dimension $k + i$), the kernel of $B$ (at least $n - k$) and the span of the right singular vectors of the matrix $A - B$ with numbers from $i$ to $n$ (dimension $n - i + 1$). By the estimate the intersection of the first two has dimension at least $(k + i) + (n - k) - n = i$, and its intersection with the third at least $i + (n - i + 1) - n = 1$. Take a unit vector $\mathbf x$ in the common intersection. As in the main proof, $|A\mathbf x| \ge \sigma_{k+i}(A)$. On the other hand, $A\mathbf x = (A - B)\mathbf x$, and $\mathbf x$ expands in the singular vectors of $A - B$ with numbers from $i$ on, which are stretched by at most $\sigma_i(A - B)$, so $|(A - B)\mathbf x| \le \sigma_i(A - B)$. It remains to add up the squares: $\|A - B\|_F^2 = \sum_i \sigma_i^2(A - B) \ge \sum_i \sigma_{k+i}^2(A) = \|A - A_k\|_F^2$.

The theorem was published in 1936 by the physicist Carl Eckart and the mathematician Gale Young, in the journal Psychometrika, of all places: psychologists needed to describe large tables of test results by a few common “abilities”. In computer compression of photographs SVD is hardly used: the vectors $\mathbf u_i$, $\mathbf v_i$ themselves have to be stored, different for each picture. The JPEG format expands $8 \times 8$ blocks of pixels in the cosines from the chapter on Fourier series, the same for all images. But the idea is the same: keep the loud terms and throw out the quiet ones, and the Eckart–Young theorem says that for one particular matrix nothing does better than SVD.

A picture is $400 \times 300$ pixels. How many numbers have to be stored for its rank-$20$ approximation?

By the formula $k(m + n + 1) = 20 \cdot (400 + 300 + 1) = 14{,}020$. That’s about $11.7\,\%$ of the $120{,}000$ pixels of the original picture.

Experiment 5. The cloud and its principal axis

Back to the cloud from the end of the last chapter: the heights and weights of a thousand people. Here both coordinates are inaccurate, and neither can be declared “correct”. It’s natural to measure the misses along the perpendicular to the line. This problem was posed by Karl Pearson in 1901 in the paper On Lines and Planes of Closest Fit to Systems of Points in Space.

Move the origin to the centre of mass of the cloud, subtracting from each coordinate its average, and write the points as the rows of an $n \times 2$ matrix $X$ (or $n \times d$, if there are $d$ features).

For each point Pythagoras’ theorem works: the square of its distance to the centre splits into the shadow along the line and the miss across it. The total doesn’t depend on the line, so the misses can only be reduced by lengthening the shadows.

Let $\mathbf w$ be the unit vector of the line’s direction, and $\mathbf x_1, \dots, \mathbf x_n$ the points after moving to the centre of mass. A point’s shadow on the line has length $|\langle \mathbf x_i, \mathbf w\rangle|$, and the miss $d_i$ is the distance from the point to the line. By Pythagoras’ theorem $|\mathbf x_i|^2 = \langle \mathbf x_i, \mathbf w\rangle^2 + d_i^2$. Add over all the points: $\sum |\mathbf x_i|^2 = \sum \langle \mathbf x_i, \mathbf w\rangle^2 + \sum d_i^2$. The left side doesn’t depend on the line; the bars in the picture always have the same height. The misses are smaller exactly when the sum of the squared shadows is bigger. The sum of the squared shadows is a quadratic form: $\sum \langle \mathbf x_i, \mathbf w\rangle^2 = |X\mathbf w|^2 = \mathbf w^{\mathsf T}(X^{\mathsf T}X)\mathbf w$. The matrix $S = X^{\mathsf T}X$ is symmetric, and by the spectral theorem it has orthonormal eigenvectors $\mathbf v_1, \mathbf v_2$ with $\lambda_1 \ge \lambda_2 \ge 0$. Expand $\mathbf w = c_1\mathbf v_1 + c_2\mathbf v_2$, where $c_1^2 + c_2^2 = 1$. Then $\mathbf w^{\mathsf T}S\mathbf w = \lambda_1c_1^2 + \lambda_2c_2^2 \le \lambda_1(c_1^2 + c_2^2) = \lambda_1$, with equality at $\mathbf w = \pm\mathbf v_1$. In $d$ dimensions the calculation is the same, only with $d$ terms. The best line runs along $\mathbf v_1$; the sum of the squared shadows on it is $\lambda_1 = \sigma_1^2$, and the sum of the squared misses is $\lambda_2 = \sigma_2^2$ (in $d$ dimensions, the sum of the remaining $\lambda_i$).

Through a height–weight cloud one least-squares line is drawn predicting weight from height, and a second predicting height from weight. How are they placed?

They’re the same line Both pass through the centre of mass, and the principal axis lies between them They’re parallel They’re perpendicular

Each least-squares line passes through the centre of mass. The line “weight from height” is flatter than the principal axis, and “height from weight”, drawn in the same axes, is steeper. Vertical misses “press” the line towards the horizontal, horizontal ones towards the vertical, while perpendicular ones favour no axis. The closer the cloud is to a line, the closer all three are.

The term “principal components” was introduced by Harold Hotelling in 1933, and since then the method has become the first thing done with multidimensional data. If there are a hundred features but the principal axes that hold almost all the spread number three or four, the data can be viewed in three or four dimensions. In 2008 a group of geneticists led by John Novembre studied about half a million genetic markers in a sample of three thousand Europeans and drew the people as points by their first two principal components. A recognisable map of Europe emerged: neighbours on the map turned out to be neighbours in their genes too.

Here too it’s easy to go wrong. The principal axes depend on the units of measurement: convert height from metres to centimetres, and the spread in height grows $10{,}000$-fold, and the principal axis turns towards it. So the features are usually brought to a common scale first. And if the cloud is almost round, $\lambda_1 \approx \lambda_2$, the principal axis is unstable: move one point, and it swings round. In the widget you can see this at once.

The method is directly linked to the singular value decomposition. If $X = U\Sigma V^{\mathsf T}$, then $X^{\mathsf T}X = V\Sigma^{\mathsf T}\Sigma V^{\mathsf T}$: the principal axes are the columns of $V$, and the spread along the $i$th axis is $\sigma_i^2 / n$. The Eckart–Young approximation $X_k$ is the cloud projected onto the first $k$ principal axes.

Experiment 6. Who will watch what

The last experiment is with a table in which almost nothing is filled in. In October 2006 Netflix released $100{,}480{,}507$ ratings that $480{,}189$ viewers had given to $17{,}770$ films, and offered a million-dollar prize to whoever could predict the ratings $10\,\%$ more accurately than its own program. The “viewers × films” table is a little more than $1\,\%$ filled, and the task is to guess the remaining cells.

What helps is an assumption of low rank. Suppose a viewer’s taste is described by a few numbers, how much they like comedies, horror, art-house films, and a film is described by the same number of numbers. A rating is then roughly the dot product of the taste vector and the film vector, and the whole table is the product of two narrow matrices, that is, a matrix of low rank. The singular value decomposition of an empty table can’t be computed, but the narrow matrices can be fitted by least squares using only the known cells. Fix the film vectors, and each viewer’s vector is found by an ordinary least-squares problem from experiment 2. Fix the viewers, and the films are found the same way. Alternating, we quickly arrive at a good approximation.

In December 2006 the programmer Brandyn Webb, under the pseudonym Simon Funk, described such a model on his blog: the narrow factors are fitted using only the known ratings (he fitted them by gradient descent rather than by alternation). He called it SVD, and the method quickly spread among the contestants. The prize was awarded in September 2009 to the team BellKor’s Pragmatic Chaos, which improved the result by $10.06\,\%$. The team The Ensemble achieved the same result but submitted its solution $20$ minutes later. The winning program blended more than a hundred models, and factorising the matrix into narrow factors was among the main ones.

Where next

The orthogonal $2 \times 2$ matrices turned out to be rotations and reflections. Take from among them the eight that take a square centred at the origin to itself: four rotations (through $0^\circ$, $90^\circ$, $180^\circ$, $270^\circ$) and four reflections. Multiply any two, and you get one of the eight again. These matrices have their own $8 \times 8$ multiplication table in which nothing foreign ever appears. Similar tables exist for the turns of the faces of a Rubik’s cube, for shuffles of a deck of cards and for adding hours on a clock face, although there are neither numbers nor matrices there. What they all have in common, and what follows from it, is the subject of chapter 40.