Orthonormal Bases, Gram-Schmidt & QR

The idea in one sentence

When the columns of a matrix are perpendicular and have length 1, the expensive part of every formula in this course disappears: the inverse becomes the transpose.

The wish

Lecture 12 ended at a formula that works and costs too much. To fit a model to data that no model fits exactly, you solve

Everything about that line is fine except one thing. The matrix $A\T A$ has to be built and then undone. Its entry in row $i$, column $j$ is the dot product of column $i$ of $A$ with column $j$. So $A\T A$ is a table of how much each column of $A$ overlaps every other column. A good name for it is the coupling matrix: the columns are tangled together, and solving means untangling them.

Now make a wish. What if $A\T A$ were the identity? Then $\hat{x} = A\T b$, and there is nothing left to solve. One matrix multiplication and you are done.

The wish is really a request about the columns of $A$. Written out, $A\T A = I$ says two things at once: every column has length 1, and every pair of different columns has dot product 0. Perpendicular unit vectors. Nothing more.

So the rest of this lecture builds columns like that, and then spends them.

Q, and the one identity that matters

Vectors $q_1, \ldots, q_n$ are orthonormal when

Put them in the columns of a matrix and the matrix gets the letter $Q$. Here is one, and every entry is a third:

Check the first column against itself: $(1 + 4 + 4)/9 = 1$. Check it against the second: $(2 - 4 + 2)/9 = 0$. All nine checks work, and all nine checks are the entries of $Q\T Q$. Row $i$ of $Q\T$ is column $i$ of $Q$, so

That single identity is the whole lecture, and it asks for less than you might expect. $Q$ can be tall and thin. If $Q$ is $5 \times 2$ with two orthonormal columns, $Q\T Q$ is the $2 \times 2$ identity, and it still holds. What fails for a tall $Q$ is the other order: $QQ\T$ is $5 \times 5$ with rank 2, so it cannot be $I_5$. You will meet that matrix again in a few paragraphs, doing useful work.

When $Q$ is square, $Q\T Q = I$ makes $Q\T$ a genuine inverse, so $QQ\T = I$ too and

The expensive step in every formula so far becomes a transpose.

The orthogonal family

You have already dragged some of these. In Lecture 2 you moved four numbers around and watched the plane move with them. Three of those machines were orthogonal, and now you can say why.

A rotation turns the plane by an angle $\theta$:

Its columns are the landing spots of $e_1$ and $e_2$, both of length 1, and perpendicular because they started perpendicular. At $\theta = 30°$ the columns are $(0.866, 0.5)$ and $(-0.5, 0.866)$.

A reflection flips the plane across a line. Across the line through the origin with direction $(0.8, 0.6)$, the matrix is

and $H\T H = I$ holds because $0.28^2 + 0.96^2 = 1$ and $0.28 (0.96) + 0.96(-0.28) = 0$. The next section builds this matrix instead of quoting it.

A permutation matrix shuffles the entries of a vector. It has exactly one 1 in every row and exactly one 1 in every column, and nothing else:

Its columns are the standard basis vectors in a new order, so they are already orthonormal, and every permutation matrix is orthogonal. Reordering a list cannot change its length, which is the same fact said in words.

A mirror built from a shadow

Lecture 11 derived the shadow of one vector on another. For a unit vector $u$ the shadow of $v$ is short:

Now walk from $v$ to $p$, and keep walking the same distance again. You arrive on the far side of the line through the origin perpendicular to $u$, at the same distance from it. That is a reflection, and it is one subtraction:

The matrix $H = I - 2uu\T$ is the Householder reflection. It is symmetric, and applying it twice returns you home, so $H\T H = H^2 = I$: orthogonal, with no work.

Nothing stretches

Here is the property that makes numerical people reach for $Q$ before anything else. Take any vector $x$ and hit it with $Q$. Its length is unchanged:

Two steps, and the middle one is the only fact in play. Angles survive too, because dot products do:

Try it with the $Q$ from earlier. The vector $(1, -2, 3)$ has length $\sqrt{14} = 3.7416573\ldots$, and so does $Q(1, -2, 3)$, to the last digit your computer keeps.

Projection, with the inverse deleted

Lecture 11 built the projection onto the column space of $A$:

Set $A = Q$ and watch the middle collapse. Since $Q\T Q = I$, its inverse is $I$, and

There is nothing to invert. Take the first two columns of the $Q$ above, call that tall matrix $Q_2$, and project $b = (1,1,1)$ onto the plane they span:

Build $A(A\T A)^{-1}A\T$ from the same two columns before they were orthonormalized and you get the identical matrix, after inverting a $2 \times 2$. Same answer, more work.

There is a second way to read $QQ\T b$. Multiply the other way round and it is a sum:

Each term is one shadow. So the projection is the sum of the shadows of $b$ on each $q$ separately, with no correction between them. That is what “no coupling” means, and it is exactly what fails for columns that are not perpendicular: their shadows overlap, and the overlap is $A\T A$.

Wishful thinking: suppose $A = QR$

Real data does not arrive with orthonormal columns. So make a wish and follow it. Suppose the matrix $A$ could be written as

with $Q$ carrying orthonormal columns and $R$ carrying whatever is left over. Put that into the normal equations and turn the handle:

The middle of the left side is $Q\T Q$, which is $I$, so it goes:

Now $R\T$ sits in front of both sides. If $A$ has independent columns then $R$ is invertible, so $R\T$ is too, and it cancels:

Read what happened. The coupling matrix $A\T A$ is gone. In its place is one matrix multiplication $Q\T b$ and a system with a triangular $R$, which back substitution solves top down in $n^2/2$ steps. No inverse is ever formed.

Everything now rests on one question. Does $A = QR$ exist, and can you compute it?

Gram-Schmidt, on real numbers

Yes, and the recipe is one idea repeated: take the next column, subtract its shadows on everything you have built so far, and divide by what is left.

Take three columns and run the whole thing. Nothing here is symbolic:

Step 1. There is nothing to subtract from the first vector, so only normalize it. $\norm{a} = \sqrt{1 + 4 + 4} = 3$, and

Record the length that came out: $r_{11} = 3$.

Step 2. Take $b$ and remove its shadow on $q_1$. The shadow coefficient is

Subtract:

Check that this is perpendicular to $q_1$: $\tfrac{1}{3}\left(\tfrac{4}{3} - \tfrac{8}{3} + \tfrac{4}{3}\right) = 0$. It is, by construction, because subtracting the shadow is exactly what kills the overlap. Its length is $\tfrac{2}{3}\sqrt{4 + 4 + 1} = 2$, so $r_{22} = 2$ and

Drag the two starting vectors and run the steps yourself. What happens when they are nearly parallel is the whole story of the last section of this lecture:

Step 3. Now $c$, with two shadows to remove:

Its length is $\tfrac{1}{9}\sqrt{4 + 1 + 4} = \tfrac{1}{3}$, so $r_{33} = \tfrac{1}{3}$ and $q_3 = \tfrac{1}{3}(2, 1, -2)$.

Now collect the numbers you wrote down along the way. Each original vector is built from the $q$‘s using exactly the coefficients Gram-Schmidt produced:

Stack those coefficient lists as columns and you have $R$. Stack the $q$‘s as columns and you have $Q$:

Multiply it out to check. Column 2 of the product is $2q_1 + 2q_2 = \tfrac{2}{3}(1,2,2) + \tfrac{2}{3}(2,-2,1) = (2, 0, 2)$, which is $b$. The wish came true, and it came true for the plain reason that every step of the process was a subtraction and a division.

Why is $R$ triangular? Because $q_2$ was born after $a$ was already finished. The vector $a$ is a multiple of $q_1$ alone, so its column of $R$ has one entry and two zeros. In general $R = Q\T A$, and the entry below the diagonal is $q_i\T a_j$ with $i > j$: a later $q$ against an earlier $a$, and those are perpendicular.

import numpy as np

A = np.array([[1., 2., 1.],
              [2., 0., 1.],
              [2., 2., 1.]])         # columns a, b, c

qs, R = [], [[0.0] * 3 for _ in range(3)]
for j in range(3):
    v = A[:, j]
    for i in range(j):
        R[i][j] = float(qs[i] @ v)   # how much of q_i is still in v
        v = v - R[i][j] * qs[i]      # take the shadow away
    R[j][j] = float(np.linalg.norm(v))       # the length that is left
    qs.append(v / R[j][j])
Q = np.stack(qs, axis=1)

print(3 * Q)   # [[ 1.  2.  2.] [ 2. -2.  1.] [ 2.  1. -2.]]
print(R)       # [[3.0, 2.0, 1.6666666666666665],
               #  [0.0, 2.0, 0.33333333333333337],
               #  [0.0, 0.0, 0.3333333333333333]]
print(np.round(Q.T @ Q, 12))              # the identity
print(np.abs(Q @ np.array(R) - A).max())  # 0.0   exact here
x = np.array([1., -2., 3.])
print(np.linalg.norm(Q @ x))   # 3.741657386773942
print(np.linalg.norm(x))       # 3.7416573867739413

import torch

A = torch.tensor([[1., 2., 1.],
                  [2., 0., 1.],
                  [2., 2., 1.]], dtype=torch.float64)   # columns a, b, c

qs, R = [], [[0.0] * 3 for _ in range(3)]
for j in range(3):
    v = A[:, j]
    for i in range(j):
        R[i][j] = float(qs[i] @ v)   # how much of q_i is still in v
        v = v - R[i][j] * qs[i]      # take the shadow away
    R[j][j] = float(torch.linalg.norm(v))    # the length that is left
    qs.append(v / R[j][j])
Q = torch.stack(qs, dim=1)           # stack takes dim, not axis

print(3 * Q)   # [[ 1.  2.  2.] [ 2. -2.  1.] [ 2.  1. -2.]]
print(R)       # [[3.0, 2.0, 1.6666666666666665],
               #  [0.0, 2.0, 0.33333333333333337],
               #  [0.0, 0.0, 0.3333333333333333]]
Rt = torch.tensor(R, dtype=torch.float64)
print(torch.round(Q.T @ Q, decimals=12))   # the identity
print((Q @ Rt - A).abs().max())            # 0.0   exact here
x = torch.tensor([1., -2., 3.], dtype=torch.float64)
print(torch.linalg.norm(Q @ x))   # 3.741657386773942
print(torch.linalg.norm(x))       # 3.7416573867739413

import jax
jax.config.update("jax_enable_x64", True)   # float32 hides the last digits
import jax.numpy as jnp

A = jnp.array([[1., 2., 1.],
               [2., 0., 1.],
               [2., 2., 1.]])          # columns a, b, c

qs, R = [], [[0.0] * 3 for _ in range(3)]   # jax arrays are immutable
for j in range(3):
    v = A[:, j]
    for i in range(j):
        R[i][j] = float(qs[i] @ v)   # how much of q_i is still in v
        v = v - R[i][j] * qs[i]      # take the shadow away
    R[j][j] = float(jnp.linalg.norm(v))      # the length that is left
    qs.append(v / R[j][j])
Q = jnp.stack(qs, axis=1)

print(3 * Q)   # [[ 1.  2.  2.] [ 2. -2.  1.] [ 2.  1. -2.]]
print(R)       # [[3.0, 2.0, 1.6666666666666665],
               #  [0.0, 2.0, 0.33333333333333337],
               #  [0.0, 0.0, 0.3333333333333333]]
print(jnp.round(Q.T @ Q, 12))              # the identity
print(jnp.abs(Q @ jnp.array(R) - A).max()) # 0.0   exact here
x = jnp.array([1., -2., 3.])
print(jnp.linalg.norm(Q @ x))   # 3.741657386773942
print(jnp.linalg.norm(x))       # 3.7416573867739413

import tensorflow as tf

A = tf.constant([[1., 2., 1.],
                 [2., 0., 1.],
                 [2., 2., 1.]], dtype=tf.float64)   # columns a, b, c

qs, R = [], [[0.0] * 3 for _ in range(3)]
for j in range(3):
    v = A[:, j]
    for i in range(j):
        R[i][j] = float(tf.tensordot(qs[i], v, 1))  # tf: no 1-D @ operator
        v = v - R[i][j] * qs[i]                     # take the shadow away
    R[j][j] = float(tf.norm(v))                     # the length that is left
    qs.append(v / R[j][j])
Q = tf.stack(qs, axis=1)

print(3 * Q)   # [[ 1.  2.  2.] [ 2. -2.  1.] [ 2.  1. -2.]]
print(R)       # [[3.0, 2.0, 1.6666666666666665],
               #  [0.0, 2.0, 0.33333333333333337],
               #  [0.0, 0.0, 0.3333333333333333]]
Rt = tf.constant(R, tf.float64)
print(tf.round(tf.transpose(Q) @ Q * 1e12) / 1e12)   # the identity
print(tf.reduce_max(tf.abs(Q @ Rt - A)))             # 0.0   exact here
x = tf.constant([1., -2., 3.], dtype=tf.float64)
print(tf.norm(tf.linalg.matvec(Q, x)))   # 3.741657386773942
print(tf.norm(x))                        # 3.7416573867739413

The same process, done stably

Written as an algorithm, Gram-Schmidt is six lines:

for j = 1 to n
    v = column j of A
    for i = 1 to j-1
        R[i,j] = q_i dot v
        v      = v - R[i,j] * q_i
    R[j,j] = length of v
    q_j    = v / R[j,j]

Look closely at line 4. The dot product is taken against v, the vector as it stands right now, with earlier shadows already removed. The classical version takes it against the original column instead: R[i,j] = q_i dot A[:,j]. In exact arithmetic the two agree, since the pieces already removed were perpendicular to $q_i$ anyway. In floating point they do not, because the first version measures the small leftover directly while the second measures a big vector and hopes the large parts cancel.

The version above is called modified Gram-Schmidt, and the difference is not academic. Run both on a matrix whose columns nearly coincide and check how far $Q\T Q$ drifts from the identity:

import numpy as np

def gram_schmidt(A, classical):
    qs = []
    for j in range(A.shape[1]):
        v = A[:, j].copy()
        for i in range(j):
            r = qs[i] @ (A[:, j] if classical else v)   # the line that differs
            v = v - r * qs[i]
        qs.append(v / np.linalg.norm(v))
    return np.stack(qs, axis=1)

eps = 1e-8
A = np.array([[1., 1., 1.],
              [eps, 0., 0.],
              [0., eps, 0.],
              [0., 0., eps]])         # three columns that almost agree

for name, Q in [("classical", gram_schmidt(A, True)),
                ("modified ", gram_schmidt(A, False)),
                ("library  ", np.linalg.qr(A)[0])]:
    print(name, f"{np.abs(Q.T @ Q - np.eye(3)).max():.3e}")
# classical 5.000e-01     the columns are not perpendicular at all
# modified  7.071e-09     eight digits kept
# library   6.661e-16     Householder mirrors, full precision

Half. The classical version returns columns that are not remotely perpendicular, and it does so silently. Modified Gram-Schmidt keeps eight digits. Libraries do better still by never building $Q$ this way at all: they apply a chain of Householder mirrors, each one an exact orthogonal matrix, and lengths cannot drift when nothing stretches.

Two roads to least squares

Now the payoff, priced honestly. You have two ways to solve a least squares problem.

The first road builds $A\T A$ and solves. The second road factors $A = QR$ and solves $R\hat{x} = Q\T b$ by back substitution. On clean data they agree to the last digit. Take the measurements table from Lecture 12, put a column of ones beside the two feature columns, and both roads return $\hat{x} = (14.855866, -2.079330, -0.541899)$.

The roads part when the columns of $A$ nearly line up, which is what real features do. Here is a matrix built to lean that way, with a knob $d$ that squeezes two columns together:

The condition number of a matrix is the ratio of its largest stretch to its smallest, and it predicts how many digits a solve can lose. For this $A$ it grows like $1/d$. For $A\T A$ it grows like $1/d^2$, because forming $A\T A$ squares every stretch factor. Squaring a condition number means throwing away half your digits before you start.

import numpy as np

def make(d):
    A = np.array([[1., 1.], [1., 1. + d], [1., 1. - d]])
    return A, A @ np.array([1., 1.])   # so the answer is exactly (1, 1)

for d in [1e-2, 1e-6, 1e-8]:
    A, b = make(d)
    s = np.linalg.svd(A, compute_uv=False)
    t = np.linalg.svd(A.T @ A, compute_uv=False)
    print(d, s[0] / s[-1], t[0] / t[-1])
# 0.01  2.45e+02  6.00e+04
# 1e-06 2.45e+06  6.00e+12
# 1e-08 2.45e+08  8.83e+16   squaring the condition costs half the digits

A, b = make(1e-6)
print(np.linalg.solve(A.T @ A, A.T @ b))   # [1.00022207 0.99977793]
Q, R = np.linalg.qr(A)
print(np.linalg.solve(R, Q.T @ b))         # [1. 1.]        exact

A, b = make(1e-8)
print(np.linalg.det(A.T @ A))       # 0.0   AtA is now singular
Q, R = np.linalg.qr(A)
print(np.linalg.solve(R, Q.T @ b))  # [1.00000001 0.99999999]   QR is fine

import torch

def make(d):
    rows = [[1., 1.], [1., 1. + d], [1., 1. - d]]
    A = torch.tensor(rows, dtype=torch.float64)
    return A, A @ torch.ones(2, dtype=torch.float64)

for d in [1e-2, 1e-6, 1e-8]:
    A, b = make(d)
    s = torch.linalg.svdvals(A)      # svdvals here, not svd(...)
    t = torch.linalg.svdvals(A.T @ A)
    print(d, (s[0] / s[-1]).item(), (t[0] / t[-1]).item())
# 0.01  2.45e+02  6.00e+04
# 1e-06 2.45e+06  6.00e+12
# 1e-08 2.45e+08  8.83e+16   squaring the condition costs half the digits

A, b = make(1e-6)
print(torch.linalg.solve(A.T @ A, A.T @ b))   # [1.00022207 0.99977793]
Q, R = torch.linalg.qr(A)
print(torch.linalg.solve(R, Q.T @ b))         # [1. 1.]      exact

A, b = make(1e-8)
print(torch.linalg.det(A.T @ A))     # 0.0   AtA is now singular
Q, R = torch.linalg.qr(A)
print(torch.linalg.solve(R, Q.T @ b))  # [1.00000001 0.99999999]  QR is fine

import jax
jax.config.update("jax_enable_x64", True)
import jax.numpy as jnp

def make(d):
    A = jnp.array([[1., 1.], [1., 1. + d], [1., 1. - d]])
    return A, A @ jnp.array([1., 1.])

for d in [1e-2, 1e-6, 1e-8]:
    A, b = make(d)
    s = jnp.linalg.svd(A, compute_uv=False)
    t = jnp.linalg.svd(A.T @ A, compute_uv=False)
    print(d, s[0] / s[-1], t[0] / t[-1])
# 0.01  2.45e+02  6.00e+04
# 1e-06 2.45e+06  6.00e+12
# 1e-08 2.45e+08  8.83e+16   squaring the condition costs half the digits

A, b = make(1e-6)
print(jnp.linalg.solve(A.T @ A, A.T @ b))   # [1.00022207 0.99977793]
Q, R = jnp.linalg.qr(A)
print(jnp.linalg.solve(R, Q.T @ b))        # [1. 1.]        exact

A, b = make(1e-8)
print(jnp.linalg.det(A.T @ A))      # 0.0   AtA is now singular
Q, R = jnp.linalg.qr(A)
print(jnp.linalg.solve(R, Q.T @ b)) # [1.00000001 0.99999999]   QR is fine

import tensorflow as tf

def make(d):
    A = tf.constant([[1., 1.], [1., 1. + d], [1., 1. - d]], dtype=tf.float64)
    return A, tf.linalg.matvec(A, tf.constant([1., 1.], tf.float64))

for d in [1e-2, 1e-6, 1e-8]:
    A, b = make(d)
    s = tf.linalg.svd(A, compute_uv=False)   # sorted, largest first
    t = tf.linalg.svd(tf.transpose(A) @ A, compute_uv=False)
    print(d, (s[0] / s[-1]).numpy(), (t[0] / t[-1]).numpy())
# 0.01  2.45e+02  6.00e+04
# 1e-06 2.45e+06  6.00e+12
# 1e-08 2.45e+08  8.83e+16   squaring the condition costs half the digits

A, b = make(1e-6)
At = tf.transpose(A)
rhs = tf.reshape(tf.linalg.matvec(At, b), [2, 1])   # tf wants a 2-D rhs
print(tf.linalg.solve(At @ A, rhs))   # [[1.00022207] [0.99977793]]
Q, R = tf.linalg.qr(A)
qb = tf.reshape(tf.linalg.matvec(Q, b, transpose_a=True), [2, 1])
print(tf.linalg.solve(R, qb))         # [[1.] [1.]]        exact

A, b = make(1e-8)
print(tf.linalg.det(tf.transpose(A) @ A))   # 0.0   AtA is singular

At $d = 10^{-6}$ the normal equations lose four digits and QR loses none. At $d = 10^{-8}$ the normal equations stop working: all four entries of $A\T A$ round to the same number 3, its determinant is exactly zero, and the solver refuses. The matrix $A$ itself is untouched. It still has two independent columns and an ordinary condition number. All the damage was done by the step that formed $A\T A$, and QR never takes that step.

Why deep networks are initialised with Q

A deep network multiplies its input by one matrix per layer. Sixty layers is sixty multiplications, and whatever length the signal has at the end is whatever those sixty matrices left of it. The callout above claimed an orthogonal matrix can do neither harm. Here is that claim run for six steps and measured.

Start from $x = (1, 0)$ and multiply repeatedly by three different matrices. The first is the rotation $Q$ with columns $(0.6, 0.8)$ and $(-0.8, 0.6)$. The other two are ordinary matrices with no special structure.

Six steps already spread the two ordinary matrices by a factor of nearly $200$. At sixty layers the same arithmetic gives numbers your machine cannot hold at one end and cannot distinguish from zero at the other. That is the exploding and vanishing signal problem, and it is a fact about $\norm{Ax}$ before it is a fact about training.

This is why orthogonal initialisation exists. Give every layer a $Q$ and the signal arrives at the last layer the same length it left the first, whatever the depth. The gradient runs the same gauntlet backwards through the transposes, and $Q\T$ is orthogonal too, so it is protected by the identical two-line proof.

Where this is going

You now have a factorization that turns any independent columns into perpendicular unit ones, and a reason to want it. Two things follow.

The first is that $Q$ shows up again as the answer to a different question. Lecture 15 asks which directions a matrix leaves alone, and for symmetric matrices the answer is a full orthonormal set: $S = Q\Lambda Q\T$, a rotation, a stretch along the axes, and the rotation undone. Every fact in this lecture is waiting there.

The second is that Gram-Schmidt needs independent columns, and dies when they are not. When two columns coincide exactly, $r_{22}$ comes out zero and you divide by it. The tool that handles rank honestly, and handles nearly-dependent columns better than anything here, is the singular value decomposition, which Part VI reaches by combining this lecture with the next two.