c08 · Section 3

Orthonormal bases and QR

Why is an orthogonal basis easier to use?

8. Projection and the Best FitL1 · L2 · L3E1 E2 E3

Terms: basis, condition number, matrix, orthonormal basis, QR factorisation, transpose

Symbols: q

L1Core

The idea, one concrete case and enough to use it. Every simplification is marked and linked to the precise statement.

A ship reports its position as 5 km east and 1 km north. Your chart is rotated: its first axis runs 3 east and 4 north for every 5 km, and its second axis runs 4 west and 3 north. Reading the ship's two chart coordinates normally means solving two equations in two unknowns, and the arithmetic grows with every extra axis you add.

There is a shortcut when the axes are perpendicular to one another and one unit long, and it replaces the solve with a pair of dot products. The chart is a basis: a set of axes that reaches every point of the space exactly once, so every list is a combination of the axes and no axis is a combination of the others. The position is one list written in two different bases, and the shortcut is what makes the second description cheap to compute. !Left: the axes q1=(0.6,0.8) and q2=(−0.8,0.6) are perpendicular unit vectors, so the ship at (5,1) has coordinates 3.8 and −3.4, each read by one dot product. Right: with the skewed axes (1,0) and (1,1) the same list has true coordinates 4 and 1, but its dot product with the second axis is 6, because the first axis leans into the second. Status: re-simulated.

In an orthonormal basis the coordinates of a list are its dot products with the basis lists [S:25:137]. Write the basis lists as $q_{1}, \dots, q_{n}$ and the list as $x$ . Then the coordinates are the numbers $x \cdot q_{1}$ through $x \cdot q_{n}$ , and no system of equations needs solving [S:21:106].

x = (x \cdot q_{1}) q_{1} + (x \cdot q_{2}) q_{2} + \dots + (x \cdot q_{n}) q_{n}

The identity says that a list is rebuilt by taking its dot product with each axis and laying that many copies of the axis end to end. Each axis is invisible to the others, so the dot product with one axis contains that axis's coefficient and nothing else. Written as a table, the axes become the columns of a matrix, a rectangular array whose columns are those lists.

The transpose of a matrix is the same table flipped about its diagonal, so its columns become rows, and multiplying by the transpose performs all the dot products with the columns in one multiplication. In a general basis that multiplication still runs, but it no longer returns the coefficients. Each dot product picks up contributions from every axis, so the $n$ numbers are $n$ equations in the $n$ unknown coefficients, and you have to solve them.

Worked example. The chart axes are $q_{1} = (3/5, 4/5)$ and $q_{2} = (- 4/5, 3/5)$ , which have length 1 and dot product 0. Against the first axis the ship's dot product is 3 + 0.8 = 3.8, and against the second it is −4 + 0.6 = −3.4. Rebuilding gives 3.8 × (0.6, 0.8) − 3.4 × (−0.8, 0.6), which is (2.28 + 2.72, 3.04 − 2.04) = (5, 1). The two coordinates are 3.8 km along the first axis and 3.4 km backwards along the second, and the rebuild returns the ship exactly.

Lengths survive the change of chart in a clean way. The squared length of the ship's position is 26, because 25 + 1 = 26, and the squared coordinates add to 14.44 + 11.56 = 26 [S:24:214]. Squaring each coordinate and adding gives the squared length, with no cross terms to carry. That is Parseval's identity, and it holds for every orthonormal basis, in any number of dimensions. In a general basis the same sum needs extra terms mixing pairs of coordinates, and those terms are exactly the bookkeeping that an orthonormal basis removes.

Stop and try

use the same two axes to write the list (1, 7). Compute the two coordinates, then square each and add. What does the total tell you about the length of (1, 7)?

Key idea

In an orthonormal basis each coordinate is one dot product, so finding coordinates costs no linear solve. Perpendicular unit axes do not leak into one another, and squared lengths add without cross terms. That is why formulas are written in orthonormal coordinates.

The standard axes are the simplest orthonormal basis, with $q_{1} = (1, 0)$ and $q_{2} = (0, 1)$ , and the formula returns the entries of the list unchanged, which is a consistency check worth doing once. A basis that is not perpendicular fails the same test in a visible way.

Take the axes (1, 0) and (1, 1) with the list (5, 1): the true coordinates are 4 and 1, because 4 × (1, 0) + 1 × (1, 1) = (5, 1). The raw dot product with the second axis is 6, and the arithmetic shows the two separate defects: $6 = 1 \times 2 + 4 \times 1$ . The first term is the coefficient times the squared length of that axis and the second is the leak from the first axis. Dividing by the squared length, 2, removes the length defect and leaves 3. Only the remaining excess 3 − 1 = 2 measures how far the first axis leans into the second, and the two defects have to be corrected in that order.

Most bases you meet are not orthonormal, and turning one into an orthonormal basis takes work. The loop that does it is Gram–Schmidt: divide the first axis by its length, then subtract from each later axis its projection onto the axes already built, and divide what remains by its length. Applied to the columns of a data table $A$ , the loop produces an orthonormal basis $Q$ for the same column space. It also produces a triangular table $R$ holding the coefficients that express each original column in that basis, so that $A = QR$ .

That is the QR factorisation: Q for the orthonormal axes, R for the triangular bookkeeping, and it exists whenever the columns are independent. Least squares then becomes $R c = Q^{⊤} b$ , which back-substitution solves one coefficient at a time, from the last row upward. The Deepen carries out the construction.

The loop can lose accuracy when two axes nearly coincide, and the loss shows up in a number that measures the problem rather than the algorithm. The condition number of a table is the ratio of the largest stretch it applies to a list to the smallest. A large value means that a small change in the data can move the answer a long way. Near-coincident columns make that number large, and forming the Gram matrix squares it [S:28:109], a consequence Section 8.4 states in full.

A QR route avoids the squaring, at the price of more arithmetic than forming $A^{⊤} A$ , and Solomon presents it as the alternative for poorly conditioned problems [S:28:110]. Section 8.4 takes the normal equations, the route through $A^{⊤} A$ , and turns them into the least-squares formula; the QR route above is the alternative when the columns are nearly dependent.

L2Deepen

The formal statement, the derivation and the edge cases. It refines L1 and never contradicts it.

The formal statement. Let $a_{1}, \dots, a_{n}$ be linearly independent lists of length $m$ . Gram–Schmidt produces lists $q_{1}, \dots, q_{n}$ that have unit length and are mutually orthogonal, with the same span at every step. Write $r_{j k} = q_{j} \cdot a_{k}$ and $r_{k k} = ∥ a_{k} - \sum_{j < k} r_{j k} q_{j} ∥$ . Then $A = QR$ , where the columns of $Q$ are the $q_{k}$ and $R$ is upper triangular with those entries. The hypothesis is linear independence: it is what keeps every $r_{k k}$ away from zero.

The derivation. Take $a_{1}$ and divide by its length. That gives $q_{1}$ , and $a_{1} = r_{11} q_{1}$ . Suppose $q_{1}, \dots, q_{k - 1}$ are built. Set

$v_{k} = a_{k} - j < k \sum (a_{k} \cdot q_{j}) q_{j} .$

For each $j < k$ , the dot product of $v_{k}$ with $q_{j}$ is $a_{k} \cdot q_{j} - a_{k} \cdot q_{j} = 0$ , because every other term vanishes by orthogonality. So $v_{k}$ is orthogonal to all earlier axes, and $q_{k} = v_{k} /∥ v_{k} ∥$ is the next unit axis. Rearranging gives $a_{k} = \sum_{j < k} r_{j k} q_{j} + r_{k k} q_{k}$ , which is column $k$ of $A = QR$ . The induction covers every column, so the factorisation follows. The omitted step is the stable way to compute it: the algebra above is exact and the arithmetic is not.

The edge case. Suppose two columns are nearly parallel. Then $v_{k}$ is the difference of two nearly equal lists, and its computed length is dominated by rounding. The computed $q_{k}$ can fail to be orthogonal to the earlier axes by far more than the machine precision suggests. L1 said the construction "can lose accuracy when two axes nearly coincide" [S:28:110]. Householder reflections avoid the subtraction of nearly equal lists altogether and give the same factorisation through orthogonal row operations [S:28:117]. A reorthogonalisation pass is the cheaper repair when you keep the Gram–Schmidt loop.

What L1 hid. L1 required axes that are perpendicular and one unit long, then read each coordinate as a single dot product. With perpendicular axes alone, the coordinate is the dot product divided by the squared length of the axis [S:25:137]. The unit-length condition is doing work, and dropping it changes the formula. L1 also said the shortcut needs "no system of equations". That is exact in exact arithmetic. In floating point the dot products carry rounding, and the quality of the answer depends on how orthogonal the computed axes really are. The simplification was safe because the worked example used axes that were exactly orthonormal, so both claims held to the last digit.

Where the analogy breaks. A rotated chart is a change of description that leaves lengths and angles alone. Gram–Schmidt is not a rotation: it builds new axes out of the old columns, in the order you supply them. Reorder the columns and you get a different $Q$ and a different $R$ , though the same column space. A chart has no preferred first axis; the algorithm does.

L3Research

Full rigour, the original sources and the open questions. This is the research frontier.

Original sources. Treil gives the coordinate formula for an orthogonal basis and states Parseval's identity in the exercises that follow it [S:25:137]. Boyd states the expansion formula for an orthonormal basis and shows how the coefficients come out as dot products [S:21:106]. Solomon develops Gram–Schmidt, shows why the numerical version is unstable, and introduces Householder's reflection of 1958 as the stable route to the same factorisation [S:28:117].

An open question. The field has not settled how much loss of orthogonality is acceptable in a computed basis. The answer depends on what the basis feeds: a least-squares solve, an eigenvalue iteration, or a downstream model. A bound on the departure of the computed $Q^{T} Q$ from the identity does not by itself bound the error in the quantity you care about. Evidence that would settle a case: a task-level error bound derived from the measured orthogonality defect, tested on matrices across a range of condition numbers, rather than the backward error of the factorisation alone.

A frontier method. Randomized range finders compute an approximate basis for the column space by multiplying the matrix by a random projection and then taking a QR factorisation of the smaller result [S:52:63]. Randomized algorithms now exist for the pivoted QR decomposition as well [S:52:62]. Both depend on this section's object at the last step: the output is an orthonormal basis, and the quality of the approximation is measured in the coordinates that basis supplies.

What to read next. Read Solomon's chapter 5 [S:28:109] for the normal equations, Gram–Schmidt and Householder in one place, Treil's section on orthogonal bases [S:25:137] for the formulas and Parseval, and Boyd's chapter 5 [S:21:106] for the expansion identity. Then read Brunton and Kutz's section 1.8 [S:52:62] for how a randomized QR sits inside a modern matrix decomposition, and [S:28:117] for the reflection that replaced the unstable loop.