Orthogonality of the Four Subspaces

The idea in one sentence

The four subspaces come in two perpendicular pairs, and each pair covers its whole space with nothing left over, which is why every reachable $b$ comes from one and only one $x$.

The matrix and its four bases

Everything in this chapter happens inside one $3 \times 4$ matrix, the same one Lecture 9 took apart. Its four bases are already in hand:

In $\R^4$: the row space $C(A\T)$ has basis $(1, 0, 2, 0)$ and $(0, 1, 0, 2)$, and the null space $N(A)$ has basis $s_1 = (-2, 0, 1, 0)$ and $s_2 = (0, -2, 0, 1)$.

In $\R^3$: the column space $C(A)$ has basis $(1, 3, 4)$ and $(2, 8, 10)$, and the left null space $N(A\T)$ has basis $(-1, -1, 1)$.

Four sets of vectors, four dimensions that add to 4 and to 3. What the last chapter did not say is how they sit relative to each other. That is one line of algebra away.

The proof is one line

Take any $x$ in the null space. The statement $Ax = 0$ looks like a fact about a matrix. Read it by rows and it is a fact about angles. Each entry of $Ax$ is a row of $A$ dotted with $x$, so $Ax = 0$ says every one of those dot products is zero:

A zero dot product means a right angle. So $x$ is perpendicular to every row of $A$. It is then perpendicular to every combination of the rows, because dotting distributes over sums and scalar multiples. The combinations of the rows are the row space. Done:

Now turn the whole argument around by writing $A\T$ where $A$ was. Take any $y$ in the left null space. Then $A\T y = 0$ says every row of $A\T$, which is to say every column of $A$, is perpendicular to $y$. So $y$ is perpendicular to every combination of the columns:

Check it on the numbers. The left null space vector is $(-1, -1, 1)$, and the two column space basis vectors give $-1 - 3 + 4 = 0$ and $-2 - 8 + 10 = 0$. Both zero, both right angles.

Two subspaces are called orthogonal when every vector in one is perpendicular to every vector in the other. That is a strong condition, and it has an immediate consequence. If some vector $v$ belonged to both, it would have to be perpendicular to itself, so $v \cdot v = 0$, so $v = 0$. Orthogonal subspaces meet at the origin and nowhere else.

The right angle you can draw

Abstract right angles are easy to nod at and hard to feel. Here is one with numbers. Take $x = (1, 1, 1, 1)$ and split it into the row space piece and the null space piece from Lecture 9:

Their dot product is $\tfrac{6}{25} + \tfrac{6}{25} - \tfrac{6}{25} - \tfrac{6}{25} = 0$. So the two pieces of $x$ sit at a right angle, and the three vectors $x_r$, $x_n$ and $x$ form a right triangle inside $\R^4$. That triangle lives in an ordinary flat plane, so you can draw it exactly, and Pythagoras has to hold:

Move $x$ around and the triangle changes shape, but the corner at the origin stays square. Try it.

Nothing perpendicular is missing

Perpendicular is less than it sounds. In $\R^3$ the $x$ axis and the $y$ axis are perpendicular, yet together they reach only a plane. Everything in the $z$ direction is missing. Two subspaces can meet at right angles and still leave most of the space uncovered.

The four subspaces do better than that. Nothing is left over, and the proof runs both ways in two short steps.

If $x$ is in $N(A)$ then $Ax = 0$, so $x$ is perpendicular to every row, so $x$ is perpendicular to the row space. That is the direction already proved. Now the other one. Suppose $x$ is perpendicular to every row of $A$. Then every entry of $Ax$, being a row dotted with $x$, is zero, so $Ax = 0$ and $x$ is in $N(A)$.

Put those together and the null space holds every vector perpendicular to the row space, and nothing else. The name for that is the orthogonal complement, and the dimensions are the receipt:

For our matrix, $2 + 2 = 4$ in $\R^4$ and $2 + 1 = 3$ in $\R^3$. If the null space had come out one dimension short, some direction perpendicular to the row space would have had nowhere to live, and the count would have caught it.

The $\R^3$ side is small enough to draw honestly. The column space is a plane through the origin. The left null space is the single line through $(-1, -1, 1)$. The plane and the line are perpendicular, and a plane plus its normal line is all of $\R^3$.

Every reachable b comes from exactly one x

Here is the payoff.

Take any $b$ in the column space. Since $b$ is a combination of columns, there is some $x$ with $Ax = b$. Split it: $x = x_r + x_n$, and $Ax_n = 0$, so $Ax_r = b$ as well. So $b$ comes from at least one vector in the row space.

Could it come from two? Suppose $x_r$ and $x_r'$ both sit in the row space and both satisfy $Ax_r = Ax_r' = b$. Then $A(x_r - x_r') = 0$, so the difference is in the null space. It is also in the row space, since the row space is closed under subtraction. A vector in both spaces is perpendicular to itself, and the only such vector is zero. So $x_r = x_r'$.

The matrix may crush a whole two-dimensional null space to zero. Underneath the crushing there is still a perfect one-to-one correspondence between two spaces of dimension $r$, and $A$ restricted to the row space is invertible.

Watch it happen. The vectors $x = (1, 1, 1, 1)$ and $x' = (-1, 1, 2, 1)$ differ by the null space vector $s_1$, so both give the same output:

Feed each through the split and both hand back the same row space vector, $x_r = \left(\tfrac{3}{5}, \tfrac{3}{5}, \tfrac{6}{5}, \tfrac{6}{5}\right)$. Two inputs, one output, one row space representative. And that representative is the shortest of all the solutions, since $\norm{x_r} \approx 1.897$ while $\norm{x} = 2$ and $\norm{x'} \approx 2.646$. Pythagoras again: any other solution carries an extra perpendicular leg, and a leg can only make a vector longer.

The four rank cases, seen through the right angle

Lecture 8 sorted every matrix into four cases by comparing the rank $r$ with the shape $m$ and $n$, and gave a worked example of each. The Big Picture says what those four cases really are: which of the two null spaces has collapsed to a point.

Read the panels off the two null spaces. A null space that is only the zero vector means no two inputs collide, so a solution is unique when it exists. A left null space that is only the zero vector means the column space fills $\R^m$, so every $b$ is reachable. Case 1 has both, which is what invertible means. Case 4 has neither, and every unpleasant thing that can happen to a linear system happens there.

The part that survives all four cases is the one-to-one link across the middle. Even in case 4, where $A$ crushes a null space and misses part of $\R^m$, the row space and the column space are still perfectly matched.

The pseudoinverse, defined honestly

A matrix with a nonzero null space has no inverse. There is no undoing the crush: once $s_1$ and $0$ have both been sent to $0$, nothing downstream can tell them apart. That looks like the end of the story. The one-to-one link says otherwise.

Here is what $A^{+}$, the pseudoinverse, does with a vector $b$ in $\R^m$.

1. Split $b$ into the part inside $C(A)$ and the part inside $N(A\T)$. Those two spaces are complements, so the split exists and is unique.
2. Throw away the $N(A\T)$ part. It was never reachable, so no $x$ was ever going to produce it.
3. Undo the one-to-one link on what is left, landing on the single $x_r$ in the row space with $A x_r$ equal to the kept part.

In one line: invert the part that is invertible and ignore the null spaces. Every property of $A^{+}$ follows from those three steps. $A^{+}A$ is the projection onto the row space, because it sends $x$ to $x_r$. $AA^{+}$ is the projection onto the column space, because it sends $b$ to its reachable part. And $AA^{+}A = A$, since $A$ ignores the discarded pieces anyway.

On the running matrix, numpy.linalg.pinv reproduces all of this exactly. Feeding it $b = (9, 33, 42)$ returns $\left(\tfrac{3}{5}, \tfrac{3}{5}, \tfrac{6}{5}, \tfrac{6}{5}\right)$, the row space representative. Feed it the unreachable $b' = (8, 32, 43)$, which is $b$ plus the left null space vector, and it returns the same answer, because step 2 discarded the difference.

import numpy as np
np.set_printoptions(precision=4, suppress=True)

A = np.array([[1., 2., 2., 4.], [3., 8., 6., 16.], [4., 10., 8., 20.]])
row = np.array([[1., 0., 2., 0.], [0., 1., 0., 2.]])     # row space basis
nul = np.array([[-2., 0., 1., 0.], [0., -2., 0., 1.]])   # null space basis
lnul = np.array([[-1., -1., 1.]])                        # left null basis

print(row @ nul.T)          # [[0. 0.] [0. 0.]]   every pair perpendicular
print(A[:, :2].T @ lnul.T)  # [[0.] [0.]]         and on the other side too

x = np.array([1., 1., 1., 1.])
xr = (x @ row[0]) / 5 * row[0] + (x @ row[1]) / 5 * row[1]
xn = x - xr
print(xr, xn)               # [0.6 0.6 1.2 1.2] [ 0.4  0.4 -0.2 -0.2]
print(xr @ xn)              # 1.1102230246251565e-16   zero, in floating point
print(np.array([xr @ xr, xn @ xn, x @ x]))     # [3.6 0.4 4. ]  Pythagoras
print(A @ x, A @ xr, A @ xn)  # [ 9. 33. 42.] [ 9. 33. 42.] [0. 0. 0.]

import torch

A = torch.tensor([[1., 2., 2., 4.], [3., 8., 6., 16.], [4., 10., 8., 20.]])
row = torch.tensor([[1., 0., 2., 0.], [0., 1., 0., 2.]])     # row space basis
nul = torch.tensor([[-2., 0., 1., 0.], [0., -2., 0., 1.]])   # null space basis
lnul = torch.tensor([[-1., -1., 1.]])                        # left null basis

print(row @ nul.T)          # tensor([[0., 0.], [0., 0.]])
print(A[:, :2].T @ lnul.T)  # tensor([[0.], [0.]])

x = torch.tensor([1., 1., 1., 1.])
xr = (x @ row[0]) / 5 * row[0] + (x @ row[1]) / 5 * row[1]
xn = x - xr
print(xr, xn)   # tensor([0.6, 0.6, 1.2, 1.2]) tensor([0.4, 0.4, -0.2, -0.2])
print(xr @ xn)              # tensor(0.)   zero, in floating point
print(xr @ xr, xn @ xn, x @ x)   # tensor(3.6) tensor(0.4) tensor(4.)
print(A @ x, A @ xr, A @ xn)     # [9., 33., 42.] [9., 33., 42.] [0., 0., 0.]

import jax.numpy as jnp

A = jnp.array([[1., 2., 2., 4.], [3., 8., 6., 16.], [4., 10., 8., 20.]])
row = jnp.array([[1., 0., 2., 0.], [0., 1., 0., 2.]])     # row space basis
nul = jnp.array([[-2., 0., 1., 0.], [0., -2., 0., 1.]])   # null space basis
lnul = jnp.array([[-1., -1., 1.]])                        # left null basis

print(row @ nul.T)          # [[0. 0.] [0. 0.]]
print(A[:, :2].T @ lnul.T)  # [[0.] [0.]]

x = jnp.array([1., 1., 1., 1.])
xr = (x @ row[0]) / 5 * row[0] + (x @ row[1]) / 5 * row[1]
xn = x - xr
print(xr, xn)               # [0.6 0.6 1.2 1.2] [ 0.4  0.4 -0.2 -0.2]
print(xr @ xn)              # 0.0   zero, in floating point
print(jnp.array([xr @ xr, xn @ xn, x @ x]))    # [3.6 0.4 4. ]  Pythagoras
print(A @ x, A @ xr, A @ xn)  # [ 9. 33. 42.] [ 9. 33. 42.] [0. 0. 0.]

import tensorflow as tf

A = tf.constant([[1., 2., 2., 4.], [3., 8., 6., 16.], [4., 10., 8., 20.]])
row = tf.constant([[1., 0., 2., 0.], [0., 1., 0., 2.]])     # row space basis
nul = tf.constant([[-2., 0., 1., 0.], [0., -2., 0., 1.]])   # null space basis
lnul = tf.constant([[-1., -1., 1.]])                        # left null basis

print(row @ tf.transpose(nul))              # [[0. 0.] [0. 0.]]
print(tf.transpose(A[:, :2]) @ tf.transpose(lnul))   # [[0.] [0.]]

x = tf.constant([1., 1., 1., 1.])
xr = tf.tensordot(x, row[0], 1) / 5 * row[0] + tf.tensordot(x, row[1], 1) / 5 * row[1]
xn = x - xr
print(xr, xn)               # [0.6 0.6 1.2 1.2]  [ 0.4  0.4 -0.2 -0.2]
print(tf.tensordot(xr, xn, 1))              # 0.0   zero, in floating point
print(tf.reduce_sum(xr * xr), tf.reduce_sum(xn * xn))   # 3.6  0.4
print(tf.linalg.matvec(A, x), tf.linalg.matvec(A, xr))  # [9. 33. 42.] twice

The second block is the pseudoinverse doing the three steps.

import numpy as np
np.set_printoptions(precision=4, suppress=True)

A = np.array([[1., 2., 2., 4.], [3., 8., 6., 16.], [4., 10., 8., 20.]])
Aplus = np.linalg.pinv(A)

b = np.array([9., 33., 42.])          # reachable: b is in C(A)
print(Aplus @ b)                      # [0.6 0.6 1.2 1.2]   the row space x
print(A @ (Aplus @ b))                # [ 9. 33. 42.]       reached exactly

x = np.array([1., 1., 1., 1.])
s1 = np.array([-2., 0., 1., 0.])
print(Aplus @ A @ x)                  # [0.6 0.6 1.2 1.2]   x projected onto the row space
print(Aplus @ A @ s1)                 # [0. 0. 0. 0.]       the null part dies

b_off = np.array([8., 32., 43.])      # b plus the left null space vector
print(Aplus @ b_off)                  # [0.6 0.6 1.2 1.2]   same answer
print(A @ (Aplus @ b_off) - b_off)    # [ 1.  1. -1.]  the part it threw away
print(np.allclose(A @ Aplus @ A, A))  # True

import torch

A = torch.tensor([[1., 2., 2., 4.], [3., 8., 6., 16.], [4., 10., 8., 20.]])
Aplus = torch.linalg.pinv(A)

b = torch.tensor([9., 33., 42.])      # reachable: b is in C(A)
print(Aplus @ b)                      # tensor([0.6, 0.6, 1.2, 1.2])
print(A @ (Aplus @ b))                # tensor([ 9., 33., 42.])

x = torch.tensor([1., 1., 1., 1.])
s1 = torch.tensor([-2., 0., 1., 0.])
print(Aplus @ A @ x)                  # tensor([0.6, 0.6, 1.2, 1.2])
print(Aplus @ A @ s1)                 # tensor([0., 0., 0., 0.])

b_off = torch.tensor([8., 32., 43.])  # b plus the left null space vector
print(Aplus @ b_off)                  # tensor([0.6, 0.6, 1.2, 1.2])   same answer
print(A @ (Aplus @ b_off) - b_off)    # tensor([ 1.,  1., -1.])
print(torch.allclose(A @ Aplus @ A, A))   # True

import jax.numpy as jnp

A = jnp.array([[1., 2., 2., 4.], [3., 8., 6., 16.], [4., 10., 8., 20.]])
Aplus = jnp.linalg.pinv(A)

b = jnp.array([9., 33., 42.])         # reachable: b is in C(A)
print(Aplus @ b)                      # [0.6 0.6 1.2 1.2]   the row space x
print(A @ (Aplus @ b))                # [ 9. 33. 42.]       reached exactly

x = jnp.array([1., 1., 1., 1.])
s1 = jnp.array([-2., 0., 1., 0.])
print(Aplus @ A @ x)                  # [0.6 0.6 1.2 1.2]
print(Aplus @ A @ s1)                 # [0. 0. 0. 0.]

b_off = jnp.array([8., 32., 43.])     # b plus the left null space vector
print(Aplus @ b_off)                  # [0.6 0.6 1.2 1.2]   same answer
print(A @ (Aplus @ b_off) - b_off)    # [ 1.  1. -1.]
print(jnp.allclose(A @ Aplus @ A, A))     # True

import tensorflow as tf

A = tf.constant([[1., 2., 2., 4.], [3., 8., 6., 16.], [4., 10., 8., 20.]])
Aplus = tf.linalg.pinv(A)

b = tf.constant([[9.], [33.], [42.]])   # reachable: b is in C(A)
print(Aplus @ b)                        # [[0.6] [0.6] [1.2] [1.2]]
print(A @ (Aplus @ b))                  # [[ 9.] [33.] [42.]]

x = tf.constant([[1.], [1.], [1.], [1.]])
s1 = tf.constant([[-2.], [0.], [1.], [0.]])
print(Aplus @ (A @ x))                  # [[0.6] [0.6] [1.2] [1.2]]
print(Aplus @ (A @ s1))                 # [[0.] [0.] [0.] [0.]]

b_off = tf.constant([[8.], [32.], [43.]])   # b plus the left null space vector
print(Aplus @ b_off)                    # [[0.6] [0.6] [1.2] [1.2]]  same answer
print(A @ (Aplus @ b_off) - b_off)      # [[ 1.] [ 1.] [-1.]]
print(tf.reduce_all(tf.abs(A @ Aplus @ A - A) < 1e-5))   # True

Where this is going

Two things in this chapter were used before they were earned. The split $x = x_r + x_n$ was described as though computing it were easy, and step 1 of the pseudoinverse asked for the part of $b$ that lies inside the column space. Both are the same operation: dropping a vector perpendicularly onto a subspace. That operation is called projection, it has a short formula, and it is the next chapter. Once you can project, the whole of least squares is one substitution away, and least squares is how every linear model in machine learning gets fitted.