Every orthogonality tool

Existence. Let \{\mathbf{u}_1, \dots, \mathbf{u}_p\} be an orthogonal basis for W (one always exists, by Gram–Schmidt if no other is in hand). Define

Each term of \hat{\mathbf y} is a scalar multiple of some \mathbf{u}_i \in W, and W is closed under linear combinations, so \hat{\mathbf y} \in W outright — no further check needed there. The work is showing \mathbf{z} \in W^\perp, and for that it suffices to show \mathbf{z} is orthogonal to every basis vector \mathbf{u}_i, since orthogonality to a spanning set of W forces orthogonality to all of W. Fix i and compute \mathbf{z}\cdot\mathbf{u}_i, using that \mathbf{u}_i\cdot\mathbf{u}_j = 0 for j \ne i (orthogonality of the basis) so every term of \hat{\mathbf y} except the i-th drops out of the dot product:

So \mathbf{z} \perp \mathbf{u}_i for every i, hence \mathbf{z} \in W^\perp, and \mathbf{y} = \hat{\mathbf y} + \mathbf{z} is a decomposition of the required form. That \hat{\mathbf y} is exactly \operatorname{proj}_W\mathbf{y} is just the definition of the projection onto W built from an orthogonal basis of W.

Uniqueness. Suppose \mathbf{y} also decomposes as \mathbf{y} = \hat{\mathbf y}' + \mathbf{z}' with \hat{\mathbf y}' \in W and \mathbf{z}' \in W^\perp. Subtracting the two decompositions of \mathbf{y} gives

The left side lies in W, since W is a subspace and both \hat{\mathbf y}, \hat{\mathbf y}' \in W; the right side lies in W^\perp for the same reason applied to \mathbf{z}, \mathbf{z}' \in W^\perp. So this common vector belongs to W \cap W^\perp. But a vector orthogonal to every vector in W — including itself, if it lies in W — must have zero length: taking \mathbf{v} = \hat{\mathbf y} - \hat{\mathbf y}' as that vector, \mathbf{v} \in W and \mathbf{v} \in W^\perp together force \mathbf{v}\cdot\mathbf{v} = 0, so \mathbf{v} = \mathbf{0}. Hence \hat{\mathbf y} = \hat{\mathbf y}' and, from the subtracted equation, \mathbf{z} = \mathbf{z}' as well. The decomposition is unique.

Sanity check: take W = \operatorname{span}\{(1,0,0),(0,1,0)\}, the xy-plane inside \mathbb{R}^3, and \mathbf{y} = (2,3,5). The orthogonal basis \{(1,0,0),(0,1,0)\} gives \hat{\mathbf y} = (2,3,0) and \mathbf{z} = (0,0,5), which is indeed orthogonal to both basis vectors — exactly the geometrically obvious split of \mathbf{y} into its in-plane part and its part along the perpendicular z-axis.

The claim has two parts: each \mathbf{u}_k built this way is orthogonal to every \mathbf{u}_i that came before it, and the running set \{\mathbf{u}_1, \ldots, \mathbf{u}_k\} spans the same subspace as \{\mathbf{v}_1, \ldots, \mathbf{v}_k\} at every stage. Both are proved together by induction on k.

Base case. \mathbf{u}_1 = \mathbf{v}_1, so trivially \operatorname{span}\{\mathbf{u}_1\} = \operatorname{span}\{\mathbf{v}_1\}, and there is no earlier vector to be orthogonal to.

Inductive step. Suppose \{\mathbf{u}_1, \ldots, \mathbf{u}_{k-1}\} is already orthogonal and spans the same set as \{\mathbf{v}_1, \ldots, \mathbf{v}_{k-1}\}. Form \mathbf{u}_k = \mathbf{v}_k - \sum_{i=1}^{k-1} \operatorname{proj}_{\mathbf{u}_i}(\mathbf{v}_k) = \mathbf{v}_k - \sum_{i=1}^{k-1} \frac{\mathbf{v}_k \cdot \mathbf{u}_i}{\mathbf{u}_i \cdot \mathbf{u}_i}\,\mathbf{u}_i.

Orthogonality. Fix any j with 1 \le j \le k-1 and dot both sides with \mathbf{u}_j: \mathbf{u}_k \cdot \mathbf{u}_j = \mathbf{v}_k \cdot \mathbf{u}_j - \sum_{i=1}^{k-1} \frac{\mathbf{v}_k \cdot \mathbf{u}_i}{\mathbf{u}_i \cdot \mathbf{u}_i}\,(\mathbf{u}_i \cdot \mathbf{u}_j). By the inductive hypothesis \{\mathbf{u}_1, \ldots, \mathbf{u}_{k-1}\} is orthogonal, so \mathbf{u}_i \cdot \mathbf{u}_j = 0 for every i \ne j in that sum, and every term except i = j drops out. What's left is \mathbf{u}_k \cdot \mathbf{u}_j = \mathbf{v}_k \cdot \mathbf{u}_j - \frac{\mathbf{v}_k \cdot \mathbf{u}_j}{\mathbf{u}_j \cdot \mathbf{u}_j}\,(\mathbf{u}_j \cdot \mathbf{u}_j) = \mathbf{v}_k \cdot \mathbf{u}_j - \mathbf{v}_k \cdot \mathbf{u}_j = 0. Since j was arbitrary, \mathbf{u}_k is orthogonal to every earlier \mathbf{u}_i, and since \{\mathbf{u}_1, \ldots, \mathbf{u}_{k-1}\} is already pairwise orthogonal, \{\mathbf{u}_1, \ldots, \mathbf{u}_k\} is too. (This step needs \mathbf{u}_j \ne \mathbf{0} for each earlier j, which the spanning argument below guarantees: a zero \mathbf{u}_j would mean \mathbf{v}_j already lay in the span of the vectors before it, contradicting linear independence of the original basis.)

Same span. The formula for \mathbf{u}_k writes \mathbf{v}_k as \mathbf{u}_k plus a linear combination of \mathbf{u}_1, \ldots, \mathbf{u}_{k-1}, i.e. \mathbf{v}_k = \mathbf{u}_k + \sum_{i=1}^{k-1} \frac{\mathbf{v}_k \cdot \mathbf{u}_i}{\mathbf{u}_i \cdot \mathbf{u}_i}\,\mathbf{u}_i. By the inductive hypothesis each \mathbf{u}_i on the right is itself a linear combination of \mathbf{v}_1, \ldots, \mathbf{v}_{k-1}, so \mathbf{v}_k \in \operatorname{span}\{\mathbf{v}_1, \ldots, \mathbf{v}_{k-1}, \mathbf{u}_k\}, and conversely \mathbf{u}_k \in \operatorname{span}\{\mathbf{v}_1, \ldots, \mathbf{v}_k\} directly from its defining formula. Combined with the inductive hypothesis that \operatorname{span}\{\mathbf{u}_1, \ldots, \mathbf{u}_{k-1}\} = \operatorname{span}\{\mathbf{v}_1, \ldots, \mathbf{v}_{k-1}\}, this gives \operatorname{span}\{\mathbf{u}_1, \ldots, \mathbf{u}_k\} = \operatorname{span}\{\mathbf{v}_1, \ldots, \mathbf{v}_k\}.

By induction, both properties hold for every k up to n, so \{\mathbf{u}_1, \ldots, \mathbf{u}_n\} is an orthogonal basis for the same subspace spanned by \{\mathbf{v}_1, \ldots, \mathbf{v}_n\}. None of the \mathbf{u}_k can be \mathbf{0}: that would force \mathbf{v}_k into the span of \mathbf{v}_1, \ldots, \mathbf{v}_{k-1}, contradicting linear independence, so every division by \mathbf{u}_i \cdot \mathbf{u}_i above is division by a positive number.

Normalizing. Once \{\mathbf{u}_1, \ldots, \mathbf{u}_n\} is orthogonal with no zero vector among them, setting \mathbf{q}_k = \mathbf{u}_k / \lVert\mathbf{u}_k\rVert scales each vector to unit length without disturbing the directions, so \{\mathbf{q}_1, \ldots, \mathbf{q}_n\} is orthonormal and spans the same subspace.

Run the Gram–Schmidt process on the columns \mathbf{a}_1, \dots, \mathbf{a}_n, producing orthonormal vectors \mathbf{q}_1, \dots, \mathbf{q}_n one at a time. At each step, build an intermediate vector \mathbf{w}_k by stripping off \mathbf{a}_k's component along every \mathbf{q}_i already found, \mathbf{w}_k = \mathbf{a}_k - \sum_{i=1}^{k-1} (\mathbf{a}_k \cdot \mathbf{q}_i)\,\mathbf{q}_i, \qquad \mathbf{q}_k = \frac{\mathbf{w}_k}{\lVert\mathbf{w}_k\rVert}. Linear independence of \mathbf{a}_1, \dots, \mathbf{a}_n guarantees \mathbf{w}_k \neq \mathbf{0} at every step (a zero \mathbf{w}_k would make \mathbf{a}_k a combination of \mathbf{a}_1, \dots, \mathbf{a}_{k-1}), so each division is legal and \lVert\mathbf{w}_k\rVert > 0.

Solve the defining relation for \mathbf{a}_k itself by moving the sum to the other side and substituting \mathbf{w}_k = \lVert\mathbf{w}_k\rVert\, \mathbf{q}_k: \mathbf{a}_k = \sum_{i=1}^{k-1} (\mathbf{a}_k \cdot \mathbf{q}_i)\, \mathbf{q}_i \;+\; \lVert\mathbf{w}_k\rVert\,\mathbf{q}_k. Every vector on the right is one of \mathbf{q}_1, \dots, \mathbf{q}_k, so this says each column of A is a linear combination of the \mathbf{q}_i's with an explicit, computable set of coefficients — exactly the data needed to assemble a matrix equation.

Let Q be the m \times n matrix with columns \mathbf{q}_1, \dots, \mathbf{q}_n, and define R entry by entry from the coefficients just found: r_{ik} = \mathbf{a}_k \cdot \mathbf{q}_i \ \ (i < k), \qquad r_{kk} = \lVert\mathbf{w}_k\rVert, \qquad r_{ik} = 0 \ \ (i > k). The zero pattern is exactly upper triangular, since \mathbf{a}_k was only ever built from \mathbf{q}_1, \dots, \mathbf{q}_k — it carries no component along \mathbf{q}_{k+1}, \dots, \mathbf{q}_n, which is what r_{ik} = 0 for i > k records. With this R, column k of the product QR is \sum_{i=1}^{n} r_{ik}\,\mathbf{q}_i = \sum_{i=1}^{k} r_{ik}\,\mathbf{q}_i, which is precisely the expression for \mathbf{a}_k derived above. Since this holds for every column, A = QR.

Finally, Q has orthonormal columns by construction of Gram–Schmidt, so Q^{\mathsf T}Q = I_n, and R is upper triangular with diagonal entries r_{kk} = \lVert\mathbf{w}_k\rVert > 0, as required.

Sanity check: take n = 1, so A is a single nonzero column \mathbf{a}_1. Then \mathbf{w}_1 = \mathbf{a}_1, \mathbf{q}_1 = \mathbf{a}_1 / \lVert\mathbf{a}_1\rVert, and R = [\lVert\mathbf{a}_1\rVert] — a single positive number. Indeed QR = \mathbf{q}_1 \cdot \lVert\mathbf{a}_1\rVert = \mathbf{a}_1, recovering A exactly, as the factorization must in the smallest nontrivial case.

Since the columns of A are linearly independent, A^{\mathsf T}A is invertible and the least-squares solution is the unique solution of the normal equations A^{\mathsf T}A\,\mathbf{x} = A^{\mathsf T}\mathbf{b}.

Substitute A = QR into both sides. On the left, A^{\mathsf T}A = (QR)^{\mathsf T}(QR) = R^{\mathsf T}Q^{\mathsf T}QR. The columns of Q are orthonormal, which is exactly the statement Q^{\mathsf T}Q = I_n, so this collapses to A^{\mathsf T}A = R^{\mathsf T}R. On the right, A^{\mathsf T}\mathbf{b} = (QR)^{\mathsf T}\mathbf{b} = R^{\mathsf T}Q^{\mathsf T}\mathbf{b}. The normal equations become R^{\mathsf T}R\,\mathbf{x} = R^{\mathsf T}Q^{\mathsf T}\mathbf{b}.

R is invertible — it is a factor of the full-column-rank matrix A in a factorization where Q has full column rank — so R^{\mathsf T} is also invertible, and left-multiplying both sides by \left(R^{\mathsf T}\right)^{-1} cancels it entirely: R\mathbf{x} = Q^{\mathsf T}\mathbf{b}. Since R is invertible this system has the unique solution \hat{\mathbf{x}} = R^{-1}Q^{\mathsf T}\mathbf{b}, which must therefore be the least-squares solution guaranteed by the normal equations.

Why this is the practical version of the normal equations. Forming A^{\mathsf T}A squares the conditioning of the problem — round-off error that was small in A gets amplified when A is multiplied against itself. The QR route never forms A^{\mathsf T}A at all: it computes Q^{\mathsf T}\mathbf{b} directly and then solves the triangular system R\hat{\mathbf x} = Q^{\mathsf T}\mathbf{b} by back substitution, which is both cheaper and numerically safer. The two formulas agree exactly — R^{\mathsf T}R above is nothing but a triangular factorization of A^{\mathsf T}A itself — they just reach the same \hat{\mathbf{x}} by different arithmetic paths.