Subspaces with Python

The previous lesson defined subspaces and established the test for recognising them. This one is about operating with them: deciding whether a vector belongs to one, finding the nearest point when it does not, and choosing the subspace that best represents a set of data.

Membership

Whether $\vec{b}$ belongs to the column space of $A$ needs no new machinery. By definition, $\operatorname{col}(A)$ is the set of linear combinations of the columns, so

\vec{b} \in \operatorname{col}(A) \iff \exists\,\vec{x} : A\vec{x} = \vec{b}

Solving a system and asking about membership are the same operation stated two ways.

import numpy as np

V = np.column_stack([[1., 0., 0.],
                     [0., 1., 0.]])        # span = the xy plane
w = np.array([3., 2., 0.])

c, *_ = np.linalg.lstsq(V, w, rcond=None)
np.allclose(V @ c, w)                      # True: it belongs

The rank criterion of the first lesson restates itself in this language. The system is consistent when adding $\vec{b}$ to the columns of $A$ does not enlarge the subspace they span:

\operatorname{rank}(A) = \operatorname{rank}([A \mid \vec{b}]) \iff \vec{b} \in \operatorname{col}(A)

This is the Rouché-Capelli theorem, which the lesson on linear systems stated in terms of ranks and which turns out here to be a statement about membership.

Orthogonal projection

When $\vec{b}$ does not belong to the subspace, the useful question stops being whether a solution exists and becomes which point of the subspace is nearest to $\vec{b}$ . That point is the orthogonal projection.

For the subspace spanned by a single vector $\vec{v}$ , the projection is given by a matrix:

P = \frac{\vec{v}\vec{v}^\top}{\vec{v}^\top\vec{v}}, \qquad P\vec{w} = \frac{\vec{v}^\top\vec{w}}{\vec{v}^\top\vec{v}}\,\vec{v}

Drag the direction of the line, or the point being projected.

Pw = (1.65, 0.55)

‖w − Pw‖ = 2.06

0.900.300.300.10

P² = P: projecting an already projected point does not move it. ✓

The dashed segment is the residual $\vec{w} - P\vec{w}$ , and it is perpendicular to the line wherever the point sits. That perpendicularity is not a visual observation but the condition defining the projection: the residual lies in the orthogonal complement of the subspace.

From it the normal equation of the lesson on the inverse and the transpose follows directly. Requiring the residual to be orthogonal to every column of $A$ is to write

A^\top(\vec{b} - A\vec{x}) = \vec{0} \quad\Longrightarrow\quad A^\top A\,\vec{x} = A^\top\vec{b}

The normal equations are not a computational recipe: they are the orthogonality condition on the residual.

For a higher-dimensional subspace spanned by the columns of $V$ , the projection matrix generalises the expression above:

P = V(V^\top V)^{-1}V^\top

and reduces to $P = QQ^\top$ when the columns of $Q$ are orthonormal, which is the form returned by the QR factorisation or by null_space. Three properties characterise it:

Q, _ = np.linalg.qr(V)
P = Q @ Q.T

np.allclose(P @ P, P)          # idempotent: projecting twice is projecting once
np.allclose(P, P.T)            # symmetric
np.trace(P)                    # the dimension of the subspace

The trace deserves attention: for a projection it equals the rank, so it gives the dimension of the subspace without computing it any other way.

The four fundamental subspaces

A matrix $A \in \mathbb{R}^{m\times n}$ determines four subspaces, two in each space:

Subspace	Lives in	Dimension
column space, $\operatorname{col}(A)$	$\mathbb{R}^m$	$r$
null space, $\ker(A)$	$\mathbb{R}^n$	$n - r$
row space, $\operatorname{col}(A^\top)$	$\mathbb{R}^n$	$r$
left null space, $\ker(A^\top)$	$\mathbb{R}^m$	$m - r$

where $r = \operatorname{rank}(A)$ . The two pairs are orthogonal complements:

\ker(A) \perp \operatorname{col}(A^\top) \quad\text{in } \mathbb{R}^n, \qquad \ker(A^\top) \perp \operatorname{col}(A) \quad\text{in } \mathbb{R}^m

The first relation was already used in the lesson on particular and general solutions: the minimum-norm solution is the component in the row space, precisely because that is the orthogonal complement of the null space.

The second explains what happens when a system is inconsistent. The vector $\vec{b}$ decomposes uniquely into a part inside $\operatorname{col}(A)$ , which is what the system can reach, and another in $\ker(A^\top)$ , which is the residual no choice of $\vec{x}$ removes.

Application: the best subspace

Principal component analysis solves a problem that can now be stated precisely: given a set of centred points, find the subspace of dimension $k$ minimising the sum of the squared distances from the points to it.

The line turns; the segments are the distances being squared and summed.

θ0°

residual = 9.25 · explained 72.7 %

minimum at θ = 30.8°

For $k = 1$ the subspace is a line through the origin and the problem reduces to choosing an angle. The total residual varies between $9.25$ at $0°$ and $24.68$ at $90°$ , and reaches its minimum, $0.76$ , at $30.8°$ . That direction retains $97.8\,\%$ of the variance of the data.

The equivalence that makes the problem tractable appears in the figure: since the residual and the projection are orthogonal, the Pythagorean theorem gives

\lVert \vec{x} \rVert^2 = \lVert P\vec{x} \rVert^2 + \lVert \vec{x} - P\vec{x} \rVert^2

and the sum over all points is constant. Minimising the residual and maximising the projected variance are therefore the same problem. The first formulation is geometric and the second statistical, and they coincide by orthogonality.

The solution requires no search. The optimal direction is the eigenvector of largest eigenvalue of the covariance matrix $\tfrac{1}{n}X^\top X$ , which the lesson on the inverse and the transpose identified as symmetric and positive semidefinite. That property is what guarantees the eigenvalues are real and an orthogonal basis of eigenvectors exists, without which the problem would not be well posed.

X = X - X.mean(axis=0)                       # centre
C = X.T @ X / len(X)
values, vectors = np.linalg.eigh(C)          # eigh: exploits the symmetry
direction = vectors[:, -1]                   # largest eigenvalue

np.linalg.eigh is used rather than eig because it exploits the symmetry: it is faster and returns real eigenvalues by construction, rather than real up to rounding error.

What a model can represent

The vocabulary of this lesson describes the limits of a linear model precisely. The column space of the design matrix is the set of reachable predictions: if the target vector does not belong to it, no adjustment of the weights will reach it, and the best available is its projection. The null space is the set of directions the model does not distinguish, which is the collinearity described in the lesson on the null space.

Increasing the capacity of a linear model means enlarging its column space by adding columns — interactions, non-linear transformations of the features — until the target falls inside it or close enough.

Exercise. Build the projection matrix onto the plane spanned by $(1,0,0)$ and $(0,1,0)$ and check that its trace is $2$ . Then verify that $P$ and $I - P$ are both idempotent, that their sum is the identity, and explain which subspace $I - P$ projects onto.