One subtraction clears a direction

The previous rung ended with a vector nobody was looking for. Proving Cauchy–Schwarz produced, as the thing whose square is the gap between the two sides, the part of $a$ that is left when its shadow on $b$ has been taken away. That leftover was measured and discarded.

This essay keeps it.

That is the entire construction. It has a name, two inventors and a numerical hazard, and applied to the powers of $x$ it manufactures a family of special functions that nobody was aiming at.

Why a basis wants to be perpendicular #

A basis is a set of directions such that every vector is a combination of them in exactly one way. Any independent set that is large enough will do, and most of them are unpleasant to work with for one reason: finding the coefficients means solving equations.

Given $v = c_1u_1 + c_2u_2 + \cdots$, the numbers $c_i$ are what a coordinate system is for, and extracting them from a general basis means inverting a matrix. The work grows as the cube of the dimension, the answer is sensitive to how nearly dependent the $u_i$ are, and nothing about the picture suggests where to start.

For a perpendicular basis there are no equations. Take the inner product of both sides with $u_1$: every term except the first contributes nothing, because $u_1$ casts no shadow on anything perpendicular to it, and what is left is $v\cdot u_1 = c_1(u_1\cdot u_1)$. So

$$c_1 = \frac{v\cdot u_1}{u_1\cdot u_1},$$

and if the basis vectors have been scaled to length one the denominator is one and the coefficient is simply $v\cdot u_1$ — a shadow, computed in one pass and independently of every other coefficient.

That independence is the whole prize. Each coordinate can be found without knowing any of the others, which means the calculation parallelises, that adding a new direction does not disturb the coefficients already computed, and that dropping one leaves the rest still correct. Fourier coefficients are computed one at a time for exactly this reason and no other.

It also makes truncation meaningful. Keep the first few coefficients of a general basis and the result is not the best approximation using those directions — it is not even close, because the discarded terms were propping up the kept ones. Keep the first few coefficients of a perpendicular basis and the result is the best approximation available from those directions, exactly, with the error being the leftover and nothing else. Every expansion that is truncated in practice — a Fourier series cut off at the twentieth harmonic, a signal compressed by dropping small coefficients, a polynomial fit stopped at degree four — depends on that, and none of it works in a basis whose members overlap.

There is a third property, less often stated and more often relied on. In a perpendicular basis the squared length of a vector is the sum of the squared coefficients, which is Pythagoras with as many terms as there are directions. So a coefficient’s size is a direct measure of how much of the vector it accounts for, and comparing two coefficients is meaningful. In a general basis it is not: a large coefficient may be cancelling another large one, and the two together may describe a small vector.

The step #

Given $v_1$ and $v_2$, keep $v_1$ and replace $v_2$ by what is left of it after its shadow on $v_1$ is removed:

$$u_2 = v_2 - \frac{v_2\cdot v_1}{v_1\cdot v_1}\,v_1 .$$

Two things have to be true for that to be worth doing, and both are checked in the figures rather than asserted in the prose.

The result is perpendicular to $v_1$. Take the inner product with $v_1$ and the two terms cancel by construction — which is the definition of the coefficient, run backwards.

The pair still spans the same plane. This is the half that is easy to lose and impossible to see. An orthogonal pair that happens to span something else is not a basis for the space that was handed over, and no drawing of two perpendicular arrows says which plane they came from. The check the generator makes is a determinant: $u_2$ differs from $v_2$ by a multiple of $v_1$, and adding a multiple of one column to another leaves the area they span unchanged, so the two parallelograms have exactly equal area and the spans coincide.

With more vectors the step repeats. At each stage, subtract the shadows on everything already built, and what is left is perpendicular to all of them:

$$u_k = v_k - \sum_{j<k}\frac{v_k\cdot u_j}{u_j\cdot u_j}\,u_j .$$

What the coefficients do when they are written down #

The numbers thrown away at each step are not waste. Collected, they say exactly how each original vector sits in the basis built from it, and they have a shape.

The $k$-th original vector is a combination of $u_1, \ldots, u_k$ and no others, because that is the order the construction ran in. So the table of coefficients is triangular — zeros above the diagonal — and the triangularity is a fact about the order of the algorithm rather than about the vectors.

Written as matrices, this says $A = QR$, with $Q$ the orthonormal directions as columns and $R$ the triangular table. It is one of the two or three decompositions the whole of computational linear algebra is organised around, and it arrives here as bookkeeping: nothing was computed for it that the construction was not computing anyway.

The diagonal entries are worth a moment. The $k$-th of them is the length of the leftover at step $k$ — the same quantity the previous rung’s parabola bottomed out at. A small diagonal entry means the $k$-th vector was nearly a combination of the ones before it, and a zero one means it was exactly that, which is how the construction reports dependence: it does not fail with a message, it hands back a vector of length nought.

Multiplying the diagonal entries together gives something with a name of its own. Each is the distance from a vector to the space its predecessors span, which is a height; the product of the heights is the volume of the box the original vectors span, in exactly the way a rectangular box’s volume is a product of its sides. So the construction computes the determinant on its way past, as a product of $n$ perpendicular distances rather than as a signed sum over $n!$ permutations — and the two agree, because the shear that removes a shadow is precisely the operation the determinant is blind to.

That reading also explains the inequality named after Hadamard: the volume of a box is at most the product of the lengths of its sides, since each height is at most the vector it is measured from. Equality holds exactly when every leftover is the whole vector, which is to say exactly when the vectors were already perpendicular. It is Cauchy–Schwarz in $n$ vectors rather than two, and it is the same sentence — a shadow is never longer than what casts it — counted $n$ times.

What happens when the vectors are not independent #

The construction was described for an independent set and the description did not say what it does otherwise, which is worth settling, because the answer is the useful one.

Suppose $v_3$ is a combination of $v_1$ and $v_2$. At step three, everything of $v_3$ is shadow: the subtraction removes all of it, and what is left is the zero vector. There is nothing to normalise, the diagonal entry of the table is nought, and the procedure has detected the dependence without having been asked to look for it.

That is the honest behaviour and it is why the construction doubles as a test. Deciding whether a set of vectors is independent is otherwise a question about a determinant or a rank, both of which are answers with no evidence attached; Gram–Schmidt answers it and hands over the combination as well, since the coefficients it computed at that step are exactly how $v_3$ is built from the two before it.

In floating-point arithmetic the leftover is never exactly zero, and this is where the diagnosis becomes a judgement. A leftover of $10^{-17}$ against inputs of size one means dependence; a leftover of $10^{-3}$ means the vectors are nearly dependent, which is a different and more awkward situation. Nothing in the arithmetic distinguishes “dependent, with rounding” from “independent, but barely”, and every practical implementation draws the line with a tolerance somebody chose.

The diagonal of the table is therefore the most informative thing the construction produces, and it is the part usually discarded. It is a list saying how much genuinely new direction each vector contributed, in order, and reading it says whether the basis is a basis in any useful sense or a set of nearly the same arrow written down several times.

What it costs #

Orthogonalising $n$ vectors in $m$ dimensions costs about $2mn^2$ multiplications: each of the $n$ steps takes inner products against everything built so far, and each inner product is $m$ multiplications. For a square problem that is $2n^3$, which is the same order as solving the equations directly and about twice the constant.

So Gram–Schmidt is not a way of making a calculation cheaper. What it buys is that the calculation is done once and every subsequent coordinate extraction is free — one pass per coordinate instead of a solve. A basis used once is not worth orthogonalising; a basis used repeatedly is worth it after the second use.

The other thing it buys is stability, or would if it were done in the right order, which is the subject of the next section. An orthonormal matrix does not amplify anything — it takes vectors to vectors of the same length, so an error present in the input comes out the same size rather than magnified. That single property is why so much numerical linear algebra factors matrices into orthogonal pieces, and it is the shadow argument again: a shadow is never longer than the arrow casting it, so nothing anywhere in the process can grow.

The comparison worth making is with the elimination everybody meets first. Gaussian elimination does the same thing with a weaker requirement: it subtracts multiples of one vector from the others to clear entries rather than to clear directions, so the coefficient it uses is a ratio of coordinates instead of a shadow, and no inner product appears. It is cheaper — about $\tfrac23 n^3$ against $2n^3$ — and it does not produce a perpendicular anything, so the matrix it leaves behind can amplify an error by an amount nothing bounds. The extra factor of three is the price of the guarantee, and which of the two is right depends entirely on whether the guarantee is needed.

Both are the same manoeuvre at bottom: take one vector as fixed, and remove from every other whatever part of it the fixed one accounts for. What differs is the sense of “accounts for”, and the whole difference in behaviour comes from that one choice.

The order the subtractions are done in #

Here is a fact that no amount of algebra will produce, because in exact arithmetic there is nothing to say.

Two orders give the same answer on paper. The classical order computes every coefficient against the original vector $v_k$ and subtracts all the shadows at the end. The modified order subtracts each shadow as it is computed, and takes the next coefficient against what is left rather than against the original. The two differ by rearranging terms that are exactly equal, and any student asked to prove them equal would do so in a line.

The test matrix is Läuchli’s: a row of ones over a diagonal of $\varepsilon$. Its columns are independent for every $\varepsilon$ above nought and are within $\varepsilon$ of being the same vector, which is exactly the situation the construction is meant to handle and exactly where it is fragile.

The reason for the difference is one subtraction. When $v_2$ is nearly parallel to $v_1$, the leftover $u_2$ is tiny — it is a difference of two nearly equal quantities, so most of its significant digits are lost. The classical order then computes the third coefficient against $v_3$, which knows nothing about that loss, so the shadow it subtracts is a shadow on a direction that is no longer quite the direction it was. The modified order takes the coefficient against the partially-reduced vector, so each subtraction is made against the vector that actually survives, and the errors do not compound.

The two orders are the same algebra and different computations, and nothing in the mathematics distinguishes them. This is the sharpest instance in this ladder of a distinction that runs through the whole subject: a derivation establishes what a quantity is, and says nothing whatever about what happens when a machine tries to produce it.

For serious work neither order is used. Householder reflections build the same $Q$ and $R$ by reflecting one vector onto an axis at a time, which never forms a small difference at all, and the orthogonality they produce is good to the arithmetic’s own precision whatever the input. Gram–Schmidt survives because it is the version that explains itself, and because in the modified order it is good enough for most purposes.

The powers of x, straightened #

Nothing in the construction mentioned arrows, so it runs wherever an inner product does. Give it the functions $1, x, x^2, x^3$ on the interval $[-1,1]$, with

$$\langle f,g\rangle = \int_{-1}^{1} f(x)g(x)\,dx,$$

and it will straighten them.

What comes out is the Legendre polynomials, scaled so that the leading coefficient is one. They were introduced in 1785 in a study of gravitational potential, they are the standard basis for expanding a function on an interval, and they are the nodes a Gaussian quadrature rule uses. None of that was being aimed at here: the input was the four dullest functions available and the procedure was the one that makes two arrows perpendicular.

The reason this is worth more than a curiosity is that the monomials are a terrible basis and the fact is invisible until they are orthogonalised. On $[-1,1]$ the functions $x^{10}$ and $x^{12}$ are almost the same function — both are nearly zero across the middle and shoot up at the ends — so the diagonal entry of the table at that step is minute, and fitting a polynomial by solving in the monomial basis is the numerically hopeless calculation that Runge’s phenomenon is usually blamed for. Orthogonalise first and the same fit is stable. The construction’s real output is not the perpendicular vectors; it is a diagnosis of the basis it was handed.

Two names and neither of them first #

Jørgen Pedersen Gram published the process in 1883, in a paper on least squares. Erhard Schmidt published it in 1907, in a paper on integral equations, and cited Gram. Neither was first: Laplace used it in 1812 and Cauchy in 1836, both without remarking on it, because to them it was an obvious way to arrange a calculation rather than a procedure worth naming.

That pattern is worth noticing because it recurs and because it says something about what a name is for. The construction is short enough that anybody needing it invents it, so being first is not the achievement; the achievement is isolating it — separating the manoeuvre from the problem it was used on, so that it can be handed to something else. Schmidt’s paper is the one that did that, which is why his name is on it and Laplace’s is not.

Gram’s own paper is worth a sentence for what it was about, because the subject was already the one at the top of this ladder. He was studying least squares — fitting a function to data by minimising a sum of squares — and the orthogonalisation was the device that made the fit computable one coefficient at a time. The construction and the application it exists for arrived together, and the next rung of this anchor is that application.

The numerical hazard was not noticed for a long time after either paper, for the plain reason that neither author was working with a machine. Rice described the difference between the two orders in 1966 and Björck analysed it in 1967, a hundred and thirty years after Cauchy first used the procedure and about fifteen years after anybody had a computer to lose digits on. The mathematics was finished long before the question the last section asks could even be posed, which is the usual relationship between a construction and its numerical behaviour.

What the pictures cannot show #

Every drawing here is of two vectors in a plane or three in space, and the construction’s value is entirely at sizes no drawing reaches. A basis of four is the largest table on the page; the calculations that matter run to thousands, where the triangle is the whole picture and the arrows are gone.

More sharply, the figures cannot show the failure they are about. The loss-of-orthogonality plot is a graph of numbers a computer produced, and every point on it would be on the axis if the arithmetic were exact — so what is drawn is not a property of Gram–Schmidt at all but a property of Gram–Schmidt plus a particular arithmetic, and no picture of the mathematics could contain it. A reader looking at the two-vector panels for the reason the classical order fails will not find it, because the reason is in the digits.

And the Legendre figure draws four curves and shows no orthogonality. That every pair integrates to nought is a statement about areas that cancel, and the areas are on both sides of the axis and are not drawn; the assertion behind the figure is a rational arithmetic exactly equal to nought, which is stronger than anything the eye could confirm and completely invisible.

The ladder from here #

Rungs above: projection onto a subspace, where the leftover is what a fit leaves behind and the normal equations are read off the picture. The adjoint, defined by moving a matrix across an inner product, and why a symmetric matrix has perpendicular eigenvectors — the same definition twice. Bessel’s inequality, which says the shadows on an orthonormal set can never add to more than the vector has, and Parseval’s identity, which says when they add to exactly that. The QR algorithm, where repeatedly factoring and multiplying back in the other order drives a matrix towards its own eigenvalues. And the cross product, which keeps the perpendicular part rather than measuring it, and exists in three dimensions and almost nowhere else.

Keeping what an argument threw away #

The habit is worth naming because it is one of the cheapest sources of new mathematics there is: when a proof produces an object on the way to a number, look at the object.

Cauchy–Schwarz needs only the size of the leftover, and the standard presentation discards the leftover the moment its size is known. Keeping it gives Gram–Schmidt, which gives orthonormal bases, which gives the QR factorisation, which gives an eigenvalue algorithm and a family of special functions. All of it was sitting inside a proof of an inequality, being thrown away at the last line.

The same move accounts for a good deal of the rest of this collection. The remainder in a division is what is left when a quotient has been taken, and keeping it is the whole of Euclid’s algorithm. The residual in a fit is what a model does not explain, and keeping it is the whole of the theory of errors. In each case the discarded object is the part of the input the procedure could not account for, which is exactly the part worth looking at next.