An inner product on a real vector space \(V\) is a map \(\langle \cdot, \cdot \rangle : V \times V \to \mathbb{R}\) satisfying, for all \(\mathbf{x}, \mathbf{y}, \mathbf{z} \in V\) and scalars \(a, b\),
On \(\mathbb{R}^{n}\) the standard choice is \(\langle \mathbf{x}, \mathbf{y}\rangle = \mathbf{x}'\mathbf{y} = \sum_i x_i y_i\), and then \(\|\mathbf{x}\| = \sqrt{\mathbf{x}'\mathbf{x}}\) is the ordinary length. Two vectors are orthogonal when \(\mathbf{x}'\mathbf{y} = 0\).
Why statistics cares. With the inner product \(\langle X, Y \rangle = E(XY)\) on random variables, orthogonality is zero correlation for centred variables, the length \(\|X\|\) is the standard deviation, and the angle between two variables has cosine equal to their correlation coefficient. Least squares, developed in Unit 4, is then a statement about orthogonal projection and nothing more.
Given linearly independent \(\mathbf{v}_1, \ldots, \mathbf{v}_k\), set
\[ \begin{aligned} \mathbf{u}_1 &= \mathbf{v}_1, \\ \mathbf{u}_2 &= \mathbf{v}_2 - \frac{\langle \mathbf{v}_2, \mathbf{u}_1\rangle}{\langle \mathbf{u}_1, \mathbf{u}_1\rangle}\,\mathbf{u}_1, \\ \mathbf{u}_r &= \mathbf{v}_r - \sum_{j=1}^{r-1} \frac{\langle \mathbf{v}_r, \mathbf{u}_j\rangle}{\langle \mathbf{u}_j, \mathbf{u}_j\rangle}\,\mathbf{u}_j, \end{aligned} \]and finally normalise, \(\mathbf{e}_r = \mathbf{u}_r / \|\mathbf{u}_r\|\). The \(\mathbf{u}_r\) are orthogonal and the \(\mathbf{e}_r\) orthonormal, and at every stage \(\operatorname{span}\{\mathbf{u}_1, \ldots, \mathbf{u}_r\} = \operatorname{span}\{\mathbf{v}_1, \ldots, \mathbf{v}_r\}\).
What each step does. The subtracted term is the orthogonal projection of \(\mathbf{v}_r\) onto the space already built. Removing it leaves exactly the part of \(\mathbf{v}_r\) that the earlier vectors could not explain — which is the same operation as taking a residual in regression.
Given. \(\mathbf{v}_1 = (1,1,0)'\), \(\mathbf{v}_2 = (1,0,1)'\), \(\mathbf{v}_3 = (0,1,1)'\).
Asked. Construct an orthonormal basis by Gram–Schmidt.
Step 1 — take \(\mathbf{u}_1 = \mathbf{v}_1 = (1,1,0)'\). Its squared length is \(\mathbf{u}_1'\mathbf{u}_1 = 1 + 1 + 0 = 2\).
Step 2 — the coefficient for \(\mathbf{v}_2\).
\[ \frac{\mathbf{v}_2'\mathbf{u}_1}{\mathbf{u}_1'\mathbf{u}_1} = \frac{(1)(1) + (0)(1) + (1)(0)}{2} = \frac{1}{2}. \]Step 3 — subtract.
\[ \mathbf{u}_2 = (1,0,1)' - \tfrac12 (1,1,0)' = \left(\tfrac12, -\tfrac12, 1\right)', \qquad \mathbf{u}_2'\mathbf{u}_2 = \tfrac14 + \tfrac14 + 1 = \tfrac32. \]Check orthogonality: \(\mathbf{u}_1'\mathbf{u}_2 = (1)\tfrac12 + (1)(-\tfrac12) + 0 = 0\). \(\checkmark\)
Step 4 — the two coefficients for \(\mathbf{v}_3\).
\[ \frac{\mathbf{v}_3'\mathbf{u}_1}{\mathbf{u}_1'\mathbf{u}_1} = \frac{0 + 1 + 0}{2} = \frac{1}{2}, \qquad \frac{\mathbf{v}_3'\mathbf{u}_2}{\mathbf{u}_2'\mathbf{u}_2} = \frac{0\left(\tfrac12\right) + 1\left(-\tfrac12\right) + 1(1)}{\tfrac32} = \frac{\tfrac12}{\tfrac32} = \frac{1}{3}. \]Step 5 — subtract both.
\[ \mathbf{u}_3 = (0,1,1)' - \tfrac12(1,1,0)' - \tfrac13\left(\tfrac12,-\tfrac12,1\right)' \]Component by component:
\[ u_{31} = 0 - \tfrac12 - \tfrac16 = -\tfrac23, \quad u_{32} = 1 - \tfrac12 + \tfrac16 = \tfrac23, \quad u_{33} = 1 - 0 - \tfrac13 = \tfrac23, \] \[ \mathbf{u}_3 = \left(-\tfrac23, \tfrac23, \tfrac23\right)', \qquad \mathbf{u}_3'\mathbf{u}_3 = 3 \times \tfrac49 = \tfrac43. \]Checks: \(\mathbf{u}_1'\mathbf{u}_3 = -\tfrac23 + \tfrac23 + 0 = 0\) and \(\mathbf{u}_2'\mathbf{u}_3 = -\tfrac13 - \tfrac13 + \tfrac23 = 0\). \(\checkmark\)
Step 6 — normalise. Dividing each by its length,
\[ \mathbf{e}_1 = \frac{1}{\sqrt{2}}(1,1,0)', \qquad \mathbf{e}_2 = \frac{1}{\sqrt{6}}(1,-1,2)', \qquad \mathbf{e}_3 = \frac{1}{\sqrt{3}}(-1,1,1)', \]where \(\mathbf{e}_2\) used \(\|\mathbf{u}_2\| = \sqrt{3/2} = 1.224745\) and the factor \(\tfrac12\) was cleared, and \(\mathbf{e}_3\) used \(\|\mathbf{u}_3\| = \sqrt{4/3} = 1.154701\).
Interpretation. The three starting vectors were independent but at \(60^{\circ}\) to one another; the output is a genuine set of coordinate axes for the same space. In Unit 3 of Distribution Theory exactly this construction — there in the form of the Helmert matrix — is what makes the sample mean and the sample variance independent.
Given. \(\mathbf{v}_1 = (3,1)'\) and \(\mathbf{v}_2 = (1,2)'\), the vectors drawn in Fig 1.1.
Step 1 — the coefficient.
\[ c = \frac{\mathbf{v}_2'\mathbf{v}_1}{\mathbf{v}_1'\mathbf{v}_1} = \frac{(1)(3) + (2)(1)}{3^{2} + 1^{2}} = \frac{5}{10} = 0.5. \]Step 2 — the projection. \(c\,\mathbf{v}_1 = 0.5\,(3,1)' = (1.5,\, 0.5)'\).
Step 3 — the residual.
\[ \mathbf{u}_2 = (1,2)' - (1.5, 0.5)' = (-0.5,\ 1.5)'. \]Step 4 — verify orthogonality. \(\mathbf{v}_1'\mathbf{u}_2 = 3(-0.5) + 1(1.5) = -1.5 + 1.5 = 0\). \(\checkmark\)
Let \(S = \mathcal{C}(A)\) be the column space of an \(n \times p\) matrix \(A\) of full column rank. The orthogonal projection of \(\mathbf{y}\) onto \(S\) is \(P\mathbf{y}\) with
\[ P = A\left(A'A\right)^{-1}A'. \]Its three defining properties, each checkable directly:
\[ P' = P \ \text{(symmetric)}, \qquad P^{2} = P \ \text{(idempotent)}, \qquad \operatorname{trace}(P) = \operatorname{rank}(P) = p. \]The residual \(\mathbf{y} - P\mathbf{y} = (I - P)\mathbf{y}\) is orthogonal to every column of \(A\), because \(A'(I-P) = A' - A'A(A'A)^{-1}A' = A' - A' = 0\). That single line is the normal equations, the Pythagorean decomposition of the sum of squares, and the geometry of least squares all at once.
In the one-way layout of Unit 4 the design matrix has more columns than its rank, so \(X'X\) is singular and \((X'X)^{-1}\) does not exist. The normal equations \(X'X\hat{\boldsymbol\beta} = X'\mathbf{y}\) are still consistent and still have solutions — infinitely many of them. A generalized inverse is the tool that produces one.
For any \(m \times n\) matrix \(A\), a matrix \(G\) is
\[ \begin{aligned} \text{a \textbf{generalized inverse}} \ A^{-}: \quad & \text{(i)}\ AGA = A \\ \text{a \textbf{reflexive} g-inverse}: \quad & \text{(i) and (ii)}\ GAG = G \\ \text{the \textbf{Moore–Penrose} inverse} \ A^{+}: \quad & \text{(i), (ii), (iii)}\ (AG)' = AG, \ \text{(iv)}\ (GA)' = GA \end{aligned} \]A generalized inverse always exists and is not unique. The Moore–Penrose inverse always exists and is unique — the extra two symmetry conditions are exactly what pins it down. When \(A\) is square and non-singular, all of these collapse to \(A^{-1}\).
Given.
\[ A = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ 1 & 1 \end{pmatrix}, \qquad 3 \times 2. \]Asked. Find \(A^{+}\) and verify the four Penrose conditions.
Step 1 — check the rank. The two columns are not multiples of one another, so \(\operatorname{rank}(A) = 2\), which is full column rank. Then \(A^{+} = (A'A)^{-1}A'\).
Step 2 — form \(A'A\).
\[ A'A = \begin{pmatrix} 1 & 0 & 1 \\ 0 & 1 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ 1 & 1 \end{pmatrix} = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}. \]Step 3 — invert it. \(\det(A'A) = 4 - 1 = 3\), so
\[ \left(A'A\right)^{-1} = \frac{1}{3}\begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix}. \]Step 4 — multiply by \(A'\).
\[ A^{+} = \frac{1}{3}\begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix} \begin{pmatrix} 1 & 0 & 1 \\ 0 & 1 & 1 \end{pmatrix} = \frac{1}{3}\begin{pmatrix} 2 & -1 & 1 \\ -1 & 2 & 1 \end{pmatrix}. \]Step 5 — condition (ii)'s ingredient, \(A^{+}A\).
\[ A^{+}A = \frac{1}{3}\begin{pmatrix} 2 & -1 & 1 \\ -1 & 2 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ 1 & 1 \end{pmatrix} = \frac{1}{3}\begin{pmatrix} 3 & 0 \\ 0 & 3 \end{pmatrix} = I_2. \]Step 6 — the other product, \(AA^{+}\).
\[ AA^{+} = \frac{1}{3}\begin{pmatrix} 2 & -1 & 1 \\ -1 & 2 & 1 \\ 1 & 1 & 2 \end{pmatrix}. \]This is not \(I_3\), and could not be: its rank is \(\operatorname{rank}(A) = 2 < 3\). Its trace is \(\tfrac{2 + 2 + 2}{3} = 2\), which equals the rank as an idempotent matrix's must.
Step 7 — the four conditions.
| condition | check | holds? |
|---|---|---|
| (i) \(AA^{+}A = A\) | \(A(A^{+}A) = A I_2 = A\) | yes |
| (ii) \(A^{+}AA^{+} = A^{+}\) | \(I_2 A^{+} = A^{+}\) | yes |
| (iii) \((AA^{+})' = AA^{+}\) | the matrix in Step 6 is symmetric | yes |
| (iv) \((A^{+}A)' = A^{+}A\) | \(I_2\) is symmetric | yes |
Interpretation. \(AA^{+}\) is the projection matrix onto the column space of \(A\) — compare section 3, where \(P = A(A'A)^{-1}A'\) is precisely this product. So the Moore–Penrose inverse is not an abstract construction: applied to \(\mathbf{y}\) it returns the least squares coefficients, and \(AA^{+}\mathbf{y}\) returns the fitted values.
The system \(A\mathbf{x} = \mathbf{b}\) with \(A\) of order \(m \times n\) is
\[ \begin{aligned} \textbf{consistent} \iff & \ \operatorname{rank}(A) = \operatorname{rank}\left([A \mid \mathbf{b}]\right) \iff AA^{-}\mathbf{b} = \mathbf{b} \\ \textbf{unique} \iff & \ \text{consistent and } \operatorname{rank}(A) = n \\ \textbf{many solutions} \iff & \ \text{consistent and } \operatorname{rank}(A) = r < n, \ \text{with } n - r \text{ free parameters} \end{aligned} \]and when it is consistent the general solution is
\[ \mathbf{x} = A^{-}\mathbf{b} + \left(I - A^{-}A\right)\mathbf{z}, \qquad \mathbf{z} \in \mathbb{R}^{n} \text{ arbitrary}. \]The second term runs over the null space of \(A\) as \(\mathbf{z}\) varies, so it generates every solution and nothing else.
Homogeneous systems. \(A\mathbf{x} = \mathbf{0}\) is always consistent, since \(\mathbf{x} = \mathbf{0}\) solves it. It has a non-trivial solution if and only if \(\operatorname{rank}(A) < n\), and the solution space then has dimension \(n - \operatorname{rank}(A)\).
Given.
\[ 2x_1 + x_2 + x_3 = 7, \qquad x_1 + 2x_2 + x_3 = 8, \qquad x_1 + x_2 + 2x_3 = 9. \]Step 1 — the coefficient matrix and its determinant.
\[ A = \begin{pmatrix} 2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{pmatrix}, \qquad \det(A) = 4 \ne 0, \]so \(\operatorname{rank}(A) = 3 = n\) and the solution is unique.
Step 2 — the inverse. Computed in Example 2.3 of the next unit by two different routes,
\[ A^{-1} = \frac{1}{4}\begin{pmatrix} 3 & -1 & -1 \\ -1 & 3 & -1 \\ -1 & -1 & 3 \end{pmatrix}. \]Step 3 — solve.
\[ \mathbf{x} = A^{-1}\mathbf{b} = \frac{1}{4} \begin{pmatrix} 3 & -1 & -1 \\ -1 & 3 & -1 \\ -1 & -1 & 3 \end{pmatrix} \begin{pmatrix} 7 \\ 8 \\ 9 \end{pmatrix}. \]Row by row:
\[ x_1 = \tfrac{1}{4}\left(21 - 8 - 9\right) = \tfrac{4}{4} = 1, \quad x_2 = \tfrac{1}{4}\left(-7 + 24 - 9\right) = \tfrac{8}{4} = 2, \quad x_3 = \tfrac{1}{4}\left(-7 - 8 + 27\right) = \tfrac{12}{4} = 3. \]Step 4 — check by substitution, not by trusting the algebra.
\[ 2(1) + 2 + 3 = 7 \checkmark, \qquad 1 + 2(2) + 3 = 8 \checkmark, \qquad 1 + 2 + 2(3) = 9 \checkmark. \]Interpretation. The solution is exactly \((1, 2, 3)'\). In Unit 4 the same machinery is applied to the normal equations, where the matrix is \(X'X\) and the right-hand side is \(X'\mathbf{y}\) — and where the determinant is often zero, which is why the generalized inverse of section 4 is needed.