7 Vector Spaces

So far vectors have been lists of numbers and matrices have been machines that act on them. This chapter asks a different kind of question: given a matrix, what set of vectors can it reach, and what does it destroy?

The answers are the most reusable ideas in the subject. Span, independence, basis, rank, and orthogonality are the vocabulary you need to say precisely why a regression is unidentifiable, what PCA is choosing, and what “effective number of parameters” means.

7.1 What makes a vector space

A vector space is a set of objects you can add together and scale, where the results stay in the set and the usual arithmetic works.

Formally, a set $V$ is a vector space if for any $\mathbf{u}, \mathbf{v} \in V$ and any scalars $c, d$:

$\mathbf{u} + \mathbf{v} \in V$ — closed under addition
$c\mathbf{u} \in V$ — closed under scalar multiplication

plus eight bookkeeping axioms (associativity, commutativity, a zero vector, additive inverses, and so on) that hold automatically for everything in this book. The two closure conditions are the ones that ever fail in practice, so they are the two to check.

$\mathbb{R}^n$ is the vector space we care about. But the definition is deliberately abstract, and that generality pays: polynomials of degree $\leq 3$ form a vector space, as do continuous functions on $[0,1]$, and as do $m \times n$ matrices. Anything proved about vector spaces in general applies to all of them.

In machine learning

The abstraction is not decoration. Kernel methods work by treating functions as vectors in an infinite-dimensional space, where “inner product” and “projection” still mean what they mean in $\mathbb{R}^n$. Everything in this chapter transfers.

7.2 Subspaces

A subspace is a subset of a vector space that is itself a vector space — a set closed under addition and scaling.

The closure requirement has an immediate consequence that catches people out: taking $c = 0$ shows every subspace must contain the zero vector. A subspace of $\mathbb{R}^2$ therefore cannot be an arbitrary line; it must be a line through the origin.

xr <- c(-3, 3)
draw_plane(
  curves = list(
    `through the origin` = cbind(xr, 0.5 * xr),
    `not through the origin` = cbind(xr, 0.5 * xr + 1.5)
  ),
  points = cbind(0, 0),
  notes = list("0" = c(0.3, -0.45))
)

Figure 7.1: A line through the origin is a subspace: scale or add any two of its points and you stay on it. A line that misses the origin is not — it fails at the very first test, since it does not contain $\mathbf{0}$.

The subspaces of $\mathbb{R}^3$ are exactly: the origin alone, every line through the origin, every plane through the origin, and $\mathbb{R}^3$ itself. Nothing else qualifies. That short list is worth holding onto — it means “subspace” always describes something flat, unbounded, and anchored at zero.

7.3 Span

The span of a set of vectors is the set of all their linear combinations:

\[ \operatorname{span}\{\mathbf{v}_1, \dots, \mathbf{v}_k\} = \left\{ c_1\mathbf{v}_1 + \cdots + c_k\mathbf{v}_k \;:\; c_i \in \mathbb{R} \right\} \]

It is everything you can reach from those vectors using only the two operations you have. A span is always a subspace — it contains $\mathbf{0}$ (take all $c_i = 0$) and is closed by construction.

The span of one non-zero vector is a line through the origin. The span of two independent vectors in $\mathbb{R}^2$ is the whole plane, as we saw in Figure 4.4.

But two vectors do not always give you a plane.

v1 <- c(2, 1)
v2 <- c(4, 2)
tt <- c(-1.6, 1.9)
draw_plane(
  vectors = list(v1 = v1, v2 = v2),
  curves = list(cbind(2 * tt, 1 * tt)),
  curve_line_type = "dashed",
  curve_color = amds_gray,
  label_at = c("mid", "tip")
)

Figure 7.2: Two vectors, one line. Because $\mathbf{v}_2 = 2\mathbf{v}_1$, every combination $c_1\mathbf{v}_1 + c_2\mathbf{v}_2$ is still a multiple of $\mathbf{v}_1$, so their span is only a line.

The second vector added nothing. This is the situation the next section names.

7.4 Linear independence

Vectors $\mathbf{v}_1, \dots, \mathbf{v}_k$ are linearly independent when the only way to combine them into the zero vector is to use all-zero coefficients:

\[ c_1\mathbf{v}_1 + \cdots + c_k\mathbf{v}_k = \mathbf{0} \quad\Longrightarrow\quad c_1 = \cdots = c_k = 0 \tag{7.1}\]

If some non-trivial combination gives $\mathbf{0}$, they are linearly dependent, and that equation lets you write at least one of them in terms of the others. A dependent vector is redundant: removing it does not shrink the span.

Worked example. For $\mathbf{v}_1 = (2,1)$ and $\mathbf{v}_2 = (4,2)$, we have $2\mathbf{v}_1 - \mathbf{v}_2 = \mathbf{0}$ with coefficients $(2, -1) \neq (0,0)$. Dependent.

For $\mathbf{a} = (3,1)$ and $\mathbf{b} = (1,2)$, suppose $c_1\mathbf{a} + c_2\mathbf{b} = \mathbf{0}$. That is two equations, $3c_1 + c_2 = 0$ and $c_1 + 2c_2 = 0$, whose only solution is $c_1 = c_2 = 0$. Independent.

Notice that Equation 7.1 is exactly the statement $\mathbf{V}\mathbf{c} = \mathbf{0}$ for the matrix $\mathbf{V}$ whose columns are the vectors. So:

The columns of a matrix are independent exactly when $\mathbf{V}\mathbf{c} = \mathbf{0}$ has only the zero solution.

For a square matrix that is the same as being invertible, and the same as having a non-zero determinant.

dependent <- cbind(c(2, 1), c(4, 2))
independent <- cbind(c(3, 1), c(1, 2))
c(det(dependent), det(independent))

[1] 0 5

In machine learning

Linearly dependent columns in a design matrix make a model unidentifiable: infinitely many coefficient vectors fit the data equally well, so no unique answer exists. The classic case is the dummy-variable trap — one indicator per category plus an intercept, where the indicators sum to the intercept column. R silently drops the redundant column and reports NA for its coefficient.

7.5 Basis and dimension

A basis for a subspace is a set of vectors that is both

spanning — every vector in the subspace is some combination of them, and
independent — none of them is redundant.

A basis is a minimal complete description: enough vectors to reach everything, with nothing wasted.

The dimension of a subspace is the number of vectors in a basis. Any two bases for the same subspace have the same size, which is what makes dimension well defined.

The standard basis for $\mathbb{R}^n$ is $\mathbf{e}_1, \dots, \mathbf{e}_n$, where $\mathbf{e}_i$ has a 1 in position $i$ and zeros elsewhere.

diag(3) # columns are e1, e2, e3

     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

But it is not the only basis, and often not the best one. $\{(2,1), (-1,2)\}$ is also a basis for $\mathbb{R}^2$: two independent vectors in a 2-dimensional space always are.

B <- cbind(c(2, 1), c(-1, 2))
det(B) # non-zero, so independent, so a basis

[1] 5

In machine learning

Choosing a good basis is much of unsupervised learning. PCA finds the basis whose first few directions capture the most variance. The Fourier basis makes periodic signals sparse; wavelets do it for localized ones. Same data, different description, and the right description makes the structure visible.

7.6 Coordinates

Once you fix a basis, every vector has exactly one set of coordinates — the coefficients that build it.

This is worth pausing on. When you write $\mathbf{w} = (3, 4)$, you are not stating something intrinsic about $\mathbf{w}$; you are stating its coordinates in the standard basis, because $\mathbf{w} = 3\mathbf{e}_1 + 4\mathbf{e}_2$. The vector is the arrow. The numbers are a description of the arrow relative to a choice.

Worked example. Express $\mathbf{w} = (3,4)$ in the basis $\{\mathbf{v}_1, \mathbf{v}_2\}$ with $\mathbf{v}_1 = (2,1)$, $\mathbf{v}_2 = (-1,2)$. We need $c_1, c_2$ with

\[ c_1\begin{bmatrix}2\\1\end{bmatrix} + c_2\begin{bmatrix}-1\\2\end{bmatrix} = \begin{bmatrix}3\\4\end{bmatrix} \]

which is the linear system $\mathbf{B}\mathbf{c} = \mathbf{w}$.

w <- c(3, 4)
coords <- solve(B, w)
coords

[1] 2 1

So $\mathbf{w} = 2\mathbf{v}_1 + 1\mathbf{v}_2$. Same arrow, different numbers: $(3,4)$ in the standard basis, $(2,1)$ in this one.

draw_plane(
  vectors = list(
    `2v1` = 2 * B[, 1],
    v2 = B[, 2],
    w = w
  ),
  origin = list(c(0, 0), 2 * B[, 1], c(0, 0)),
  line_type = c("dashed", "dashed", "solid"),
  label_at = c("mid", "mid", "tip")
)

Figure 7.3: The vector $\mathbf{w}$ decomposed in the basis $\{\mathbf{v}_1, \mathbf{v}_2\}$: two steps along $\mathbf{v}_1$, then one along $\mathbf{v}_2$. Its coordinates in this basis are $(2, 1)$, not $(3, 4)$.

7.7 Change of basis

The previous section is a general recipe. Let $\mathbf{B}$ be the matrix whose columns are the new basis vectors. Then

\[ \mathbf{w} = \mathbf{B}\mathbf{c} \qquad\Longleftrightarrow\qquad \mathbf{c} = \mathbf{B}^{-1}\mathbf{w} \tag{7.2}\]

Reading it in the useful direction: $\mathbf{B}$ converts new-basis coordinates into standard ones, and $\mathbf{B}^{-1}$ goes the other way.

drop(B %*% coords) # coordinates -> standard

[1] 3 4

drop(solve(B, w)) # standard -> coordinates

[1] 2 1

Because $\mathbf{B}$ has independent columns it is invertible, so the conversion never loses information. A change of basis is a relabeling, not a transformation of the underlying object.

When the basis is orthonormal (Section 7.11) this gets much cheaper: $\mathbf{B}^{-1} = \mathbf{B}^\top$, and finding coordinates costs a few dot products instead of a linear solve.

7.8 Rank

The rank of a matrix is the dimension of the span of its columns — the number of genuinely independent directions it contains.

Equivalently, and not obviously, it is the number of pivots in row echelon form (Section 6.4), and also the dimension of the span of its rows. Row rank always equals column rank, which is one of the small miracles of the subject.

A matrix is full rank when its rank is as large as its shape allows, $\min(m, n)$, and rank deficient otherwise.

A <- matrix(
  c(1, 2, 3, 2, 4, 6, 1, 1, 2),
  nrow = 3, byrow = TRUE
)
A

     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    4    6
[3,]    1    1    2

qr(A)$rank

[1] 2

Rank 2, not 3. Row 2 is twice row 1, and the third column is the sum of the first two — the same deficiency showing up in both directions, as row rank equals column rank requires.

A[, 1] + A[, 2] == A[, 3]

[1] TRUE TRUE TRUE

Watch out

Rank is a discrete quantity computed in floating-point arithmetic, which makes it fragile. A matrix that is rank deficient in exact arithmetic will usually come back full-rank numerically, with one singular value at $10^{-16}$ instead of 0.

For real data, ask “how close to rank deficient” rather than “what is the rank” — the condition number (Section 6.10) or the singular values (Section 9.6) answer that, and qr()$rank does not.

7.9 The four fundamental subspaces

Every $m \times n$ matrix determines four subspaces, and together they say everything about what the matrix does.

Subspace	Definition	Lives in	Dimension
Column space $C(\mathbf{A})$	span of the columns	$\mathbb{R}^m$	$r$
Null space $N(\mathbf{A})$	all $\mathbf{x}$ with $\mathbf{A}\mathbf{x} = \mathbf{0}$	$\mathbb{R}^n$	$n - r$
Row space $C(\mathbf{A}^\top)$	span of the rows	$\mathbb{R}^n$	$r$
Left null space $N(\mathbf{A}^\top)$	all $\mathbf{y}$ with $\mathbf{A}^\top\mathbf{y} = \mathbf{0}$	$\mathbb{R}^m$	$m - r$

where $r = \operatorname{rank}(\mathbf{A})$.

Two of them you have already met under other names. The column space is the set of reachable outputs — $\mathbf{A}\mathbf{x} = \mathbf{b}$ is solvable exactly when $\mathbf{b} \in C(\mathbf{A})$. The null space is what the matrix destroys: every vector in it is mapped to zero, so information about it cannot be recovered.

The dimension column contains a theorem worth stating on its own. The rank–nullity theorem:

\[ \operatorname{rank}(\mathbf{A}) + \dim N(\mathbf{A}) = n \tag{7.3}\]

Every input dimension is either preserved or annihilated; there is no third option.

Worked example. Our rank-2 matrix $\mathbf{A}$ is $3 \times 3$, so its null space must have dimension $3 - 2 = 1$. Since $\mathbf{a}_1 + \mathbf{a}_2 - \mathbf{a}_3 = \mathbf{0}$, the vector $(1, 1, -1)$ is in it:

nullvec <- c(1, 1, -1)
drop(A %*% nullvec)

[1] 0 0 0

And the left null space has dimension $3 - 2 = 1$ too. Row 2 is twice row 1, so $2\mathbf{r}_1 - \mathbf{r}_2 = \mathbf{0}$ and $(2, -1, 0)$ works:

drop(crossprod(A, c(2, -1, 0)))

[1] 0 0 0

In machine learning

A non-trivial null space in a design matrix is exactly non-identifiability: if $\mathbf{n} \in N(\mathbf{X})$ then $\boldsymbol{\beta}$ and $\boldsymbol{\beta} + \mathbf{n}$ produce identical predictions for every observation, so the data cannot distinguish them. Ridge regression (Section 6.11) resolves the tie by preferring the smaller coefficient vector.

7.10 Orthogonality

Two subspaces are orthogonal when every vector in one is orthogonal to every vector in the other.

The four fundamental subspaces pair up this way, and the pairing is the structural heart of the subject:

\[ N(\mathbf{A}) \perp C(\mathbf{A}^\top) \qquad\text{and}\qquad N(\mathbf{A}^\top) \perp C(\mathbf{A}) \]

The first is easy to see. If $\mathbf{A}\mathbf{x} = \mathbf{0}$, then every row of $\mathbf{A}$ dotted with $\mathbf{x}$ gives zero — which is precisely what $\mathbf{A}\mathbf{x} = \mathbf{0}$ says, read one row at a time. So $\mathbf{x}$ is orthogonal to every row, hence to their whole span.

A %*% nullvec # each entry is a row dotted with nullvec

     [,1]
[1,]    0
[2,]    0
[3,]    0

More than orthogonal, these pairs are orthogonal complements: together they fill the whole space, and their dimensions add to it. $\mathbb{R}^n$ splits cleanly into the row space and the null space, with every vector having exactly one piece in each.

That decomposition is the engine of least squares. The second pairing says the residual, which lives in $N(\mathbf{A}^\top)$, is orthogonal to the column space — the geometric fact behind the normal equations (Section 6.9).

7.11 Orthonormal bases

A basis is orthonormal when its vectors are mutually orthogonal and each has length 1:

\[ \mathbf{u}_i \cdot \mathbf{u}_j = \begin{cases} 1 & i = j \\ 0 & i \neq j \end{cases} \]

Orthonormal bases are worth going out of your way for. Collect them as columns of a matrix $\mathbf{Q}$ and the definition above says exactly $\mathbf{Q}^\top\mathbf{Q} = \mathbf{I}$ — so for a square $\mathbf{Q}$,

\[ \mathbf{Q}^{-1} = \mathbf{Q}^\top \]

An inverse for the price of a transpose. Three consequences:

Coordinates are dot products. $\mathbf{c} = \mathbf{Q}^\top\mathbf{w}$, so $c_i = \mathbf{u}_i \cdot \mathbf{w}$. No solve needed.
Lengths and angles are preserved, since $\|\mathbf{Q}\mathbf{x}\|^2 = \mathbf{x}^\top\mathbf{Q}^\top\mathbf{Q}\mathbf{x} = \|\mathbf{x}\|^2$. Orthonormal matrices are exactly the rotations and reflections.
Numerically excellent. $\kappa(\mathbf{Q}) = 1$, the best possible, so orthonormal matrices never amplify error. This is why serious algorithms are built from them.

7.12 Gram–Schmidt

Given any independent set, Gram–Schmidt manufactures an orthonormal basis for the same span. The idea is one step of Section 4.14, repeated:

Take the next vector, subtract off its projection onto everything you have already orthonormalized, and normalize what remains.

Formally, for $j = 1, \dots, k$:

\[ \mathbf{w}_j = \mathbf{v}_j - \sum_{i<j}(\mathbf{u}_i \cdot \mathbf{v}_j)\,\mathbf{u}_i \qquad \mathbf{u}_j = \frac{\mathbf{w}_j}{\|\mathbf{w}_j\|} \]

Worked example. Take $\mathbf{v}_1 = (3,1)$ and $\mathbf{v}_2 = (1,2)$ from Chapter 4.

$\mathbf{u}_1 = \mathbf{v}_1/\sqrt{10} \approx (0.949, 0.316)$.

For the second: $\mathbf{u}_1 \cdot \mathbf{v}_2 = 5/\sqrt{10}$, so

\[ \mathbf{w}_2 = \begin{bmatrix}1\\2\end{bmatrix} - \frac{5}{\sqrt{10}}\cdot\frac{1}{\sqrt{10}}\begin{bmatrix}3\\1\end{bmatrix} = \begin{bmatrix}1\\2\end{bmatrix} - \begin{bmatrix}1.5\\0.5\end{bmatrix} = \begin{bmatrix}-0.5\\1.5\end{bmatrix} \]

That is exactly the residual computed in Figure 4.8 — Gram–Schmidt is the projection you already know, applied in sequence. Normalizing gives $\mathbf{u}_2 \approx (-0.316, 0.949)$.

gram_schmidt <- function(V) {
  U <- matrix(0, nrow(V), ncol(V))
  for (j in seq_len(ncol(V))) {
    w <- V[, j]
    for (i in seq_len(j - 1L)) {
      w <- w - drop(crossprod(U[, i], V[, j])) * U[, i]
    }
    U[, j] <- w / sqrt(sum(w^2))
  }
  U
}
V <- cbind(c(3, 1), c(1, 2))
U <- gram_schmidt(V)
round(U, 4)

       [,1]    [,2]
[1,] 0.9487 -0.3162
[2,] 0.3162  0.9487

round(crossprod(U), 10) # should be the identity

     [,1] [,2]
[1,]    1    0
[2,]    0    1

draw_plane(
  vectors = list(
    v1 = V[, 1],
    v2 = V[, 2],
    u1 = U[, 1],
    u2 = U[, 2]
  ),
  color = c(
    amds_gray, amds_gray,
    amds_colors[1], amds_colors[2]
  ),
  line_type = c("dashed", "dashed", "solid", "solid"),
  label_at = c("tip", "tip", "mid", "mid")
)

Figure 7.4: Gram–Schmidt in the plane. $\mathbf{u}_1$ points along $\mathbf{v}_1$; $\mathbf{u}_2$ is what is left of $\mathbf{v}_2$ after removing its component along $\mathbf{u}_1$.

It works in any dimension:

V3 <- cbind(c(1, 1, 0), c(1, 0, 1), c(0, 1, 1))
U3 <- gram_schmidt(V3)
round(U3, 4)

       [,1]    [,2]    [,3]
[1,] 0.7071  0.4082 -0.5774
[2,] 0.7071 -0.4082  0.5774
[3,] 0.0000  0.8165  0.5774

round(crossprod(U3), 10)

     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

Watch out

Classical Gram–Schmidt as written above is numerically unstable: rounding errors accumulate and the later vectors drift away from orthogonality. Use it to understand the idea, not to compute. R’s qr() uses Householder reflections instead, which are stable, and gives you the same orthonormal basis:

round(qr.Q(qr(V3)), 4)

        [,1]    [,2]    [,3]
[1,] -0.7071  0.4082 -0.5774
[2,] -0.7071 -0.4082  0.5774
[3,]  0.0000  0.8165  0.5774

Signs may differ — an orthonormal basis is only determined up to the direction of each vector.

7.13 Projection onto a subspace

Section 4.14 projected one vector onto another. The general case projects a vector onto a whole subspace, and it is the last piece of machinery least squares needs.

Let $\mathbf{A}$ have independent columns spanning the subspace. The projection of $\mathbf{b}$ onto that subspace is

\[ \operatorname{proj}(\mathbf{b}) = \mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top\mathbf{b} = \mathbf{P}\mathbf{b} \tag{7.4}\]

where $\mathbf{P} = \mathbf{A}(\mathbf{A}^\top\mathbf{A})^{-1}\mathbf{A}^\top$ is the projection matrix. Compare it with the single-vector case $\frac{\mathbf{a}\cdot\mathbf{b}}{\mathbf{a}\cdot\mathbf{a}}\mathbf{a}$: the same shape, with $\mathbf{A}^\top\mathbf{A}$ playing the role of $\mathbf{a}\cdot\mathbf{a}$.

Worked example. Project $\mathbf{b} = (1,2,3)$ onto the plane spanned by $(1,1,0)$ and $(0,1,1)$.

Ap <- cbind(c(1, 1, 0), c(0, 1, 1))
b <- c(1, 2, 3)
coefs <- solve(crossprod(Ap), crossprod(Ap, b))
proj <- drop(Ap %*% coefs)
resid <- b - proj
round(rbind(proj = proj, resid = resid), 4)

        [,1]    [,2]   [,3]
proj  0.3333  2.6667 2.3333
resid 0.6667 -0.6667 0.6667

The residual is orthogonal to both spanning vectors, so it is orthogonal to the entire plane:

round(drop(crossprod(Ap, resid)), 10)

[1] 0 0

The projection matrix itself has two defining properties.

P <- Ap %*% solve(crossprod(Ap), t(Ap))
all.equal(P %*% P, P) # idempotent: projecting twice

[1] TRUE

all.equal(P, t(P)) # symmetric

[1] TRUE

sum(diag(P)) # trace = dim of the subspace

[1] 2

Idempotent because once you are in the subspace, projecting again does nothing. Symmetric because orthogonal projection is its own adjoint. And the trace equals the dimension of the subspace being projected onto — a fact that reappears as the effective degrees of freedom of a linear model.

In machine learning

In least squares, $\hat{\mathbf{y}} = \mathbf{P}\mathbf{y}$ where $\mathbf{P} = \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top$ is called the hat matrix — it puts the hat on $\mathbf{y}$. Its diagonal entries are the leverages, measuring how much each observation pulls on its own fitted value, and its trace is the number of parameters. Regression diagnostics are projection geometry.

7.14 Summary

Concept	Question it answers	R
Span	what can these vectors reach?	—
Independence	is any of them redundant?	`det()`, `qr()$rank`
Basis	a minimal complete description	—
Dimension	how many directions?	—
Rank	how many independent columns?	`qr(A)$rank`
Null space	what does the matrix destroy?	solve $\mathbf{A}\mathbf{x}=\mathbf{0}$
Orthonormal basis	the numerically ideal basis	`qr.Q(qr(A))`
Projection	closest point in a subspace	`A %*% solve(crossprod(A), crossprod(A, b))`

7.15 Exercises

1. Is the set of vectors in $\mathbb{R}^2$ with $x_1 + x_2 = 0$ a subspace? What about $x_1 + x_2 = 1$?

Solution

The first is a subspace. It contains $\mathbf{0}$; if $u_1 + u_2 = 0$ and $v_1 + v_2 = 0$ then $(u_1+v_1) + (u_2+v_2) = 0$, and scaling preserves the equation too. Geometrically it is the line $x_2 = -x_1$ through the origin.

The second is not. It fails immediately: $\mathbf{0}$ does not satisfy $0 + 0 = 1$. It is a line parallel to the first but shifted off the origin — an affine set, not a subspace.

Any set defined by $\mathbf{A}\mathbf{x} = \mathbf{0}$ is a subspace (it is a null space); any set defined by $\mathbf{A}\mathbf{x} = \mathbf{b}$ with $\mathbf{b} \neq \mathbf{0}$ is not.

2. Are $(1,2,3)$, $(2,4,6)$, and $(1,0,1)$ linearly independent? What is the dimension of their span?

Solution

No — the second is twice the first.

M <- cbind(c(1, 2, 3), c(2, 4, 6), c(1, 0, 1))
qr(M)$rank

[1] 2

Rank 2, so the span is a 2-dimensional plane in $\mathbb{R}^3$, not all of it. Dropping the redundant second vector leaves a basis for that plane.

3. For the matrix A in Section 7.8, verify the rank–nullity theorem and find a basis for the null space.

Solution

Rank is 2 and $n = 3$, so the null space has dimension 1 — a single basis vector.

r <- qr(A)$rank
c(rank = r, nullity = ncol(A) - r)

   rank nullity 
      2       1

drop(A %*% c(1, 1, -1))

[1] 0 0 0

Any non-zero multiple of $(1,1,-1)$ is an equally valid basis. The null space is a line through the origin in $\mathbb{R}^3$.

4. Show that if $\mathbf{Q}$ has orthonormal columns then $\|\mathbf{Q}\mathbf{x}\| = \|\mathbf{x}\|$ for every $\mathbf{x}$. Verify numerically.

Solution

\[ \|\mathbf{Q}\mathbf{x}\|^2 = (\mathbf{Q}\mathbf{x})^\top(\mathbf{Q}\mathbf{x}) = \mathbf{x}^\top\mathbf{Q}^\top\mathbf{Q}\mathbf{x} = \mathbf{x}^\top\mathbf{I}\mathbf{x} = \|\mathbf{x}\|^2 \]

xv <- c(2, -3)
c(
  before = sqrt(sum(xv^2)),
  after = sqrt(sum((U %*% xv)^2))
)

  before    after 
3.605551 3.605551

Orthonormal matrices move vectors around without changing their size — they are rigid motions. This is why they are the safe building block for numerical algorithms.

5. Project $\mathbf{b} = (4, 0, 0)$ onto the plane spanned by $(1,1,0)$ and $(0,1,1)$. Confirm the residual is orthogonal to both.

Solution

b5 <- c(4, 0, 0)
c5 <- solve(crossprod(Ap), crossprod(Ap, b5))
p5 <- drop(Ap %*% c5)
round(p5, 4)

[1]  2.6667  1.3333 -1.3333

round(drop(crossprod(Ap, b5 - p5)), 10)

[1] 0 0

The projection is the closest point of the plane to $\mathbf{b}$, and the error points straight out of the plane — the shortest way back.

6. Run gram_schmidt() on the Läuchli matrix — three columns that are nearly identical, differing only in a tiny entry:

\[ \mathbf{V} = \begin{bmatrix} 1 & 1 & 1 \\ \varepsilon & 0 & 0 \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix} \qquad \varepsilon = 10^{-8} \]

Check how orthogonal the result really is, and compare against qr().

Solution

eps <- 1e-8
Vbad <- rbind(
  c(1, 1, 1),
  c(eps, 0, 0),
  c(0, eps, 0),
  c(0, 0, eps)
)
Ubad <- gram_schmidt(Vbad)
round(crossprod(Ubad), 8)

       [,1]   [,2]   [,3]
[1,]  1e+00 -1e-08 -1e-08
[2,] -1e-08  1e+00  5e-01
[3,] -1e-08  5e-01  1e+00

The last two columns have a dot product of $0.5$. They are supposed to be orthogonal, and instead they are at $60°$ — not a rounding-level error but a total failure.

Here is why. Each $\mathbf{v}_j$ is dominated by its first component, so after subtracting the projection onto $\mathbf{u}_1$ almost nothing survives: $\mathbf{w}_2$ and $\mathbf{w}_3$ both have norm $\varepsilon\sqrt{2}$. Dividing by that tiny norm magnifies whatever rounding error was in them by a factor of $1/\varepsilon$, and the errors in $\mathbf{w}_2$ and $\mathbf{w}_3$ are correlated, because both came from cancelling against the same $\mathbf{u}_1$. The algorithm then treats the amplified noise as signal.

A stable method does not lose the orthogonality:

round(crossprod(qr.Q(qr(Vbad))), 8)

     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

This is why Section 7.12 says to use Gram–Schmidt to understand and qr() to compute. Note the distinction from the collinear regression in Section 6.11: there, the data had genuinely lost the information. Here the information is present and a good algorithm recovers it — the failure belongs entirely to the method.

Subspace	Definition	Lives in	Dimension
Column space \(C(\mathbf{A})\)	span of the columns	\(\mathbb{R}^m\)	\(r\)
Null space \(N(\mathbf{A})\)	all \(\mathbf{x}\) with \(\mathbf{A}\mathbf{x} = \mathbf{0}\)	\(\mathbb{R}^n\)	\(n - r\)
Row space \(C(\mathbf{A}^\top)\)	span of the rows	\(\mathbb{R}^n\)	\(r\)
Left null space \(N(\mathbf{A}^\top)\)	all \(\mathbf{y}\) with \(\mathbf{A}^\top\mathbf{y} = \mathbf{0}\)	\(\mathbb{R}^m\)	\(m - r\)