Vector Projection and Gram-Schmidt Orthogonalization
Given vectors $\mathbf{a}$ and $\mathbf{b}$ in $\mathbb{R}^n$:
- Derive the formula for the projection of $\mathbf{b}$ onto $\mathbf{a}$, and compute the residual (the component of $\mathbf{b}$ orthogonal to $\mathbf{a}$). Verify that the residual is indeed orthogonal to $\mathbf{a}$.
- Explain how this projection extends to the Gram-Schmidt procedure: given a set of linearly independent vectors $\{\mathbf{v}_1, \ldots, \mathbf{v}_k\}$, produce an orthogonal set $\{\mathbf{u}_1, \ldots, \mathbf{u}_k\}$ spanning the same subspace.
- How does Gram-Schmidt relate to the QR decomposition?
Hints
- The projection of $\mathbf{b}$ onto $\mathbf{a}$ is the vector in the direction of $\mathbf{a}$ that minimizes $\|\mathbf{b} - c\mathbf{a}\|$. Set the derivative with respect to $c$ to zero.
- Gram-Schmidt works by subtracting from each new vector its projections onto all previously computed orthogonal vectors. The residual is the new orthogonal direction.
- The projection coefficients $\mathbf{e}_i^T \mathbf{v}_j$ form the entries of the upper-triangular matrix $R$ in the QR factorization.
Worked Solution
How to Think About It: Projection is the most fundamental operation in linear algebra -- it decomposes a vector into "the part along a direction" and "the part perpendicular to it." If you think of $\mathbf{a}$ as a signal and $\mathbf{b}$ as noisy data, the projection extracts the signal component and the residual is the noise. Gram-Schmidt just applies this idea repeatedly: project out each previously found direction, and what remains is the new orthogonal component. This is how you build an orthonormal basis from any linearly independent set.
Quick Estimate: The projection of $\mathbf{b}$ onto $\mathbf{a}$ should be a scalar multiple of $\mathbf{a}$. The scalar is the "shadow length" -- $\cos\theta \cdot \|\mathbf{b}\| / \|\mathbf{a}\|$ -- which works out to $\mathbf{a}^T \mathbf{b} / \mathbf{a}^T \mathbf{a}$.
Formal Solution:
Part 1 -- Projection and residual:
We want the vector $\hat{\mathbf{b}}$ in the span of $\mathbf{a}$ that is closest to $\mathbf{b}$. Write $\hat{\mathbf{b}} = c \, \mathbf{a}$ and minimize $\|\mathbf{b} - c \, \mathbf{a}\|^2$:
$$\frac{d}{dc} \|\mathbf{b} - c\mathbf{a}\|^2 = -2\mathbf{a}^T(\mathbf{b} - c\mathbf{a}) = 0$$
$$c = \frac{\mathbf{a}^T \mathbf{b}}{\mathbf{a}^T \mathbf{a}}$$
So the projection is:
$$\text{proj}_{\mathbf{a}}(\mathbf{b}) = \frac{\mathbf{a}^T \mathbf{b}}{\mathbf{a}^T \mathbf{a}} \, \mathbf{a}$$
The residual is:
$$\mathbf{b}_{\perp} = \mathbf{b} - \text{proj}_{\mathbf{a}}(\mathbf{b}) = \mathbf{b} - \frac{\mathbf{a}^T \mathbf{b}}{\mathbf{a}^T \mathbf{a}} \, \mathbf{a}$$
Verification: Check $\mathbf{a}^T \mathbf{b}_{\perp} = 0$:
$$\mathbf{a}^T \mathbf{b}_{\perp} = \mathbf{a}^T \mathbf{b} - \frac{\mathbf{a}^T \mathbf{b}}{\mathbf{a}^T \mathbf{a}} \, \mathbf{a}^T \mathbf{a} = \mathbf{a}^T \mathbf{b} - \mathbf{a}^T \mathbf{b} = 0 \quad \checkmark$$
Part 2 -- Gram-Schmidt orthogonalization:
Given linearly independent $\{\mathbf{v}_1, \ldots, \mathbf{v}_k\}$, produce orthogonal $\{\mathbf{u}_1, \ldots, \mathbf{u}_k\}$:
$$\mathbf{u}_1 = \mathbf{v}_1$$
For $j = 2, \ldots, k$:
$$\mathbf{u}_j = \mathbf{v}_j - \sum_{i=1}^{j-1} \frac{\mathbf{u}_i^T \mathbf{v}_j}{\mathbf{u}_i^T \mathbf{u}_i} \, \mathbf{u}_i$$
Each step subtracts from $\mathbf{v}_j$ its projections onto all previously computed orthogonal vectors. What remains is the component of $\mathbf{v}_j$ orthogonal to $\text{span}(\mathbf{u}_1, \ldots, \mathbf{u}_{j-1})$.
Why it works: By induction. After step $j$, $\mathbf{u}_j$ is orthogonal to $\mathbf{u}_1, \ldots, \mathbf{u}_{j-1}$ by construction (we subtracted all the projections). And $\text{span}(\mathbf{u}_1, \ldots, \mathbf{u}_j) = \text{span}(\mathbf{v}_1, \ldots, \mathbf{v}_j)$ because each $\mathbf{u}_j$ is a linear combination of $\mathbf{v}_1, \ldots, \mathbf{v}_j$ and vice versa.
To get an orthonormal set, normalize: $\mathbf{e}_j = \mathbf{u}_j / \|\mathbf{u}_j\|$.
Part 3 -- Connection to QR decomposition:
Gram-Schmidt on the columns of a matrix $V = [\mathbf{v}_1 | \cdots | \mathbf{v}_k]$ produces $Q = [\mathbf{e}_1 | \cdots | \mathbf{e}_k]$ with orthonormal columns, and the coefficients of the projection steps form an upper triangular matrix $R$:
$$V = QR$$
where $R_{ij} = \mathbf{e}_i^T \mathbf{v}_j$ for $i \leq j$, and $R_{ij} = 0$ for $i > j$. The upper-triangular structure follows because $\mathbf{v}_j$ is expressed as a combination of $\mathbf{e}_1, \ldots, \mathbf{e}_j$ only (never using later basis vectors).
In practice, the classical Gram-Schmidt algorithm can suffer from numerical instability. The modified Gram-Schmidt algorithm reorthogonalizes at each step, and Householder reflections provide an even more stable alternative for computing QR.
Answer: $\text{proj}_{\mathbf{a}}(\mathbf{b}) = \frac{\mathbf{a}^T \mathbf{b}}{\mathbf{a}^T \mathbf{a}} \mathbf{a}$, with residual $\mathbf{b} - \text{proj}_{\mathbf{a}}(\mathbf{b})$ orthogonal to $\mathbf{a}$. Gram-Schmidt applies this iteratively to build an orthogonal basis, and the process is equivalent to computing the QR decomposition $V = QR$.
Intuition
Projection and Gram-Schmidt are the workhorses behind half of applied linear algebra. OLS regression is projection of $Y$ onto the column space of $X$. PCA finds the directions of maximum variance by orthogonalizing the covariance eigenvectors. Signal processing decomposes signals into orthogonal frequency components. In every case, the logic is the same: decompose a vector into "the part we care about" (the projection) and "the rest" (the residual).
The connection to QR decomposition is practically important. QR is how most numerical software actually solves least-squares problems -- it is more stable than forming $X^T X$ and inverting it (the normal equations approach). When you call numpy.linalg.lstsq, it uses QR (or SVD) under the hood. Understanding that QR is just Gram-Schmidt in matrix form helps you see why it works and when it might fail (e.g., when columns of $X$ are nearly collinear, the $R$ matrix becomes nearly singular).