Projections

The idea in one sentence

When a target is out of reach, the closest thing you can reach is its shadow, and you can find that shadow because the miss it leaves behind stands at a right angle to everything reachable.

A target with no answer

Five rows of data, two knobs to turn, and no setting of the knobs that works. That is the measurements table from Lecture 8: each row carries a weight $w$ and a height $h$, and each row wants a target number built from them with the same two coefficients $u_1$ and $u_2$:

Two knobs, five demands. Take the first two rows on their own and they pin the knobs down: $5u_1 + 2u_2 = 3$ and $3u_1 + 4u_2 = 8$ give $u_1 = -2/7$ and $u_2 = 31/14$. Now look at row three. It wants $2u_1 + 8u_2 = 7$, and those knob settings deliver $120/7$, which is about $17.1$. The first two rows and the third row want different things, and no setting of two knobs can satisfy all five.

Write the two columns side by side as a matrix $A$, and the target as $b$. The reachable outputs are the combinations $Ax$, which is the column space $C(A)$: a flat two-dimensional plane sitting inside five-dimensional space. The target $b$ is not on that plane. The code at the end of this chapter checks it in one line: $A$ has rank 2, and $A$ with $b$ glued on as a third column has rank 3, so $b$ adds a direction the columns never had.

That is the ordinary situation. Real data almost never lands in the column space. Lecture 10 showed which directions are perpendicular to which. This chapter uses that to answer the question the failure raises: if you cannot reach $b$, what is the closest point you can reach?

Start with one column

Shrink the problem until it cannot get smaller. One column $a$, one target $b$. Every reachable point is a multiple $xa$, so the reachable set is the line through $a$. Which multiple sits closest to $b$?

You already know the answer as a picture. Drop a perpendicular from $b$ onto the line, and the foot of that perpendicular is the nearest point. Every other point on the line is the far corner of a right triangle whose short side is that perpendicular, so it has to be farther away.

Now turn the picture into algebra. Call the closest point $p$ and write it as $p = \hat{x}a$, where $\hat{x}$ is the multiple you are looking for. Call the miss $e = b - p$. The picture says $e$ is perpendicular to $a$, and perpendicular means the dot product is zero:

That is the whole derivation. Two numbers go in, $a\T b$ and $a\T a$, and the closest point comes out. With the numbers in the figure, $a = (5, 1)$ and $b = (2, 3)$, so $a\T b = 13$ and $a\T a = 26$, giving $\hat{x} = 1/2$ and $p = (2.5, 0.5)$. The miss is $e = (-0.5, 2.5)$, and $a\T e = -2.5 + 2.5 = 0$, exactly as promised.

One question is still open. The picture asserted that the foot of the perpendicular is the closest point, and the algebra only ever used the right angle. Take any other point $q$ on the line. Then $b$, $p$, and $q$ sit at the corners of a right triangle whose right angle is at $p$, so $\norm{b - q}^2 = \norm{b - p}^2 + \norm{p - q}^2$. The last term is never negative, and it is zero only when $q$ is $p$. No other point on the line can tie.

Two cases you can check in your head

A formula you cannot sanity check is a formula you do not own. There are two cases where you already know the answer.

If $b$ lies along $a$, say $b = a$, then $a\T b = a\T a$ and $\hat{x} = 1$, so $p = b$. The shadow of something already on the line is itself. If $b$ is perpendicular to $a$, then $a\T b = 0$ and $\hat{x} = 0$, so $p = 0$. Nothing of $b$ lies along the line, so the closest reachable point is the origin.

Look at the second picture again next to the first. The vector $(-1, 5)$ is perpendicular to $(5, 1)$, and the miss computed a moment ago was $e = (-0.5, 2.5)$, which is exactly half of $(-1, 5)$. So $b = (2, 3)$ splits as $\tfrac12 (5,1) + \tfrac12 (-1,5)$: the part the line can see, and the part it is blind to. Projection is that split.

The matrix that casts the shadow

The formula for $p$ hides a matrix. In $p = a \left( a\T b / a\T a \right)$ the quantity $a\T b$ is a single number, so it can be moved around freely. Slide it to the right of $a$ and group differently:

Read the two products carefully, because they look alike and are nothing alike. A row times a column, $a\T a$, is one number. A column times a row, $a a\T$, is a full square matrix.

For $a = (5, 1)$ the outer product $a a\T$ is the matrix with rows $(25, 5)$ and $(5, 1)$, and $P$ is that matrix divided by 26. Every row of $a a\T$ is a multiple of the same row $(5, 1)$, so all rows sit on one line and the rank is 1. That has to be true: $P$ maps everything onto a one-dimensional line, so its column space is that line.

Three properties fall out.

$P$ is symmetric, since $(a a\T)\T = a a\T$.

$P$ applied twice is $P$ applied once:

The middle step is the whole trick: the two inner factors $a\T$ and $a$ meet and collapse into the number $a\T a$, which then cancels one copy from the denominator. The geometry says the same thing. $Pb$ already sits on the line, and the shadow of a point on the line is itself, so pressing the button again changes nothing.

And the trace of $P$ is $a\T a / a\T a = 1$, which matches the rank. For projections that is a rule: the trace counts how many dimensions survive.

The simplest projection matrix in the world

Here is a projection you have used without calling it one. Take a picture six pixels wide, two rows of three, and stack the pixels into a vector in $\R^6$. Now redact the middle column by multiplying with a diagonal matrix of ones and zeros.

That diagonal matrix $P$ is symmetric because diagonal matrices are. It satisfies $P^2 = P$ because $1^2 = 1$ and $0^2 = 0$, and because redacting an already redacted picture does nothing new. Its rank is 4, its trace is 4, and its column space is the four-dimensional set of pictures whose middle column is blank. It is a projection onto that subspace, and it is not rank one.

So rank one came from projecting onto a line, and the real definition of a projection matrix is the pair of properties: $P^2 = P$ and $P = P\T$. The matrix $I - P$ is a projection too, the one that keeps exactly what $P$ throws away. In this picture $I - P$ keeps the middle column and blanks the rest.

Many columns at once

Back to the real problem. Now $A$ has $n$ independent columns instead of one, so the reachable set grows from a single line to the whole column space. Everything else survives the upgrade.

Write the closest reachable point as $p = A\hat{x}$, where $\hat{x}$ now holds one coefficient per column, and write the miss as $e = b - A\hat{x}$. Closest still means the miss is perpendicular to everything reachable, and everything reachable is built from the columns, so it is enough for $e$ to be perpendicular to each column separately. Stacking those conditions is exactly what $A\T$ does, since the rows of $A\T$ are the columns of $A$:

One rearrangement finishes it:

The picture uses columns $a_1 = (1, 1, 0)$ and $a_2 = (0, 1, 1)$ and target $b = (2, 3, 4)$. Then $A\T A$ has rows $(2, 1)$ and $(1, 2)$, and $A\T b = (5, 7)$, so $\hat{x} = (1, 3)$ and $p = 1 a_1 + 3 a_2 = (1, 4, 3)$. The miss is $e = (1, -1, 1)$, and it dots to zero against both columns: $1 - 1 + 0 = 0$ and $0 - 1 + 1 = 0$.

Two checks tell you the formula is the right one.

Set $n = 1$, so $A$ is the single column $a$. Then $A\T A$ is the number $a\T a$, its inverse is division, and $P = a (a\T a)^{-1} a\T = a a\T / a\T a$. The subspace formula contains the line formula.

Now let $A$ be square and invertible. Then every $b$ is reachable and the projection should do nothing, and indeed $P = A A^{-1} (A\T)^{-1} A\T = I$.

Why you are allowed to invert that

The formula leans on $(A\T A)^{-1}$, so it owes you a reason. Here it is, in two directions.

If $Ax = 0$ then $A\T A x = A\T 0 = 0$. That direction is free.

Going back, suppose $A\T A x = 0$. Multiply on the left by $x\T$:

A squared length is zero only for the zero vector, so $Ax = 0$. The two null spaces are the same set: $N(A\T A) = N(A)$.

Now the conclusion writes itself. If the columns of $A$ are independent, then $N(A)$ holds only the zero vector, so $N(A\T A)$ holds only the zero vector, and $A\T A$ is a square matrix whose null space is trivial. That is exactly what invertible means. Independent columns are the only condition the projection formula ever needs.

One vector, two pieces

Put the pieces together. Any $b$ splits as $b = p + e$, where $p = Pb$ lives in $C(A)$ and $e$ satisfies $A\T e = 0$, which is the definition of the left null space $N(A\T)$. Lecture 10 proved those two subspaces are perpendicular complements inside $\R^m$, so every $b$ splits this way, and only this way.

Since $p$ and $e$ meet at a right angle, the lengths obey Pythagoras: $\norm{b}^2 = \norm{p}^2 + \norm{e}^2$. For the plane above that reads $29 = 26 + 3$.

The two machines are $P$ and $I - P$. For that same plane, $I - P$ turns out to be the matrix $n n\T / n\T n$ with $n = (1, -1, 1)$, the direction perpendicular to the plane. That is the rank-one formula from the start of the chapter, doing the complementary job. Projecting onto a plane in $\R^3$ and projecting onto the line perpendicular to it are the same calculation seen from two sides.

import numpy as np

a = np.array([5., 1.])
b = np.array([2., 3.])

x_hat = (a @ b) / (a @ a)
p = x_hat * a
e = b - p
print(x_hat)                   # 0.5
print(p)                       # [2.5 0.5]
print(e)                       # [-0.5  2.5]
print(a @ e)                   # 0.0        the miss is perpendicular to a
print(b @ b, p @ p + e @ e)    # 13.0 13.0  Pythagoras

P = np.outer(a, a) / (a @ a)   # a column times a row, over a number
print(bool(abs(P @ b - p).max() < 1e-6))    # True   P finds the same shadow
print(bool(abs(P @ P - P).max() < 1e-6))    # True   P applied twice is P
print(bool(abs(P - P.T).max() < 1e-6))      # True   P is symmetric
print(round(float(np.trace(P)), 6))     # 1.0    the trace counts what is kept
print(int(np.linalg.matrix_rank(P)))        # 1
print(bool(abs(P @ a - a).max() < 1e-6))    # True   a projects to itself
print(bool(abs(P @ np.array([-1., 5.])).max() < 1e-6))  # True

import torch

a = torch.tensor([5., 1.])
b = torch.tensor([2., 3.])

x_hat = (a @ b) / (a @ a)
p = x_hat * a
e = b - p
print(x_hat)                   # 0.5
print(p)                       # [2.5 0.5]
print(e)                       # [-0.5  2.5]
print(a @ e)                   # 0.0        the miss is perpendicular to a
print(b @ b, p @ p + e @ e)    # 13.0 13.0  Pythagoras

P = torch.outer(a, a) / (a @ a)   # a column times a row, over a number
print(bool(abs(P @ b - p).max() < 1e-6))    # True   P finds the same shadow
print(bool(abs(P @ P - P).max() < 1e-6))    # True   P applied twice is P
print(bool(abs(P - P.T).max() < 1e-6))      # True   P is symmetric
print(round(float(torch.trace(P)), 6))  # 1.0    the trace counts what is kept
print(int(torch.linalg.matrix_rank(P)))     # 1
print(bool(abs(P @ a - a).max() < 1e-6))    # True   a projects to itself
print(bool(abs(P @ torch.tensor([-1., 5.])).max() < 1e-6))  # True

import jax.numpy as jnp

a = jnp.array([5., 1.])
b = jnp.array([2., 3.])

x_hat = (a @ b) / (a @ a)
p = x_hat * a
e = b - p
print(x_hat)                   # 0.5
print(p)                       # [2.5 0.5]
print(e)                       # [-0.5  2.5]
print(a @ e)                   # 0.0        the miss is perpendicular to a
print(b @ b, p @ p + e @ e)    # 13.0 13.0  Pythagoras

P = jnp.outer(a, a) / (a @ a)  # a column times a row, over a number
print(bool(abs(P @ b - p).max() < 1e-6))    # True   P finds the same shadow
print(bool(abs(P @ P - P).max() < 1e-6))    # True   P applied twice is P
print(bool(abs(P - P.T).max() < 1e-6))      # True   P is symmetric
print(round(float(jnp.trace(P)), 6))    # 1.0    the trace counts what is kept
print(int(jnp.linalg.matrix_rank(P)))       # 1
print(bool(abs(P @ a - a).max() < 1e-6))    # True   a projects to itself
print(bool(abs(P @ jnp.array([-1., 5.])).max() < 1e-6))  # True

import tensorflow as tf

a = tf.constant([5., 1.])
b = tf.constant([2., 3.])
dot = lambda u, v: tf.tensordot(u, v, 1)
mv = tf.linalg.matvec
gap = lambda z: float(tf.reduce_max(tf.abs(z)))

x_hat = dot(a, b) / dot(a, a)
p = x_hat * a
e = b - p
print(x_hat)                   # 0.5
print(p)                       # [2.5 0.5]
print(e)                       # [-0.5  2.5]
print(dot(a, e))               # 0.0        the miss is perpendicular to a
print(dot(b, b), dot(p, p) + dot(e, e))     # 13.0 13.0  Pythagoras

P = tf.einsum('i,j->ij', a, a) / dot(a, a)  # a column times a row
print(gap(mv(P, b) - p) < 1e-6)             # True   P finds the same shadow
print(gap(P @ P - P) < 1e-6)                # True   P applied twice is P
print(gap(P - tf.transpose(P)) < 1e-6)      # True   P is symmetric
print(round(float(tf.linalg.trace(P)), 6))  # 1.0   the trace counts
print(int(tf.linalg.matrix_rank(P)))        # 1
print(gap(mv(P, a) - a) < 1e-6)             # True   a projects to itself
print(gap(mv(P, tf.constant([-1., 5.]))) < 1e-6)   # True  blind direction

The second block moves up to a subspace. It confirms that the measurements table really is unreachable, then projects onto a plane and checks every claim this chapter made about $P$.

import numpy as np

w = np.array([5., 3., 2., 1., 4.])
h = np.array([2., 4., 8., 5., 7.])
target = np.array([3., 8., 7., 9., 2.])
A = np.column_stack([w, h])
print(int(np.linalg.matrix_rank(A)))                              # 2
print(int(np.linalg.matrix_rank(np.column_stack([A, target]))))   # 3

B = np.array([[1., 0.], [1., 1.], [0., 1.]])   # columns (1,1,0) and (0,1,1)
b = np.array([2., 3., 4.])
print(B.T @ B)                      # [[2. 1.] [1. 2.]]
x_hat = np.linalg.solve(B.T @ B, B.T @ b)
print(x_hat)                        # [1. 3.]
p = B @ x_hat
e = b - p
print(p)                            # [1. 4. 3.]
print(e)                            # [ 1. -1.  1.]
print(B.T @ e)                      # [0. 0.]  perpendicular to both columns
print(b @ b, p @ p + e @ e)         # 29.0 29.0  Pythagoras again

P = B @ np.linalg.inv(B.T @ B) @ B.T
print(bool(abs(P @ b - p).max() < 1e-6))    # True
print(bool(abs(P @ P - P).max() < 1e-6))    # True   still idempotent
print(bool(abs(P - P.T).max() < 1e-6))      # True   still symmetric
print(round(float(np.trace(P)), 6))     # 2.0    two dimensions kept this time
n = np.array([1., -1., 1.])         # the direction perpendicular to the plane
print(bool(abs(np.eye(3) - P - np.outer(n, n) / (n @ n)).max() < 1e-6))

import torch

w = torch.tensor([5., 3., 2., 1., 4.])
h = torch.tensor([2., 4., 8., 5., 7.])
target = torch.tensor([3., 8., 7., 9., 2.])
A = torch.stack([w, h], dim=1)
print(int(torch.linalg.matrix_rank(A)))                                    # 2
print(int(torch.linalg.matrix_rank(torch.cat([A, target[:, None]], 1))))   # 3

B = torch.tensor([[1., 0.], [1., 1.], [0., 1.]])
b = torch.tensor([2., 3., 4.])
print(B.T @ B)                      # [[2. 1.] [1. 2.]]
x_hat = torch.linalg.solve(B.T @ B, B.T @ b)
print(x_hat)                        # [1. 3.]
p = B @ x_hat
e = b - p
print(p)                            # [1. 4. 3.]
print(e)                            # [ 1. -1.  1.]
print(B.T @ e)                      # [0. 0.]  perpendicular to both columns
print(b @ b, p @ p + e @ e)         # 29.0 29.0  Pythagoras again

P = B @ torch.linalg.inv(B.T @ B) @ B.T
print(bool(abs(P @ b - p).max() < 1e-6))    # True
print(bool(abs(P @ P - P).max() < 1e-6))    # True   still idempotent
print(bool(abs(P - P.T).max() < 1e-6))      # True   still symmetric
print(round(float(torch.trace(P)), 6))  # 2.0    two dimensions kept this time
n = torch.tensor([1., -1., 1.])     # the direction perpendicular to the plane
print(bool(abs(torch.eye(3) - P - torch.outer(n, n) / (n @ n)).max() < 1e-6))

import jax.numpy as jnp

w = jnp.array([5., 3., 2., 1., 4.])
h = jnp.array([2., 4., 8., 5., 7.])
target = jnp.array([3., 8., 7., 9., 2.])
A = jnp.column_stack([w, h])
print(int(jnp.linalg.matrix_rank(A)))                              # 2
print(int(jnp.linalg.matrix_rank(jnp.column_stack([A, target]))))  # 3

B = jnp.array([[1., 0.], [1., 1.], [0., 1.]])
b = jnp.array([2., 3., 4.])
print(B.T @ B)                      # [[2. 1.] [1. 2.]]
x_hat = jnp.linalg.solve(B.T @ B, B.T @ b)
print(x_hat)                        # [1. 3.]
p = B @ x_hat
e = b - p
print(p)                            # [1. 4. 3.]
print(e)                            # [ 1. -1.  1.]
print(B.T @ e)                      # [0. 0.]  perpendicular to both columns
print(b @ b, p @ p + e @ e)         # 29.0 29.0  Pythagoras again

P = B @ jnp.linalg.inv(B.T @ B) @ B.T
print(bool(abs(P @ b - p).max() < 1e-6))    # True
print(bool(abs(P @ P - P).max() < 1e-6))    # True   still idempotent
print(bool(abs(P - P.T).max() < 1e-6))      # True   still symmetric
print(round(float(jnp.trace(P)), 6))    # 2.0    two dimensions kept this time
n = jnp.array([1., -1., 1.])        # the direction perpendicular to the plane
print(bool(abs(jnp.eye(3) - P - jnp.outer(n, n) / (n @ n)).max() < 1e-6))

import tensorflow as tf

w = tf.constant([5., 3., 2., 1., 4.])
h = tf.constant([2., 4., 8., 5., 7.])
target = tf.constant([3., 8., 7., 9., 2.])
mv, dot = tf.linalg.matvec, lambda u, v: tf.tensordot(u, v, 1)
gap = lambda z: float(tf.reduce_max(tf.abs(z)))
A = tf.stack([w, h], axis=1)
print(int(tf.linalg.matrix_rank(A)))                                       # 2
print(int(tf.linalg.matrix_rank(tf.concat([A, target[:, None]], 1))))  # 3

B = tf.constant([[1., 0.], [1., 1.], [0., 1.]])
b = tf.constant([2., 3., 4.])
Bt = tf.transpose(B)
print(Bt @ B)                       # [[2. 1.] [1. 2.]]
x_hat = tf.linalg.solve(Bt @ B, mv(Bt, b)[:, None])[:, 0]
print(x_hat)                        # [1. 3.]
p = mv(B, x_hat)
e = b - p
print(p)                            # [1. 4. 3.]
print(e)                            # [ 1. -1.  1.]
print(mv(Bt, e))                    # [0. 0.]  perpendicular to both columns
print(dot(b, b), dot(p, p) + dot(e, e))     # 29.0 29.0  Pythagoras again

P = B @ tf.linalg.inv(Bt @ B) @ Bt
print(gap(mv(P, b) - p) < 1e-6)             # True
print(gap(P @ P - P) < 1e-6)                # True   still idempotent
print(gap(P - tf.transpose(P)) < 1e-6)      # True   still symmetric
print(round(float(tf.linalg.trace(P)), 6))  # 2.0   two dimensions kept
n = tf.constant([1., -1., 1.])      # the direction perpendicular to the plane
print(gap(tf.eye(3) - P - tf.einsum('i,j->ij', n, n) / dot(n, n)) < 1e-6)

Deleting a direction on purpose

The redaction matrix threw away coordinates. A direction a model has learned is almost never one of your coordinates, and you can throw it away anyway.

Take the word embeddings inside a language model and train a small linear classifier to read some property off them. The classifier’s weight vector is a direction $a$ in embedding space, and how much of the property a word carries is how far that word reaches along $a$. To remove the property, project every embedding onto the perpendicular complement of $a$. That is $I - P$ with $P = a a\T / a\T a$, the same rank-one matrix from the start of this chapter, now used for what it destroys instead of what it keeps.

The readings are not shrunk. They are exactly zero, because $a\T(I - P) = a\T - a\T = 0$, and a classifier whose weight vector is $a$ now reads the same number off every word in the language. What the words keep is everything perpendicular to $a$, which in a real model is all but one of several hundred directions.

One pass rarely finishes the job. Train a fresh classifier on what survived and it often finds another direction that still carries the property, so the method is run in a loop: learn a direction, project it away, learn the next one. Each pass subtracts one more rank-one $P$, and each pass costs the model one more direction of the space it had to work in.

Where this is going

You now have the closest reachable point and the matrix that finds it. What you have not been told is what the coefficients $\hat{x}$ mean. For the measurements table they are the numbers a model would learn, and the miss $e$ is the error it has to live with. Naming that properly, and fitting a line to real data with it, is the next chapter: least squares.

One warning to carry forward. The formula $P = A(A\T A)^{-1}A\T$ is how you think about projection. A good program computes it another way, because forming $A\T A$ squares how badly conditioned the problem is. Lecture 13 replaces the inverse with orthonormal columns and turns the whole calculation into back substitution.