9. Inner Product Spaces and Orthogonality

The vector space axioms of Chapter 6 describe addition and scaling, but say nothing about length, distance, angle, or perpendicularity. This chapter adds that geometry through an inner product, a generalized dot product that works for polynomials, functions, matrices, and complex vectors. It then turns perpendicularity into algorithms: orthonormal bases, projections, Gram-Schmidt, the QR factorization, orthogonal and unitary matrices, and least squares, which is the standard way to fit models to noisy data.

<a id="91-inner-products"></a>

9.1 Inner Products

This subchapter abstracts the dot product into a set of axioms, adapts them to complex scalars, and shows that length, distance, and angle follow from the axioms alone in any space where they hold.

Inner Product Axioms

Why it matters. In $\mathbb{R}^n$ the dot product gives lengths, angles, perpendicularity, and projections (Section 2.3). Many other spaces need the same tools. Examples are the best quadratic approximation to $e^x$, two audio signals that do not interfere, and a distance on data whose features use different units. These spaces either have no obvious dot product, or have one that gives the wrong geometry. The fix is to list the few properties that made the dot product work and accept any function with those properties. Every proof that uses only those properties then works in every such space.

Analogy. A vector space is a box of arrows that can be added and scaled, but it comes with no ruler and no protractor. An inner product is a measuring kit installed on that box. Different kits give different geometries on the same arrows, just as a road map and a travel-time map of the same city disagree about which places are "close".

Definition (real case). An inner product on a real vector space $V$ is a function that assigns to each pair $\mathbf{u}, \mathbf{v} \in V$ a real number $\langle \mathbf{u}, \mathbf{v} \rangle$, such that for all $\mathbf{u}, \mathbf{v}, \mathbf{w} \in V$ and all scalars $a, b$:

Axiom	Statement	What it guarantees
Symmetry	$\langle \mathbf{u}, \mathbf{v} \rangle = \langle \mathbf{v}, \mathbf{u} \rangle$	the angle from $\mathbf{u}$ to $\mathbf{v}$ equals the angle from $\mathbf{v}$ to $\mathbf{u}$
Linearity in the first argument	$\langle a\mathbf{u} + b\mathbf{w}, \mathbf{v} \rangle = a\langle \mathbf{u}, \mathbf{v} \rangle + b\langle \mathbf{w}, \mathbf{v} \rangle$	compatibility with the vector space operations
Positive definiteness	$\langle \mathbf{v}, \mathbf{v} \rangle \ge 0$, with equality only for $\mathbf{v} = \mathbf{0}$	every nonzero vector has a positive length

A vector space together with an inner product is an inner product space. A finite-dimensional real inner product space is also called a Euclidean space.

Immediate consequences.

• Bilinearity. Symmetry turns linearity in the first argument into linearity in the second: $\langle \mathbf{u}, a\mathbf{v} + b\mathbf{w} \rangle = a\langle \mathbf{u}, \mathbf{v} \rangle + b\langle \mathbf{u}, \mathbf{w} \rangle$.
• The zero vector. $\langle \mathbf{0}, \mathbf{v} \rangle = \langle 0 \cdot \mathbf{0}, \mathbf{v} \rangle = 0 \cdot \langle \mathbf{0}, \mathbf{v} \rangle = 0$ for every $\mathbf{v}$.
• Vectors are determined by their inner products. If $\langle \mathbf{u}, \mathbf{v} \rangle = \langle \mathbf{w}, \mathbf{v} \rangle$ for every $\mathbf{v}$, then $\mathbf{u} = \mathbf{w}$. Take $\mathbf{v} = \mathbf{u} - \mathbf{w}$: then $\langle \mathbf{u} - \mathbf{w}, \mathbf{u} - \mathbf{w} \rangle = 0$, so $\mathbf{u} - \mathbf{w} = \mathbf{0}$ by positive definiteness. This argument is used constantly to prove that two vectors are equal.

Examples.

Space	Inner product	Notes
$\mathbb{R}^n$	$\langle \mathbf{u}, \mathbf{v} \rangle = \mathbf{u} \cdot \mathbf{v} = \mathbf{u}^T\mathbf{v}$	the standard (Euclidean) inner product
$\mathbb{R}^n$	$\sum_i w_i u_i v_i$ with all $w_i > 0$	weighted dot product (later in this subchapter)
$\mathbb{R}^{m \times n}$	$\langle A, B \rangle = \operatorname{tr}(A^TB) = \sum_{i,j} a_{ij}b_{ij}$	Frobenius inner product: the dot product of the matrices read as long vectors
$C[a, b]$ (continuous functions)	$\langle f, g \rangle = \int_a^b f(x)g(x),dx$	the $L^2$ inner product
$\mathbb{P}_n$	$\sum_{i=0}^{n} p(x_i)q(x_i)$ for distinct $x_0, \dots, x_n$	the dot product of the value vectors

Non-examples. Each of these fails one axiom, and the failure is instructive.

Candidate	Fails	Witness
$u_1v_1$ on $\mathbb{R}^2$	positive definiteness	$\mathbf{v} = (0, 1) \ne \mathbf{0}$, yet $\langle \mathbf{v}, \mathbf{v} \rangle = 0$
$u_1v_1 - u_2v_2$ on $\mathbb{R}^2$	positive definiteness	$\langle (0, 1), (0, 1) \rangle = -1$; this is the Minkowski form of special relativity
$\lvert \mathbf{u} \cdot \mathbf{v} \rvert$ on $\mathbb{R}^n$	linearity	$\langle -\mathbf{u}, \mathbf{u} \rangle \ne -\langle \mathbf{u}, \mathbf{u} \rangle$
$p(0)q(0) + p(1)q(1)$ on $\mathbb{P}_2$	positive definiteness	$p(x) = x(x - 1) \ne 0$ has $\langle p, p \rangle = 0$: two points cannot pin down a quadratic

Every inner product on $\mathbb{R}^n$ is a matrix. Write $\mathbf{u} = \sum_i u_i\mathbf{e}_i$ and $\mathbf{v} = \sum_j v_j\mathbf{e}_j$ and expand using bilinearity:

$$\langle \mathbf{u}, \mathbf{v} \rangle = \sum_{i,j} u_i v_j \langle \mathbf{e}_i, \mathbf{e}j \rangle = \mathbf{u}^T G \mathbf{v}, \qquad G{ij} = \langle \mathbf{e}_i, \mathbf{e}_j \rangle.$$

$G$ is the Gram matrix of the standard basis. The axioms translate into conditions on $G$: symmetry of $\langle\cdot,\cdot\rangle$ means $G = G^T$, and positive definiteness means $\mathbf{x}^TG\mathbf{x} > 0$ for all $\mathbf{x} \ne \mathbf{0}$. A matrix with both properties is called symmetric positive definite (SPD). So inner products on $\mathbb{R}^n$ correspond one-to-one to SPD matrices. Section 10.3 gives practical tests for positive definiteness. The simplest is that all eigenvalues of the symmetric matrix $G$ are positive. The same construction works in any finite-dimensional space once a basis is chosen.

import numpy as np

rng = np.random.default_rng(0)

def check_inner_product(G, trials=1000):
    """Test the axioms for <u, v> = u^T G v on R^n (random tests find violations, eigenvalues prove positivity)."""
    ip = lambda u, v: u @ G @ v
    n = G.shape[0]
    symmetric = linear = True
    for _ in range(trials):
        u, v, w = rng.standard_normal((3, n))
        a, b = rng.standard_normal(2)
        symmetric &= np.isclose(ip(u, v), ip(v, u))
        linear &= np.isclose(ip(a * u + b * w, v), a * ip(u, v) + b * ip(w, v))
    # For symmetric G, x^T G x > 0 for all x != 0 exactly when every eigenvalue is positive (Section 10.3)
    positive = bool(symmetric and np.all(np.linalg.eigvalsh(G) > 0))
    return bool(symmetric), bool(linear), positive

print(check_inner_product(np.eye(3)))                          # (True, True, True): the dot product
print(check_inner_product(np.array([[2., 1.], [1., 2.]])))     # (True, True, True)
print(check_inner_product(np.array([[1., 2.], [2., 1.]])))     # (True, True, False)
print(check_inner_product(np.array([[1., 1.], [0., 1.]])))     # (False, True, False)

# The failing candidate gives the nonzero vector (1, -1) a negative "squared length"
x = np.array([1., -1.])
print(x @ np.array([[1., 2.], [2., 1.]]) @ x)                  # -2.0

# Recover G from any inner product: G_ij = <e_i, e_j>
ip = lambda u, v: 2 * u[0] * v[0] + u[0] * v[1] + u[1] * v[0] + 2 * u[1] * v[1]
E = np.eye(2)
print(np.array([[ip(E[i], E[j]) for j in range(2)] for i in range(2)]))   # [[2. 1.]
                                                                         #  [1. 2.]]

Random testing can expose a violated axiom but can never prove that all of them hold. Linearity is automatic for any $\mathbf{u}^TG\mathbf{v}$, symmetry is the exact check $G = G^T$, and positivity is the eigenvalue test.

Key takeaways

• An inner product is a symmetric, bilinear, positive definite function of two vectors; an inner product space is a vector space together with one.
• Positive definiteness is the demanding axiom: it is what makes "length zero" mean "the zero vector".
• If $\langle \mathbf{u}, \mathbf{v} \rangle = \langle \mathbf{w}, \mathbf{v} \rangle$ for all $\mathbf{v}$, then $\mathbf{u} = \mathbf{w}$.
• Inner products on $\mathbb{R}^n$ are exactly $\mathbf{u}^TG\mathbf{v}$ with $G$ symmetric positive definite, where $G_{ij} = \langle \mathbf{e}_i, \mathbf{e}_j \rangle$.
• One vector space can carry many inner products, and each one defines a different geometry.

🧪 Practice

1. Decide which of these are inner products on $\mathbb{R}^2$, naming the failing axiom when one fails: (a) $3u_1v_1 + 5u_2v_2$; (b) $u_1v_2 + u_2v_1$; (c) $u_1v_1 + 2u_1v_2 + 2u_2v_1 + 5u_2v_2$; (d) $u_1^2v_1^2 + u_2^2v_2^2$.
2. Show that $\langle p, q \rangle = p(-1)q(-1) + p(0)q(0) + p(1)q(1)$ is an inner product on $\mathbb{P}_2$ but not on $\mathbb{P}_3$.
3. Prove that the Frobenius inner product satisfies $\operatorname{tr}(A^TB) = \sum_{i,j} a_{ij}b_{ij}$, and use it to show that $\langle A, A \rangle = 0$ only for $A = 0$.
4. Interview: A kernel method in machine learning (an SVM, or a Gaussian process) uses a similarity function $k(\mathbf{x}, \mathbf{y})$ as if it were an inner product in some feature space. What property must the matrix $K_{ij} = k(\mathbf{x}_i, \mathbf{x}_j)$ have for every finite set of points, and what breaks if it does not? Hint: compare $K$ with the Gram matrix $G$ above; the requirement is symmetric positive semidefinite, and without it some "squared distances" in feature space become negative and the optimization problem is no longer convex.

Real and Complex Inner Products

Why it matters. Quantum states, alternating currents, and Fourier coefficients are complex vectors. Copying the real dot product $\sum u_iv_i$ to $\mathbb{C}^n$ fails immediately: for $\mathbf{z} = (1, i)$ it gives $1 \cdot 1 + i \cdot i = 0$, so a nonzero vector would have length zero. The fix is the same conjugate that makes $\lvert z \rvert^2 = z\bar{z}$ real and nonnegative (Section 1.2). It changes symmetry into conjugate symmetry and gives the complex inner product a slightly different algebra.

Definition (complex case). An inner product on a complex vector space $V$ assigns a complex number $\langle \mathbf{u}, \mathbf{v} \rangle$ to each pair of vectors, such that

Axiom	Statement
Conjugate symmetry	$\langle \mathbf{u}, \mathbf{v} \rangle = \overline{\langle \mathbf{v}, \mathbf{u} \rangle}$
Linearity in the first argument	$\langle a\mathbf{u} + b\mathbf{w}, \mathbf{v} \rangle = a\langle \mathbf{u}, \mathbf{v} \rangle + b\langle \mathbf{w}, \mathbf{v} \rangle$
Positive definiteness	$\langle \mathbf{v}, \mathbf{v} \rangle \ge 0$, with equality only for $\mathbf{v} = \mathbf{0}$

Conjugate symmetry with $\mathbf{u} = \mathbf{v}$ gives $\langle \mathbf{v}, \mathbf{v} \rangle = \overline{\langle \mathbf{v}, \mathbf{v} \rangle}$. So $\langle \mathbf{v}, \mathbf{v} \rangle$ is automatically real, and the inequality in the third axiom makes sense. For real vectors and real scalars, the complex axioms reduce to the real ones.

Conjugate-linear in the second argument. Combining the first two axioms,

$$\langle \mathbf{u}, c\mathbf{v} \rangle = \overline{\langle c\mathbf{v}, \mathbf{u} \rangle} = \overline{c,\langle \mathbf{v}, \mathbf{u} \rangle} = \bar{c},\langle \mathbf{u}, \mathbf{v} \rangle.$$

A scalar comes out of the second slot conjugated. A function that is linear in one argument and conjugate-linear in the other is called sesquilinear ("one and a half times linear").

The standard inner product on $\mathbb{C}^n$.

$$\langle \mathbf{u}, \mathbf{v} \rangle = \sum_{i=1}^n u_i\overline{v_i} = \mathbf{v}^*\mathbf{u},$$

where $\mathbf{v}^* = \overline{\mathbf{v}}^T$ is the conjugate transpose (Section 4.2). Then $\langle \mathbf{v}, \mathbf{v} \rangle = \sum \lvert v_i \rvert^2 > 0$ for $\mathbf{v} \ne \mathbf{0}$. As in the real case, every inner product on $\mathbb{C}^n$ has the form $\mathbf{v}^G\mathbf{u}$, with $G$ Hermitian positive definite ($G^ = G$ and $\mathbf{x}^*G\mathbf{x} > 0$ for $\mathbf{x} \ne \mathbf{0}$).

Worked example. Let $\mathbf{u} = (1 + i, 2)$ and $\mathbf{v} = (i, 1 - i)$.

$$\langle \mathbf{u}, \mathbf{v} \rangle = (1 + i)\overline{i} + 2,\overline{(1 - i)} = (1 + i)(-i) + 2(1 + i) = (1 - i) + (2 + 2i) = 3 + i.$$

Swapping the arguments gives $\langle \mathbf{v}, \mathbf{u} \rangle = i(1 - i) + (1 - i) \cdot 2 = (1 + i) + (2 - 2i) = 3 - i$, the conjugate. Also $\langle \mathbf{u}, \mathbf{u} \rangle = \lvert 1 + i \rvert^2 + \lvert 2 \rvert^2 = 2 + 4 = 6$.

Two conventions. Texts disagree about which argument is conjugated, and libraries follow one or the other. Check the convention before copying a formula.

	Mathematics (this guide)	Physics (bra-ket notation)
Linear in	first argument	second argument
Standard inner product	$\langle \mathbf{u}, \mathbf{v} \rangle = \mathbf{v}^*\mathbf{u}$	$\langle \phi \vert \psi \rangle = \boldsymbol{\phi}^*\boldsymbol{\psi}$
Scalar out of the "other" slot	$\langle \mathbf{u}, c\mathbf{v} \rangle = \bar{c}\langle \mathbf{u}, \mathbf{v} \rangle$	$\langle c\phi \vert \psi \rangle = \bar{c}\langle \phi \vert \psi \rangle$
NumPy	np.vdot(v, u)	np.vdot(phi, psi)

np.vdot(a, b) computes $\sum \overline{a_i},b_i$, conjugating its first argument. The operator @ and np.dot conjugate nothing.

Why not just drop the conjugate, or keep only the real part? Without the conjugate, positivity fails, as $(1, i)$ shows. Keeping only $\operatorname{Re}\langle \mathbf{u}, \mathbf{v} \rangle$ gives a valid real inner product: identifying $\mathbb{C}^n$ with $\mathbb{R}^{2n}$, it is the ordinary dot product of the real and imaginary parts. But it discards the imaginary part, which carries phase information. For example, $\mathbf{u} = (1)$ and $\mathbf{v} = (i)$ have $\operatorname{Re}\langle \mathbf{u}, \mathbf{v} \rangle = 0$, yet $\mathbf{v} = i\mathbf{u}$ is a complex multiple of $\mathbf{u}$, so they span the same complex line.

import numpy as np

def inner(u, v):
    """Standard complex inner product <u, v> = sum u_i conj(v_i) = v^* u (linear in u)."""
    return np.vdot(v, u)              # np.vdot conjugates its FIRST argument

u = np.array([1 + 1j, 2])
v = np.array([1j, 1 - 1j])
print(inner(u, v), inner(v, u))      # (3+1j) (3-1j): conjugate symmetry
print(inner(u, u))                   # (6+0j): always real and >= 0

# Linear in the first slot, conjugate-linear in the second
c = 2 - 3j
print(np.isclose(inner(c * u, v), c * inner(u, v)))            # True
print(np.isclose(inner(u, c * v), np.conj(c) * inner(u, v)))   # True

# Pitfall: @ and np.dot do NOT conjugate
z = np.array([1, 1j])
print(z @ z, inner(z, z))            # 0j (2+0j): without the conjugate, a nonzero vector has "length" 0

# A Hermitian positive definite G gives another inner product <u, v>_G = v^* G u
G = np.array([[2, 1j], [-1j, 2]])
print(np.allclose(G, G.conj().T), np.linalg.eigvalsh(G))       # True [1. 3.]
ipG = lambda a, b: b.conj() @ G @ a
print(np.isclose(ipG(u, v), np.conj(ipG(v, u))), ipG(z, z).real > 0)   # True True

Key takeaways

• Complex inner products replace symmetry with conjugate symmetry, $\langle \mathbf{u}, \mathbf{v} \rangle = \overline{\langle \mathbf{v}, \mathbf{u} \rangle}$, which forces $\langle \mathbf{v}, \mathbf{v} \rangle$ to be real.
• They are linear in one argument and conjugate-linear in the other (sesquilinear); this guide is linear in the first argument.
• The standard inner product on $\mathbb{C}^n$ is $\langle \mathbf{u}, \mathbf{v} \rangle = \mathbf{v}^*\mathbf{u}$; general ones are $\mathbf{v}^*G\mathbf{u}$ with $G$ Hermitian positive definite.
• Mathematics and physics conjugate opposite arguments; in NumPy, np.vdot conjugates the first argument and @ conjugates nothing.

🧪 Practice

1. For $\mathbf{u} = (2i, 1)$ and $\mathbf{v} = (1 - i, 3i)$, compute $\langle \mathbf{u}, \mathbf{v} \rangle$, $\langle \mathbf{v}, \mathbf{u} \rangle$, $\lVert \mathbf{u} \rVert$, and $\lVert \mathbf{v} \rVert$ by hand, then check them with np.vdot.
2. Prove that $\langle \mathbf{u}, a\mathbf{v} + b\mathbf{w} \rangle = \bar{a}\langle \mathbf{u}, \mathbf{v} \rangle + \bar{b}\langle \mathbf{u}, \mathbf{w} \rangle$ using only the three axioms.
3. Show that $\langle \mathbf{u}, \mathbf{v} \rangle = \mathbf{v}^G\mathbf{u}$ satisfies conjugate symmetry exactly when $G^ = G$.
4. Interview: A colleague's signal-processing code computes the "correlation" of two complex baseband signals as np.sum(x * y) and reports a value near zero for two identical signals. What is wrong, and what should the code compute instead? Hint: compare $\sum x_iy_i$ with $\sum x_i\overline{y_i}$ for $\mathbf{y} = \mathbf{x} = (1, i, -1, -i)$.

Induced Norms and Distance

Why it matters. Once an inner product is available, length and distance follow from it, and the familiar theorems of plane geometry carry over with the same proofs. These include Pythagoras, the Cauchy-Schwarz inequality, the triangle inequality, and the parallelogram law. The distance between two functions then means something precise, for example the root-mean-square error of an approximation.

Definitions. In an inner product space, the induced norm and distance are

$$\lVert \mathbf{v} \rVert = \sqrt{\langle \mathbf{v}, \mathbf{v} \rangle}, \qquad d(\mathbf{u}, \mathbf{v}) = \lVert \mathbf{u} - \mathbf{v} \rVert.$$

A vector with $\lVert \mathbf{u} \rVert = 1$ is a unit vector, and $\mathbf{v}/\lVert \mathbf{v} \rVert$ is the normalization of $\mathbf{v} \ne \mathbf{0}$.

Expanding a squared norm. Everything below rests on one identity. It follows from linearity in both arguments, with conjugates in the complex case:

$$\lVert \mathbf{u} + \mathbf{v} \rVert^2 = \langle \mathbf{u} + \mathbf{v}, \mathbf{u} + \mathbf{v} \rangle = \lVert \mathbf{u} \rVert^2 + 2\operatorname{Re}\langle \mathbf{u}, \mathbf{v} \rangle + \lVert \mathbf{v} \rVert^2.$$

In the real case $\operatorname{Re}\langle \mathbf{u}, \mathbf{v} \rangle$ is just $\langle \mathbf{u}, \mathbf{v} \rangle$. In particular, if $\langle \mathbf{u}, \mathbf{v} \rangle = 0$, then $\lVert \mathbf{u} + \mathbf{v} \rVert^2 = \lVert \mathbf{u} \rVert^2 + \lVert \mathbf{v} \rVert^2$. This is the Pythagorean theorem in any inner product space.

Cauchy-Schwarz inequality. For all $\mathbf{u}, \mathbf{v}$,

$$\lvert \langle \mathbf{u}, \mathbf{v} \rangle \rvert \le \lVert \mathbf{u} \rVert,\lVert \mathbf{v} \rVert,$$

with equality exactly when $\mathbf{u}$ and $\mathbf{v}$ are linearly dependent.

Proof. If $\mathbf{v} = \mathbf{0}$ both sides are $0$. Otherwise let $\mathbf{p} = \dfrac{\langle \mathbf{u}, \mathbf{v} \rangle}{\langle \mathbf{v}, \mathbf{v} \rangle}\mathbf{v}$, the "shadow" of $\mathbf{u}$ on $\mathbf{v}$. Then

$$\langle \mathbf{u} - \mathbf{p}, \mathbf{v} \rangle = \langle \mathbf{u}, \mathbf{v} \rangle - \frac{\langle \mathbf{u}, \mathbf{v} \rangle}{\langle \mathbf{v}, \mathbf{v} \rangle}\langle \mathbf{v}, \mathbf{v} \rangle = 0,$$

so $\mathbf{u} - \mathbf{p}$ is orthogonal to $\mathbf{v}$, and hence to its multiple $\mathbf{p}$. Pythagoras applied to $\mathbf{u} = \mathbf{p} + (\mathbf{u} - \mathbf{p})$ gives

$$\lVert \mathbf{u} \rVert^2 = \lVert \mathbf{p} \rVert^2 + \lVert \mathbf{u} - \mathbf{p} \rVert^2 \ge \lVert \mathbf{p} \rVert^2 = \frac{\lvert \langle \mathbf{u}, \mathbf{v} \rangle \rvert^2}{\lVert \mathbf{v} \rVert^2}.$$

Multiplying by $\lVert \mathbf{v} \rVert^2$ gives the inequality. Equality holds exactly when $\mathbf{u} - \mathbf{p} = \mathbf{0}$, that is, when $\mathbf{u}$ is a multiple of $\mathbf{v}$. $\blacksquare$

The proof never used coordinates, so it applies to every inner product space at once. For functions on $[a, b]$, for instance, it says

$$\left(\int_a^b fg,dx\right)^2 \le \int_a^b f^2,dx \int_a^b g^2,dx.$$

Norm properties. The induced norm satisfies the three defining properties of a norm (Section 11.3):

1. $\lVert \mathbf{v} \rVert \ge 0$, with equality only for $\mathbf{v} = \mathbf{0}$;
2. $\lVert c\mathbf{v} \rVert = \lvert c \rvert,\lVert \mathbf{v} \rVert$, because $\langle c\mathbf{v}, c\mathbf{v} \rangle = c\bar{c}\langle \mathbf{v}, \mathbf{v} \rangle$;
3. triangle inequality $\lVert \mathbf{u} + \mathbf{v} \rVert \le \lVert \mathbf{u} \rVert + \lVert \mathbf{v} \rVert$, because $\lVert \mathbf{u} + \mathbf{v} \rVert^2 \le \lVert \mathbf{u} \rVert^2 + 2\lvert \langle \mathbf{u}, \mathbf{v} \rangle \rvert + \lVert \mathbf{v} \rVert^2 \le (\lVert \mathbf{u} \rVert + \lVert \mathbf{v} \rVert)^2$ by Cauchy-Schwarz.

The distance therefore behaves as a distance should: $d(\mathbf{u}, \mathbf{v}) = d(\mathbf{v}, \mathbf{u})$, $d(\mathbf{u}, \mathbf{v}) = 0$ only when $\mathbf{u} = \mathbf{v}$, and $d(\mathbf{u}, \mathbf{w}) \le d(\mathbf{u}, \mathbf{v}) + d(\mathbf{v}, \mathbf{w})$.

Parallelogram law and polarization. Adding the expansions of $\lVert \mathbf{u} + \mathbf{v} \rVert^2$ and $\lVert \mathbf{u} - \mathbf{v} \rVert^2$ cancels the cross terms:

$$\lVert \mathbf{u} + \mathbf{v} \rVert^2 + \lVert \mathbf{u} - \mathbf{v} \rVert^2 = 2\lVert \mathbf{u} \rVert^2 + 2\lVert \mathbf{v} \rVert^2.$$

Geometrically, the squared diagonals of a parallelogram add up to the squared sides. Subtracting the expansions instead isolates the inner product. This is the polarization identity:

$$\text{real: } \langle \mathbf{u}, \mathbf{v} \rangle = \tfrac{1}{4}\left(\lVert \mathbf{u} + \mathbf{v} \rVert^2 - \lVert \mathbf{u} - \mathbf{v} \rVert^2\right), \qquad \text{complex: } \langle \mathbf{u}, \mathbf{v} \rangle = \tfrac{1}{4}\sum_{k=0}^{3} i^k,\lVert \mathbf{u} + i^k\mathbf{v} \rVert^2.$$

So the norm alone determines the inner product. A theorem of Jordan and von Neumann says that a norm comes from some inner product exactly when it satisfies the parallelogram law. This test settles which of the common norms have an underlying geometry of angles:

Norm on $\mathbb{R}^n$	Formula	Parallelogram law	Induced by an inner product?
Euclidean $\ell^2$	$\sqrt{\sum v_i^2}$	holds	yes, the dot product
$G$-norm	$\sqrt{\mathbf{v}^TG\mathbf{v}}$, $G$ SPD	holds	yes, $\mathbf{u}^TG\mathbf{v}$
Manhattan $\ell^1$	$\sum \lvert v_i \rvert$	fails for $\mathbf{e}_1, \mathbf{e}_2$: $8 \ne 4$	no
Maximum $\ell^\infty$	$\max \lvert v_i \rvert$	fails for $\mathbf{e}_1, \mathbf{e}_2$: $2 \ne 4$	no

   Parallelogram law in the l2 norm:            The same vectors in the l1 norm:

          u + v                                   ||e1 + e2||_1 = 2,  ||e1 - e2||_1 = 2
       .--------->                                ||e1||_1 = ||e2||_1 = 1
   v  /         /                                 2^2 + 2^2 = 8,  but  2(1) + 2(1) = 4
     /  u - v  /                                  -> no inner product produces the l1 norm,
    /  <----  /                                      so l1 has no notion of angle
   .--------->
        u
   |u+v|^2 + |u-v|^2 = 2 + 2 = 2|u|^2 + 2|v|^2    (u = e1, v = e2)

import numpy as np

def parallelogram_gap(norm, u, v):
    """||u+v||^2 + ||u-v||^2 - 2||u||^2 - 2||v||^2: zero for every u, v iff the norm comes from an inner product."""
    return norm(u + v)**2 + norm(u - v)**2 - 2 * norm(u)**2 - 2 * norm(v)**2

u, v = np.array([1., 0.]), np.array([0., 1.])
G = np.array([[2., 1.], [1., 2.]])
norms = {
    "l2":      lambda x: np.linalg.norm(x),
    "G-norm":  lambda x: np.sqrt(x @ G @ x),
    "l1":      lambda x: np.linalg.norm(x, 1),
    "l_inf":   lambda x: np.linalg.norm(x, np.inf),
}
for name, nrm in norms.items():
    print(f"{name:7s} {np.round(parallelogram_gap(nrm, u, v), 10) + 0.0:+.1f}")   # + 0.0 turns -0.0 into 0.0
# l2      +0.0
# G-norm  +0.0
# l1      +4.0
# l_inf   -2.0

# Polarization recovers the inner product from the norm alone (real case)
rng = np.random.default_rng(1)
x, y = rng.standard_normal((2, 2))
nG = norms["G-norm"]
print(np.isclose((nG(x + y)**2 - nG(x - y)**2) / 4, x @ G @ y))     # True

# Complex polarization: <u, v> = (1/4) sum_k i^k ||u + i^k v||^2  (inner product linear in u)
a, b = rng.standard_normal((2, 3)) + 1j * rng.standard_normal((2, 3))
pol = sum(1j**k * np.linalg.norm(a + 1j**k * b)**2 for k in range(4)) / 4
print(np.isclose(pol, np.vdot(b, a)))                                # True

# Cauchy-Schwarz and the triangle inequality hold for the G-norm too
print(abs(x @ G @ y) <= nG(x) * nG(y), nG(x + y) <= nG(x) + nG(y))    # True True

Key takeaways

• Every inner product induces a norm $\lVert \mathbf{v} \rVert = \sqrt{\langle \mathbf{v}, \mathbf{v} \rangle}$ and a distance $d(\mathbf{u}, \mathbf{v}) = \lVert \mathbf{u} - \mathbf{v} \rVert$.
• Cauchy-Schwarz, $\lvert \langle \mathbf{u}, \mathbf{v} \rangle \rvert \le \lVert \mathbf{u} \rVert\lVert \mathbf{v} \rVert$, holds in every inner product space; the proof uses only the axioms and Pythagoras.
• The triangle inequality follows from Cauchy-Schwarz.
• The parallelogram law characterizes the norms that come from inner products; $\ell^1$ and $\ell^\infty$ fail it.
• Polarization recovers the inner product from the norm, so a map that preserves norms also preserves inner products (Section 9.4).

🧪 Practice

1. In $C[0, 1]$ with $\langle f, g \rangle = \int_0^1 fg,dx$, compute $\lVert x \rVert$, $\lVert x^2 \rVert$, and the distance $d(x, x^2)$.
2. Use Cauchy-Schwarz in $C[0, 1]$ to prove $\left(\int_0^1 f,dx\right)^2 \le \int_0^1 f^2,dx$ for every continuous $f$. When does equality hold?
3. Verify the parallelogram law for the Frobenius norm on $2 \times 2$ matrices, and show that it fails for the matrix norm $\max_{i,j} \lvert a_{ij} \rvert$.
4. Interview: You are given a black-box function that returns $\lVert \mathbf{x} \rVert$ for an unknown inner product on $\mathbb{R}^n$. How many calls do you need to reconstruct the full inner product, and how? Hint: polarization recovers each $G_{ij} = \langle \mathbf{e}_i, \mathbf{e}_j \rangle$ from norms of $\mathbf{e}_i$, $\mathbf{e}_j$, and $\mathbf{e}_i + \mathbf{e}_j$; count the distinct norms you need.

Weighted and Function Space Inner Products

Why it matters. The standard dot product treats every coordinate equally, which is often wrong. A precise sensor should count more than a noisy one. Height in meters and income in dollars should not be compared on the same raw scale. A continuous signal has infinitely many "coordinates". Choosing the inner product means choosing what "close" means, and that choice determines which approximation is "best" in every projection and least squares problem later in this chapter.

Weighted inner products on $\mathbb{R}^n$. For positive weights $w_1, \dots, w_n$,

$$\langle \mathbf{u}, \mathbf{v} \rangle_w = \sum_{i=1}^n w_i u_i v_i = \mathbf{u}^TW\mathbf{v}, \qquad W = \operatorname{diag}(w_1, \dots, w_n).$$

All weights must be positive: a zero weight lets a nonzero vector have length $0$, and a negative one lets a squared length be negative. More generally, any SPD matrix $G$ gives $\langle \mathbf{u}, \mathbf{v} \rangle_G = \mathbf{u}^TG\mathbf{v}$ (previous topics). Changing the inner product reshapes the unit circle into an ellipse:

   standard: x^2 + y^2 = 1            weighted, W = diag(4, 1): 4x^2 + y^2 = 1

            y                                     y
            |                                     |
        .-'''''-.                               .-+-.  (0, 1)
      /     |     \                            /  |  \
     |      |      |                          |   |   |
   --+------+------+-- x               -------+---+---+------- x
     |      |      |                   (-1/2, 0)  |   (1/2, 0)
      \     |     /                           |   |   |
        '-.....-'                              \  |  /
            |                                   '-+-'
                                                  |
   all unit vectors equally long         the x-direction now counts double:
                                         (1/2, 0) already has length 1

Typical choices of weights:

Application	Inner product	Effect
Measurements with variances $\sigma_i^2$	$w_i = 1/\sigma_i^2$	precise measurements count more (weighted least squares, Section 9.5)
Features on different scales	$w_i = 1/s_i^2$, $s_i$ = feature spread	same as standardizing the features first
Correlated data with covariance $\Sigma$	$G = \Sigma^{-1}$	Mahalanobis distance $\sqrt{(\mathbf{x} - \mathbf{y})^T\Sigma^{-1}(\mathbf{x} - \mathbf{y})}$, used for outlier detection
Finite element methods	$G$ = mass matrix	discrete version of a function inner product

Function space inner products. On $C[a, b]$, the space of continuous functions on $[a, b]$,

$$\langle f, g \rangle = \int_a^b f(x)g(x),dx, \qquad \lVert f - g \rVert = \sqrt{\int_a^b (f - g)^2,dx}.$$

Symmetry and linearity follow from properties of the integral. Positive definiteness needs continuity: if $f(x_0) \ne 0$ at some point, then $f^2 > 0$ on a small interval around $x_0$, so $\int f^2 > 0$. Without continuity the axiom fails: a function that is $1$ at a single point and $0$ elsewhere is nonzero but has "length" $0$. The space $L^2$ of Section 13.4 handles such functions by treating functions that agree except on a negligible set as equal.

Why this is the natural extension of the dot product. Sample $f$ and $g$ at $N$ equally spaced points of $[a, b]$, with spacing $\Delta x = (b - a)/N$. The weighted dot product $\sum_i f(x_i)g(x_i),\Delta x$ of the sample vectors is a Riemann sum, so as $N \to \infty$ it tends to $\int_a^b fg,dx$. A function is a vector with infinitely many coordinates, and the integral is the dot product of those coordinates.

   f sampled:  (f(x_1), f(x_2), ..., f(x_N))          <f, g>  ~  sum f(x_i) g(x_i) dx
   g sampled:  (g(x_1), g(x_2), ..., g(x_N))                     --> integral of f g  as N grows

Weight functions. A positive weight $w(x)$ gives $\langle f, g \rangle_w = \int_a^b f(x)g(x)w(x),dx$. It makes some regions of the interval count more than others. Orthogonalizing $1, x, x^2, \dots$ under different weights (Section 9.3) produces the classical orthogonal polynomial families:

Interval	Weight $w(x)$	Orthogonal polynomials	Typical use
$[-1, 1]$	$1$	Legendre	Gaussian quadrature, multipole expansions
$[-1, 1]$	$1/\sqrt{1 - x^2}$	Chebyshev	near-minimax approximation, spectral methods
$(-\infty, \infty)$	$e^{-x^2}$	Hermite	quantum harmonic oscillator, probability
$[0, \infty)$	$e^{-x}$	Laguerre	hydrogen atom, queueing

The interval changes the geometry. On $[-1, 1]$ the functions $1$ and $x$ are orthogonal, since $\int_{-1}^1 x,dx = 0$. On $[0, 1]$ they are not: $\langle 1, x \rangle = \tfrac{1}{2}$. The Gram matrix of $1, x, \dots, x^{n-1}$ on $[0, 1]$ is the Hilbert matrix $H_{ij} = \int_0^1 x^{i+j-2},dx = 1/(i + j - 1)$. Its huge condition number (Section 11.3) shows that high powers of $x$ are nearly parallel in this geometry. This is why stable polynomial fitting uses orthogonal polynomials rather than monomials.

import numpy as np
import sympy as sp

x = sp.symbols("x")

def ip(f, g, a=0, b=1, w=1):
    """Function inner product <f, g> = integral_a^b f(x) g(x) w(x) dx, computed exactly."""
    return sp.integrate(f * g * w, (x, a, b))

print(ip(1, x), ip(x, x), ip(1, x, -1, 1))          # 1/2 1/3 0: the interval changes the geometry
print(sp.sqrt(ip(x**2, x**2)))                      # sqrt(5)/5: the length of x^2 on [0, 1]

# Chebyshev weight on [-1, 1]: T_1 = x and T_2 = 2x^2 - 1 are orthogonal
w_cheb = 1 / sp.sqrt(1 - x**2)
print(ip(x, 2 * x**2 - 1, -1, 1, w_cheb), ip(1, 2 * x**2 - 1, -1, 1, w_cheb))    # 0 0

# The Gram matrix of 1, x, ..., x^(n-1) on [0, 1] is the Hilbert matrix H_ij = 1/(i + j - 1)
n = 4
H = sp.Matrix(n, n, lambda i, j: ip(x**i, x**j))
print(H)        # Matrix([[1, 1/2, 1/3, 1/4], [1/2, 1/3, 1/4, 1/5], [1/3, 1/4, 1/5, 1/6], [1/4, 1/5, 1/6, 1/7]])
for n in [4, 8, 12]:
    Hn = np.array([[1 / (i + j + 1) for j in range(n)] for i in range(n)])
    print(n, f"{np.linalg.cond(Hn):.1e}")
# 4 1.6e+04
# 8 1.5e+10
# 12 1.8e+16     the monomial basis is nearly dependent in this geometry

# A function inner product is the limit of a weighted dot product of samples
f, g = np.exp, np.cos
for N in [10, 100, 1000]:
    t = (np.arange(N) + 0.5) / N                   # midpoints of N equal subintervals of [0, 1]
    print(N, f(t) @ g(t) / N)                      # dot product of sample vectors, weight 1/N
# 10 1.3787817190027822
# 100 1.3780321912359899
# 1000 1.3780246893249128
print(sp.N(ip(sp.exp(x), sp.cos(x))))              # 1.37802461354736 (exact value (e(cos 1 + sin 1) - 1)/2)

Key takeaways

• Weighted inner products $\sum w_iu_iv_i$ (all $w_i > 0$) and $\mathbf{u}^TG\mathbf{v}$ ($G$ SPD) let some directions count more than others; unit circles become ellipses.
• Inverse variances and inverse covariances are the standard weights; the latter gives the Mahalanobis distance.
• $\langle f, g \rangle = \int_a^b fg,dx$ is the limit of a dot product of sample vectors; continuity is what makes it positive definite on $C[a, b]$.
• Weight functions and intervals change which functions are orthogonal and produce the classical orthogonal polynomials.
• Monomials are badly conditioned as a basis (Hilbert matrix), which motivates orthogonal bases.

🧪 Practice

1. With $\langle \mathbf{u}, \mathbf{v} \rangle = 4u_1v_1 + u_2v_2$ on $\mathbb{R}^2$, find all unit vectors, and find a nonzero vector orthogonal to $(1, 1)$.
2. In $C[-1, 1]$ with $\langle f, g \rangle = \int_{-1}^1 fg,dx$, show that every even function is orthogonal to every odd function.
3. Two points have the same Euclidean distance from the mean of a dataset with covariance $\Sigma = \begin{bmatrix} 4 & 0 \ 0 & 1 \end{bmatrix}$: $\mathbf{x}_1 = (2, 0)$ and $\mathbf{x}_2 = (0, 2)$, both measured from the mean. Compute their Mahalanobis distances and explain which one is more unusual.
4. Interview: You run $k$-nearest neighbors on features measured in very different units, such as age in years and income in dollars. Why do the results depend on the units, and how does choosing a weighted inner product fix it? Hint: in the standard dot product the income axis dominates every distance; weighting by $1/s_i^2$ is the same as standardizing, and $\Sigma^{-1}$ also removes correlations.

Angles and Orthogonality in General Spaces

Why it matters. "The angle between two functions" sounds meaningless, yet it is exactly what correlation coefficients, interference between signals, and the independence of Fourier modes measure. With Cauchy-Schwarz available, angles and perpendicularity extend to every inner product space. Orthogonality, the case $\langle \mathbf{u}, \mathbf{v} \rangle = 0$, is the most important case for the rest of the chapter.

Angle (real spaces). For nonzero $\mathbf{u}, \mathbf{v}$ in a real inner product space, Cauchy-Schwarz guarantees $-1 \le \dfrac{\langle \mathbf{u}, \mathbf{v} \rangle}{\lVert \mathbf{u} \rVert\lVert \mathbf{v} \rVert} \le 1$. So there is a unique $\theta \in [0, \pi]$ with

$$\cos\theta = \frac{\langle \mathbf{u}, \mathbf{v} \rangle}{\lVert \mathbf{u} \rVert,\lVert \mathbf{v} \rVert}.$$

This $\theta$ is the angle between $\mathbf{u}$ and $\mathbf{v}$. In complex spaces, $\langle \mathbf{u}, \mathbf{v} \rangle$ is a complex number, and the usual substitute is $\cos\theta = \lvert \langle \mathbf{u}, \mathbf{v} \rangle \rvert / (\lVert \mathbf{u} \rVert\lVert \mathbf{v} \rVert) \in [0, 1]$. It measures how close $\mathbf{u}$ is to the complex line through $\mathbf{v}$. In quantum mechanics, its square is the transition probability between two normalized states.

Orthogonality. Vectors $\mathbf{u}$ and $\mathbf{v}$ are orthogonal, written $\mathbf{u} \perp \mathbf{v}$, if $\langle \mathbf{u}, \mathbf{v} \rangle = 0$. This definition works over $\mathbb{R}$ and over $\mathbb{C}$, and it is symmetric because $\overline{0} = 0$. Two facts follow at once:

• $\mathbf{0}$ is orthogonal to every vector.
• $\mathbf{0}$ is the only vector orthogonal to itself, by positive definiteness. So a vector orthogonal to everything, or orthogonal to a spanning set, is $\mathbf{0}$.

Orthogonality depends on the inner product. The same pair of vectors can be perpendicular in one geometry and not in another:

Pair	Inner product	Inner product value	Angle
$(1, 1)$, $(1, -1)$	dot product	$0$	$90°$
$(1, 1)$, $(1, -1)$	$4u_1v_1 + u_2v_2$	$3$	$\arccos(3/5) \approx 53.13°$
$1$, $x$	$\int_{-1}^1$	$0$	$90°$
$1$, $x$	$\int_0^1$	$1/2$	$30°$

Important orthogonal families.

• On $[-\pi, \pi]$ with $\int_{-\pi}^{\pi}$: the functions $1, \cos x, \sin x, \cos 2x, \sin 2x, \dots$ are mutually orthogonal. For example, $\int_{-\pi}^{\pi}\sin(mx)\sin(nx),dx = 0$ for $m \ne n$. This fact underlies Fourier series (Section 9.5).
• On a symmetric interval $[-a, a]$: every even function is orthogonal to every odd function, because their product is odd.
• In $\mathbb{R}^{n \times n}$ with the Frobenius inner product: every symmetric matrix is orthogonal to every skew-symmetric matrix.
• For random variables with mean zero, $\langle X, Y \rangle = \mathbb{E}[XY]$ is the covariance. "Uncorrelated" means "orthogonal", and the correlation coefficient is $\cos\theta$. For data vectors, this is the cosine of the angle between the mean-centered vectors (Section 2.3).

Pythagoras, and a complex surprise. If $\mathbf{u} \perp \mathbf{v}$, then $\lVert \mathbf{u} + \mathbf{v} \rVert^2 = \lVert \mathbf{u} \rVert^2 + \lVert \mathbf{v} \rVert^2$ in every inner product space. In real spaces the converse also holds, since the expansion $\lVert \mathbf{u} + \mathbf{v} \rVert^2 = \lVert \mathbf{u} \rVert^2 + 2\langle \mathbf{u}, \mathbf{v} \rangle + \lVert \mathbf{v} \rVert^2$ forces $\langle \mathbf{u}, \mathbf{v} \rangle = 0$. In complex spaces the converse fails: the expansion only forces $\operatorname{Re}\langle \mathbf{u}, \mathbf{v} \rangle = 0$. In $\mathbb{C}^1$, $\mathbf{u} = 1$ and $\mathbf{v} = i$ satisfy Pythagoras, yet $\langle \mathbf{u}, \mathbf{v} \rangle = 1 \cdot \bar{i} = -i \ne 0$.

import numpy as np
import sympy as sp

def angle_deg(ip, u, v):
    """Angle between u and v for ANY real inner product ip(u, v)."""
    c = ip(u, v) / sp.sqrt(ip(u, u) * ip(v, v))
    return sp.nsimplify(sp.deg(sp.acos(c)))

x = sp.symbols("x")
L2 = lambda a, b: (lambda f, g: sp.integrate(f * g, (x, a, b)))   # L2 inner product on [a, b]

print(angle_deg(L2(0, 1), 1, x))                       # 30
print(angle_deg(L2(-1, 1), 1, x))                      # 90
print(angle_deg(L2(-sp.pi, sp.pi), sp.sin(x), sp.cos(x)))           # 90
print(angle_deg(L2(-sp.pi, sp.pi), sp.sin(2*x), sp.sin(3*x)))       # 90
print(sp.N(angle_deg(L2(-1, 1), x, x**3), 4))          # 23.58: cos(theta) = sqrt(21)/5

# Orthogonality depends on the inner product
G = sp.diag(4, 1)
dot = lambda a, b: (a.T * b)[0]
dotG = lambda a, b: (a.T * G * b)[0]
u, v = sp.Matrix([1, 1]), sp.Matrix([1, -1])
print(angle_deg(dot, u, v), sp.N(angle_deg(dotG, u, v), 4))         # 90 53.13

# Frobenius inner product <A, B> = trace(A^T B): symmetric and skew-symmetric matrices are orthogonal
rng = np.random.default_rng(0)
M = rng.standard_normal((3, 3))
S, K = (M + M.T) / 2, (M - M.T) / 2
print(np.isclose(np.trace(S.T @ K), 0))                # True

# Pythagoras: the converse fails over C
a, b = np.array([1.0 + 0j]), np.array([1j])
print(round(np.linalg.norm(a + b)**2, 12), np.linalg.norm(a)**2 + np.linalg.norm(b)**2, np.vdot(b, a))
# 2.0 2.0 -1j     ||a+b||^2 = ||a||^2 + ||b||^2, yet <a, b> = -i != 0

The angle of about $23.6°$ between $x$ and $x^3$ on $[-1, 1]$ quantifies something visible in their graphs: the two functions are nearly parallel. This is the function-space version of the ill-conditioning seen with the Hilbert matrix.

Key takeaways

• In a real inner product space, $\cos\theta = \langle \mathbf{u}, \mathbf{v} \rangle / (\lVert \mathbf{u} \rVert\lVert \mathbf{v} \rVert)$ defines an angle in $[0, \pi]$, because Cauchy-Schwarz keeps the ratio in $[-1, 1]$.
• $\mathbf{u} \perp \mathbf{v}$ means $\langle \mathbf{u}, \mathbf{v} \rangle = 0$; only $\mathbf{0}$ is orthogonal to itself.
• Whether two vectors are orthogonal depends on the chosen inner product, including the interval and weight for functions.
• Sines and cosines of different frequencies, even and odd functions, and symmetric and skew-symmetric matrices are standard orthogonal pairs.
• Pythagoras holds for orthogonal vectors in every inner product space; its converse holds over $\mathbb{R}$ but not over $\mathbb{C}$.

🧪 Practice

1. In $C[0, 1]$, find the angle between $x$ and $x^2$, and find a constant $c$ such that $x - c$ is orthogonal to $1$.
2. Prove that $\int_{-\pi}^{\pi}\cos(mx)\cos(nx),dx = 0$ for nonnegative integers $m \ne n$, using $\cos A\cos B = \tfrac{1}{2}[\cos(A - B) + \cos(A + B)]$.
3. Show that symmetric and skew-symmetric matrices are orthogonal under $\langle A, B \rangle = \operatorname{tr}(A^TB)$, using only properties of the trace and the transpose.
4. Interview: Two features in a dataset have a Pearson correlation of exactly $0$. Does that mean they are unrelated? Give an example. Hint: correlation is the cosine of the angle between centered vectors, so it detects only linear relationships; try $y = x^2$ with $x$ sampled symmetrically around $0$.

<a id="92-orthogonal-sets-and-bases"></a>

9.2 Orthogonal Sets and Bases

This subchapter shows why perpendicular vectors are the best building blocks: they are automatically independent, coordinates in them are single inner products, and they turn "find the closest point in a subspace" into a formula.

Orthogonal and Orthonormal Sets

Why it matters. Perpendicular directions do not interfere with each other. Walking north does not change how far east you are, so the two distances can be measured separately. A set of mutually perpendicular vectors has the same property: each vector carries information that none of the others carries, and computations break into independent one-dimensional pieces. Fourier analysis, wavelets, principal components, and the QR factorization all rely on this.

Analogy. An orchestra in which each instrument plays in its own frequency band can be recorded and mixed one instrument at a time, because turning one up does not change what you hear from the others. Orthogonal vectors are channels that do not leak into each other.

Definitions. A set ${\mathbf{v}_1, \dots, \mathbf{v}_k}$ in an inner product space is

• orthogonal if $\langle \mathbf{v}_i, \mathbf{v}_j \rangle = 0$ whenever $i \ne j$;
• orthonormal if, in addition, every $\lVert \mathbf{v}_i \rVert = 1$. Equivalently, $\langle \mathbf{v}_i, \mathbf{v}j \rangle = \delta{ij}$, which is $1$ if $i = j$ and $0$ otherwise (the Kronecker delta).

Any orthogonal set of nonzero vectors becomes orthonormal after dividing each vector by its norm. Normalizing does not change any of the right angles.

Theorem. An orthogonal set of nonzero vectors is linearly independent.

Proof. Suppose $c_1\mathbf{v}_1 + \cdots + c_k\mathbf{v}_k = \mathbf{0}$. Take the inner product of both sides with $\mathbf{v}_j$. Every term with $i \ne j$ vanishes, leaving $c_j\langle \mathbf{v}_j, \mathbf{v}_j \rangle = 0$. Since $\mathbf{v}_j \ne \mathbf{0}$, $\langle \mathbf{v}_j, \mathbf{v}_j \rangle > 0$, so $c_j = 0$. This holds for every $j$. $\blacksquare$

So in an $n$-dimensional space, an orthogonal set has at most $n$ nonzero vectors, and $n$ of them form a basis. Independence can be read off from inner products, with no row reduction. The word "nonzero" matters: ${\mathbf{e}_1, \mathbf{0}}$ is orthogonal but dependent.

Generalized Pythagoras. For an orthogonal set, all cross terms in the expansion of $\lVert \sum c_i\mathbf{v}_i \rVert^2$ vanish:

$$\left\lVert \sum_{i=1}^k c_i\mathbf{v}i \right\rVert^2 = \sum{i=1}^k \lvert c_i \rvert^2,\lVert \mathbf{v}_i \rVert^2.$$

Matrix test. Put vectors $\mathbf{q}_1, \dots, \mathbf{q}_k \in \mathbb{R}^m$ into the columns of $Q$ ($m \times k$). The entries of $Q^TQ$ are the dot products $\mathbf{q}_i^T\mathbf{q}_j$. This $Q^TQ$ is the Gram matrix of the columns. Therefore

$$\text{columns orthogonal} \iff Q^TQ \text{ is diagonal}, \qquad \text{columns orthonormal} \iff Q^TQ = I_k.$$

Over $\mathbb{C}$, replace $Q^T$ by $Q^*$. When $Q$ is tall ($k < m$), $Q^TQ = I_k$ does not imply $QQ^T = I_m$. Instead, $QQ^T$ is the projection onto the column space (later in this subchapter). Only a square $Q$ with orthonormal columns satisfies both equations. Section 9.4 studies those matrices.

Examples.

Space and inner product	Orthogonal set
$\mathbb{R}^n$, dot product	standard basis $\mathbf{e}_1, \dots, \mathbf{e}_n$ (orthonormal)
$\mathbb{R}^3$, dot product	$(1, 1, 1)$, $(1, -1, 0)$, $(1, 1, -2)$
$\mathbb{R}^4$, dot product	Haar vectors $(1, 1, 1, 1)$, $(1, 1, -1, -1)$, $(1, -1, 0, 0)$, $(0, 0, 1, -1)$
$C[-\pi, \pi]$, $\int_{-\pi}^{\pi}$	$1, \cos x, \sin x, \cos 2x, \sin 2x, \dots$
$\mathbb{P}2$, $\int{-1}^{1}$	$1$, $x$, $x^2 - \tfrac{1}{3}$ (Legendre, up to scaling)

   Haar vectors in R^4: averages and differences at two scales

   v1 = ( 1,  1,  1,  1)    overall level           ####  ####  ####  ####
   v2 = ( 1,  1, -1, -1)    left half vs right      ####  ####  ....  ....
   v3 = ( 1, -1,  0,  0)    detail in left half     ####  ....
   v4 = ( 0,  0,  1, -1)    detail in right half                ####  ....

   any two of them: matching entries multiply and cancel  ->  inner product 0

import numpy as np

# Haar vectors in R^4: averages and differences at two scales
V = np.array([[1,  1,  1,  0],
              [1,  1, -1,  0],
              [1, -1,  0,  1],
              [1, -1,  0, -1]], dtype=float)     # columns v1, v2, v3, v4

print(V.T @ V)                                   # Gram matrix: diagonal, so the columns are orthogonal
# [[4. 0. 0. 0.]
#  [0. 4. 0. 0.]
#  [0. 0. 2. 0.]
#  [0. 0. 0. 2.]]

Q = V / np.linalg.norm(V, axis=0)                # normalize each column
print(np.allclose(Q.T @ Q, np.eye(4)))           # True: orthonormal

# Orthogonal nonzero vectors are independent
print(np.linalg.matrix_rank(V))                  # 4

# Generalized Pythagoras: ||sum c_i v_i||^2 = sum c_i^2 ||v_i||^2
c = np.array([3., -1., 2., 5.])
print(np.linalg.norm(V @ c)**2, np.sum(c**2 * np.sum(V**2, axis=0)))   # 98.0 98.0

# A tall matrix with orthonormal columns: Q^T Q = I, but Q Q^T is NOT I
Q2 = Q[:, :2]
print(np.allclose(Q2.T @ Q2, np.eye(2)), np.allclose(Q2 @ Q2.T, np.eye(4)))    # True False
print(Q2 @ Q2.T)                                 # a projection matrix (Orthogonal Projection, below)
# [[0.5 0.5 0.  0. ]
#  [0.5 0.5 0.  0. ]
#  [0.  0.  0.5 0.5]
#  [0.  0.  0.5 0.5]]

Key takeaways

• Orthogonal: pairwise inner products are zero. Orthonormal: additionally, every vector has length 1, so $\langle \mathbf{v}_i, \mathbf{v}j \rangle = \delta{ij}$.
• Orthogonal sets of nonzero vectors are automatically linearly independent.
• For orthogonal vectors, $\lVert \sum c_i\mathbf{v}_i \rVert^2 = \sum \lvert c_i \rvert^2\lVert \mathbf{v}_i \rVert^2$.
• The columns of $Q$ are orthonormal exactly when $Q^TQ = I$ (or $Q^*Q = I$); for tall $Q$, $QQ^T \ne I$.

🧪 Practice

1. Check that $(1, 2, 2)$, $(2, 1, -2)$, $(2, -2, 1)$ are mutually orthogonal, normalize them, and verify $Q^TQ = I$ for the resulting matrix.
2. Find all vectors $(a, b, c)$ that are orthogonal to both $(1, 1, 1)$ and $(1, -1, 0)$, and explain why the answer is a line.
3. Show that $1$, $x$, and $x^2 - \tfrac{1}{3}$ are mutually orthogonal in $C[-1, 1]$ with $\langle f, g \rangle = \int_{-1}^1 fg,dx$. Are they orthogonal in $C[0, 1]$?
4. Interview: You compute $Q^TQ$ for a matrix whose columns should be orthonormal and find off-diagonal entries of size $10^{-6}$ in double precision. Is that acceptable, and what might have caused it? Hint: compare $10^{-6}$ with machine epsilon ($\approx 10^{-16}$); loss of orthogonality on that scale points to classical Gram-Schmidt on ill-conditioned input (Section 9.3).

Orthonormal Bases and Coordinates

Why it matters. In a general basis, finding the coordinates of a vector means solving a linear system (Section 6.5), which costs $O(n^3)$ operations. In an orthonormal basis, each coordinate is a single inner product. Lengths and angles can then be computed directly from coordinates, as if the basis were the standard one. This is why Fourier coefficients, wavelet coefficients, and principal component scores are all computed as inner products.

Analogy. To locate a point on a map with a perpendicular grid, you drop a perpendicular to each axis and read off a number. Nothing about one axis affects the reading on the other. On a map with a skewed grid, the two readings depend on each other and must be solved for together.

Theorem (coordinates in an orthonormal basis). If ${\mathbf{u}_1, \dots, \mathbf{u}_n}$ is an orthonormal basis of $V$, then every $\mathbf{v} \in V$ satisfies

$$\mathbf{v} = \sum_{i=1}^n \langle \mathbf{v}, \mathbf{u}_i \rangle,\mathbf{u}_i.$$

Proof. Write $\mathbf{v} = \sum_j c_j\mathbf{u}_j$ (possible since it is a basis) and take the inner product with $\mathbf{u}_i$: $\langle \mathbf{v}, \mathbf{u}_i \rangle = \sum_j c_j\langle \mathbf{u}_j, \mathbf{u}_i \rangle = c_i$. $\blacksquare$

For an orthogonal basis that is not normalized, the same computation gives $c_i = \dfrac{\langle \mathbf{v}, \mathbf{v}_i \rangle}{\langle \mathbf{v}_i, \mathbf{v}_i \rangle}$. Over $\mathbb{C}$ the order of the arguments matters: with the convention of Section 9.1, the coefficient is $\langle \mathbf{v}, \mathbf{u}_i \rangle$, not $\langle \mathbf{u}_i, \mathbf{v} \rangle$ (which is its conjugate).

   Orthonormal basis: each coordinate is a perpendicular "shadow"

            u2
            ^
            |      v
   c2 ------+-----*
            |    /:
            |   / :
            |  /  :         c1 = <v, u1>    (shadow on u1)
            | /   :         c2 = <v, u2>    (shadow on u2)
            |/    :
            O-----+-------> u1
                  c1

Matrix form. In $\mathbb{R}^n$, let $Q$ have the orthonormal basis vectors as columns. The coordinate vector solves $Q\mathbf{c} = \mathbf{v}$, and since $Q^TQ = I$,

$$[\mathbf{v}]_{\mathcal{B}} = \mathbf{c} = Q^T\mathbf{v}, \qquad \mathbf{v} = Q\mathbf{c}.$$

The change-of-basis matrix is inverted by transposing it, with no elimination required.

Parseval's identity. In an orthonormal basis, inner products of vectors equal dot products of their coordinate vectors. If $\mathbf{v} = \sum c_i\mathbf{u}_i$ and $\mathbf{w} = \sum d_i\mathbf{u}_i$, then

$$\langle \mathbf{v}, \mathbf{w} \rangle = \sum_{i=1}^n c_i\overline{d_i}, \qquad \lVert \mathbf{v} \rVert^2 = \sum_{i=1}^n \lvert \langle \mathbf{v}, \mathbf{u}_i \rangle \rvert^2.$$

In a non-orthonormal basis, the formula needs the Gram matrix: $\langle \mathbf{v}, \mathbf{w} \rangle = \mathbf{d}^*G\mathbf{c}$ with $G_{ji} = \langle \mathbf{b}_i, \mathbf{b}_j \rangle$. A consequence is that every $n$-dimensional real inner product space looks exactly like $\mathbb{R}^n$ with the dot product once an orthonormal basis is chosen, and every complex one looks like $\mathbb{C}^n$. Gram-Schmidt (Section 9.3) shows that an orthonormal basis always exists.

Bessel's inequality. If ${\mathbf{u}_1, \dots, \mathbf{u}_k}$ is orthonormal but not necessarily a basis, then

$$\sum_{i=1}^k \lvert \langle \mathbf{v}, \mathbf{u}_i \rangle \rvert^2 \le \lVert \mathbf{v} \rVert^2.$$

The coefficients capture part of the "energy" of $\mathbf{v}$, and never more than all of it. Equality holds exactly when $\mathbf{v}$ lies in the span of the $\mathbf{u}_i$ (proof in Orthogonal Decomposition Theorem below).

Worked example. Take the orthonormal basis $\mathbf{q}_1 = \tfrac{1}{\sqrt{3}}(1, 1, 1)$, $\mathbf{q}_2 = \tfrac{1}{\sqrt{2}}(1, -1, 0)$, $\mathbf{q}_3 = \tfrac{1}{\sqrt{6}}(1, 1, -2)$ and $\mathbf{v} = (3, 1, 2)$:

$$c_1 = \frac{3 + 1 + 2}{\sqrt{3}} = 2\sqrt{3}, \qquad c_2 = \frac{3 - 1}{\sqrt{2}} = \sqrt{2}, \qquad c_3 = \frac{3 + 1 - 4}{\sqrt{6}} = 0.$$

Parseval check: $c_1^2 + c_2^2 + c_3^2 = 12 + 2 + 0 = 14 = 3^2 + 1^2 + 2^2$.

Application: wavelet compression. The orthonormal Haar basis splits a signal into averages and differences at several scales. Smooth regions produce many zero or tiny coefficients. Keeping only the largest coefficients compresses the signal, and Parseval gives the exact squared error: it equals the sum of the squares of the discarded coefficients.

import numpy as np

Q = np.column_stack([np.array([1, 1, 1]) / np.sqrt(3),
                     np.array([1, -1, 0]) / np.sqrt(2),
                     np.array([1, 1, -2]) / np.sqrt(6)])
v = np.array([3., 1., 2.])

c_solve = np.linalg.solve(Q, v)        # general basis: solve Q c = v, O(n^3)
c_dot = Q.T @ v                        # orthonormal basis: c_i = <v, q_i>, O(n^2)
print(np.allclose(c_solve, c_dot), np.round(c_dot, 4))     # True [3.4641 1.4142 0.    ]
print(np.allclose(Q @ c_dot, v))                           # True: v = sum c_i q_i

# Parseval: lengths and inner products can be computed from coordinates
w = np.array([0., 2., -1.])
d = Q.T @ w
print(np.isclose(v @ v, c_dot @ c_dot), np.isclose(v @ w, c_dot @ d))    # True True

# Haar wavelet basis of R^8 (orthonormal), built recursively from averages and differences
def haar_rows(n):
    if n == 1:
        return np.array([[1.0]])
    coarse = np.kron(haar_rows(n // 2), [1, 1])      # each coarser pattern, stretched to length n
    fine = np.kron(np.eye(n // 2), [1, -1])          # differences of neighboring pairs
    return np.vstack([coarse, fine])

H = haar_rows(8)
Q8 = (H / np.linalg.norm(H, axis=1, keepdims=True)).T    # columns = orthonormal basis vectors
print(np.allclose(Q8.T @ Q8, np.eye(8)))         # True

signal = np.array([4., 4., 5., 5., 9., 9., 8., 0.])
coef = Q8.T @ signal                             # coordinates: one inner product each
print(np.round(coef, 3))                         # [15.556 -2.828 -1.     5.     0.     0.     0.     5.657]

# Compression: keep the 4 largest coefficients; the squared error is the energy discarded
keep = np.argsort(-np.abs(coef))[:4]
approx = Q8[:, keep] @ coef[keep]
print(approx)                                    # [4.5 4.5 4.5 4.5 9.  9.  8.  0. ]
err2 = np.sum((signal - approx)**2)
dropped = np.sum(np.delete(coef, keep)**2)
print(np.isclose(err2, dropped), round(err2, 10))   # True 1.0: error^2 = sum of the dropped coefficients^2

Three of the eight coefficients are exactly zero because the signal is constant on the pairs $(4, 4)$, $(5, 5)$, and $(9, 9)$. Dropping the smallest nonzero coefficient as well halves the storage at a squared error of exactly $1$.

Key takeaways

• In an orthonormal basis, coordinates are inner products: $\mathbf{v} = \sum \langle \mathbf{v}, \mathbf{u}_i \rangle\mathbf{u}_i$, or $\mathbf{c} = Q^T\mathbf{v}$ in matrix form.
• For an orthogonal (not normalized) basis, divide by $\langle \mathbf{v}_i, \mathbf{v}_i \rangle$.
• Parseval: inner products and norms can be computed from orthonormal coordinates as ordinary dot products.
• Bessel: an orthonormal set that is not a basis captures at most $\lVert \mathbf{v} \rVert^2$ of the energy.
• Truncating an orthonormal expansion leaves a squared error equal to the sum of the squares of the dropped coefficients, which is the basis of transform compression.

🧪 Practice

1. Find the coordinates of $\mathbf{v} = (1, 0, 0)$ in the orthonormal basis of the worked example, and check Parseval.
2. In $\mathbb{R}^2$, the basis $\mathbf{b}_1 = (1, 0)$, $\mathbf{b}_2 = (1, 1)$ is not orthogonal. Show that the "shadow" numbers $\mathbf{v} \cdot \mathbf{b}_i / \lVert \mathbf{b}_i \rVert^2$ for $\mathbf{v} = (0, 1)$ are not its coordinates.
3. Use Bessel's inequality with the orthonormal functions $\tfrac{1}{\sqrt{\pi}}\sin(kx)$, $k = 1, \dots, N$, in $C[-\pi, \pi]$ to show that $\sum_{k=1}^N \left(\int_{-\pi}^{\pi} f(x)\sin(kx),dx\right)^2 \le \pi\int_{-\pi}^{\pi} f^2,dx$.
4. Interview: Why do signal and image compression schemes (JPEG, MP3) transform the data into an orthonormal basis before discarding information, instead of discarding raw samples? Hint: in a well-chosen basis most of the energy sits in a few coefficients, and Parseval turns "drop the small coefficients" into "make the smallest possible squared error".

Orthogonal Projection onto a Subspace

Why it matters. "What is the closest point of a subspace to a given vector?" may be the most frequently solved question in applied linear algebra. It is the best fit of a model to data (least squares), the best approximation of a function by a polynomial or a trigonometric sum, the denoised version of a signal, and the shadow of a point on a plane. The answer is the orthogonal projection, and an orthogonal basis gives it by a direct formula.

Analogy. The shadow of a pole on flat ground at noon, with the sun directly overhead, lands at the point of the ground closest to the pole's tip. The light rays are perpendicular to the ground, and perpendicular is what "closest" means.

From lines to subspaces. Section 2.3 projected onto a single vector $\mathbf{u}$: $\operatorname{proj}_{\mathbf{u}}\mathbf{v} = \dfrac{\langle \mathbf{v}, \mathbf{u} \rangle}{\langle \mathbf{u}, \mathbf{u} \rangle}\mathbf{u}$. If $W$ has an orthogonal basis $\mathbf{w}_1, \dots, \mathbf{w}_k$, then the projection onto $W$ is the sum of the projections onto the basis vectors:

$$\operatorname{proj}W \mathbf{v} = \sum{i=1}^k \frac{\langle \mathbf{v}, \mathbf{w}_i \rangle}{\langle \mathbf{w}_i, \mathbf{w}_i \rangle}\mathbf{w}_i, \qquad \text{or, with an orthonormal basis,} \qquad \operatorname{proj}W \mathbf{v} = \sum{i=1}^k \langle \mathbf{v}, \mathbf{q}_i \rangle\mathbf{q}_i.$$

The residual is perpendicular to $W$. Let $\mathbf{p} = \operatorname{proj}_W\mathbf{v}$. For each basis vector,

$$\langle \mathbf{v} - \mathbf{p}, \mathbf{w}_j \rangle = \langle \mathbf{v}, \mathbf{w}_j \rangle - \frac{\langle \mathbf{v}, \mathbf{w}_j \rangle}{\langle \mathbf{w}_j, \mathbf{w}_j \rangle}\langle \mathbf{w}_j, \mathbf{w}_j \rangle = 0,$$

because every other term contains $\langle \mathbf{w}_i, \mathbf{w}_j \rangle = 0$. Being orthogonal to a basis of $W$ means being orthogonal to all of $W$. This calculation is also why the basis must be orthogonal. With a non-orthogonal basis, the cross terms survive and the sum of line projections is not the projection (see the pitfall in the code below).

Best Approximation Theorem. $\mathbf{p} = \operatorname{proj}_W\mathbf{v}$ is the unique point of $W$ closest to $\mathbf{v}$: $\lVert \mathbf{v} - \mathbf{p} \rVert < \lVert \mathbf{v} - \mathbf{w} \rVert$ for every $\mathbf{w} \in W$ with $\mathbf{w} \ne \mathbf{p}$.

Proof. Write $\mathbf{v} - \mathbf{w} = (\mathbf{v} - \mathbf{p}) + (\mathbf{p} - \mathbf{w})$. The first part is orthogonal to $W$, and the second lies in $W$, so by Pythagoras

$$\lVert \mathbf{v} - \mathbf{w} \rVert^2 = \lVert \mathbf{v} - \mathbf{p} \rVert^2 + \lVert \mathbf{p} - \mathbf{w} \rVert^2 > \lVert \mathbf{v} - \mathbf{p} \rVert^2 \quad \text{when } \mathbf{w} \ne \mathbf{p}. \qquad \blacksquare$$

                          v
                          *
                         /|
                        / |   v - p   (perpendicular to W: the shortest route to W)
                       /  |
         ____________ /___|__________________
        /           /     * p = proj_W v      /
       /           O       \                 /     W (a plane through O)
      /                     \               /
     /                       * w  (any other point of W is farther:
    /_________________________ |v - w|^2 = |v - p|^2 + |p - w|^2)

Worked example. Let $W = \operatorname{span}{\mathbf{w}_1, \mathbf{w}_2}$ with $\mathbf{w}_1 = (1, 1, 0)$ and $\mathbf{w}_2 = (1, -1, 1)$, which are orthogonal, and let $\mathbf{v} = (1, 2, 3)$.

$$\mathbf{p} = \frac{3}{2}(1, 1, 0) + \frac{2}{3}(1, -1, 1) = \left(\tfrac{13}{6}, \tfrac{5}{6}, \tfrac{2}{3}\right), \qquad \mathbf{v} - \mathbf{p} = \tfrac{7}{6}(-1, 1, 2).$$

Check: $(-1, 1, 2) \cdot (1, 1, 0) = 0$ and $(-1, 1, 2) \cdot (1, -1, 1) = 0$. The distance from $\mathbf{v}$ to $W$ is $\tfrac{7}{6}\sqrt{6} = 7/\sqrt{6} \approx 2.858$.

A function example. In $C[-1, 1]$ with $\int_{-1}^{1}$, the closest linear function to $x^2$ is its projection onto $W = \operatorname{span}{1, x}$. The basis ${1, x}$ is orthogonal on this interval, so

$$\operatorname{proj}_W x^2 = \frac{\langle x^2, 1 \rangle}{\langle 1, 1 \rangle},1 + \frac{\langle x^2, x \rangle}{\langle x, x \rangle},x = \frac{2/3}{2} + 0 = \frac{1}{3}.$$

The best linear approximation of $x^2$ on $[-1, 1]$, in the mean-square sense, is the constant $\tfrac{1}{3}$. The residual $x^2 - \tfrac{1}{3}$ is the third Legendre polynomial up to scaling.

The projection matrix. In $\mathbb{R}^n$ with an orthonormal basis of $W$ in the columns of $Q$ ($n \times k$),

$$\operatorname{proj}_W\mathbf{v} = \sum_i \mathbf{q}_i(\mathbf{q}_i^T\mathbf{v}) = QQ^T\mathbf{v}, \qquad P = QQ^T.$$

Its properties characterize orthogonal projections:

Property	Reason
$P^2 = P$ (idempotent)	projecting a second time changes nothing: $QQ^TQQ^T = Q(Q^TQ)Q^T = QQ^T$
$P^T = P$ (symmetric)	this is what makes the projection orthogonal rather than oblique
eigenvalues $1$ (on $W$) and $0$ (on $W^\perp$)	$P\mathbf{w} = \mathbf{w}$ for $\mathbf{w} \in W$; $P\mathbf{x} = \mathbf{0}$ for $\mathbf{x} \perp W$
$\operatorname{tr} P = \operatorname{rank} P = k = \dim W$	the trace is the sum of the eigenvalues
$\lVert P\mathbf{v} \rVert \le \lVert \mathbf{v} \rVert$	Pythagoras: $\lVert \mathbf{v} \rVert^2 = \lVert P\mathbf{v} \rVert^2 + \lVert \mathbf{v} - P\mathbf{v} \rVert^2$

Conversely, any matrix with $P^2 = P$ and $P^T = P$ is the orthogonal projection onto its column space. A projection with $P^2 = P$ but $P^T \ne P$ is oblique: it projects along a direction that is not perpendicular (Section 7.1). If the basis of $W$ is not orthonormal but sits in the columns of a matrix $A$ with independent columns, the same projection is $P = A(A^TA)^{-1}A^T$, derived from the normal equations in Section 9.5.

import numpy as np

w1, w2 = np.array([1., 1., 0.]), np.array([1., -1., 1.])     # orthogonal basis of a plane W
v = np.array([1., 2., 3.])

p = (v @ w1) / (w1 @ w1) * w1 + (v @ w2) / (w2 @ w2) * w2   # sum of projections onto each basis vector
r = v - p
print(p * 6, r * 6)                       # [13.  5.  4.] [-7.  7. 14.]: p = (13, 5, 4)/6, r = 7(-1, 1, 2)/6
print(np.isclose(r @ w1, 0), np.isclose(r @ w2, 0))         # True True: residual is orthogonal to W
print(np.linalg.norm(r))                  # 2.8577380332470415 = 7/sqrt(6), the distance from v to W

# Matrix form with an orthonormal basis: P = Q Q^T
Q = np.column_stack([w1 / np.linalg.norm(w1), w2 / np.linalg.norm(w2)])
P = Q @ Q.T
print(np.allclose(P @ v, p), np.allclose(P @ P, P), np.allclose(P, P.T))   # True True True
print(np.round(np.linalg.eigvalsh(P), 10) + 0.0, round(np.trace(P), 10))   # [0. 1. 1.] 2.0

# Best approximation: no point of W is closer to v than p
rng = np.random.default_rng(0)
others = rng.standard_normal((10000, 2)) @ np.vstack([w1, w2])     # 10000 random points of W
print(np.min(np.linalg.norm(others - v, axis=1)) >= np.linalg.norm(r))   # True

# Pitfall: the sum-of-projections formula needs an ORTHOGONAL basis
u1, u2 = np.array([1., 1., 0.]), np.array([1., 0., 1.])     # also spans a plane, but u1 . u2 = 1
wrong = (v @ u1) / (u1 @ u1) * u1 + (v @ u2) / (u2 @ u2) * u2
print(np.round((v - wrong) @ u1, 4), np.round((v - wrong) @ u2, 4))   # -2.0 -1.5: not orthogonal, so not the projection

Key takeaways

• With an orthogonal basis of $W$, $\operatorname{proj}_W\mathbf{v} = \sum \frac{\langle \mathbf{v}, \mathbf{w}_i \rangle}{\langle \mathbf{w}_i, \mathbf{w}_i \rangle}\mathbf{w}_i$; with an orthonormal basis in $Q$, $P = QQ^T$.
• The residual $\mathbf{v} - \operatorname{proj}_W\mathbf{v}$ is orthogonal to $W$, and the projection is the unique closest point of $W$ (Best Approximation Theorem).
• Orthogonal projection matrices are exactly the matrices with $P^2 = P$ and $P^T = P$; their eigenvalues are $0$ and $1$, and $\operatorname{tr} P = \dim W$.
• The sum-of-line-projections formula is wrong for a non-orthogonal basis; orthogonalize first (Section 9.3) or use $A(A^TA)^{-1}A^T$.

🧪 Practice

1. Project $\mathbf{v} = (2, 0, 4)$ onto the plane spanned by $(1, 2, 2)$ and $(2, 1, -2)$, and find the distance from $\mathbf{v}$ to the plane.
2. Find the projection matrix onto the line spanned by $(1, 2)$ in $\mathbb{R}^2$, and verify $P^2 = P$, $P^T = P$, and $\operatorname{tr} P = 1$.
3. In $C[-1, 1]$, find the best approximation of $f(x) = x^3$ by polynomials of degree at most $1$, using the orthogonal basis ${1, x}$. Compute the mean-square error $\lVert x^3 - \operatorname{proj},x^3 \rVert^2$.
4. Interview: You are given a matrix $P$ and told it is a projection. How would you check numerically whether it is an orthogonal projection, and what is its rank? Hint: check $P^2 \approx P$ and $P^T \approx P$; for a projection, the rank equals the trace, which is far cheaper than an SVD.

Orthogonal Complements

Why it matters. A subspace can be described in two ways: by vectors that span it, or by the conditions it satisfies. A plane in $\mathbb{R}^3$ is "all combinations of two vectors" or "all vectors perpendicular to a normal vector". The orthogonal complement converts one description into the other. It is the canonical partner of a subspace: of all subspaces that together with $W$ fill the space, it is the only one perpendicular to $W$. Section 6.6 already used it to describe the four fundamental subspaces.

Definition. For a subspace $W$ of an inner product space $V$, the orthogonal complement is

$$W^\perp = {\mathbf{v} \in V : \langle \mathbf{v}, \mathbf{w} \rangle = 0 \text{ for every } \mathbf{w} \in W}.$$

To test membership it is enough to check orthogonality to a spanning set of $W$. If $\mathbf{v} \perp \mathbf{w}_i$ for each $i$, then $\langle \mathbf{v}, \sum c_i\mathbf{w}_i \rangle = \sum \overline{c_i}\langle \mathbf{v}, \mathbf{w}_i \rangle = 0$.

Properties. For subspaces $U, W$ of a finite-dimensional inner product space $V$ with $\dim V = n$:

Property	Statement
Subspace	$W^\perp$ is a subspace (it is defined by linear equations)
Trivial intersection	$W \cap W^\perp = {\mathbf{0}}$ (a vector in both is orthogonal to itself)
Dimensions	$\dim W + \dim W^\perp = n$
Double complement	$(W^\perp)^\perp = W$
Reverses inclusions	$U \subseteq W \implies W^\perp \subseteq U^\perp$
Sums and intersections	$(U + W)^\perp = U^\perp \cap W^\perp$ and $(U \cap W)^\perp = U^\perp + W^\perp$
Extremes	${\mathbf{0}}^\perp = V$ and $V^\perp = {\mathbf{0}}$

Proof of the dimension formula and the double complement. Choose an orthonormal basis $\mathbf{q}_1, \dots, \mathbf{q}_k$ of $W$ and extend it to an orthonormal basis $\mathbf{q}_1, \dots, \mathbf{q}_n$ of $V$ (possible by Gram-Schmidt, Section 9.3). A vector $\mathbf{v} = \sum c_i\mathbf{q}_i$ is in $W^\perp$ exactly when $c_i = \langle \mathbf{v}, \mathbf{q}i \rangle = 0$ for $i \le k$. So $W^\perp = \operatorname{span}{\mathbf{q}{k+1}, \dots, \mathbf{q}_n}$, which has dimension $n - k$. Applying the same argument to $W^\perp$ gives $(W^\perp)^\perp = \operatorname{span}{\mathbf{q}_1, \dots, \mathbf{q}_k} = W$. $\blacksquare$

Computing complements with matrices. If $W = \operatorname{Col} A$ in $\mathbb{R}^m$ with the dot product, then $\mathbf{v} \in W^\perp$ means $\mathbf{a}_j^T\mathbf{v} = 0$ for every column, that is, $A^T\mathbf{v} = \mathbf{0}$. This recovers the relations of Section 6.6:

$$(\operatorname{Col} A)^\perp = \operatorname{Nul} A^T, \qquad (\operatorname{Row} A)^\perp = \operatorname{Nul} A.$$

Over $\mathbb{C}$ the transposes become conjugate transposes: $(\operatorname{Col} A)^\perp = \operatorname{Nul} A^*$. With a weighted inner product $\langle \mathbf{u}, \mathbf{v} \rangle_G = \mathbf{u}^TG\mathbf{v}$, the condition becomes $A^TG\mathbf{v} = \mathbf{0}$, so the complement changes when the geometry changes.

Examples.

$W$	$W^\perp$
plane $x + y + z = 0$ in $\mathbb{R}^3$	line spanned by the normal vector $(1, 1, 1)$
line spanned by $(1, 2)$ in $\mathbb{R}^2$	line spanned by $(-2, 1)$
even functions in $C[-1, 1]$	odd functions
constants in $\mathbb{P}2$, with $\int{-1}^{1}$	polynomials with zero mean: $\operatorname{span}{x, x^2 - \tfrac{1}{3}}$
symmetric matrices in $\mathbb{R}^{n \times n}$ (Frobenius)	skew-symmetric matrices

Link to annihilators. The annihilator $W^0 \subseteq V^*$ of Section 7.4 is the set of linear functionals that vanish on $W$. An inner product turns each vector $\mathbf{v}$ into the functional $\langle \cdot, \mathbf{v} \rangle$, and in finite dimensions every functional arises this way (the Riesz representation theorem). Under this identification, $W^\perp$ is the annihilator $W^0$. The inner product makes the abstract annihilator into a concrete subspace of $V$ itself.

A warning about infinite dimensions. In infinite-dimensional spaces, $W^\perp$ is still a subspace and $W \cap W^\perp = {\mathbf{0}}$, but $(W^\perp)^\perp = W$ can fail. Polynomials are dense in the space of square-integrable functions on $[-1, 1]$, so their complement there is ${0}$, and the double complement is the whole space rather than the polynomials. The correct statement is $(W^\perp)^\perp = \overline{W}$, the closure of $W$ (Section 13.4).

import numpy as np
from scipy.linalg import null_space

# W = span of the columns of A in R^4
A = np.array([[1., 0.],
              [1., 1.],
              [0., 1.],
              [1., 2.]])
Wperp = null_space(A.T)                  # W-perp = Nul(A^T): vectors orthogonal to every column
print(Wperp.shape[1])                    # 2: dim W + dim W-perp = 2 + 2 = 4
print(np.allclose(A.T @ Wperp, 0))       # True

# Double complement: (W-perp)-perp = W
Wpp = null_space(Wperp.T)
same_span = np.linalg.matrix_rank(np.hstack([A, Wpp])) == 2
print(same_span)                         # True: the same 2-dimensional subspace

# The complement depends on the inner product: with <u, v>_G = u^T G v, W-perp = Nul(A^T G)
G = np.diag([1., 2., 3., 4.])
WperpG = null_space(A.T @ G)
print(np.allclose(A.T @ G @ WperpG, 0), np.allclose(A.T @ WperpG, 0))   # True False

# A plane in R^3 given by an equation: its complement is the line through the normal vector
plane = null_space(np.array([[1., 1., 1.]]))     # x + y + z = 0
normal = null_space(plane.T)
print(np.round(normal.ravel() / normal[0, 0], 6))     # [1. 1. 1.]

Key takeaways

• $W^\perp$ is the set of vectors orthogonal to all of $W$; checking a spanning set is enough.
• In finite dimensions, $W \cap W^\perp = {\mathbf{0}}$, $\dim W + \dim W^\perp = n$, and $(W^\perp)^\perp = W$.
• $(\operatorname{Col} A)^\perp = \operatorname{Nul} A^T$ and $(\operatorname{Row} A)^\perp = \operatorname{Nul} A$ (with $A^*$ over $\mathbb{C}$), so complements are computed as null spaces.
• The complement depends on the inner product; with weight $G$ it is $\operatorname{Nul}(A^TG)$.
• Complements convert between "spanned by" and "cut out by equations" descriptions of a subspace.

🧪 Practice

1. Find a basis for $W^\perp$ where $W = \operatorname{span}{(1, 0, 2, 1), (0, 1, 1, -1)}$ in $\mathbb{R}^4$, and check the dimension formula.
2. Prove that $(U + W)^\perp = U^\perp \cap W^\perp$ directly from the definition.
3. In $\mathbb{P}2$ with $\langle p, q \rangle = \int_0^1 pq,dx$, find $W^\perp$ for $W = \operatorname{span}{1}$. Compare with the answer for $\int{-1}^{1}$ in the table above.
4. Interview: A system $A\mathbf{x} = \mathbf{b}$ arising from a physical model is singular. How can you decide whether it has a solution without solving it, and what is the physical meaning of the vectors you check against? Hint: $\mathbf{b}$ must lie in $(\operatorname{Nul} A^T)^\perp$; for a pure-Neumann heat or pressure problem, $\operatorname{Nul} A^T$ is spanned by the constant vector, so the sources must sum to zero (a conservation law).

Orthogonal Decomposition Theorem

Why it matters. Splitting a vector into "the part that lies in $W$" and "the rest" is meaningless for a general complement, because different complements give different splits (Section 6.2). With the perpendicular complement the split is unique and respects lengths: the squared lengths of the two parts add up to the squared length of the whole. Signal plus noise, explained plus unexplained variation in statistics, and the row-space and null-space parts of a solution are all instances of this one theorem.

Theorem (orthogonal decomposition). Let $W$ be a finite-dimensional subspace of an inner product space $V$. Every $\mathbf{v} \in V$ can be written uniquely as

$$\mathbf{v} = \hat{\mathbf{v}} + \mathbf{z}, \qquad \hat{\mathbf{v}} \in W, \quad \mathbf{z} \in W^\perp,$$

and $\hat{\mathbf{v}} = \operatorname{proj}_W\mathbf{v}$. In other words, $V = W \oplus W^\perp$.

Proof. Existence: $W$ has an orthonormal basis (Gram-Schmidt, Section 9.3). Set $\hat{\mathbf{v}} = \operatorname{proj}_W\mathbf{v}$ and $\mathbf{z} = \mathbf{v} - \hat{\mathbf{v}}$. The previous topics showed $\hat{\mathbf{v}} \in W$ and $\mathbf{z} \in W^\perp$. Uniqueness: if $\hat{\mathbf{v}} + \mathbf{z} = \hat{\mathbf{v}}' + \mathbf{z}'$, then $\hat{\mathbf{v}} - \hat{\mathbf{v}}' = \mathbf{z}' - \mathbf{z}$ lies in $W \cap W^\perp = {\mathbf{0}}$. $\blacksquare$

Only $W$ needs to be finite-dimensional. $V$ can be a function space, which is what makes best polynomial and Fourier approximation work (Section 9.5).

                 W-perp
                   ^
                   |
                   |  z = v - proj_W v
            z  +   |- - - - - - - * v
                   |              :
                   |              :    ||v||^2 = ||v_hat||^2 + ||z||^2
                   |              :
   ----------------O--------------+--------------->  W
                                v_hat = proj_W v

   the split along W and W-perp is unique; any other complement gives a different,
   non-perpendicular split, and the squared lengths no longer add up

Consequences.

1. Complementary projections. $\operatorname{proj}{W^\perp}\mathbf{v} = \mathbf{v} - \operatorname{proj}W\mathbf{v}$. In matrix form, $P{W^\perp} = I - P_W$, and $P_WP{W^\perp} = 0$.
2. Energy splits. $\lVert \mathbf{v} \rVert^2 = \lVert \operatorname{proj}W\mathbf{v} \rVert^2 + \lVert \operatorname{proj}{W^\perp}\mathbf{v} \rVert^2$. This proves Bessel's inequality, with equality exactly when $\mathbf{v} \in W$.
3. Distance to a subspace. $\operatorname{dist}(\mathbf{v}, W) = \lVert \operatorname{proj}_{W^\perp}\mathbf{v} \rVert$. Sometimes the complement is the easier space to project onto. The distance from $\mathbf{v}$ to a hyperplane with unit normal $\mathbf{n}$ is just $\lvert \langle \mathbf{v}, \mathbf{n} \rangle \rvert$.
4. More than two pieces. If $V = W_1 \oplus \cdots \oplus W_k$ with the $W_i$ mutually orthogonal, then $I = P_1 + \cdots + P_k$ with $P_iP_j = 0$ for $i \ne j$, and $\lVert \mathbf{v} \rVert^2 = \sum \lVert P_i\mathbf{v} \rVert^2$. The spectral theorem (Chapter 10) writes every symmetric matrix as $\sum \lambda_iP_i$ for such a decomposition.

Instances of the theorem.

$V$	$W$	$W^\perp$	Decomposition
$\mathbb{R}^n$	$\operatorname{Row} A$	$\operatorname{Nul} A$	$\mathbf{x} = \mathbf{x}_r + \mathbf{x}_n$ (Section 6.6)
$\mathbb{R}^n$ (data)	constant vectors	vectors with zero mean	$\mathbf{y} = \bar{y}\mathbf{1} + (\mathbf{y} - \bar{y}\mathbf{1})$
$\mathbb{R}^{n \times n}$ (Frobenius)	symmetric	skew-symmetric	$M = \tfrac{1}{2}(M + M^T) + \tfrac{1}{2}(M - M^T)$
$C[-a, a]$	even functions	odd functions	$f(x) = \tfrac{f(x) + f(-x)}{2} + \tfrac{f(x) - f(-x)}{2}$
$\mathbb{R}^m$ (regression)	$\operatorname{Col} X$	residual space	$\mathbf{y} = \hat{\mathbf{y}} + \mathbf{e}$ (Section 9.5)

The second row is the identity behind the variance formula. $\sum y_i^2 = n\bar{y}^2 + \sum(y_i - \bar{y})^2$ is Pythagoras for the split into mean and deviations. In regression with an intercept, the same reasoning gives the decomposition $\text{SST} = \text{SSR} + \text{SSE}$ of the total sum of squares into explained and residual parts.

import numpy as np

# 1) Data vector = mean part (in W = span{1}) + deviations (in W-perp)
y = np.array([3., 5., 4., 8.])
ones = np.ones_like(y)
mean_part = (y @ ones) / (ones @ ones) * ones      # projection onto span{(1, 1, 1, 1)}
dev = y - mean_part
print(mean_part, dev)                              # [5. 5. 5. 5.] [-2.  0. -1.  3.]
print(y @ y, mean_part @ mean_part + dev @ dev)    # 114.0 114.0: sum y^2 = n*ybar^2 + sum (y - ybar)^2

# 2) Matrix = symmetric part + skew-symmetric part, orthogonal in the Frobenius inner product
M = np.array([[1., 4.], [2., 3.]])
S, K = (M + M.T) / 2, (M - M.T) / 2
fro = lambda A, B: np.trace(A.T @ B)
print(fro(S, K), fro(M, M), fro(S, S) + fro(K, K))  # 0.0 30.0 30.0

# 3) Projections onto complementary subspaces add up to the identity
Q = np.linalg.qr(np.array([[1., 0.], [1., 1.], [0., 1.], [1., 2.]]))[0]   # orthonormal basis of W (Section 9.3)
P = Q @ Q.T                                        # projection onto W
P_perp = np.eye(4) - P                             # projection onto W-perp
print(np.allclose(P @ P_perp, 0), np.allclose(P_perp @ P_perp, P_perp))   # True True
v = np.array([1., 2., 3., 4.])
print(np.isclose(v @ v, (P @ v) @ (P @ v) + (P_perp @ v) @ (P_perp @ v)))  # True

Key takeaways

• $V = W \oplus W^\perp$: every vector splits uniquely into $\operatorname{proj}_W\mathbf{v}$ plus a part orthogonal to $W$; only $W$ needs to be finite-dimensional.
• The two parts satisfy Pythagoras, $\lVert \mathbf{v} \rVert^2 = \lVert \hat{\mathbf{v}} \rVert^2 + \lVert \mathbf{z} \rVert^2$, which gives Bessel's inequality and distance formulas.
• $P_{W^\perp} = I - P_W$; mutually orthogonal pieces give projections that sum to $I$ and annihilate each other.
• Mean plus deviations, symmetric plus skew-symmetric, even plus odd, and row space plus null space are all orthogonal decompositions.

🧪 Practice

1. Decompose $\mathbf{v} = (1, 2, 3)$ as a vector in the line spanned by $(1, 1, 1)$ plus a vector orthogonal to it, and verify Pythagoras.
2. Write $f(x) = e^x$ on $[-1, 1]$ as an even part plus an odd part, and identify them as familiar functions. Why are the two parts orthogonal?
3. Show that if $P$ is an orthogonal projection, then $I - P$ is also an orthogonal projection, and that $\operatorname{Col}(I - P) = (\operatorname{Col} P)^\perp$.
4. Interview: In ordinary least squares with an intercept, why does $\text{SST} = \text{SSR} + \text{SSE}$ hold exactly, and why can $R^2 = 1 - \text{SSE}/\text{SST}$ become negative if you fit a model without an intercept but still compute SST around the mean? Hint: the identity is Pythagoras for $\mathbf{y} - \bar{y}\mathbf{1} = (\hat{\mathbf{y}} - \bar{y}\mathbf{1}) + \mathbf{e}$, which needs $\mathbf{1} \in \operatorname{Col} X$ so that $\mathbf{e} \perp (\hat{\mathbf{y}} - \bar{y}\mathbf{1})$.

<a id="93-gram-schmidt-and-qr"></a>

9.3 Gram-Schmidt and QR

This subchapter turns any basis into an orthonormal one, records the process as the matrix factorization $A = QR$, and compares the three ways to compute it (Gram-Schmidt, Householder reflections, and Givens rotations) for numerical stability and cost.

Gram-Schmidt Process

Why it matters. Orthonormal bases make coordinates, projections, and least squares easy, but the bases that problems provide are rarely orthonormal. Examples are the columns of a data matrix, the monomials $1, x, x^2, \dots$, and a set of measured directions. Gram-Schmidt converts any linearly independent list into an orthonormal list with the same span, one vector at a time. It also proves that every finite-dimensional inner product space has an orthonormal basis, a fact that Section 9.2 relied on.

Analogy. Imagine setting up perpendicular axes from a set of slanted sticks. Keep the first stick as it is. For the second stick, remove the part that points along the first, so that only the part sticking out perpendicularly remains. For the third, remove the parts along both earlier axes. Each new axis keeps only what is genuinely new in the corresponding stick.

The algorithm. Given linearly independent $\mathbf{x}_1, \dots, \mathbf{x}_k$:

$$\begin{aligned} \mathbf{v}_1 &= \mathbf{x}_1, \ \mathbf{v}_2 &= \mathbf{x}_2 - \frac{\langle \mathbf{x}_2, \mathbf{v}_1 \rangle}{\langle \mathbf{v}_1, \mathbf{v}_1 \rangle}\mathbf{v}_1, \ \mathbf{v}_j &= \mathbf{x}j - \sum{i=1}^{j-1} \frac{\langle \mathbf{x}_j, \mathbf{v}_i \rangle}{\langle \mathbf{v}_i, \mathbf{v}_i \rangle}\mathbf{v}i = \mathbf{x}j - \operatorname{proj}{W{j-1}}\mathbf{x}_j, \end{aligned}$$

where $W_{j-1} = \operatorname{span}{\mathbf{v}1, \dots, \mathbf{v}{j-1}}$. Finally normalize: $\mathbf{q}_j = \mathbf{v}_j / \lVert \mathbf{v}_j \rVert$. Equivalently, normalize as you go and use $\mathbf{v}_j = \mathbf{x}j - \sum{i<j}\langle \mathbf{x}_j, \mathbf{q}_i \rangle\mathbf{q}_i$.

Why it works. By induction on $j$:

• $\mathbf{v}j$ is the residual of projecting $\mathbf{x}j$ onto $W{j-1}$, so it is orthogonal to $W{j-1}$ (Section 9.2), and in particular to all earlier $\mathbf{v}_i$.
• $\mathbf{v}_j \ne \mathbf{0}$: otherwise $\mathbf{x}j \in W{j-1} = \operatorname{span}{\mathbf{x}1, \dots, \mathbf{x}{j-1}}$, contradicting independence.
• Nested spans: $\operatorname{span}{\mathbf{q}_1, \dots, \mathbf{q}_j} = \operatorname{span}{\mathbf{x}_1, \dots, \mathbf{x}_j}$ for every $j$. Each $\mathbf{v}_j$ is $\mathbf{x}_j$ minus a combination of earlier vectors, and conversely each $\mathbf{x}_j$ is $\mathbf{v}_j$ plus such a combination.

The nested-span property is stronger than "same final span". It is the reason the factorization in the QR topic below has a triangular factor.

Consequences.

• Every finite-dimensional inner product space has an orthonormal basis: apply Gram-Schmidt to any basis.
• Every orthonormal set extends to an orthonormal basis: extend it to a basis (Section 6.4), then apply Gram-Schmidt. The given vectors pass through unchanged, because they are already orthonormal.
• If the input is dependent, some $\mathbf{v}_j$ comes out as $\mathbf{0}$. Skipping such vectors yields an orthonormal basis of the span, and the number kept is the rank. In floating point, "$\mathbf{0}$" must be replaced by "small relative to $\lVert \mathbf{x}_j \rVert$".

Worked example. Orthonormalize $\mathbf{x}_1 = (1, 1, 0)$, $\mathbf{x}_2 = (1, 0, 1)$, $\mathbf{x}_3 = (0, 1, 1)$.

1. $\mathbf{v}_1 = (1, 1, 0)$, with $\lVert \mathbf{v}_1 \rVert^2 = 2$.
2. $\langle \mathbf{x}_2, \mathbf{v}_1 \rangle = 1$, so $\mathbf{v}_2 = (1, 0, 1) - \tfrac{1}{2}(1, 1, 0) = (\tfrac{1}{2}, -\tfrac{1}{2}, 1)$. Scaling does not affect orthogonality, so use $\mathbf{v}_2 = (1, -1, 2)$, with $\lVert \mathbf{v}_2 \rVert^2 = 6$.
3. $\langle \mathbf{x}_3, \mathbf{v}_1 \rangle = 1$ and $\langle \mathbf{x}_3, \mathbf{v}_2 \rangle = 1$, so $\mathbf{v}_3 = (0, 1, 1) - \tfrac{1}{2}(1, 1, 0) - \tfrac{1}{6}(1, -1, 2) = (-\tfrac{2}{3}, \tfrac{2}{3}, \tfrac{2}{3})$, which scales to $(-1, 1, 1)$.

Check: $(-1, 1, 1) \cdot (1, 1, 0) = 0$ and $(-1, 1, 1) \cdot (1, -1, 2) = 0$. The orthonormal basis is

$$\mathbf{q}_1 = \tfrac{1}{\sqrt{2}}(1, 1, 0), \qquad \mathbf{q}_2 = \tfrac{1}{\sqrt{6}}(1, -1, 2), \qquad \mathbf{q}_3 = \tfrac{1}{\sqrt{3}}(-1, 1, 1).$$

   Step 2 of Gram-Schmidt: remove from x2 its shadow on v1

              x2
              *
             /|
            / |   v2 = x2 - proj_v1(x2)     (what is new in x2: perpendicular to v1)
           /  |
          /   |
         O----+----------------> v1 = x1
          proj_v1(x2) = (<x2, v1> / <v1, v1>) v1

In function spaces. Applied to $1, x, x^2, x^3$ in $C[-1, 1]$ with $\int_{-1}^{1}$, Gram-Schmidt produces $1$, $x$, $x^2 - \tfrac{1}{3}$, $x^3 - \tfrac{3}{5}x$. These are the Legendre polynomials, scaled to leading coefficient $1$. With other weights and intervals, the same procedure produces the other families in Section 9.1.

import numpy as np
import sympy as sp

def gram_schmidt(X, tol=1e-12):
    """Classical Gram-Schmidt on the columns of X. Returns orthonormal columns spanning the same space."""
    Q = []
    for x in X.T:                                   # process columns one at a time
        v = x - sum((x @ q) * q for q in Q)         # subtract the projection onto span(previous q's)
        norm = np.linalg.norm(v)
        if norm > tol * np.linalg.norm(x):          # v ~ 0 means x depends on the previous columns
            Q.append(v / norm)
    return np.column_stack(Q)

X = np.array([[1., 1., 0.],
              [1., 0., 1.],
              [0., 1., 1.]])
Q = gram_schmidt(X)
print(np.round(Q * [np.sqrt(2), np.sqrt(6), np.sqrt(3)], 6))    # columns (1,1,0), (1,-1,2), (-1,1,1)
# [[ 1.  1. -1.]
#  [ 1. -1.  1.]
#  [ 0.  2.  1.]]
print(np.allclose(Q.T @ Q, np.eye(3)))                          # True

# Nested spans: span{q1..qk} = span{x1..xk} for every k
for k in range(1, 4):
    print(k, np.linalg.matrix_rank(np.hstack([X[:, :k], Q[:, :k]])))    # 1 1, 2 2, 3 3

# A dependent column is detected and skipped
Xdep = np.column_stack([X[:, 0], X[:, 1], X[:, 0] + 2 * X[:, 1]])
print(gram_schmidt(Xdep).shape)                                 # (3, 2)

# The same algorithm in a function space: Gram-Schmidt on 1, x, x^2, x^3 in C[-1, 1]
x = sp.symbols("x")
ip = lambda f, g: sp.integrate(f * g, (x, -1, 1))
basis = []
for f in [1, x, x**2, x**3]:
    v = f - sum(ip(f, b) / ip(b, b) * b for b in basis)        # orthogonal (not normalized) version
    basis.append(sp.expand(v))
print(basis)                                                    # [1, x, x**2 - 1/3, x**3 - 3*x/5]
print([sp.expand(sp.legendre(n, x) / sp.Poly(sp.legendre(n, x), x).LC()) for n in range(4)])
# [1, x, x**2 - 1/3, x**3 - 3*x/5]: monic Legendre polynomials

Cost. For $n$ vectors in $\mathbb{R}^m$, step $j$ computes $j - 1$ inner products and subtracts $j - 1$ multiples of vectors, about $4m(j - 1)$ flops. The total is about $2mn^2$ flops.

Key takeaways

• Gram-Schmidt subtracts from each vector its projection onto the span of the previous ones, then normalizes.
• It preserves nested spans: $\operatorname{span}{\mathbf{q}_1, \dots, \mathbf{q}_j} = \operatorname{span}{\mathbf{x}_1, \dots, \mathbf{x}_j}$ for every $j$.
• It proves that orthonormal bases exist and that orthonormal sets extend to orthonormal bases.
• A zero (or tiny) $\mathbf{v}_j$ signals a dependent input vector.
• In function spaces it generates orthogonal polynomials, such as Legendre from $1, x, x^2, \dots$.

🧪 Practice

1. Apply Gram-Schmidt to $(1, 0, 1, 0)$, $(1, 1, 1, 1)$, $(0, 1, 2, 1)$ in $\mathbb{R}^4$, and check that the results are orthonormal.
2. Apply Gram-Schmidt to $1, x, x^2$ in $C[0, 1]$ with $\int_0^1$. Compare the result with the Legendre polynomials on $[-1, 1]$ after the substitution $x \mapsto 2x - 1$.
3. Show that if $\mathbf{x}_1, \dots, \mathbf{x}_k$ are already orthogonal, Gram-Schmidt returns them unchanged up to normalization.
4. Interview: Implement Gram-Schmidt for the columns of a NumPy array and test it on a matrix with two nearly parallel columns, such as $(1, 1, 1)$ and $(1, 1, 1 + 10^{-10})$. What do you observe about the orthogonality of the output, and why? Hint: the second vector is the difference of two nearly equal vectors; cancellation destroys most significant digits, and the error is then amplified by normalization (next topic).

Modified Gram-Schmidt for Stability

Why it matters. In exact arithmetic, the Gram-Schmidt formula of the previous topic is correct. In floating point, the classical version can return "orthonormal" vectors that are far from orthogonal when the input columns are nearly dependent. Nearly dependent columns are common in practice: polynomial features, correlated measurements, and Krylov vectors (Section 12.3) all produce them. A small reordering of the same arithmetic, modified Gram-Schmidt, is far more robust. The difference shows how two algorithms that are identical on paper can behave very differently on a computer.

Classical vs. modified. Both compute $\mathbf{v}_k = \mathbf{x}k - \sum{j<k}\langle \cdot, \mathbf{q}_j \rangle\mathbf{q}_j$. They differ in which vector goes into each inner product:

	Classical Gram-Schmidt (CGS)	Modified Gram-Schmidt (MGS)
Coefficient $r_{jk}$	$\langle \mathbf{x}_k, \mathbf{q}_j \rangle$, using the original vector	$\langle \mathbf{v}^{(j)}, \mathbf{q}_j \rangle$, using the vector already cleaned of $\mathbf{q}1, \dots, \mathbf{q}{j-1}$
Update	all projections subtracted at once	$\mathbf{v}^{(j+1)} = \mathbf{v}^{(j)} - r_{jk}\mathbf{q}_j$, one at a time
Exact arithmetic	same result	same result
Loss of orthogonality $\lVert Q^TQ - I \rVert$	up to about $\kappa(A)^2,\varepsilon$, can reach $O(1)$	about $\kappa(A),\varepsilon$
Parallelism	the $r_{jk}$ for fixed $k$ are one matrix-vector product	sequential inner loop

Here $\varepsilon \approx 1.1 \times 10^{-16}$ is the double-precision unit roundoff, and $\kappa(A)$ is the condition number of $A$ (Section 11.3), which is large when the columns are nearly dependent.

Why are the two the same in exact arithmetic? $\mathbf{v}^{(j)}$ differs from $\mathbf{x}_k$ by multiples of $\mathbf{q}1, \dots, \mathbf{q}{j-1}$, all of which are orthogonal to $\mathbf{q}_j$, so $\langle \mathbf{v}^{(j)}, \mathbf{q}_j \rangle = \langle \mathbf{x}_k, \mathbf{q}_j \rangle$. In floating point, the computed $\mathbf{q}_i$ are only nearly orthogonal. CGS measures every component against the raw $\mathbf{x}_k$, so errors in the earlier $\mathbf{q}$'s leak into every coefficient. MGS measures each component of what is actually left, so each step also cleans up the component that rounding errors introduced in the previous step.

   CGS, step k:  r_1k = <x_k, q1>,  r_2k = <x_k, q2>,  r_3k = <x_k, q3>   (all from the raw x_k)
                 v = x_k - r_1k q1 - r_2k q2 - r_3k q3

   MGS, step k:  v = x_k
                 r_1k = <v, q1>;  v = v - r_1k q1      <- v is now cleaned of q1
                 r_2k = <v, q2>;  v = v - r_2k q2      <- measured on the cleaned v
                 r_3k = <v, q3>;  v = v - r_3k q3

The Läuchli example. Take $\varepsilon = 10^{-8}$, so that $\varepsilon^2 = 10^{-16}$ is below double precision ($1 + \varepsilon^2$ rounds to $1$), and

$$A = \begin{bmatrix} 1 & 1 & 1 \ \varepsilon & 0 & 0 \ 0 & \varepsilon & 0 \ 0 & 0 & \varepsilon \end{bmatrix}.$$

The columns are independent but nearly parallel. In CGS, $\mathbf{q}_2$ and $\mathbf{q}_3$ come out at $60°$ to each other instead of $90°$: their inner product is $0.5$. MGS keeps the error near $10^{-8}$.

Doing even better. Two further remedies are standard:

• Reorthogonalization (CGS2). Run the classical projection step twice. The rule "twice is enough" (Kahan, Parlett) says that the second pass restores orthogonality to about $\varepsilon$, provided the matrix is not numerically rank deficient. CGS2 costs twice as much as CGS, but it keeps CGS's matrix-vector structure, which runs well on parallel hardware. It is common in Krylov methods (GMRES, Arnoldi).
• Householder QR (two topics ahead) achieves $\lVert Q^TQ - I \rVert \approx \varepsilon$ regardless of conditioning. It is the default in LAPACK and in np.linalg.qr.

import numpy as np

def cgs(A):
    """Classical Gram-Schmidt: coefficients use the ORIGINAL column a_k."""
    m, n = A.shape
    Q, R = np.zeros((m, n)), np.zeros((n, n))
    for k in range(n):
        R[:k, k] = Q[:, :k].T @ A[:, k]            # all projections computed from a_k at once
        v = A[:, k] - Q[:, :k] @ R[:k, k]
        R[k, k] = np.linalg.norm(v)
        Q[:, k] = v / R[k, k]
    return Q, R

def mgs(A):
    """Modified Gram-Schmidt: each projection is removed from the CURRENT, partly cleaned vector."""
    m, n = A.shape
    Q, R = np.zeros((m, n)), np.zeros((n, n))
    for k in range(n):
        v = A[:, k].copy()
        for j in range(k):
            R[j, k] = Q[:, j] @ v                  # uses v, which is already orthogonal to q_1..q_{j-1}
            v -= R[j, k] * Q[:, j]
        R[k, k] = np.linalg.norm(v)
        Q[:, k] = v / R[k, k]
    return Q, R

def cgs2(A):
    """Classical Gram-Schmidt with one full reorthogonalization pass ("twice is enough")."""
    m, n = A.shape
    Q, R = np.zeros((m, n)), np.zeros((n, n))
    for k in range(n):
        v = A[:, k].copy()
        for _ in range(2):
            c = Q[:, :k].T @ v
            v -= Q[:, :k] @ c
            R[:k, k] += c
        R[k, k] = np.linalg.norm(v)
        Q[:, k] = v / R[k, k]
    return Q, R

loss = lambda Q: np.linalg.norm(Q.T @ Q - np.eye(Q.shape[1]))     # 0 for perfectly orthonormal columns

# Lauchli matrix: nearly dependent columns; eps^2 is below machine precision
eps = 1e-8
L = np.array([[1, 1, 1], [eps, 0, 0], [0, eps, 0], [0, 0, eps]])
for name, f in [("CGS", cgs), ("MGS", mgs), ("CGS2", cgs2), ("Householder", np.linalg.qr)]:
    print(f"{name:12s} {loss(f(L)[0]):.1e}")
# CGS          7.1e-01      q2 and q3 end up at 60 degrees instead of 90
# MGS          1.2e-08      about kappa * machine epsilon (kappa ~ 1.7e8)
# CGS2         5.2e-17
# Householder  2.7e-16

# Loss of orthogonality versus condition number (Hilbert-like test matrices)
print(f"{'kappa':>8} {'CGS':>9} {'MGS':>9} {'CGS2':>9} {'Householder':>12}")
for n in [4, 6, 8, 10]:
    A = np.array([[1 / (i + j + 1) for j in range(n)] for i in range(n + 2)])   # (n+2) x n
    k = np.linalg.cond(A)
    print(f"{k:8.0e} {loss(cgs(A)[0]):9.0e} {loss(mgs(A)[0]):9.0e} "
          f"{loss(cgs2(A)[0]):9.0e} {loss(np.linalg.qr(A)[0]):12.0e}")
#    kappa       CGS       MGS      CGS2  Householder
#    6e+03     5e-11     4e-13     3e-16        4e-16
#    4e+06     3e-05     9e-11     3e-16        8e-16
#    4e+09     1e+00     6e-08     6e-16        1e-15
#    3e+12     3e+00     9e-05     6e-16        1e-15

The table shows the predicted scaling: the MGS column grows in proportion to $\kappa$, while the CGS column grows roughly like $\kappa^2$ and reaches $O(1)$, at which point the output is useless. Exact digits vary slightly between machines and BLAS libraries.

An important subtlety. Even when MGS's $Q$ is far from orthogonal, its $R$ is accurate, and $QR$ reproduces $A$ to about $\varepsilon$. MGS applied to the augmented matrix $[A ;; \mathbf{b}]$ also solves least squares problems stably (Björck). The loss of orthogonality matters when $Q$ itself is used, for example as a basis for projections.

Key takeaways

• CGS and MGS are identical in exact arithmetic; MGS projects each component out of the partly cleaned vector instead of the original one.
• Loss of orthogonality: CGS up to $\kappa^2\varepsilon$ (it can fail completely), MGS about $\kappa\varepsilon$, CGS2 and Householder about $\varepsilon$.
• Nearly dependent columns (large $\kappa$) are what trigger the problem, so always prefer MGS over CGS.
• When orthogonality to working precision is required, use reorthogonalization (CGS2) or Householder QR.

🧪 Practice

1. Run CGS and MGS by hand, keeping only 3 significant digits after each operation, on the columns $(1, 0.001, 0)$, $(1, 0, 0.001)$. Compare the inner products of the results.
2. Rewrite mgs in the "row-oriented" form: as soon as $\mathbf{q}_k$ is computed, subtract its component from all remaining columns. Check that it gives the same $Q$ and $R$.
3. Extend the experiment table to n = 12 and explain what happens to CGS2. Hint: "twice is enough" assumes $\kappa(A)\varepsilon < 1$.
4. Interview: You inherit an Arnoldi or Lanczos implementation (Section 12.3) whose computed eigenvalues show spurious duplicate copies after a few hundred iterations. How does this relate to Gram-Schmidt, and what are your options? Hint: the Krylov basis loses orthogonality, so converged directions reappear ("ghost" eigenvalues); fixes are full or selective reorthogonalization (CGS2 or MGS), or restarting.

QR Decomposition

Why it matters. Gram-Schmidt is more than a way to produce orthonormal vectors. Its bookkeeping, recording how each original column is built from the new orthonormal ones, is a matrix factorization $A = QR$. This factorization is the standard tool for least squares (Section 9.5). It is the engine of the QR algorithm, which computes eigenvalues (Section 12.3), and it gives orthonormal bases for column spaces as well as rank estimates. Software rarely forms an orthonormal basis without also returning $R$.

Theorem (QR decomposition). Let $A$ be an $m \times n$ real matrix with linearly independent columns (so $m \ge n$). Then

$$A = QR,$$

where $Q$ is $m \times n$ with orthonormal columns ($Q^TQ = I_n$) and $R$ is $n \times n$, upper triangular, with positive diagonal entries. This factorization is unique. Over $\mathbb{C}$, $Q^*Q = I$ and the diagonal of $R$ is real and positive.

Derivation from Gram-Schmidt. Gram-Schmidt gives $\mathbf{v}k = \mathbf{x}k - \sum{j<k} r{jk}\mathbf{q}j$ with $r{jk} = \langle \mathbf{x}_k, \mathbf{q}_j \rangle$, and $\mathbf{q}k = \mathbf{v}k / r{kk}$ with $r{kk} = \lVert \mathbf{v}_k \rVert > 0$. Solving for $\mathbf{x}_k$:

$$\mathbf{x}k = r{1k}\mathbf{q}1 + r{2k}\mathbf{q}2 + \cdots + r{kk}\mathbf{q}_k.$$

Column $k$ of $A$ uses only $\mathbf{q}_1, \dots, \mathbf{q}_k$. This is the nested-span property, and it means the coefficient matrix is upper triangular:

$$\underbrace{\begin{bmatrix} | & & | \ \mathbf{x}1 & \cdots & \mathbf{x}n \ | & & | \end{bmatrix}}{A} = \underbrace{\begin{bmatrix} | & & | \ \mathbf{q}1 & \cdots & \mathbf{q}n \ | & & | \end{bmatrix}}{Q}\underbrace{\begin{bmatrix} r{11} & r{12} & \cdots & r_{1n} \ 0 & r_{22} & \cdots & r_{2n} \ \vdots & & \ddots & \vdots \ 0 & 0 & \cdots & r_{nn} \end{bmatrix}}_{R}.$$

Since $Q^TQ = I$, the factor $R$ can also be computed directly as $R = Q^TA$, whose entries are $r_{jk} = \mathbf{q}_j^T\mathbf{x}_k$.

Uniqueness. If $QR = Q'R'$, then $Q'^TQ = R'R^{-1}$. The right side is upper triangular. The left side, $M = Q'^TQ$, satisfies $M^TM = I$ when both $Q$ and $Q'$ span $\operatorname{Col} A$. An upper triangular matrix $M$ with $M^TM = I$ is diagonal with entries $\pm 1$: its first column has a single nonzero entry, which must be $\pm 1$, and induction finishes the argument. Positive diagonals in $R$ and $R'$ force $M = I$. $\blacksquare$

Worked example. For the matrix whose columns were orthonormalized in the Gram-Schmidt topic,

$$\begin{bmatrix} 1 & 1 & 0 \ 1 & 0 & 1 \ 0 & 1 & 1 \end{bmatrix} = \begin{bmatrix} \tfrac{1}{\sqrt{2}} & \tfrac{1}{\sqrt{6}} & -\tfrac{1}{\sqrt{3}} \ \tfrac{1}{\sqrt{2}} & -\tfrac{1}{\sqrt{6}} & \tfrac{1}{\sqrt{3}} \ 0 & \tfrac{2}{\sqrt{6}} & \tfrac{1}{\sqrt{3}} \end{bmatrix}\begin{bmatrix} \sqrt{2} & \tfrac{1}{\sqrt{2}} & \tfrac{1}{\sqrt{2}} \ 0 & \tfrac{3}{\sqrt{6}} & \tfrac{1}{\sqrt{6}} \ 0 & 0 & \tfrac{2}{\sqrt{3}} \end{bmatrix}.$$

For instance, $r_{23} = \mathbf{q}_2 \cdot \mathbf{x}_3 = \tfrac{1}{\sqrt{6}}(0 - 1 + 2) = \tfrac{1}{\sqrt{6}}$.

Reduced vs. complete QR.

	Reduced (thin) QR	Complete (full) QR
Shapes	$Q$: $m \times n$, $R$: $n \times n$	$Q$: $m \times m$, $R$: $m \times n$ (bottom $m - n$ rows zero)
$Q$	orthonormal basis of $\operatorname{Col} A$	square orthogonal matrix; its last $m - n$ columns span $(\operatorname{Col} A)^\perp = \operatorname{Nul} A^T$
Typical use	least squares, projections	complements, residual norms, theory
NumPy	np.linalg.qr(A) (default mode="reduced")	np.linalg.qr(A, mode="complete")

What QR is used for.

Task	How QR helps
Orthonormal basis of $\operatorname{Col} A$	the columns of $Q$; the projection onto $\operatorname{Col} A$ is $QQ^T$
Solve square $A\mathbf{x} = \mathbf{b}$	$R\mathbf{x} = Q^T\mathbf{b}$ by back substitution; more stable than LU but costs about twice as much
Least squares	$R\hat{\mathbf{x}} = Q^T\mathbf{b}$, without forming $A^TA$ (Section 9.5)
Determinant	$\lvert\det A\rvert = \prod_i r_{ii}$ for square $A$
Rank and near-dependence	QR with column pivoting, $AP = QR$, orders $\lvert r_{11} \rvert \ge \lvert r_{22} \rvert \ge \cdots$; a sharp drop marks the numerical rank
Eigenvalues	the QR algorithm repeats $A_k = Q_kR_k$, $A_{k+1} = R_kQ_k$ (Section 12.3)

The sign convention in software. LAPACK, and therefore NumPy and SciPy, computes QR with Householder reflections. It does not force the diagonal of $R$ to be positive, so some $r_{kk}$ may be negative and the corresponding columns of $Q$ flipped. Multiplying by $D = \operatorname{diag}(\pm 1)$ restores the convention: $A = (QD)(DR)$ with $D^2 = I$.

import numpy as np
from scipy.linalg import qr, solve_triangular

A = np.array([[1., 1., 0.],
              [1., 0., 1.],
              [0., 1., 1.]])
Q, R = np.linalg.qr(A)
print(np.round(R, 4))
# [[-1.4142 -0.7071 -0.7071]
#  [ 0.      1.2247  0.4082]
#  [ 0.      0.      1.1547]]       NumPy (LAPACK) does not force a positive diagonal

# Normalize signs so that diag(R) > 0: this is the unique QR, the one Gram-Schmidt produces
s = np.sign(np.diag(R))
Q, R = Q * s, s[:, None] * R           # A = (Q D)(D R) with D = diag(s), D^2 = I
print(np.round(R, 4))
# [[1.4142 0.7071 0.7071]              r11 = sqrt(2), r12 = r13 = 1/sqrt(2)
#  [0.     1.2247 0.4082]              r22 = 3/sqrt(6), r23 = 1/sqrt(6)
#  [0.     0.     1.1547]]             r33 = 2/sqrt(3)
print(np.allclose(Q @ R, A), np.allclose(Q.T @ Q, np.eye(3)))   # True True

# Solve A x = b: Q R x = b  =>  R x = Q^T b (back substitution, no inverse)
b = np.array([2., 3., 5.])
x = solve_triangular(R, Q.T @ b)
print(x, np.allclose(A @ x, b))       # [0. 2. 3.] True

# |det A| is the product of the diagonal of R
print(round(np.prod(np.diag(R)), 12), round(abs(np.linalg.det(A)), 12))   # 2.0 2.0

# Reduced vs complete QR of a tall matrix
T = np.array([[1., 2.], [3., 4.], [5., 6.]])
Qr, Rr = np.linalg.qr(T, mode="reduced")
Qc, Rc = np.linalg.qr(T, mode="complete")
print(Qr.shape, Rr.shape, Qc.shape, Rc.shape)                  # (3, 2) (2, 2) (3, 3) (3, 2)
print(np.round(Rc[2], 12) + 0.0)                                # [0. 0.]: the last row of R is zero

# Rank-revealing QR with column pivoting: A P = Q R with |r_11| >= |r_22| >= ...
D = np.array([[1., 2., 3.], [2., 4., 6.], [1., 0., 1.]])        # column 3 = column 1 + column 2: rank 2
Qp, Rp, perm = qr(D, pivoting=True)
print(np.round(np.abs(np.diag(Rp)), 4) + 0.0, perm)              # [6.7823 0.6594 0.    ] [2 0 1]: one zero pivot => rank 2

The signs in the first printout depend on the LAPACK build, so they may differ on your machine. After sign normalization, every machine gives the same result.

Key takeaways

• $A = QR$ with orthonormal columns in $Q$ and upper triangular $R$; it is Gram-Schmidt written as a matrix identity, with $r_{jk} = \langle \mathbf{x}_k, \mathbf{q}_j \rangle$.
• $R$ is triangular because of nested spans; with $r_{kk} > 0$ the factorization is unique.
• Reduced QR gives an orthonormal basis of $\operatorname{Col} A$; complete QR also gives one of $\operatorname{Nul} A^T$.
• QR solves square systems and least squares problems through $R\mathbf{x} = Q^T\mathbf{b}$, and column pivoting reveals numerical rank.
• Library QR may return negative diagonal entries in $R$; flip signs with $D = \operatorname{diag}(\pm 1)$ if the convention matters.

🧪 Practice

1. Compute the QR factorization of $A = \begin{bmatrix} 3 & -6 \ 4 & -8 \ 0 & 1 \end{bmatrix}$ by hand, and verify $QR = A$.
2. Show that if $A = QR$ with $A$ square and invertible, then $A^TA = R^TR$. This is the Cholesky factorization of $A^TA$ (Section 10.3).
3. Use QR to solve $\begin{bmatrix} 1 & 1 \ 1 & -1 \end{bmatrix}\mathbf{x} = \begin{bmatrix} 3 \ 1 \end{bmatrix}$, and explain why no elimination is needed beyond back substitution.
4. Interview: Why do numerical libraries compute QR with Householder reflections rather than Gram-Schmidt, even though Gram-Schmidt is easier to explain? When might you still choose Gram-Schmidt? Hint: compare loss of orthogonality (previous topic) and whether $Q$ is needed explicitly; Gram-Schmidt produces the basis one vector at a time, which suits Krylov methods that generate vectors on the fly.

Householder Reflections

Why it matters. Gram-Schmidt builds $Q$ by subtracting projections: it multiplies $A$ by triangular matrices until the columns are orthonormal ($AR^{-1} = Q$). Householder's method does the opposite. It multiplies $A$ by orthogonal matrices until it becomes triangular ($Q^TA = R$). Orthogonal matrices never amplify errors (Section 9.4), so the computed $Q$ is orthogonal to working precision regardless of how ill-conditioned $A$ is. Householder QR is the standard algorithm in LAPACK. The same reflections reduce matrices to Hessenberg and tridiagonal form, the first step of every dense eigenvalue solver.

Analogy. A Householder reflection is a mirror. Given any vector, you can always position a mirror so that the vector's reflection lands exactly on the first coordinate axis, keeping its length. Do this for column 1 of $A$, then for the rest of column 2 with a smaller mirror, and so on. After $n$ mirrors the matrix is triangular.

Definition. For a nonzero vector $\mathbf{v} \in \mathbb{R}^m$, the Householder reflection is

$$H = I - 2\frac{\mathbf{v}\mathbf{v}^T}{\mathbf{v}^T\mathbf{v}}.$$

It reflects across the hyperplane $\mathbf{v}^\perp$. The component of any vector along $\mathbf{v}$ is reversed, and the components perpendicular to $\mathbf{v}$ are unchanged:

Property	Statement
Action	$H\mathbf{v} = -\mathbf{v}$; $H\mathbf{w} = \mathbf{w}$ for $\mathbf{w} \perp \mathbf{v}$
Symmetric	$H^T = H$
Orthogonal	$H^TH = I$, so $H$ preserves lengths and angles
Involution	$H^2 = I$: reflecting twice undoes the reflection, so $H^{-1} = H$
Spectrum	eigenvalue $-1$ (once, direction $\mathbf{v}$) and $1$ (with multiplicity $m - 1$); $\det H = -1$

Since $H = I - 2P$, where $P$ is the projection onto $\operatorname{span}{\mathbf{v}}$, the reflection overshoots the projection by an equal distance on the other side.

                      x
                      *
                     /|
                    / |                   mirror: the hyperplane v-perp
                   /  |                  /
                  /   |                 /
                 /    v = x - y        /
                /     |               /
               O------+--------------*------> e1
                      |            y = H x = alpha e1
          (x and y have the same length; the mirror bisects the angle between them)

Key construction: reflect $\mathbf{x}$ onto a coordinate axis. If $\lVert \mathbf{x} \rVert = \lVert \mathbf{y} \rVert$ and $\mathbf{v} = \mathbf{x} - \mathbf{y}$, then $H\mathbf{x} = \mathbf{y}$.

Proof. $\mathbf{v}^T\mathbf{v} = \lVert \mathbf{x} \rVert^2 - 2\mathbf{y}^T\mathbf{x} + \lVert \mathbf{y} \rVert^2 = 2(\lVert \mathbf{x} \rVert^2 - \mathbf{y}^T\mathbf{x}) = 2\mathbf{v}^T\mathbf{x}$. Hence $H\mathbf{x} = \mathbf{x} - 2\mathbf{v}\dfrac{\mathbf{v}^T\mathbf{x}}{\mathbf{v}^T\mathbf{v}} = \mathbf{x} - \mathbf{v} = \mathbf{y}$. $\blacksquare$

Take $\mathbf{y} = \alpha\mathbf{e}_1$ with $\alpha = \pm\lVert \mathbf{x} \rVert$. Choose the sign opposite to $x_1$: $\alpha = -\operatorname{sign}(x_1)\lVert \mathbf{x} \rVert$. Then $v_1 = x_1 - \alpha = x_1 + \operatorname{sign}(x_1)\lVert \mathbf{x} \rVert$ adds two numbers of the same sign. The other choice would subtract two nearly equal numbers whenever $\mathbf{x}$ is close to the $\mathbf{e}_1$ axis, which causes catastrophic cancellation.

Example. $\mathbf{x} = (3, 4, 0)$: $\lVert \mathbf{x} \rVert = 5$, $\alpha = -5$, $\mathbf{v} = (8, 4, 0)$. Then $\mathbf{v}^T\mathbf{x} = 40$, $\mathbf{v}^T\mathbf{v} = 80$, and $H\mathbf{x} = \mathbf{x} - \mathbf{v} = (-5, 0, 0)$.

Householder QR. For $k = 1, \dots, n$: build $H_k$ that zeros the entries below the diagonal in column $k$, acting only on rows $k$ to $m$, and apply it to the remaining submatrix. Earlier columns are not disturbed, because $H_k$ leaves rows $1, \dots, k - 1$ alone and the entries it touches in those columns are already zero.

   A (5 x 3)          H1 A               H2 H1 A            H3 H2 H1 A = R
   x x x              * * *              * * *              * * *
   x x x              0 * *              0 * *              0 * *
   x x x     -->      0 * *     -->      0 0 *     -->      0 0 *
   x x x              0 * *              0 0 *              0 0 0
   x x x              0 * *              0 0 *              0 0 0

   Q^T = H3 H2 H1,   so   A = (H1 H2 H3) R = Q R      (each H_k is its own inverse)

Never form $H$. Applying $H$ to a vector costs $O(m)$: $H\mathbf{x} = \mathbf{x} - \mathbf{v}\left(2\mathbf{v}^T\mathbf{x}/\mathbf{v}^T\mathbf{v}\right)$. Forming the $m \times m$ matrix costs $O(m^2)$ to store and $O(m^2)$ per vector to apply. Libraries store only the vectors $\mathbf{v}_k$, compactly, in the zeroed part of $A$, and apply $Q$ or $Q^T$ on demand. Householder QR costs about $2mn^2 - \tfrac{2}{3}n^3$ flops, the same as Gram-Schmidt when $m \gg n$, or $\tfrac{4}{3}n^3$ for a square matrix.

import numpy as np

def householder_vector(x):
    """Return v such that (I - 2 v v^T / v^T v) x = alpha e1, with alpha = -sign(x1) ||x||."""
    alpha = -np.copysign(np.linalg.norm(x), x[0])   # sign chosen to avoid cancellation in x1 - alpha
    v = x.astype(float).copy()
    v[0] -= alpha
    return v, alpha

def apply_reflection(v, X):
    """Compute H X = X - 2 v (v^T X) / (v^T v) without forming H: O(m n) instead of O(m^2 n)."""
    return X - np.outer(v, 2 * (v @ X) / (v @ v))

x = np.array([3., 4., 0.])
v, alpha = householder_vector(x)
print(v, alpha, apply_reflection(v, x[:, None]).ravel())     # [8. 4. 0.] -5.0 [-5.  0.  0.]

H = np.eye(3) - 2 * np.outer(v, v) / (v @ v)
print(np.allclose(H, H.T), np.allclose(H @ H, np.eye(3)), round(np.linalg.det(H)))   # True True -1

def householder_qr(A):
    """Householder QR: returns Q (m x m) and R (m x n) with A = Q R."""
    m, n = A.shape
    R = A.astype(float).copy()
    vs = []
    for k in range(min(n, m - 1)):
        v, _ = householder_vector(R[k:, k])
        R[k:, k:] = apply_reflection(v, R[k:, k:])  # zero out column k below the diagonal
        vs.append(v)
    Q = np.eye(m)
    for k in reversed(range(len(vs))):              # Q = H_1 H_2 ... H_n, applied to I from the right end
        Q[k:, :] = apply_reflection(vs[k], Q[k:, :])
    return Q, R

A = np.array([[12., -51., 4.],
              [6., 167., -68.],
              [-4., 24., -41.]])
Q, R = householder_qr(A)
print(np.round(R, 6) + 0.0)
# [[ -14.  -21.   14.]
#  [   0. -175.   70.]
#  [   0.    0.  -35.]]
print(np.allclose(Q @ R, A), np.allclose(Q.T @ Q, np.eye(3)))   # True True

# Stability on the nearly dependent Lauchli matrix (compare with Gram-Schmidt in the previous topic)
eps = 1e-8
L = np.array([[1, 1, 1], [eps, 0, 0], [0, eps, 0], [0, 0, eps]])
QL, RL = householder_qr(L)
print(f"{np.linalg.norm(QL.T @ QL - np.eye(4)):.1e}", np.allclose(QL @ RL, L))   # 4.3e-16 True

The negative diagonal of $R$ comes from the sign rule $\alpha = -\operatorname{sign}(x_1)\lVert \mathbf{x} \rVert$. This is why library QR routines often return negative $r_{kk}$.

Key takeaways

• $H = I - 2\mathbf{v}\mathbf{v}^T/(\mathbf{v}^T\mathbf{v})$ reflects across $\mathbf{v}^\perp$; it is symmetric, orthogonal, its own inverse, and has determinant $-1$.
• With $\mathbf{v} = \mathbf{x} - \alpha\mathbf{e}_1$ and $\alpha = -\operatorname{sign}(x_1)\lVert \mathbf{x} \rVert$, $H$ maps $\mathbf{x}$ to $\alpha\mathbf{e}_1$ without cancellation.
• Householder QR triangularizes $A$ by $n$ reflections, $Q^T = H_n \cdots H_1$; it is backward stable and $Q$ is orthogonal to machine precision.
• Store the vectors $\mathbf{v}_k$, never the matrices $H_k$; applying a reflection costs $O(m)$ per vector.

🧪 Practice

1. Find the Householder reflection that maps $\mathbf{x} = (1, 2, 2)$ to a multiple of $\mathbf{e}_1$. Write out $H$ and check $H\mathbf{x}$, $H^T = H$, and $H^2 = I$.
2. Compute by hand the first Householder step for $A = \begin{bmatrix} 1 & 1 \ 2 & 0 \ 2 & 3 \end{bmatrix}$, and verify that the first column of $H_1A$ is $(-3, 0, 0)$.
3. Show that every Householder reflection has determinant $-1$, and deduce that $\det Q = (-1)^n$ for the $Q$ produced by $n$ reflections.
4. Interview: Your Householder implementation produces NaN for some inputs, but only when a column below the diagonal is already zero. What went wrong, and how do you fix it? Hint: if $\mathbf{x} = \mathbf{0}$, then $\mathbf{v} = \mathbf{0}$ and $\mathbf{v}^T\mathbf{v} = 0$; also check np.copysign when $x_1 = 0$. Skip the reflection (use $H = I$) when the column is already zero.

Givens Rotations

Why it matters. A Householder reflection zeros a whole column at once and modifies every row it touches. Sometimes only a single entry needs to be zeroed. This happens with sparse or structured matrices (Hessenberg, tridiagonal, banded), where most entries are already zero, and when a QR factorization must be updated because one data row arrives. A Givens rotation zeros one chosen entry by rotating just two rows, leaves everything else alone, and lets independent rotations run in parallel.

Analogy. Householder is a bulldozer that clears a whole column in one pass. Givens is a shovel that removes one entry at a time. On a mostly empty field, the shovel does less work.

Definition. A Givens rotation $G(i, j, \theta)$ is the identity matrix except in rows and columns $i, j$, where it contains the $2 \times 2$ rotation

$$\begin{bmatrix} c & s \ -s & c \end{bmatrix}, \qquad c = \cos\theta, \quad s = \sin\theta.$$

Multiplying a matrix by $G(i, j, \theta)$ on the left changes only rows $i$ and $j$, so it costs $O(n)$ for an $m \times n$ matrix. $G$ is orthogonal with $\det G = 1$: a rotation in the $(i, j)$ coordinate plane.

Zeroing one entry. To rotate $(a, b)$ onto $(r, 0)$, choose

$$r = \sqrt{a^2 + b^2}, \qquad c = \frac{a}{r}, \qquad s = \frac{b}{r}, \qquad \begin{bmatrix} c & s \ -s & c \end{bmatrix}\begin{bmatrix} a \ b \end{bmatrix} = \begin{bmatrix} (a^2 + b^2)/r \ (-ab + ab)/r \end{bmatrix} = \begin{bmatrix} r \ 0 \end{bmatrix}.$$

No angle or trigonometric function is ever computed. Compute $r$ with hypot, which avoids overflow when $a$ or $b$ is huge.

   Givens QR of a 4 x 3 matrix: zero one entry at a time, bottom-up in each column

   x x x      x x x      x x x      * * *      * * *      * * *      * * *
   x x x      x x x      * * *      0 * *      0 * *      0 * *      0 * *
   x x x  ->  * * *  ->  0 * *  ->  0 * *  ->  0 * *  ->  0 * *  ->  0 0 *
   x x x      0 * *      0 * *      0 * *      0 0 *      0 0 *      0 0 0
          G(3,4)     G(2,3)     G(1,2)     G(3,4)     G(2,3)     G(3,4)

   rotations that touch disjoint row pairs can run at the same time

QR by Givens rotations. Zero the subdiagonal entries column by column, working from the bottom up within each column. A dense $m \times n$ matrix needs about $mn - n^2/2$ rotations and about $3mn^2 - n^3$ flops, roughly 50% more than Householder. Structured matrices are where Givens wins:

Matrix structure	Entries to zero	Givens cost	Householder cost
dense $n \times n$	$n(n-1)/2$	about $2n^3$	about $\tfrac{4}{3}n^3$
upper Hessenberg $n \times n$	$n - 1$	$O(n^2)$	$O(n^2)$ with care, but each reflection touches more
tridiagonal	$n - 1$	$O(n)$	wasteful (fill-in)
$R$ plus one new row	$n$	$O(n^2)$	refactor from scratch: $O(mn^2)$

The Hessenberg case is the inner step of the QR eigenvalue algorithm (Section 12.3), and the new-row case is recursive least squares.

Updating QR when data arrive. Suppose $A = QR$ and a new observation row $\mathbf{a}^T$ is appended. Stacking $\begin{bmatrix} R \ \mathbf{a}^T \end{bmatrix}$ gives an upper triangular block plus one full row. $n$ Givens rotations, one per column, rotate the new row into $R$ and produce the triangular factor of the enlarged matrix. Old data never need to be revisited. The check is $R_{\text{new}}^TR_{\text{new}} = A^TA + \mathbf{a}\mathbf{a}^T$.

import numpy as np

def givens(a, b):
    """Return (c, s, r) with [[c, s], [-s, c]] @ [a, b] = [r, 0]."""
    r = np.hypot(a, b)                  # sqrt(a^2 + b^2) without overflow or underflow
    if r == 0:
        return 1.0, 0.0, 0.0
    return a / r, b / r, r

c, s, r = givens(3., 4.)
print(c, s, r, np.round(np.array([[c, s], [-s, c]]) @ [3., 4.], 12) + 0.0)    # 0.6 0.8 5.0 [5. 0.]

def apply_givens(M, i, j, c, s):
    """Rotate rows i and j of M in place: only 2 rows change, O(n) work."""
    Mi, Mj = M[i].copy(), M[j].copy()
    M[i], M[j] = c * Mi + s * Mj, -s * Mi + c * Mj

def givens_qr(A):
    """QR by Givens rotations, zeroing subdiagonal entries column by column, bottom to top."""
    m, n = A.shape
    R, Qt = A.astype(float).copy(), np.eye(m)
    count = 0
    for k in range(n):
        for i in range(m - 1, k, -1):
            if R[i, k] != 0:                        # skip entries that are already zero (sparsity)
                c, s, _ = givens(R[i - 1, k], R[i, k])
                apply_givens(R, i - 1, i, c, s)
                apply_givens(Qt, i - 1, i, c, s)    # accumulate Q^T
                count += 1
    return Qt.T, R, count

# Upper Hessenberg matrix: only n - 1 entries below the diagonal, so only n - 1 rotations
Hs = np.triu(np.arange(1., 26.).reshape(5, 5), -1)
Q, R, count = givens_qr(Hs)
print(count, np.allclose(Q @ R, Hs), np.allclose(np.tril(R, -1), 0))    # 4 True True

# Dense 5 x 5: all 10 subdiagonal entries need a rotation
_, _, count = givens_qr(np.random.default_rng(0).standard_normal((5, 5)))
print(count)                                                          # 10

# Updating a QR factorization when a new data row arrives (streaming least squares)
rng = np.random.default_rng(1)
A = rng.standard_normal((100, 3))
R = np.linalg.qr(A, mode="r")                     # 3 x 3 upper triangular
new_row = rng.standard_normal(3)
Rn = np.vstack([R, new_row])                      # 4 x 3: triangular except the last row
for k in range(3):                                # one rotation per column restores triangular form
    c, s, _ = givens(Rn[k, k], Rn[3, k])
    apply_givens(Rn, k, 3, c, s)
R_updated = Rn[:3]
A2 = np.vstack([A, new_row])
print(np.allclose(R_updated.T @ R_updated, A2.T @ A2))               # True: same R^T R = A^T A

Comparison of the QR algorithms.

Method	Orthogonality of computed $Q$	Flops (dense $m \times n$)	Best for
Classical Gram-Schmidt	about $\kappa^2\varepsilon$ (can fail)	$2mn^2$	teaching; with reorthogonalization, parallel Krylov methods
Modified Gram-Schmidt	about $\kappa\varepsilon$	$2mn^2$	building an orthonormal basis one vector at a time
Householder	about $\varepsilon$	$2mn^2 - \tfrac{2}{3}n^3$	general dense QR (LAPACK default)
Givens	about $\varepsilon$	$3mn^2 - n^3$	sparse or structured matrices, updates, parallel zeroing

Key takeaways

• A Givens rotation acts on two rows and zeros one chosen entry, using $c = a/r$ and $s = b/r$ with $r = \operatorname{hypot}(a, b)$.
• Givens QR is about 50% more expensive than Householder for dense matrices, but much cheaper for Hessenberg, tridiagonal, and banded matrices.
• Appending one row to a QR factorization takes $n$ rotations, which makes Givens the tool for recursive least squares.
• Rotations on disjoint row pairs commute and can run in parallel; like Householder, Givens QR is backward stable.

🧪 Practice

1. Find the Givens rotation that zeros the second entry of $(1, 1)$ and of $(-5, 12)$. What are $c$, $s$, and $r$ in each case?
2. Compute the QR factorization of $\begin{bmatrix} 3 & 1 \ 4 & 2 \end{bmatrix}$ with one Givens rotation, and compare it with Householder QR.
3. Count the Givens rotations needed to triangularize an $n \times n$ tridiagonal matrix, and show that $R$ has at most two nonzero superdiagonals.
4. Interview: A sensor streams one new measurement row per millisecond into a linear regression with 20 features, and you must report the current least-squares coefficients after every row. How do you avoid refitting from scratch each time? Hint: keep only the $20 \times 20$ factor $R$ and the vector $Q^T\mathbf{b}$; each new row costs 20 Givens rotations ($O(n^2)$ work) plus one triangular solve, and exponential forgetting just scales $R$ before each update.

<a id="94-orthogonal-and-unitary-matrices"></a>

9.4 Orthogonal and Unitary Matrices

This subchapter studies the square matrices with orthonormal columns: the rigid motions of linear algebra. They preserve lengths and angles, their inverse is a transpose, and they form groups that describe every rotation and reflection.

Orthogonal Matrices and Properties

Why it matters. Rotations, reflections, and permutations move objects without stretching or distorting them. In matrix form, they are exactly the square matrices with orthonormal columns. Inverting one costs only a transpose, it never magnifies errors, and products of such matrices stay in the family. These properties make orthogonal matrices the building blocks of stable numerical algorithms (Householder, Givens, QR, SVD), of 3D graphics and robotics, and of every change of basis between orthonormal bases.

Analogy. Picking up a rigid object and setting it down in a different orientation changes where its points are, but not the distances between them or the angles at its corners. An orthogonal matrix is that kind of motion, with the origin held fixed.

Definition. A real $n \times n$ matrix $Q$ is orthogonal if

$$Q^TQ = I.$$

The name is a historical accident: the columns must be orthonormal, not merely orthogonal. A matrix whose columns are orthogonal but not unit length satisfies $Q^TQ = D$ (diagonal), not $I$.

Equivalent conditions. For a square real matrix $Q$, the following are equivalent:

1. $Q^TQ = I$: the columns are orthonormal.
2. $QQ^T = I$: the rows are orthonormal.
3. $Q$ is invertible and $Q^{-1} = Q^T$.
4. $(Q\mathbf{x}) \cdot (Q\mathbf{y}) = \mathbf{x} \cdot \mathbf{y}$ for all $\mathbf{x}, \mathbf{y}$: dot products are preserved.
5. $\lVert Q\mathbf{x} \rVert = \lVert \mathbf{x} \rVert$ for all $\mathbf{x}$: lengths are preserved.

Proof sketch. (1) $\Rightarrow$ (3): a one-sided inverse of a square matrix is two-sided (Section 4.3), and (3) gives (2). (1) $\Rightarrow$ (4): $(Q\mathbf{x})^T(Q\mathbf{y}) = \mathbf{x}^TQ^TQ\mathbf{y} = \mathbf{x}^T\mathbf{y}$. (4) $\Rightarrow$ (5): take $\mathbf{y} = \mathbf{x}$. (5) $\Rightarrow$ (4): polarization (Section 9.1). (4) $\Rightarrow$ (1): take $\mathbf{x} = \mathbf{e}_i$ and $\mathbf{y} = \mathbf{e}j$, so that $(Q^TQ){ij} = \mathbf{q}_i \cdot \mathbf{q}j = \delta{ij}$. $\blacksquare$

Condition (2) is not obvious. If the columns of a square matrix are orthonormal, then so are its rows. The matrix $\tfrac{1}{3}\begin{bmatrix} 1 & 2 & 2 \ 2 & 1 & -2 \ 2 & -2 & 1 \end{bmatrix}$ illustrates it.

Properties.

Property	Reason
$\det Q = \pm 1$	$1 = \det(Q^TQ) = (\det Q)^2$
Every eigenvalue has $\lvert \lambda \rvert = 1$	$Q\mathbf{v} = \lambda\mathbf{v}$ gives $\lVert \mathbf{v} \rVert = \lVert Q\mathbf{v} \rVert = \lvert \lambda \rvert\lVert \mathbf{v} \rVert$ (with complex $\mathbf{v}$ if needed)
Real eigenvalues are $\pm 1$	the only real numbers of absolute value 1
$Q_1Q_2$, $Q^T$, $Q^{-1}$ are orthogonal	$(Q_1Q_2)^T(Q_1Q_2) = Q_2^TQ_1^TQ_1Q_2 = I$
Condition number $\kappa(Q) = 1$	$\lVert Q \rVert = \lVert Q^{-1} \rVert = 1$ in the 2-norm (Section 11.3)
Angles preserved	dot products and lengths are both preserved

A determinant of $\pm 1$ is necessary but not sufficient: $\begin{bmatrix} 2 & 0 \ 0 & 1/2 \end{bmatrix}$ has determinant $1$ and is not orthogonal.

All $2 \times 2$ orthogonal matrices. The first column is a unit vector $(\cos\theta, \sin\theta)$. The second column is a unit vector perpendicular to it, so there are two choices:

$$\underbrace{\begin{bmatrix} \cos\theta & -\sin\theta \ \sin\theta & \cos\theta \end{bmatrix}}{\text{rotation by } \theta, \ \det = +1} \qquad\qquad \underbrace{\begin{bmatrix} \cos\theta & \sin\theta \ \sin\theta & -\cos\theta \end{bmatrix}}{\text{reflection across the line at angle } \theta/2, \ \det = -1}$$

   rotation (det +1): orientation kept         reflection (det -1): orientation flipped

        e2   Q e1                                     line at angle theta/2
         ^   ^                                       /
         |  /                                       /    Q e1
         | / theta                                 /   ^
         |/                                       /  /
         O------> e1                             O------> e1
                                                  \
   counterclockwise order of (Q e1, Q e2)          v  Q e2   (clockwise order:
   is preserved                                               a mirror image)

Examples of orthogonal matrices.

Matrix	Determinant	Source
Rotation in a coordinate plane	$+1$	Givens (Section 9.3)
$I - 2\mathbf{v}\mathbf{v}^T/\mathbf{v}^T\mathbf{v}$	$-1$	Householder (Section 9.3)
Permutation matrix	$\pm 1$ (the sign of the permutation)	row swaps (Section 4.4)
Columns = orthonormal eigenvectors of a symmetric matrix	$\pm 1$	spectral theorem (Chapter 10)
$U$ and $V$ in the SVD $A = U\Sigma V^T$	$\pm 1$	Chapter 11
Normalized Hadamard matrix $\tfrac{1}{2}\begin{bmatrix} 1 & 1 & 1 & 1 \ 1 & -1 & 1 & -1 \ 1 & 1 & -1 & -1 \ 1 & -1 & -1 & 1 \end{bmatrix}$	$\pm 1$	Walsh transforms, error-correcting codes

Square vs. rectangular. A tall $m \times n$ matrix with orthonormal columns satisfies $Q^TQ = I_n$, but $QQ^T$ is the projection onto $\operatorname{Col} Q$, not $I_m$. It preserves lengths ($\lVert Q\mathbf{x} \rVert = \lVert \mathbf{x} \rVert$) but has no two-sided inverse. Calling such a matrix "orthogonal" is common but inaccurate. "Matrix with orthonormal columns" is the safe phrase.

import numpy as np

def rotation(t):
    return np.array([[np.cos(t), -np.sin(t)], [np.sin(t), np.cos(t)]])

def reflection(t):
    """Reflection across the line through the origin at angle t/2."""
    return np.array([[np.cos(t), np.sin(t)], [np.sin(t), -np.cos(t)]])

for name, Q in [("rotation", rotation(np.pi / 3)), ("reflection", reflection(np.pi / 3)),
                ("permutation", np.eye(3)[[2, 0, 1]])]:
    print(f"{name:12s} Q^TQ=I: {np.allclose(Q.T @ Q, np.eye(len(Q)))}  QQ^T=I: {np.allclose(Q @ Q.T, np.eye(len(Q)))}"
          f"  det = {np.linalg.det(Q):+.0f}  |eigenvalues| = {np.round(np.abs(np.linalg.eigvals(Q)), 12)}")
# rotation     Q^TQ=I: True  QQ^T=I: True  det = +1  |eigenvalues| = [1. 1.]
# reflection   Q^TQ=I: True  QQ^T=I: True  det = -1  |eigenvalues| = [1. 1.]
# permutation  Q^TQ=I: True  QQ^T=I: True  det = +1  |eigenvalues| = [1. 1. 1.]

# The inverse is the transpose: solving Q x = b costs one matrix-vector product
Q = rotation(0.7) @ reflection(0.2)            # products of orthogonal matrices are orthogonal
b = np.array([1., 2.])
print(np.allclose(np.linalg.solve(Q, b), Q.T @ b))                   # True

# Lengths, dot products, and angles are preserved; condition number is 1
u, v = np.array([3., 4.]), np.array([-1., 2.])
print(round(np.linalg.norm(Q @ u), 12), np.isclose((Q @ u) @ (Q @ v), u @ v))  # 5.0 True
print(round(np.linalg.cond(Q), 12))                                  # 1.0

# Repeated multiplication: orthogonal matrices never blow up or shrink vectors
x = np.array([1., 0.])
for _ in range(10000):
    x = Q @ x
print(np.linalg.norm(x))              # 1.000000000001113: only rounding drift (~1e-12) after 10000 steps

Key takeaways

• $Q$ is orthogonal when $Q^TQ = I$; for square $Q$, this is equivalent to $QQ^T = I$, $Q^{-1} = Q^T$, preservation of dot products, and preservation of lengths.
• $\det Q = \pm 1$ and every eigenvalue lies on the unit circle; $\det Q = \pm 1$ alone does not make a matrix orthogonal.
• The $2 \times 2$ orthogonal matrices are rotations ($\det = +1$) and reflections ($\det = -1$).
• Orthogonal matrices have condition number 1, so multiplying by them never amplifies errors.

🧪 Practice

1. Show that $\tfrac{1}{3}\begin{bmatrix} 1 & 2 & 2 \ 2 & 1 & -2 \ 2 & -2 & 1 \end{bmatrix}$ is orthogonal, find its determinant, and decide whether it is a rotation.
2. Prove that if $Q$ is orthogonal and upper triangular, then $Q$ is diagonal with entries $\pm 1$.
3. Find all values of $a$ and $b$ for which $\begin{bmatrix} a & b \ -b & a \end{bmatrix}$ is orthogonal, and interpret them geometrically.
4. Interview: Why is it preferable to compute $Q^T\mathbf{b}$ rather than $Q^{-1}\mathbf{b}$, or to call solve(Q, b), when $Q$ is known to be orthogonal? Is there any case where you should not trust $Q^T$ as the inverse? Hint: the cost is $O(n^2)$ instead of $O(n^3)$ with no error amplification; but if $Q$ came from an unstable process (classical Gram-Schmidt on ill-conditioned data), it may not be orthogonal to working precision.

Unitary Matrices

Why it matters. Quantum gates, the discrete Fourier transform, and the eigenvector matrices of Hermitian matrices are complex, and the complex version of "orthogonal" is "unitary". Unitary matrices preserve the complex inner product. In quantum mechanics, this means they preserve total probability, which is why every closed quantum system evolves by a unitary matrix. In signal processing, it means the energy of a signal is the same in the time domain and in the frequency domain.

Definition. A complex $n \times n$ matrix $U$ is unitary if

$$U^U = I, \qquad \text{where } U^ = \overline{U}^T \text{ (Section 4.2)}.$$

The equivalent conditions parallel the orthogonal case: the columns are orthonormal in $\mathbb{C}^n$, the rows are orthonormal, $U^{-1} = U^*$, $\langle U\mathbf{x}, U\mathbf{y} \rangle = \langle \mathbf{x}, \mathbf{y} \rangle$, and $\lVert U\mathbf{x} \rVert = \lVert \mathbf{x} \rVert$. A real matrix is unitary exactly when it is orthogonal.

Properties.

Property	Orthogonal ($\mathbb{R}$)	Unitary ($\mathbb{C}$)
Definition	$Q^TQ = I$	$U^*U = I$
Inverse	$Q^T$	$U^*$
Determinant	$\pm 1$	$\lvert\det U\rvert = 1$, so $\det U = e^{i\phi}$
Eigenvalues	on the unit circle	on the unit circle, $\lambda = e^{i\theta}$
Diagonalizable?	not always over $\mathbb{R}$ (rotations have complex eigenvalues)	always, by a unitary matrix, with orthogonal eigenvectors (normal matrices, Chapter 10)
Group	$O(n)$	$U(n)$

The pitfall. Testing U.T @ U instead of U.conj().T @ U is the most common bug with complex matrices. The plain transpose is the wrong test: the unitary DFT matrix below passes the correct test and fails the wrong one.

Examples.

• Diagonal phases. $\operatorname{diag}(e^{i\theta_1}, \dots, e^{i\theta_n})$. Every unitary matrix is unitarily similar to one of these.
• The discrete Fourier transform. With $\omega = e^{-2\pi i/n}$, the matrix $F_{jk} = \omega^{jk}/\sqrt{n}$ ($j, k = 0, \dots, n - 1$) is unitary. Its columns are sampled complex exponentials, mutually orthogonal because $\sum_{j=0}^{n-1}\omega^{j(k - l)} = 0$ for $k \ne l$ (a geometric series of roots of unity). The FFT computes $F\mathbf{x}$ in $O(n\log n)$ operations. Unitarity is Parseval's theorem for the DFT: $\lVert F\mathbf{x} \rVert = \lVert \mathbf{x} \rVert$ (Section 13.3).
• Quantum gates. The Hadamard gate $\tfrac{1}{\sqrt{2}}\begin{bmatrix} 1 & 1 \ 1 & -1 \end{bmatrix}$, the phase gate $\begin{bmatrix} 1 & 0 \ 0 & i \end{bmatrix}$, and the Pauli matrices. A quantum state is a unit vector $\boldsymbol{\psi}$ whose squared entries $\lvert \psi_k \rvert^2$ are probabilities. Unitarity keeps their sum equal to $1$.
• Exponentials of Hermitian matrices. If $H^* = H$, then $U = e^{iH}$ is unitary, because $(e^{iH})^* = e^{-iH^*} = e^{-iH}$. Conversely, every unitary matrix has this form (Chapter 10). In quantum mechanics, $e^{-iHt/\hbar}$ is time evolution under the Hamiltonian $H$. This parallels $e^{i\theta}$ lying on the unit circle for real $\theta$.

import numpy as np

n = 8
j, k = np.meshgrid(np.arange(n), np.arange(n))
omega = np.exp(-2j * np.pi / n)
F = omega ** (j * k) / np.sqrt(n)                 # unitary DFT matrix, F_jk = w^(jk) / sqrt(n)

print(np.allclose(F.conj().T @ F, np.eye(n)))    # True: unitary
print(np.allclose(F.T @ F, np.eye(n)))           # False: plain transpose is the wrong test for complex matrices

x = np.random.default_rng(0).standard_normal(n)
print(np.allclose(F @ x, np.fft.fft(x, norm="ortho")))           # True: the FFT computes F x in O(n log n)
print(np.isclose(np.linalg.norm(F @ x), np.linalg.norm(x)))      # True: Parseval, energy is preserved

ev = np.linalg.eigvals(F)
print(np.allclose(np.abs(ev), 1), np.round(abs(np.linalg.det(F)), 12))   # True 1.0
print([int(np.sum(np.isclose(ev, t))) for t in (1, -1, 1j, -1j)])   # multiplicities of the eigenvalues 1, -1, i, -i

# Quantum gates are unitary, so total probability stays 1
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)      # Hadamard gate (real, so also orthogonal)
S = np.array([[1, 0], [0, 1j]])                   # phase gate (complex)
psi = np.array([0.6, 0.8j])                       # |0.6|^2 + |0.8|^2 = 1
for U in [H, S, H @ S @ H]:
    out = U @ psi
    print(np.round(np.sum(np.abs(out)**2), 12), end=" ")
print()                                           # 1.0 1.0 1.0: total probability preserved

# A unitary matrix as the exponential of i times a Hermitian matrix (Chapter 10)
from scipy.linalg import expm
Hm = np.array([[1, 2 - 1j], [2 + 1j, -1]])       # Hermitian: Hm^* = Hm
U = expm(1j * Hm)
print(np.allclose(U.conj().T @ U, np.eye(2)))     # True

Key takeaways

• $U$ is unitary when $U^U = I$: orthonormal columns in $\mathbb{C}^n$, with $U^{-1} = U^$.
• Unitary matrices preserve complex inner products and norms; their eigenvalues and determinants have absolute value 1.
• Real unitary matrices are exactly the orthogonal matrices; unitary matrices are always diagonalizable by unitary matrices.
• The normalized DFT matrix is unitary (Parseval), quantum gates are unitary (probability is conserved), and $e^{iH}$ is unitary for Hermitian $H$.
• Always test with the conjugate transpose: U.conj().T @ U.

🧪 Practice

1. Show that $\tfrac{1}{\sqrt{2}}\begin{bmatrix} 1 & i \ i & 1 \end{bmatrix}$ is unitary, and find its eigenvalues. Check that they lie on the unit circle.
2. Write out the $4 \times 4$ unitary DFT matrix and verify by hand that its first two columns are orthogonal.
3. Prove that the product of two unitary matrices is unitary and that $\lvert\det U\rvert = 1$.
4. Interview: In a quantum circuit simulator, the total probability of a state drifts from 1 to 0.97 after many gates. List the possible causes. Hint: a gate that is not exactly unitary (rounded matrix entries, a missing $1/\sqrt{2}$), a transpose used where a conjugate transpose was needed, or accumulated floating-point error; renormalizing the state hides the bug rather than fixing it.

Isometries and Norm Preservation

Why it matters. "Which transformations preserve distance?" is the question behind rigid-body motion, the registration of medical images, the alignment of 3D scans, and the numerical stability of algorithms. The answer is short: up to a translation, exactly the orthogonal (or unitary) maps. One consequence is that registering two point clouds reduces to finding the best orthogonal matrix, which the SVD solves in a few lines.

Definitions. A linear map $T : V \to W$ between inner product spaces is a (linear) isometry if $\lVert T\mathbf{v} \rVert = \lVert \mathbf{v} \rVert$ for every $\mathbf{v}$. A general map $f$ between them is distance-preserving if $\lVert f(\mathbf{u}) - f(\mathbf{v}) \rVert = \lVert \mathbf{u} - \mathbf{v} \rVert$ for all $\mathbf{u}, \mathbf{v}$.

Theorem 1 (norms determine inner products). A linear map preserves norms if and only if it preserves inner products.

Proof. Preserving inner products clearly preserves norms. Conversely, by the polarization identity (Section 9.1), $\langle T\mathbf{u}, T\mathbf{v} \rangle$ is a fixed combination of $\lVert T\mathbf{u} + i^kT\mathbf{v} \rVert^2 = \lVert T(\mathbf{u} + i^k\mathbf{v}) \rVert^2 = \lVert \mathbf{u} + i^k\mathbf{v} \rVert^2$, and the same combination gives $\langle \mathbf{u}, \mathbf{v} \rangle$. (Over $\mathbb{R}$, use $k = 0, 2$ only.) $\blacksquare$

So a linear isometry automatically preserves angles and orthogonality, and maps orthonormal sets to orthonormal sets. Between $\mathbb{R}^n$ and $\mathbb{R}^m$, a linear map is an isometry exactly when its matrix has orthonormal columns. It is injective (its kernel is ${\mathbf{0}}$, since only $\mathbf{0}$ has length $0$), and it is surjective only when $m = n$.

Theorem 2 (rigid motions). Every distance-preserving map $f : \mathbb{R}^n \to \mathbb{R}^n$ has the form $f(\mathbf{x}) = Q\mathbf{x} + \mathbf{b}$ with $Q$ orthogonal. Linearity is not assumed: it is a consequence.

Proof. Let $g(\mathbf{x}) = f(\mathbf{x}) - f(\mathbf{0})$. Then $g$ preserves distances and $g(\mathbf{0}) = \mathbf{0}$, so it preserves norms: $\lVert g(\mathbf{x}) \rVert = \lVert g(\mathbf{x}) - g(\mathbf{0}) \rVert = \lVert \mathbf{x} \rVert$. From $\langle \mathbf{x}, \mathbf{y} \rangle = \tfrac{1}{2}\left(\lVert \mathbf{x} \rVert^2 + \lVert \mathbf{y} \rVert^2 - \lVert \mathbf{x} - \mathbf{y} \rVert^2\right)$, $g$ also preserves inner products. Hence $g(\mathbf{e}_1), \dots, g(\mathbf{e}_n)$ is an orthonormal basis. The coordinates of $g(\mathbf{x})$ in that basis are $\langle g(\mathbf{x}), g(\mathbf{e}_i) \rangle = \langle \mathbf{x}, \mathbf{e}_i \rangle = x_i$, so $g(\mathbf{x}) = \sum x_i,g(\mathbf{e}_i) = Q\mathbf{x}$, where $Q$ has columns $g(\mathbf{e}_i)$. $\blacksquare$

The hypothesis "distance-preserving" cannot be weakened to "norm-preserving". The map $\mathbf{x} \mapsto \lVert \mathbf{x} \rVert\mathbf{e}_1$ preserves every norm, but it sends $\mathbf{e}_1$ and $\mathbf{e}_2$ to the same point.

What isometries preserve.

Preserved	Not necessarily preserved
lengths, distances, inner products	orientation (reflections reverse it)
angles, orthogonality	coordinates (they change)
volumes: $\lvert\det Q\rvert = 1$	non-Euclidean norms such as $\ell^1$ and $\ell^\infty$
orthonormal bases map to orthonormal bases

The last row on the right is important. A rotation by $45°$ maps $(1, 0)$, with $\ell^1$ norm $1$, to $(\tfrac{1}{\sqrt{2}}, \tfrac{1}{\sqrt{2}})$, with $\ell^1$ norm $\sqrt{2}$. The linear isometries of the $\ell^1$ and $\ell^\infty$ norms are only the signed permutation matrices. Being an "isometry" always refers to a particular norm.

Why numerical analysts love isometries. An error vector $\mathbf{e}$ in the input of $Q\mathbf{x}$ becomes exactly $Q\mathbf{e}$ in the output, with the same length. An algorithm built from orthogonal transformations (Householder, Givens, the QR algorithm, SVD algorithms) therefore cannot amplify rounding errors from step to step. This is the core of their backward stability (Section 11.3).

Application: point cloud alignment (orthogonal Procrustes). Given points $\mathbf{p}_i$ and their measured positions $\mathbf{y}_i$ after an unknown rigid motion, find the rotation $R$ and translation $\mathbf{t}$ minimizing $\sum\lVert R\mathbf{p}_i + \mathbf{t} - \mathbf{y}_i \rVert^2$. Centering both clouds removes $\mathbf{t}$. Then, with $H = \sum \tilde{\mathbf{y}}_i\tilde{\mathbf{p}}_i^T = U\Sigma V^T$ (SVD, Chapter 11), the optimal rotation is $R = UDV^T$, where $D = \operatorname{diag}(1, \dots, 1, \det(UV^T))$ rules out reflections. This is the Kabsch algorithm, used in molecular biology, robotics (ICP scan matching), and motion capture.

import numpy as np

rng = np.random.default_rng(0)
Q, _ = np.linalg.qr(rng.standard_normal((3, 3)))      # a random orthogonal matrix
x = rng.standard_normal(3)

# An orthogonal matrix is an isometry of the Euclidean norm only
for p in [2, 1, np.inf]:
    print(p, np.isclose(np.linalg.norm(Q @ x, p), np.linalg.norm(x, p)))
# 2 True
# 1 False
# inf False

# Norm-preserving but NOT distance-preserving: nonlinear maps escape the theorem
f = lambda v: np.linalg.norm(v) * np.array([1., 0., 0.])        # sends everything onto the x-axis
a, b = np.array([1., 0., 0.]), np.array([0., 1., 0.])
print(np.linalg.norm(f(a)) == np.linalg.norm(a), np.linalg.norm(f(a) - f(b)), np.linalg.norm(a - b))
# True 0.0 1.4142135623730951

# Recover a rigid motion x -> R x + t from noisy point correspondences (Kabsch / orthogonal Procrustes)
theta = 0.4
R_true = np.array([[np.cos(theta), -np.sin(theta), 0], [np.sin(theta), np.cos(theta), 0], [0, 0, 1]])
t_true = np.array([1., -2., 0.5])
P = rng.standard_normal((50, 3))                                 # 50 points (rows)
Y = P @ R_true.T + t_true + 0.01 * rng.standard_normal((50, 3))  # moved and slightly noisy

Pc, Yc = P - P.mean(0), Y - Y.mean(0)                            # remove translations first
U, s, Vt = np.linalg.svd(Yc.T @ Pc)                              # SVD of the cross-covariance (Chapter 11)
D = np.diag([1, 1, np.sign(np.linalg.det(U @ Vt))])              # forbid reflections: keep det = +1
R = U @ D @ Vt
t = Y.mean(0) - R @ P.mean(0)
print(np.round(np.abs(R - R_true).max(), 3), np.round(t, 2))     # small error, t close to (1, -2, 0.5)
print(np.allclose(R.T @ R, np.eye(3)), round(np.linalg.det(R), 6))   # True 1.0

Key takeaways

• A linear map preserves norms if and only if it preserves inner products (polarization); in $\mathbb{R}^n$ this means its matrix has orthonormal columns.
• Every distance-preserving map of $\mathbb{R}^n$ is $\mathbf{x} \mapsto Q\mathbf{x} + \mathbf{b}$ with $Q$ orthogonal; linearity follows from the distance condition.
• Isometries preserve lengths, angles, and volumes, but not orientation, and only for the Euclidean norm they were defined with.
• Orthogonal transformations do not amplify errors, which is why stable algorithms are built from them.
• The best rotation aligning two point sets comes from an SVD (Kabsch / orthogonal Procrustes).

🧪 Practice

1. Show that the map $T : \mathbb{R}^2 \to \mathbb{R}^3$, $T(x, y) = \left(\tfrac{x + y}{\sqrt{2}}, \tfrac{x - y}{\sqrt{2}}, 0\right)$ is a linear isometry, and explain why it has no inverse.
2. Prove that the composition of two distance-preserving maps $\mathbf{x} \mapsto Q_1\mathbf{x} + \mathbf{b}_1$ and $\mathbf{x} \mapsto Q_2\mathbf{x} + \mathbf{b}_2$ is again of this form, and find its $Q$ and $\mathbf{b}$.
3. Find all $2 \times 2$ matrices that preserve the $\ell^\infty$ norm $\max(\lvert x \rvert, \lvert y \rvert)$. Hint: the unit square must map onto itself, so corners go to corners.
4. Interview: You must align two 3D scans of the same object, with known point correspondences, and the result must be a proper rotation, never a mirror image. Outline the algorithm and the one-line safeguard. Hint: center both clouds, take the SVD of the cross-covariance matrix, set $R = UV^T$, and flip the sign of the last singular direction when $\det(UV^T) = -1$.

Rotation Groups and Orthogonal Groups

Why it matters. Composing two rotations gives a rotation, every rotation can be undone, and doing nothing is a rotation. Rotations therefore form a group, and studying them as a group explains much of their behavior. It explains why 3D rotations depend on order while 2D ones do not, why every 3D rotation has an axis, why three numbers describe a 3D orientation, and why rotation matrices in a simulation slowly stop being rotations. Robotics, aerospace navigation, computer graphics, and physics all rely on this structure.

The groups. A group is a set with an associative operation, an identity element, and inverses. The orthogonal and unitary matrices form groups under multiplication, because products and inverses of such matrices stay in the set (previous topics):

Group	Definition	Geometric meaning	Dimension (free parameters)
$O(n)$, orthogonal group	$Q^TQ = I$	rotations and reflections of $\mathbb{R}^n$	$n(n-1)/2$
$SO(n)$, special orthogonal group	$Q^TQ = I$, $\det Q = 1$	rotations only	$n(n-1)/2$
$U(n)$, unitary group	$U^*U = I$	"rotations" of $\mathbb{C}^n$	$n^2$ (real parameters)
$SU(n)$, special unitary group	$U^*U = I$, $\det U = 1$	used for quantum gates and particle physics	$n^2 - 1$

Where does $n(n-1)/2$ come from? The equation $Q^TQ = I$ is symmetric, so it imposes only $n(n+1)/2$ independent conditions on the $n^2$ entries, leaving $n^2 - n(n+1)/2 = n(n-1)/2$ free parameters. $SO(2)$ needs one angle. $SO(3)$ needs three numbers: an axis (two numbers for a direction) plus an angle, or three Euler angles.

$O(n)$ splits into two pieces: $SO(n)$ ($\det = +1$) and the reflections and improper rotations ($\det = -1$). No continuous path of orthogonal matrices connects the two pieces, because the determinant cannot jump from $1$ to $-1$ without passing through values that are not $\pm 1$. A rotation can be performed gradually, but a mirror reflection cannot be reached gradually from the identity.

$SO(2)$ is commutative. Rotations of the plane add their angles: $R_\alpha R_\beta = R_{\alpha + \beta} = R_\beta R_\alpha$. $SO(2)$ behaves exactly like the unit circle of complex numbers $e^{i\theta}$ under multiplication.

$SO(3)$ is not commutative. Rotate a book $90°$ about a vertical axis and then $90°$ about a horizontal axis. Then do the same rotations in the opposite order. The book ends up in different orientations. The order of rotations matters in 3D, which is why robotics and aerospace conventions specify the order of Euler angles precisely.

Euler's rotation theorem. Every $R \in SO(3)$ has an axis: a unit vector $\mathbf{k}$ with $R\mathbf{k} = \mathbf{k}$. Every 3D rotation is therefore a rotation by some angle about some fixed line.

Proof. We show $\det(R - I) = 0$:

$$\det(R - I) = \det(R - RR^T) = \det R,\det(I - R^T) = \det(I - R)^T = \det(I - R) = (-1)^3\det(R - I),$$

using $\det R = 1$. So $\det(R - I) = -\det(R - I)$, which forces it to be $0$, and $1$ is an eigenvalue. $\blacksquare$

The proof uses the odd dimension: in $\mathbb{R}^4$, a rotation can turn two perpendicular planes independently and fix no nonzero vector.

Axis-angle and Rodrigues' formula. The rotation by angle $\theta$ about a unit axis $\mathbf{k}$ is

$$R = I + (\sin\theta)K + (1 - \cos\theta)K^2, \qquad K = \begin{bmatrix} 0 & -k_3 & k_2 \ k_3 & 0 & -k_1 \ -k_2 & k_1 & 0 \end{bmatrix},$$

where $K\mathbf{x} = \mathbf{k} \times \mathbf{x}$ is the cross-product matrix (skew-symmetric, Section 4.2). Conversely, the axis of $R$ is its eigenvector for eigenvalue $1$, and the angle comes from the trace: $\operatorname{tr} R = 1 + 2\cos\theta$, since the eigenvalues are $1$, $e^{i\theta}$, and $e^{-i\theta}$.

Rotations as exponentials. Rodrigues' formula is the matrix exponential $R = e^{\theta K}$ (Section 8.2), summed in closed form using $K^3 = -K$. In general, every rotation in $SO(n)$ is $e^{S}$ for some skew-symmetric $S$. If $S^T = -S$, then $(e^S)^Te^S = e^{-S}e^S = I$. Skew-symmetric matrices are therefore "infinitesimal rotations": angular velocity is a skew-symmetric matrix, and $\dot{R} = \Omega R$ describes a spinning body. In Lie theory, the skew-symmetric matrices form the Lie algebra $\mathfrak{so}(n)$ of the group $SO(n)$.

Reflections generate everything. By the Cartan-Dieudonné theorem, every element of $O(n)$ is a product of at most $n$ Householder reflections. A rotation in the plane is the product of two reflections across lines at angle $\theta/2$ to each other. This is why Householder QR (a product of reflections) can produce any orthogonal $Q$.

Representations of 3D orientation.

Representation	Numbers stored	Strengths	Weaknesses
Rotation matrix	9	composes by multiplication; applies directly to vectors	redundant; drifts away from $SO(3)$
Euler angles	3	intuitive (yaw, pitch, roll)	gimbal lock (a lost degree of freedom); order conventions differ
Axis-angle	3 or 4	geometric; good for interpolation	composition is awkward
Unit quaternion	4	compact; smooth interpolation (slerp); easy to renormalize	double cover: $q$ and $-q$ give the same rotation

Unit quaternions form the group $SU(2)$, which maps two-to-one onto $SO(3)$.

Drift and re-orthonormalization. Multiplying many rotation matrices in floating point slowly produces matrices that are no longer orthogonal. Objects then shear or scale slightly. The fix is to project back onto the group periodically. The closest rotation to $M$ in the Frobenius norm is $UV^T$, where $M = U\Sigma V^T$ is the SVD (the polar decomposition, Chapter 11), with the same determinant safeguard as in the Kabsch algorithm. Re-orthonormalizing the columns with QR is a cheaper alternative.

import numpy as np
from scipy.linalg import expm

def skew(k):
    """Matrix K with K x = k x x (cross product); skew-symmetric."""
    return np.array([[0, -k[2], k[1]], [k[2], 0, -k[0]], [-k[1], k[0], 0]])

def rodrigues(axis, theta):
    """Rotation by theta about a unit axis: R = I + sin(theta) K + (1 - cos(theta)) K^2."""
    k = np.asarray(axis, float) / np.linalg.norm(axis)
    K = skew(k)
    return np.eye(3) + np.sin(theta) * K + (1 - np.cos(theta)) * K @ K

R = rodrigues([1, 1, 1], 2 * np.pi / 3)          # 120 degrees about (1, 1, 1): cycles the axes
print(np.round(R, 12) + 0.0)
# [[0. 0. 1.]
#  [1. 0. 0.]
#  [0. 1. 0.]]
print(np.allclose(R.T @ R, np.eye(3)), round(np.linalg.det(R), 12))     # True 1.0: R is in SO(3)

# The rotation is the exponential of a skew-symmetric matrix (Section 8.2)
k = np.array([1., 1., 1.]) / np.sqrt(3)
print(np.allclose(expm(2 * np.pi / 3 * skew(k)), R))                    # True

# Recover axis and angle: the axis is the eigenvector for eigenvalue 1; trace R = 1 + 2 cos(theta)
vals, vecs = np.linalg.eig(R)
axis = np.real(vecs[:, np.argmin(np.abs(vals - 1))])
angle = np.degrees(np.arccos((np.trace(R) - 1) / 2))
print(np.round(axis / axis[0], 6), round(angle, 6))                     # [1. 1. 1.] 120.0

# SO(3) is not commutative; SO(2) is
Rx, Ry = rodrigues([1, 0, 0], np.pi / 2), rodrigues([0, 1, 0], np.pi / 2)
print(np.allclose(Rx @ Ry, Ry @ Rx))                                    # False
A2, B2 = expm(0.3 * np.array([[0, -1], [1, 0]])), expm(1.1 * np.array([[0, -1], [1, 0]]))
print(np.allclose(A2 @ B2, B2 @ A2))                                    # True

# Drift: composing many small rotations in float32 slowly leaves SO(3)
step = rodrigues([0.3, -0.5, 0.8], 0.01).astype(np.float32)
M = np.eye(3, dtype=np.float32)
for _ in range(100000):
    M = M @ step
print(f"{np.linalg.norm(M.T @ M - np.eye(3)):.1e}")                     # 2.3e-03: visible drift away from SO(3)

# Re-orthonormalize: the closest rotation in the Frobenius norm is U V^T from the SVD (Chapter 11)
U, _, Vt = np.linalg.svd(M.astype(np.float64))
M_fixed = U @ Vt
print(f"{np.linalg.norm(M_fixed.T @ M_fixed - np.eye(3)):.1e}", round(np.linalg.det(M_fixed), 6))   # 8.1e-16 1.0

Key takeaways

• $O(n)$, $SO(n)$, $U(n)$, and $SU(n)$ are groups under multiplication; $SO(n)$ has dimension $n(n-1)/2$, so 3D orientation needs three parameters.
• $O(n)$ has two components, $\det = \pm 1$; reflections cannot be reached continuously from the identity.
• $SO(2)$ is commutative and $SO(3)$ is not; every 3D rotation has an axis (Euler), and $\operatorname{tr} R = 1 + 2\cos\theta$.
• Rodrigues' formula $R = I + \sin\theta K + (1 - \cos\theta)K^2 = e^{\theta K}$ connects rotations with skew-symmetric matrices.
• Accumulated rotations drift out of $SO(3)$ in floating point; project back with the SVD ($UV^T$) or QR, or store orientations as quaternions.

🧪 Practice

1. Find the axis and the angle of the rotation $R = \begin{bmatrix} 0 & -1 & 0 \ 1 & 0 & 0 \ 0 & 0 & 1 \end{bmatrix}$, and verify $\operatorname{tr} R = 1 + 2\cos\theta$.
2. Show that the product of two reflections in $\mathbb{R}^2$, across lines at angles $\alpha$ and $\beta$, is the rotation by $2(\beta - \alpha)$.
3. Verify that $K^3 = -K$ for the cross-product matrix of a unit vector, and use it to sum the series $e^{\theta K}$ into Rodrigues' formula.
4. Interview: A game engine updates each object's orientation every frame by multiplying its rotation matrix by a small incremental rotation. After an hour, objects look slightly squashed. Explain why, and give two fixes with their trade-offs. Hint: floating-point drift takes $M$ out of $SO(3)$; renormalize periodically (Gram-Schmidt or QR on the columns is cheap, the SVD gives the closest rotation), or store orientations as unit quaternions and renormalize a 4-vector instead.

<a id="95-least-squares-approximation"></a>

9.5 Least Squares Approximation

This subchapter applies orthogonal projection to the most common computational problem in data analysis: finding the best approximate solution of an inconsistent system. It covers the problem and its geometry, the normal equations and the QR method, regression and curve fitting, the pseudoinverse for rank-deficient problems, and Fourier series as projection in a function space.

Least Squares Problem

Why it matters. Real measurements contain noise, so a model with a few parameters can almost never match many observations exactly. Fitting a line to 100 points gives 100 equations in 2 unknowns, and they are inconsistent (Section 3.1). Instead of giving up, we ask for the parameters that make the errors as small as possible. In 1801 Gauss used this idea to predict where the asteroid Ceres would reappear from a handful of observations. Legendre published it in 1805. Today it underlies linear regression, calibration, GPS positioning, computer vision, and the training of linear models.

The problem. Given an $m \times n$ matrix $A$ (usually $m > n$) and $\mathbf{b} \in \mathbb{R}^m$, a least squares solution of $A\mathbf{x} = \mathbf{b}$ is a vector $\hat{\mathbf{x}}$ that minimizes

$$\lVert \mathbf{b} - A\mathbf{x} \rVert^2 = \sum_{i=1}^m \left(b_i - (A\mathbf{x})_i\right)^2.$$

The vector $\mathbf{r} = \mathbf{b} - A\mathbf{x}$ is the residual, and $\lVert \mathbf{r} \rVert^2$ is the sum of squared errors. If the system is consistent, the least squares solutions are just the ordinary solutions, with residual $\mathbf{0}$.

The geometry: it is a projection problem. As $\mathbf{x}$ ranges over $\mathbb{R}^n$, the vector $A\mathbf{x}$ ranges over $\operatorname{Col} A$. Minimizing $\lVert \mathbf{b} - A\mathbf{x} \rVert$ therefore means finding the point of the subspace $\operatorname{Col} A$ closest to $\mathbf{b}$. By the Best Approximation Theorem (Section 9.2), that point is the orthogonal projection:

$$A\hat{\mathbf{x}} = \operatorname{proj}_{\operatorname{Col} A}\mathbf{b}, \qquad \mathbf{b} - A\hat{\mathbf{x}} \perp \operatorname{Col} A.$$

                       b  (the data: not in Col A, so A x = b has no solution)
                       *
                      /|
                     / |   r = b - A x_hat     (residual, perpendicular to Col A)
                    /  |
      _____________/___|__________________________
     /            /    * A x_hat = proj of b      /
    /            O                               /    Col A = all possible predictions A x
   /                    * A x  (any other choice: farther from b)
  /_____________________________________________/

  least squares = pick the prediction in Col A closest to the data b

Existence and uniqueness.

• A least squares solution always exists, because the projection always exists.
• It is unique if and only if the columns of $A$ are linearly independent. The projection $\mathbf{p} = A\hat{\mathbf{x}}$ is always unique, but if $\operatorname{Nul} A \ne {\mathbf{0}}$, every $\hat{\mathbf{x}} + \mathbf{z}$ with $\mathbf{z} \in \operatorname{Nul} A$ gives the same $\mathbf{p}$. The pseudoinverse topic below chooses one of them.

Characterization. $\hat{\mathbf{x}}$ is a least squares solution if and only if the residual is orthogonal to every column of $A$:

$$A^T(\mathbf{b} - A\hat{\mathbf{x}}) = \mathbf{0}.$$

This is the bridge to the normal equations in the next topic.

Worked example. Fit a line $y = C + Dt$ through the points $(0, 6)$, $(1, 0)$, $(2, 0)$. Each point gives one equation:

$$\begin{bmatrix} 1 & 0 \ 1 & 1 \ 1 & 2 \end{bmatrix}\begin{bmatrix} C \ D \end{bmatrix} = \begin{bmatrix} 6 \ 0 \ 0 \end{bmatrix}.$$

No line passes through all three points. The least squares line is $y = 5 - 3t$. Its predictions are $\mathbf{p} = (5, 2, -1)$ and its residual is $\mathbf{r} = (1, -2, 1)$. The residual is orthogonal to both columns, $(1, 1, 1)$ and $(0, 1, 2)$, and the minimal sum of squared errors is $\lVert \mathbf{r} \rVert^2 = 6$.

Why squares? Other error measures are possible, and they give different answers:

Criterion	Minimizes	Best line for the example	Properties
Least squares ($\ell^2$)	$\sum r_i^2$	$y = 5 - 3t$	linear equations; geometric (projection); optimal under Gaussian noise; sensitive to outliers
Least absolute deviations ($\ell^1$)	$\sum \lvert r_i \rvert$	$y = 6 - 3t$	robust to outliers; needs linear programming; the best line may not be unique
Minimax ($\ell^\infty$)	$\max \lvert r_i \rvert$	$y = 4.5 - 3t$	controls the worst case; linear programming; errors equioscillate ($-1.5, +1.5, -1.5$)

Least squares wins in practice for three reasons. Its solution comes from a linear system. Its geometry is the projection theory of this chapter. Under independent Gaussian noise of equal variance it gives the maximum-likelihood estimate, and by the Gauss-Markov theorem it is the best linear unbiased estimator even without the Gaussian assumption.

import numpy as np

# Fit a line y = C + D t through (0, 6), (1, 0), (2, 0): three equations, two unknowns
t = np.array([0., 1., 2.])
b = np.array([6., 0., 0.])
A = np.column_stack([np.ones_like(t), t])        # columns: 1, t

x_hat, res, rank, sv = np.linalg.lstsq(A, b, rcond=None)
p = A @ x_hat                                    # projection of b onto Col A
e = b - p                                        # residual
print(np.round(x_hat, 10) + 0.0, np.round(p, 10) + 0.0, np.round(e, 10) + 0.0)   # [ 5. -3.] [ 5.  2. -1.] [ 1. -2.  1.]
print(np.round(A.T @ e, 10) + 0.0, round(e @ e, 10))   # [0. 0.] 6.0: residual is orthogonal to both columns

# Any other line does worse
rng = np.random.default_rng(0)
trials = x_hat + rng.uniform(-1, 1, size=(10000, 2))
errors = np.sum((trials @ A.T - b)**2, axis=1)
print(errors.min() >= e @ e)                     # True

# The three classical error criteria give three different "best" lines
from scipy.optimize import linprog
def fit_l1(A, b):
    """min sum |A x - b|  as a linear program in (x, s) with -s <= A x - b <= s."""
    m, n = A.shape
    c = np.r_[np.zeros(n), np.ones(m)]
    A_ub = np.block([[A, -np.eye(m)], [-A, -np.eye(m)]])
    return linprog(c, A_ub=A_ub, b_ub=np.r_[b, -b], bounds=[(None, None)] * n + [(0, None)] * m).x[:n]
def fit_linf(A, b):
    """min max |A x - b|  as a linear program in (x, s) with -s <= A x - b <= s."""
    m, n = A.shape
    c = np.r_[np.zeros(n), 1.0]
    A_ub = np.block([[A, -np.ones((m, 1))], [-A, -np.ones((m, 1))]])
    return linprog(c, A_ub=A_ub, b_ub=np.r_[b, -b], bounds=[(None, None)] * n + [(0, None)]).x[:n]
print(np.round(fit_l1(A, b), 4) + 0.0, np.round(fit_linf(A, b), 4) + 0.0)   # [ 6. -3.] [ 4.5 -3. ]

Key takeaways

• A least squares solution minimizes $\lVert \mathbf{b} - A\mathbf{x} \rVert^2$; it always exists, and it is unique exactly when $A$ has independent columns.
• Geometrically, $A\hat{\mathbf{x}}$ is the orthogonal projection of $\mathbf{b}$ onto $\operatorname{Col} A$.
• $\hat{\mathbf{x}}$ is a least squares solution if and only if the residual is orthogonal to the columns: $A^T(\mathbf{b} - A\hat{\mathbf{x}}) = \mathbf{0}$.
• Squared error is chosen because it is linear to solve, geometric, and statistically optimal under Gaussian noise; $\ell^1$ and $\ell^\infty$ give different, also useful, answers.

🧪 Practice

1. Find the least squares line through $(0, 1)$, $(1, 1)$, $(2, 2)$, $(3, 2)$. Compute the residual and check that it is orthogonal to the columns of $A$.
2. Show that if $\mathbf{b}$ is orthogonal to $\operatorname{Col} A$, then $\hat{\mathbf{x}} = \mathbf{0}$ is a least squares solution. Interpret the result.
3. For a single unknown, $A = \mathbf{a}$ (one column), show that $\hat{x} = \mathbf{a}^T\mathbf{b}/\mathbf{a}^T\mathbf{a}$, and connect it with vector projection (Section 2.3). What is $\hat{x}$ when $\mathbf{a} = (1, \dots, 1)$?
4. Interview: A teammate fits a linear model to sensor data that contains a few wildly wrong readings, and the fitted line is clearly pulled toward them. Why does least squares behave this way, and what would you try instead? Hint: squaring makes one large residual cost more than many small ones; consider $\ell^1$ regression, Huber loss, RANSAC, or iteratively reweighted least squares.

Normal Equations

Why it matters. The orthogonality condition $A^T(\mathbf{b} - A\hat{\mathbf{x}}) = \mathbf{0}$ rearranges into a square system, the normal equations, which can be solved by any method from Chapter 3. They give closed-form formulas for regression coefficients and the projection matrix, and they are the starting point for theory and for very large problems. They also have a numerical weakness that every practitioner should know about.

Derivation 1 (geometry). The residual must be orthogonal to $\operatorname{Col} A$, that is, it must lie in $(\operatorname{Col} A)^\perp = \operatorname{Nul} A^T$ (Section 9.2):

$$A^T(\mathbf{b} - A\hat{\mathbf{x}}) = \mathbf{0} \quad\Longleftrightarrow\quad A^TA,\hat{\mathbf{x}} = A^T\mathbf{b}.$$

The word "normal" means "perpendicular": the equations say that the residual is normal to the column space.

Derivation 2 (calculus). Expand $f(\mathbf{x}) = \lVert A\mathbf{x} - \mathbf{b} \rVert^2 = \mathbf{x}^TA^TA\mathbf{x} - 2\mathbf{b}^TA\mathbf{x} + \mathbf{b}^T\mathbf{b}$. Its gradient is $\nabla f = 2A^TA\mathbf{x} - 2A^T\mathbf{b}$. Setting it to zero gives the same equations. The Hessian $2A^TA$ is positive semidefinite (Section 10.3), so every critical point is a global minimum.

Properties of $A^TA$.

Fact	Reason
$A^TA$ is $n \times n$ and symmetric	$(A^TA)^T = A^TA$; its entries are dot products of columns (the Gram matrix)
$\operatorname{Nul}(A^TA) = \operatorname{Nul} A$	if $A^TA\mathbf{x} = \mathbf{0}$, then $\mathbf{x}^TA^TA\mathbf{x} = \lVert A\mathbf{x} \rVert^2 = 0$, so $A\mathbf{x} = \mathbf{0}$
$A^TA$ is invertible $\iff$ the columns of $A$ are independent	from the previous row
The normal equations are always consistent	$A^T\mathbf{b} \in \operatorname{Col} A^T = \operatorname{Col}(A^TA)$
$A^TA$ is positive definite when the columns are independent	$\mathbf{x}^TA^TA\mathbf{x} = \lVert A\mathbf{x} \rVert^2 > 0$ for $\mathbf{x} \ne \mathbf{0}$

Formulas (independent columns).

$$\hat{\mathbf{x}} = (A^TA)^{-1}A^T\mathbf{b}, \qquad \mathbf{p} = A\hat{\mathbf{x}} = \underbrace{A(A^TA)^{-1}A^T}_{P},\mathbf{b}.$$

$P$ is the orthogonal projection onto $\operatorname{Col} A$ (Section 9.2): $P^2 = P$, $P^T = P$, and $\operatorname{tr} P = n$. In statistics it is the hat matrix, because it "puts the hat on $\mathbf{y}$": $\hat{\mathbf{y}} = P\mathbf{y}$. If $A = Q$ has orthonormal columns, $Q^TQ = I$ and the formulas reduce to $\hat{\mathbf{x}} = Q^T\mathbf{b}$ and $P = QQ^T$. The matrix $(A^TA)^{-1}A^T$ is a left inverse of $A$ (Section 7.3), and it is the pseudoinverse $A^+$ in this case.

Worked example. For the line through $(0, 6)$, $(1, 0)$, $(2, 0)$:

$$A^TA = \begin{bmatrix} 3 & 3 \ 3 & 5 \end{bmatrix}, \qquad A^T\mathbf{b} = \begin{bmatrix} 6 \ 0 \end{bmatrix}, \qquad \begin{cases} 3C + 3D = 6 \ 3C + 5D = 0 \end{cases} ;\Longrightarrow; D = -3, ; C = 5.$$

Weighted least squares. If measurement $i$ has variance $\sigma_i^2$, minimize the weighted error $\sum w_i r_i^2 = \mathbf{r}^TW\mathbf{r}$ with $w_i = 1/\sigma_i^2$. The same derivation, using the inner product $\langle \mathbf{u}, \mathbf{v} \rangle_W = \mathbf{u}^TW\mathbf{v}$ (Section 9.1), gives

$$A^TWA,\hat{\mathbf{x}} = A^TW\mathbf{b}.$$

This is still an orthogonal projection, but orthogonal in the $W$-inner product. With a full covariance matrix $\Sigma$ in place of the diagonal variances, $W = \Sigma^{-1}$; this is generalized least squares.

The numerical weakness: squaring the condition number. Forming $A^TA$ squares the condition number (Section 11.3):

$$\kappa(A^TA) = \kappa(A)^2.$$

In double precision ($\varepsilon \approx 10^{-16}$), solving the normal equations loses about $2\log_{10}\kappa(A)$ digits. A method that works with $A$ directly loses only about $\log_{10}\kappa(A)$ digits when the residual is small. If $\kappa(A) = 10^8$, the normal equations can lose every digit. Worse, information can be destroyed when $A^TA$ is formed, before any solving starts. For the Läuchli matrix with $\varepsilon = 10^{-8}$, the diagonal entries $1 + \varepsilon^2$ of $A^TA$ round to exactly $1$, and the computed $A^TA$ becomes exactly singular although $A$ has independent columns.

When the normal equations are fine. They are fine when $A$ is well conditioned, for example orthogonal-ish features or a small number of standardized predictors. They suit very tall problems ($m \gg n$), where forming $A^TA$ (about $mn^2$ flops, using symmetry) and a Cholesky solve ($n^3/3$) are cheaper than QR, and streaming data, where $A^TA$ and $A^T\mathbf{b}$ can be accumulated row by row in $O(n^2)$ memory.

import numpy as np

A = np.array([[1., 0.], [1., 1.], [1., 2.]])
b = np.array([6., 0., 0.])

AtA, Atb = A.T @ A, A.T @ b
print(AtA.tolist(), Atb.tolist())                # [[3.0, 3.0], [3.0, 5.0]] [6.0, 0.0]
x_hat = np.linalg.solve(AtA, Atb)
print(x_hat)                                     # [ 5. -3.]

# Projection ("hat") matrix onto Col A
P = A @ np.linalg.solve(AtA, A.T)                # A (A^T A)^{-1} A^T without forming the inverse
print(np.allclose(P @ P, P), np.allclose(P, P.T), round(np.trace(P), 10))   # True True 2.0

# Weighted least squares: trust the first measurement 10 times more
W = np.diag([10., 1., 1.])
x_w = np.linalg.solve(A.T @ W @ A, A.T @ W @ b)
print(np.round(x_w, 4), np.round(A @ x_w - b, 4))    # [ 5.8824 -3.5294] [-0.1176  2.3529 -1.1765]: much closer to (0, 6)

# Dependent columns: A^T A is singular exactly when Nul A != {0}
A_dep = np.column_stack([A[:, 0], A[:, 1], A[:, 0] + A[:, 1]])
print(np.linalg.matrix_rank(A_dep.T @ A_dep))    # 2 (not 3): infinitely many least squares solutions

# Squaring the condition number: kappa(A^T A) = kappa(A)^2
eps = 1e-8
L = np.array([[1., 1.], [eps, 0.], [0., eps]])   # Lauchli matrix, kappa(L) ~ 1.4e8
print(f"{np.linalg.cond(L):.1e}")                # 1.4e+08
print((L.T @ L).tolist())                        # [[1.0, 1.0], [1.0, 1.0]]: 1 + eps^2 rounds to 1, so A^T A is exactly singular
x_true = np.array([1., 1.])
bL = L @ x_true
try:
    print(np.linalg.solve(L.T @ L, L.T @ bL))
except np.linalg.LinAlgError as err:
    print("normal equations:", err)              # normal equations: Singular matrix
print(np.linalg.lstsq(L, bL, rcond=None)[0])     # QR/SVD-based solver: [1. 1.]

Key takeaways

• The normal equations $A^TA\hat{\mathbf{x}} = A^T\mathbf{b}$ say that the residual is orthogonal to $\operatorname{Col} A$; they also follow from setting the gradient to zero.
• They are always consistent, and they have a unique solution exactly when $A$ has independent columns ($\operatorname{Nul}(A^TA) = \operatorname{Nul} A$).
• $\hat{\mathbf{x}} = (A^TA)^{-1}A^T\mathbf{b}$ and $P = A(A^TA)^{-1}A^T$ (the hat matrix), for independent columns.
• Weighted least squares solves $A^TWA\hat{\mathbf{x}} = A^TW\mathbf{b}$, a projection in the $W$-inner product.
• $\kappa(A^TA) = \kappa(A)^2$: use the normal equations only when $A$ is well conditioned; otherwise use QR or the SVD.

🧪 Practice

1. Use the normal equations to fit $y = C + Dt + Et^2$ to the points $(-1, 2)$, $(0, 0)$, $(1, 1)$, $(2, 3)$.
2. Prove that the normal equations are always consistent, even when $A^TA$ is singular. Hint: $\operatorname{Col}(A^TA) = \operatorname{Col}(A^T)$; use the null space fact in the table and Section 6.6.
3. Show that $P = A(A^TA)^{-1}A^T$ satisfies $P^2 = P$ and $P^T = P$, and that $I - P$ projects onto $\operatorname{Nul} A^T$.
4. Interview: A data scientist computes regression coefficients with inv(X.T @ X) @ X.T @ y and gets wildly different coefficients after adding a feature that is almost a copy of an existing one. Explain what is happening, and give two fixes. Hint: collinearity makes $\kappa(X)$ large and $\kappa(X^TX) = \kappa(X)^2$ enormous; use lstsq (QR or SVD), drop or combine the duplicate feature, or add ridge regularization $(X^TX + \lambda I)$.

Least Squares via QR

Why it matters. The QR factorization solves least squares problems without forming $A^TA$, so it avoids squaring the condition number. It costs about twice as much as the normal equations for tall matrices, but it is the default for dense least squares problems in numerical libraries. It is also the cleanest view of what least squares does: rotate the problem into coordinates in which the solvable part and the unsolvable part separate.

Derivation. Let $A = QR$ be a reduced QR factorization, with $A$ having independent columns, so $R$ is invertible. Substitute into the normal equations:

$$A^TA\hat{\mathbf{x}} = A^T\mathbf{b} ;\Longrightarrow; R^TQ^TQR\hat{\mathbf{x}} = R^TQ^T\mathbf{b} ;\Longrightarrow; R^TR\hat{\mathbf{x}} = R^TQ^T\mathbf{b} ;\Longrightarrow; R\hat{\mathbf{x}} = Q^T\mathbf{b}.$$

The last step cancels the invertible $R^T$. Solve the triangular system by back substitution. In practice $Q$ is never formed explicitly: Householder reflections are applied to $\mathbf{b}$ to compute $Q^T\mathbf{b}$ directly.

The geometric view (complete QR). Write $A = [Q_1 ;; Q_2]\begin{bmatrix} R \ 0 \end{bmatrix}$, where $Q_1$ spans $\operatorname{Col} A$ and $Q_2$ spans $(\operatorname{Col} A)^\perp$. Multiplying by an orthogonal matrix does not change lengths (Section 9.4), so

$$\lVert A\mathbf{x} - \mathbf{b} \rVert^2 = \left\lVert \begin{bmatrix} R\mathbf{x} - Q_1^T\mathbf{b} \ -Q_2^T\mathbf{b} \end{bmatrix} \right\rVert^2 = \underbrace{\lVert R\mathbf{x} - Q_1^T\mathbf{b} \rVert^2}{\text{can be made } 0} + \underbrace{\lVert Q_2^T\mathbf{b} \rVert^2}{\text{cannot be changed by } \mathbf{x}}.$$

The first term vanishes when $R\hat{\mathbf{x}} = Q_1^T\mathbf{b}$, and the second term is the minimal squared residual. In rotated coordinates, the problem splits into a part that $\mathbf{x}$ controls completely and a part that no choice of $\mathbf{x}$ can touch.

   rotate by Q^T (lengths unchanged):

   Q^T A = [ R ]  n rows        Q^T b = [ Q1^T b ]  n entries   <- matched exactly by R x_hat
           [ 0 ]  m - n rows            [ Q2^T b ]  m - n       <- the residual: ||r||^2 = ||Q2^T b||^2

Accuracy. Householder QR solves least squares in a backward stable way. The computed $\hat{\mathbf{x}}$ is the exact solution of a nearby problem. Its error is governed by $\kappa(A)$ when the residual is small, and by roughly $\kappa(A)^2\tan\theta$ when the residual is large, where $\theta$ is the angle between $\mathbf{b}$ and $\operatorname{Col} A$. That sensitivity is inherent in the problem; no algorithm can avoid it. The normal equations add a $\kappa(A)^2$ error even when the residual is zero.

Choosing a method.

Method	Flops (dense $m \times n$)	Accuracy	Use when
Normal equations + Cholesky	$mn^2 + \tfrac{1}{3}n^3$	loses $2\log_{10}\kappa$ digits	well conditioned, very tall, or streaming
Householder QR	$2mn^2 - \tfrac{2}{3}n^3$	loses about $\log_{10}\kappa$ digits	the default for dense full-rank problems
QR with column pivoting	slightly more than QR	as QR, and detects rank deficiency	columns may be nearly dependent
SVD (Chapter 11)	several times QR	best; handles any rank	rank-deficient or ill-posed problems; minimum-norm solutions
Iterative (LSQR, CGLS)	per iteration: two products with $A$	depends on conditioning and preconditioning	huge sparse problems (Section 12.2)

np.linalg.lstsq uses an SVD-based LAPACK routine. scipy.linalg.lstsq can use QR- or SVD-based drivers.

import numpy as np
from scipy.linalg import solve_triangular

def lstsq_normal(A, b):
    return np.linalg.solve(A.T @ A, A.T @ b)

def lstsq_qr(A, b):
    Q, R = np.linalg.qr(A)                       # reduced QR (Householder)
    return solve_triangular(R, Q.T @ b)          # R x = Q^T b by back substitution

# The small example: same answer, and the residual norm comes from the complete QR
A = np.array([[1., 0.], [1., 1.], [1., 2.]])
b = np.array([6., 0., 0.])
print(lstsq_qr(A, b))                            # [ 5. -3.]
Qf, _ = np.linalg.qr(A, mode="complete")
print(round(np.linalg.norm(Qf[:, 2:].T @ b)**2, 10))    # 6.0 = ||b - A x_hat||^2, from the last m - n entries of Q^T b

# Ill-conditioned polynomial fit: degree-11 polynomial on [0, 1] with a known exact answer
t = np.linspace(0, 1, 50)
V = np.vander(t, 12, increasing=True)            # columns 1, t, ..., t^11
x_true = np.ones(12)
b = V @ x_true                                   # data generated exactly by the model
print(f"kappa(V) = {np.linalg.cond(V):.1e}")
for name, solver in [("normal equations", lstsq_normal), ("QR", lstsq_qr),
                     ("lstsq (SVD)", lambda A, b: np.linalg.lstsq(A, b, rcond=None)[0])]:
    x = solver(V, b)
    print(f"{name:17s} relative error {np.linalg.norm(x - x_true) / np.linalg.norm(x_true):.1e}")

The exact digits depend on the machine, but the pattern is reliable. The normal equations lose about twice as many digits as QR on the same problem, as $\kappa^2$ versus $\kappa$ predicts.

Key takeaways

• With $A = QR$, the least squares solution solves $R\hat{\mathbf{x}} = Q^T\mathbf{b}$ by back substitution, without forming $A^TA$.
• In rotated coordinates, $\lVert A\mathbf{x} - \mathbf{b} \rVert^2 = \lVert R\mathbf{x} - Q_1^T\mathbf{b} \rVert^2 + \lVert Q_2^T\mathbf{b} \rVert^2$; the second term is the minimal residual.
• QR loses about $\log_{10}\kappa$ digits while the normal equations lose about $2\log_{10}\kappa$; the QR method costs roughly twice as much for tall matrices.
• Use QR by default, the SVD for rank-deficient problems, and the normal equations only for well-conditioned or streaming problems.

🧪 Practice

1. Solve the least squares problem of the worked example (points $(0, 6)$, $(1, 0)$, $(2, 0)$) using a QR factorization computed by hand with Gram-Schmidt.
2. Show that $\lVert \mathbf{b} \rVert^2 = \lVert Q_1^T\mathbf{b} \rVert^2 + \lVert Q_2^T\mathbf{b} \rVert^2$ and interpret the two terms as the energy of $\mathbf{b}$ that the model explains and the energy it does not.
3. Repeat the polynomial experiment with degrees 5, 8, and 14, and plot the relative errors of the three methods against $\kappa(V)$ on a log-log scale. What slopes do you expect?
4. Interview: Which would you use to solve a least squares problem: the normal equations, QR, or the SVD? Give one situation where each is the right choice. Hint: think about conditioning, rank deficiency, problem size, whether data arrive in a stream, and whether a minimum-norm solution is needed.

Linear Regression and Curve Fitting

Why it matters. Linear regression is least squares with statistical vocabulary, and it is one of the most widely used tools in science, economics, and machine learning. "Linear" refers to the parameters, not to the curve. Polynomials, exponentials, sinusoids, circles, and many other shapes can be fitted with the same machinery once the problem is written as a design matrix times a coefficient vector. The geometry of this chapter also explains the standard formulas and the $R^2$ statistic.

The linear model. Observations $y_1, \dots, y_m$ are modeled as

$$y_i \approx \beta_0 + \beta_1\phi_1(x_i) + \cdots + \beta_k\phi_k(x_i) \qquad\Longleftrightarrow\qquad \mathbf{y} \approx X\boldsymbol{\beta},$$

where the design matrix $X$ has a column of ones (the intercept) and one column per feature $\phi_j(x_i)$. The features can be any fixed functions of the inputs. The model is linear in $\boldsymbol{\beta}$, so the least squares estimate is $\hat{\boldsymbol{\beta}} = \arg\min\lVert \mathbf{y} - X\boldsymbol{\beta} \rVert^2$, computed by QR or the normal equations.

Model	Columns of $X$	Linear in the parameters?
Line $y = \beta_0 + \beta_1 x$	$1, x$	yes
Polynomial of degree $d$	$1, x, \dots, x^d$	yes
Periodic $y = \beta_0 + \beta_1\cos\omega t + \beta_2\sin\omega t$ ($\omega$ known)	$1, \cos\omega t, \sin\omega t$	yes
Plane $z = \beta_0 + \beta_1 x + \beta_2 y$	$1, x, y$	yes
Circle $x^2 + y^2 + Dx + Ey + F = 0$	$x, y, 1$, with right-hand side $-(x^2 + y^2)$	yes
Exponential $y = ae^{bx}$	after taking logs: $1, x$ for $\ln y$	only after a transformation, which changes the error weighting
$y = a\sin(\omega t + \phi)$ with $\omega$ unknown	none	no: requires nonlinear least squares

Simple linear regression from the geometry. For $y \approx \beta_0 + \beta_1 x$, the columns $\mathbf{1}$ and $\mathbf{x}$ are generally not orthogonal. Replace $\mathbf{x}$ by its centered version $\mathbf{x} - \bar{x}\mathbf{1}$. This is one step of Gram-Schmidt, and it is orthogonal to $\mathbf{1}$. With orthogonal columns, each coefficient is a separate projection (Section 9.2):

$$\hat{\beta}_1 = \frac{\langle \mathbf{y}, \mathbf{x} - \bar{x}\mathbf{1} \rangle}{\lVert \mathbf{x} - \bar{x}\mathbf{1} \rVert^2} = \frac{\sum(x_i - \bar{x})(y_i - \bar{y})}{\sum(x_i - \bar{x})^2} = \frac{\operatorname{cov}(x, y)}{\operatorname{var}(x)}, \qquad \hat{\beta}_0 = \bar{y} - \hat{\beta}_1\bar{x}.$$

The numerator can use $y_i - \bar{y}$ instead of $y_i$, because $\mathbf{x} - \bar{x}\mathbf{1}$ is orthogonal to $\bar{y}\mathbf{1}$. The intercept formula says that the fitted line always passes through $(\bar{x}, \bar{y})$.

$R^2$ is a squared cosine. With an intercept in the model, the residual is orthogonal to $\mathbf{1}$ and to $\hat{\mathbf{y}}$, so Pythagoras applied to centered vectors gives

$$\underbrace{\lVert \mathbf{y} - \bar{y}\mathbf{1} \rVert^2}{\text{SST}} = \underbrace{\lVert \hat{\mathbf{y}} - \bar{y}\mathbf{1} \rVert^2}{\text{SSR (explained)}} + \underbrace{\lVert \mathbf{y} - \hat{\mathbf{y}} \rVert^2}_{\text{SSE (residual)}}, \qquad R^2 = 1 - \frac{\text{SSE}}{\text{SST}} = \cos^2\theta,$$

where $\theta$ is the angle between the centered data and the centered fit. $R^2 = 1$ means the data lie in the model's column space, and $R^2 = 0$ means the model explains nothing beyond the mean.

                   y - ybar*1   (centered data)
                   *
                  /|
                 / |    e = y - y_hat         SST = SSR + SSE   (Pythagoras)
                /  |                          R^2 = cos^2(theta)
               /   |
              / theta
             O-----*-------------------->  centered column space of X
                y_hat - ybar*1

Choosing the model's complexity. Adding columns to $X$ enlarges $\operatorname{Col} X$, so the training error can only decrease. A polynomial of degree $m - 1$ passes through all $m$ points exactly. A small training error does not mean good predictions. Too many parameters fit the noise (overfitting), and the curve swings wildly between data points. Standard remedies:

• Validation: choose the degree by the error on held-out data, not on the data used for fitting.
• Regularization: ridge regression minimizes $\lVert X\boldsymbol{\beta} - \mathbf{y} \rVert^2 + \lambda\lVert \boldsymbol{\beta} \rVert^2$, which solves $(X^TX + \lambda I)\hat{\boldsymbol{\beta}} = X^T\mathbf{y}$. This is always invertible for $\lambda > 0$ and shrinks unstable directions (Chapter 11).
• Orthogonal polynomial bases: monomials on $[0, 1]$ are nearly dependent (Section 9.1). Fitting with Legendre or Chebyshev polynomials, or after rescaling $x$ to $[-1, 1]$ as np.polynomial.Polynomial.fit does, keeps the design matrix well conditioned.

import numpy as np

rng = np.random.default_rng(42)

# 1) Simple linear regression: matrix solution vs the textbook formulas
x = np.array([1., 2., 3., 4., 5., 6.])
y = np.array([2.1, 3.9, 6.2, 7.8, 10.1, 12.2])
X = np.column_stack([np.ones_like(x), x])                 # design matrix: intercept column + feature
beta = np.linalg.lstsq(X, y, rcond=None)[0]
slope = np.sum((x - x.mean()) * (y - y.mean())) / np.sum((x - x.mean())**2)    # cov / var
intercept = y.mean() - slope * x.mean()
print(np.round(beta, 4), round(intercept, 4), round(slope, 4))   # [-0.02  2.02] -0.02 2.02

y_hat = X @ beta
SST = np.sum((y - y.mean())**2)
SSE = np.sum((y - y_hat)**2)
SSR = np.sum((y_hat - y.mean())**2)
print(np.isclose(SST, SSR + SSE), round(1 - SSE / SST, 5))      # True 0.99821: Pythagoras, R^2
yc, yhc = y - y.mean(), y_hat - y.mean()
print(round((yc @ yhc / np.linalg.norm(yc) / np.linalg.norm(yhc))**2, 5))   # 0.99821: R^2 = cos^2(angle)

# 2) "Linear" means linear in the PARAMETERS: fit a circle x^2 + y^2 + D x + E y + F = 0
theta = rng.uniform(0, 2 * np.pi, 40)
px = 3 + 2 * np.cos(theta) + 0.05 * rng.standard_normal(40)      # center (3, -1), radius 2, noisy
py = -1 + 2 * np.sin(theta) + 0.05 * rng.standard_normal(40)
M = np.column_stack([px, py, np.ones_like(px)])
D, E, F = np.linalg.lstsq(M, -(px**2 + py**2), rcond=None)[0]
center = np.array([-D / 2, -E / 2])
radius = np.sqrt(center @ center - F)
print(np.round(center, 2), round(radius, 2))                     # [ 3. -1.] 2.01

# 3) Polynomial degree: training error always falls, test error does not
f = lambda t: np.sin(2 * np.pi * t)
t_train = np.sort(rng.uniform(0, 1, 15)); y_train = f(t_train) + 0.15 * rng.standard_normal(15)
t_test = np.linspace(0, 1, 200);          y_test = f(t_test)
for deg in [1, 3, 5, 9, 14]:
    p = np.polynomial.Polynomial.fit(t_train, y_train, deg)       # well-conditioned: rescales t to [-1, 1]
    rms = lambda t, y: np.sqrt(np.mean((p(t) - y)**2))
    print(f"degree {deg:2d}: train RMS {rms(t_train, y_train):.3f}   test RMS {rms(t_test, y_test):.3f}")
# degree  1: train RMS 0.430   test RMS 0.478     underfits: a line cannot bend
# degree  3: train RMS 0.127   test RMS 0.117     about right
# degree  5: train RMS 0.107   test RMS 0.153
# degree  9: train RMS 0.075   test RMS 3.397     starts fitting the noise
# degree 14: train RMS 0.000   test RMS 10840.648 interpolates all 15 points, useless between them

Key takeaways

• Linear regression is least squares with a design matrix; "linear" means linear in the parameters, so polynomials, sinusoids, and circles can all be fitted.
• Centering $x$ (one Gram-Schmidt step) decouples the coefficients: slope $= \operatorname{cov}(x, y)/\operatorname{var}(x)$, intercept $= \bar{y} - \hat{\beta}_1\bar{x}$.
• With an intercept, $\text{SST} = \text{SSR} + \text{SSE}$ is Pythagoras, and $R^2$ is the squared cosine of the angle between the centered data and the centered fit.
• More columns always reduce training error; choose complexity with held-out data or regularization, and use well-conditioned bases.

🧪 Practice

1. Fit $y = \beta_0 + \beta_1 x$ to $(1, 1)$, $(2, 3)$, $(3, 2)$, $(4, 5)$ with the centered formulas, then check the result with np.linalg.lstsq. Compute $R^2$.
2. Fit $y = ae^{bx}$ to $(0, 2.0)$, $(1, 2.7)$, $(2, 3.7)$, $(3, 5.0)$ by regressing $\ln y$ on $x$. Explain why this is not the same as minimizing $\sum(y_i - ae^{bx_i})^2$.
3. Show that in simple linear regression the residuals sum to zero and are uncorrelated with $x$. Which columns of $X$ are responsible for each fact?
4. Interview: A model has training $R^2 = 0.99$ and test $R^2 = 0.2$. Diagnose the situation using the geometry of least squares, and list three concrete remedies. Hint: a large column space can approximate the noise in the training vector; consider fewer features, ridge or lasso regularization, cross-validation, or more data.

Pseudoinverse and Minimum Norm Solutions

Why it matters. When the columns of $A$ are dependent, the least squares problem has infinitely many solutions. This happens with duplicated features, underdetermined systems with more unknowns than equations, and models with more parameters than data points. A single, principled answer is still needed. The natural choice is the shortest solution, and the map that produces it is the Moore-Penrose pseudoinverse $A^+$. It is defined for every matrix, square or rectangular, invertible or not, and it reduces to $A^{-1}$ when $A$ is invertible.

The minimum-norm least squares solution. The set of least squares solutions is $\hat{\mathbf{x}} + \operatorname{Nul} A$: all of them have the same $A\mathbf{x} = \operatorname{proj}_{\operatorname{Col} A}\mathbf{b}$. By the orthogonal decomposition $\mathbb{R}^n = \operatorname{Row} A \oplus \operatorname{Nul} A$ (Section 6.6), exactly one of them, $\mathbf{x}^+$, lies in $\operatorname{Row} A$. Every other solution is $\mathbf{x}^+ + \mathbf{z}$ with $\mathbf{z} \in \operatorname{Nul} A$, $\mathbf{z} \ne \mathbf{0}$, and by Pythagoras

$$\lVert \mathbf{x}^+ + \mathbf{z} \rVert^2 = \lVert \mathbf{x}^+ \rVert^2 + \lVert \mathbf{z} \rVert^2 > \lVert \mathbf{x}^+ \rVert^2.$$

So $\mathbf{x}^+$ is the unique least squares solution of minimum norm. The pseudoinverse is the map $\mathbf{b} \mapsto \mathbf{x}^+$, and it is linear: $\mathbf{x}^+ = A^+\mathbf{b}$.

          x2                                   least squares solutions: a line parallel to Nul A
          ^       \                            (all give the same A x = proj of b)
          |        \
          |         *  x+  (the solution in Row A: closest to O)
          |        / \
          |   Row A   \
          |    /       \
          |   /         \  x+ + z, z in Nul A  (longer, by Pythagoras)
          |  /
          | /
          O--------------------------> x1

What $A^+$ does, in two steps. $A^+\mathbf{b}$ (1) projects $\mathbf{b}$ onto $\operatorname{Col} A$, then (2) maps the result back through the invertible correspondence between $\operatorname{Col} A$ and $\operatorname{Row} A$ (Section 6.6). Hence

$$AA^+ = \text{projection onto } \operatorname{Col} A, \qquad A^+A = \text{projection onto } \operatorname{Row} A.$$

The Penrose conditions. $A^+$ is the unique matrix $X$ satisfying

$$AXA = A, \qquad XAX = X, \qquad (AX)^T = AX, \qquad (XA)^T = XA.$$

These four conditions are a convenient test, and they are how the pseudoinverse is usually defined.

Formulas.

Case	$A^+$	Meaning
$A$ invertible	$A^{-1}$	ordinary inverse
Independent columns (full column rank)	$(A^TA)^{-1}A^T$	left inverse: the unique least squares solution
Independent rows (full row rank)	$A^T(AA^T)^{-1}$	right inverse: the shortest exact solution of an underdetermined system
Rank factorization $A = CF$ ($C$: $m \times r$ independent columns, $F$: $r \times n$ independent rows)	$F^T(FF^T)^{-1}(C^TC)^{-1}C^T$	general case, without the SVD
SVD $A = U\Sigma V^T$	$V\Sigma^+U^T$, inverting the nonzero singular values	the standard way to compute it (Section 11.2)
Limit of ridge regression	$\lim_{\lambda \to 0^+}(A^TA + \lambda I)^{-1}A^T$	regularization with vanishing strength

For example, the single equation $x_1 + 2x_2 + 2x_3 = 9$ has $A = \begin{bmatrix} 1 & 2 & 2 \end{bmatrix}$ with independent rows. Its shortest solution is $A^T(AA^T)^{-1} \cdot 9 = (1, 2, 2) \cdot \tfrac{9}{9} = (1, 2, 2)$. It points along the row of $A$, as every vector in $\operatorname{Row} A$ must.

Worked example (rank deficient). $A = \begin{bmatrix} 1 & 1 \ 1 & 1 \end{bmatrix}$ and $\mathbf{b} = (1, 3)$. $\operatorname{Col} A$ is spanned by $(1, 1)$, so $\operatorname{proj},\mathbf{b} = (2, 2)$, and the least squares solutions are all $\mathbf{x}$ with $x_1 + x_2 = 2$. The one in $\operatorname{Row} A = \operatorname{span}{(1, 1)}$ is $\mathbf{x}^+ = (1, 1)$. Indeed $A^+ = \tfrac{1}{4}\begin{bmatrix} 1 & 1 \ 1 & 1 \end{bmatrix}$ and $A^+\mathbf{b} = (1, 1)$.

A warning: $A^+$ is not continuous. $\operatorname{diag}(1, \delta)^+ = \operatorname{diag}(1, 1/\delta)$ for $\delta \ne 0$, but $\operatorname{diag}(1, 0)^+ = \operatorname{diag}(1, 0)$. As $\delta \to 0$ the pseudoinverse blows up and then suddenly jumps back. Tiny, noisy singular values produce huge, meaningless components. Numerical routines therefore treat singular values below a tolerance as exactly zero (rcond or rtol in NumPy). In inverse problems, truncating or damping small singular values is the essence of regularization (Section 11.2).

import numpy as np

# Rank-deficient least squares: A x = b with A = [[1, 1], [1, 1]] and b = (1, 3)
A = np.array([[1., 1.], [1., 1.]])
b = np.array([1., 3.])
A_pinv = np.linalg.pinv(A)
print(np.round(A_pinv, 12).tolist())             # [[0.25, 0.25], [0.25, 0.25]]
x_plus = A_pinv @ b
print(x_plus, A @ x_plus)                        # [1. 1.] [2. 2.]: best fit is x1 + x2 = 2; shortest such x is (1, 1)
x_other = x_plus + np.array([3., -3.])           # add a null space vector: still a least squares solution
print(np.allclose(A @ x_other, A @ x_plus), np.linalg.norm(x_plus) < np.linalg.norm(x_other))   # True True

# The four Penrose conditions characterize A+ uniquely
M = np.array([[1., 2., 3.], [2., 4., 6.], [1., 0., 1.]])     # rank 2
X = np.linalg.pinv(M)
checks = [np.allclose(M @ X @ M, M), np.allclose(X @ M @ X, X),
          np.allclose((M @ X).T, M @ X), np.allclose((X @ M).T, X @ M)]
print(checks)                                    # [True, True, True, True]

# Formula via a rank factorization M = C F (C: independent columns, F: independent rows)
C = M[:, [0, 1]]                                 # two independent columns of M
F = np.linalg.lstsq(C, M, rcond=None)[0]         # coefficients: M = C F exactly
print(np.allclose(C @ F, M))                     # True
X2 = F.T @ np.linalg.inv(F @ F.T) @ np.linalg.inv(C.T @ C) @ C.T
print(np.allclose(X2, X))                        # True

# Special cases
T = np.array([[1., 0.], [1., 1.], [1., 2.]])     # independent columns: A+ = (A^T A)^{-1} A^T (left inverse)
print(np.allclose(np.linalg.pinv(T), np.linalg.solve(T.T @ T, T.T)))       # True
W = np.array([[1., 2., 2.]])                     # independent rows: A+ = A^T (A A^T)^{-1} (right inverse)
print(np.linalg.pinv(W) @ [9.])                  # [1. 2. 2.]: shortest solution of x1 + 2 x2 + 2 x3 = 9

# The pseudoinverse is not continuous: tiny entries become huge, so pinv uses a cutoff
for d in [1e-3, 1e-10, 1e-17]:
    print(d, np.linalg.pinv(np.diag([1., d]))[1, 1])
# 0.001 1000.0
# 1e-10 10000000000.0
# 1e-17 0.0        below the default cutoff (rcond ~ 1e-15) the entry is treated as exactly zero

Key takeaways

• When $\operatorname{Nul} A \ne {\mathbf{0}}$, the least squares solutions form $\mathbf{x}^+ + \operatorname{Nul} A$; the unique one in $\operatorname{Row} A$ is the shortest.
• The pseudoinverse maps $\mathbf{b}$ to that minimum-norm least squares solution; $AA^+$ and $A^+A$ project onto $\operatorname{Col} A$ and $\operatorname{Row} A$.
• $A^+$ is characterized by the four Penrose conditions; it equals $(A^TA)^{-1}A^T$ for independent columns and $A^T(AA^T)^{-1}$ for independent rows.
• The pseudoinverse is discontinuous where the rank changes, so numerical routines discard tiny singular values with a tolerance.

🧪 Practice

1. Find the minimum-norm solution of $x_1 + x_2 + x_3 = 3$, $x_1 - x_3 = 0$, and verify that it lies in the row space.
2. Verify the four Penrose conditions for $A = \begin{bmatrix} 1 & 1 \ 1 & 1 \end{bmatrix}$ and $A^+ = \tfrac{1}{4}\begin{bmatrix} 1 & 1 \ 1 & 1 \end{bmatrix}$.
3. Prove that if $A$ has independent columns, then $(A^TA)^{-1}A^T$ satisfies all four Penrose conditions.
4. Interview: A linear model has more features than training examples. You train it with gradient descent on the squared loss, starting from $\mathbf{w} = \mathbf{0}$. Infinitely many weight vectors fit the data perfectly. Which one does gradient descent find, and why? Hint: each gradient $X^T(X\mathbf{w} - \mathbf{y})$ lies in $\operatorname{Row} X$, so the iterates never leave the row space; the only exact solution in $\operatorname{Row} X$ is the minimum-norm solution $X^+\mathbf{y}$. This is an example of implicit regularization.

Fourier Series as Orthogonal Projection

Why it matters. Writing a periodic signal as a sum of sines and cosines is the foundation of audio processing, image compression, the analysis of vibrations, and the solution of heat and wave equations. The Fourier coefficient formulas look like a list of integrals to memorize, but each one is the projection formula of Section 9.2 applied in a function space with an orthogonal basis. Least squares, Parseval, and Bessel all carry over unchanged.

The setting. Use $C[-\pi, \pi]$ with $\langle f, g \rangle = \int_{-\pi}^{\pi} f(x)g(x),dx$. The trigonometric system

$$1, ;\cos x, ;\sin x, ;\cos 2x, ;\sin 2x, ;\dots$$

is orthogonal (Section 9.1), with norms

$$\lVert 1 \rVert^2 = 2\pi, \qquad \lVert \cos kx \rVert^2 = \lVert \sin kx \rVert^2 = \pi \quad (k \ge 1).$$

Let $T_N = \operatorname{span}{1, \cos x, \sin x, \dots, \cos Nx, \sin Nx}$, the trigonometric polynomials of degree at most $N$, a subspace of dimension $2N + 1$.

Fourier coefficients are projection coefficients. The projection of $f$ onto $T_N$ with an orthogonal basis is $\sum \frac{\langle f, v \rangle}{\langle v, v \rangle}v$ over the basis vectors $v$. Spelling this out gives the classical formulas:

$$S_Nf(x) = a_0 + \sum_{k=1}^N\left(a_k\cos kx + b_k\sin kx\right), \qquad a_0 = \frac{1}{2\pi}\int_{-\pi}^{\pi} f, \quad a_k = \frac{1}{\pi}\int_{-\pi}^{\pi} f\cos kx, \quad b_k = \frac{1}{\pi}\int_{-\pi}^{\pi} f\sin kx.$$

The factors $\tfrac{1}{2\pi}$ and $\tfrac{1}{\pi}$ are just $1/\lVert v \rVert^2$. Many textbooks write the constant term as $a_0/2$ with $a_0 = \tfrac{1}{\pi}\int f$; that is the same projection with a different bookkeeping convention.

Consequences of the projection viewpoint.

1. Best approximation. $S_Nf$ is the trigonometric polynomial of degree $\le N$ closest to $f$ in mean square: $\lVert f - S_Nf \rVert \le \lVert f - g \rVert$ for every $g \in T_N$ (Best Approximation Theorem).
2. Adding terms never changes earlier coefficients. The basis is orthogonal, so enlarging $T_N$ to $T_{N+1}$ adds new terms and leaves the old coefficients alone. A non-orthogonal basis, such as monomials, does not have this property.
3. Error = energy not captured. Pythagoras gives $\lVert f - S_Nf \rVert^2 = \lVert f \rVert^2 - \lVert S_Nf \rVert^2$, which decreases as $N$ grows.
4. Bessel and Parseval. $2\pi a_0^2 + \pi\sum_{k=1}^N(a_k^2 + b_k^2) \le \int_{-\pi}^{\pi} f^2$ for every $N$. For square-integrable $f$ the error tends to $0$, and equality holds in the limit (Parseval): $\int_{-\pi}^{\pi} f^2 = 2\pi a_0^2 + \pi\sum_{k=1}^\infty(a_k^2 + b_k^2)$.

Worked example: the square wave. Let $f(x) = -1$ on $(-\pi, 0)$ and $f(x) = 1$ on $(0, \pi)$. It is odd, so it is orthogonal to the constant and to every cosine (Section 9.1): $a_0 = a_k = 0$. For the sines,

$$b_k = \frac{1}{\pi}\int_{-\pi}^{\pi} f(x)\sin kx,dx = \frac{2}{\pi}\int_0^\pi \sin kx,dx = \frac{2}{\pi k}\left(1 - \cos k\pi\right) = \begin{cases} \dfrac{4}{\pi k} & k \text{ odd}, \[4pt] 0 & k \text{ even}. \end{cases}$$

So $f(x) \sim \dfrac{4}{\pi}\left(\sin x + \dfrac{\sin 3x}{3} + \dfrac{\sin 5x}{5} + \cdots\right)$. Parseval with $\lVert f \rVert^2 = 2\pi$ gives $2\pi = \pi\sum_{k \text{ odd}}\frac{16}{\pi^2k^2}$, that is,

$$1 + \frac{1}{9} + \frac{1}{25} + \frac{1}{49} + \cdots = \frac{\pi^2}{8}.$$

A statement about a numerical series falls out of Pythagoras in a function space.

                                      |
 4/pi                                 |            .......
                                      |  oo    oo..   o   ..oo    oo
    1                                 |#o##o#oo.#ooooo#ooooo#.oo#o##o##
                                      |     o..               ..o
                                      |     .                   .
                                      |o   .                     .   o
                                      |  ..                       ..
                                      | .                           .
                                      |.                             .
    0 o-------------------------------o-------------------------------o
       .                             .|
        .                           . |
         ..                       ..  |
       o   .                     .   o|
            .                   .     |
            o..               ..o     |
   -1 ##o##o#oo.#ooooo#ooooo#.oo#o##o#|
         oo    oo..   o   ..oo    oo  |
                   .......            |
                                      |
      -pi                             0                              pi

   # = f (square wave)    . = S_1 = (4/pi) sin x    o = S_5 (odd harmonics up to sin 9x)
   S_5 hugs f more closely in mean square but overshoots near every jump (Gibbs phenomenon)

Mean-square convergence vs. pointwise convergence. Projection guarantees that $\lVert f - S_Nf \rVert \to 0$, which is convergence in the average squared sense. It does not guarantee that $S_Nf(x) \to f(x)$ at every point, and near a jump it does not converge uniformly. The partial sums overshoot the jump by about $9%$ of its size no matter how large $N$ is (the Gibbs phenomenon): the maximum of $S_Nf$ tends to about $1.179$, not to $1$. The overshoot gets narrower as $N$ grows, so it carries less and less energy. That is why the mean-square error still tends to $0$.

Complex form and the DFT. With the complex inner product $\langle f, g \rangle = \int_{-\pi}^{\pi} f\bar{g}$ (Section 9.1), the functions $e^{ikx}$, $k \in \mathbb{Z}$, are orthogonal with $\lVert e^{ikx} \rVert^2 = 2\pi$, and

$$f(x) \sim \sum_{k=-\infty}^{\infty} c_ke^{ikx}, \qquad c_k = \frac{\langle f, e^{ikx} \rangle}{\langle e^{ikx}, e^{ikx} \rangle} = \frac{1}{2\pi}\int_{-\pi}^{\pi} f(x)e^{-ikx},dx.$$

The conjugate in the inner product produces the minus sign in $e^{-ikx}$. Sampling $f$ at $n$ equally spaced points turns this into the discrete Fourier transform. The sampled exponentials are the orthogonal columns of the DFT matrix (Section 9.4), and the FFT computes all $n$ projections in $O(n\log n)$ operations (Section 13.3).

The general pattern. Every orthogonal expansion in a function space works the same way, including Legendre series on $[-1, 1]$, Chebyshev series, Bessel-function expansions, spherical harmonics, and wavelets: coefficient = $\langle f, v_k \rangle / \langle v_k, v_k \rangle$, truncation = best approximation, error = energy of the dropped coefficients. Section 13.4 places this in the setting of Hilbert spaces.

import numpy as np

# Midpoint grid on [-pi, pi] for numerical inner products <f, g> = integral f g dx
N = 200000
dx = 2 * np.pi / N
x = -np.pi + (np.arange(N) + 0.5) * dx           # symmetric about 0 and never exactly 0
ip = lambda f, g: np.sum(f * g) * dx

f = np.sign(x)                                   # square wave: -1 on (-pi, 0), +1 on (0, pi)

# Fourier coefficients are projection coefficients <f, v> / <v, v>
one = np.ones_like(x)
a0 = ip(f, one) / ip(one, one)
b = [ip(f, np.sin(k * x)) / ip(np.sin(k * x), np.sin(k * x)) for k in range(1, 8)]
print(round(a0, 6) + 0.0, np.round(b, 4) + 0.0)
# 0.0 [1.2732 0.     0.4244 0.     0.2546 0.     0.1819]
print(np.round([4 / (np.pi * k) for k in (1, 3, 5, 7)], 4))   # [1.2732 0.4244 0.2546 0.1819]: b_k = 4/(pi k) for odd k

def S(n_terms, x):
    """Partial sum with the first n_terms odd harmonics: the projection onto span{sin x, sin 3x, ...}."""
    return sum(4 / (np.pi * k) * np.sin(k * x) for k in range(1, 2 * n_terms, 2))

# Mean-square error = energy of f minus energy captured (Pythagoras / Parseval); overshoot does not shrink
norm_f2 = ip(f, f)                               # 2 pi
for n in [1, 5, 25, 125]:
    err2 = ip(f - S(n, x), f - S(n, x))
    captured = np.pi * sum((4 / (np.pi * k))**2 for k in range(1, 2 * n, 2))     # pi * sum b_k^2
    print(f"{n:4d} terms: ||f - S||^2 = {err2:.5f}   ||f||^2 - captured = {norm_f2 - captured:.5f}"
          f"   max S = {S(n, x).max():.4f}")
#    1 terms: ||f - S||^2 = 1.19023   ||f||^2 - captured = 1.19023   max S = 1.2732
#    5 terms: ||f - S||^2 = 0.25381   ||f||^2 - captured = 0.25381   max S = 1.1823
#   25 terms: ||f - S||^2 = 0.05092   ||f||^2 - captured = 0.05092   max S = 1.1791
#  125 terms: ||f - S||^2 = 0.01019   ||f||^2 - captured = 0.01019   max S = 1.1790

# Parseval for the full series: ||f||^2 = 2 pi = pi * sum b_k^2  =>  sum over odd k of 1/k^2 = pi^2 / 8
print(round(sum(1 / k**2 for k in range(1, 2000001, 2)), 5), round(np.pi**2 / 8, 5))   # 1.2337 1.2337

The table shows the three consequences at once. The mean-square error equals the energy not yet captured, exactly as Pythagoras predicts. It shrinks roughly like $1/N$. The overshoot settles near $1.179$ instead of disappearing.

Key takeaways

• On $[-\pi, \pi]$, ${1, \cos kx, \sin kx}$ is an orthogonal basis, and the Fourier coefficients are projection coefficients $\langle f, v \rangle / \langle v, v \rangle$.
• The partial sum $S_Nf$ is the orthogonal projection onto trigonometric polynomials of degree $\le N$, so it is the best mean-square approximation, and adding terms never changes earlier coefficients.
• The squared error equals the energy in the dropped coefficients; Parseval's identity, $\lVert f \rVert^2 = \sum(\text{coefficient})^2\lVert v_k \rVert^2$, is Pythagoras in infinitely many dimensions.
• Mean-square convergence does not imply pointwise convergence: the Gibbs overshoot of about $9%$ of the jump persists for every $N$.
• The complex form uses $e^{ikx}$ and the conjugate in the inner product; sampling it gives the DFT and the FFT.

🧪 Practice

1. Compute the Fourier series of $f(x) = x$ on $(-\pi, \pi)$, and use Parseval to show that $\sum_{k=1}^\infty 1/k^2 = \pi^2/6$.
2. Find the best approximation of $f(x) = \lvert x \rvert$ on $[-\pi, \pi]$ by a trigonometric polynomial of degree at most $1$, and compute the mean-square error.
3. Verify numerically, with a fine grid and np.sum(...) * dx, that $\int_{-\pi}^{\pi}\sin(3x)\sin(5x),dx = 0$ and $\int_{-\pi}^{\pi}\sin^2(3x),dx = \pi$. Then explain why the coefficient formula for $b_k$ divides by $\pi$.
4. Interview: An audio codec keeps only the largest frequency coefficients of each short block of samples. Use the ideas of this topic to explain why this minimizes the squared error for a given number of kept coefficients, and why codecs still need tricks to avoid audible artifacts at sharp attacks. Hint: by Parseval the squared error is the sum of the squares of the discarded coefficients, so keeping the largest is optimal; sharp transients spread energy across many coefficients and produce Gibbs-like ringing (pre-echo), which codecs handle with shorter blocks.

<a id="10-symmetric-matrices-and-quadratic-forms"></a>