Exercises — Normal equation closed-form solution

• $\mathbf{X}$ is $m\times (n+1)=100\times 5$. The $+1$ is the leading column of $1$'s that pairs with the intercept ${\theta }_0$.
• $\mathbf{\theta }$ is $(n+1)\times 1=5\times 1$ — one entry for the intercept and one per feature.
• $\mathbf{y}$ is $m\times 1=100\times 1$.
• ${\mathbf{X}}^T\mathbf{X}$: ${\mathbf{X}}^T$ is $5\times 100$, times $\mathbf{X}$ ($100\times 5$), giving $5\times 5$.

Picture it: ${\mathbf{X}}^T\mathbf{X}$ is always a small square of side $(n+1)$ no matter how many rows of data you have — that is why the normal equation is cheap when features are few even if data is huge.

$\mathbf{X}=\left[\begin{array}{cc}1& 5\\ 1& 8\\ 1& 8\end{array}\right]$ The first column is all $1$'s so that row $i$ dotted with $\mathbf{\theta }$ gives ${\theta }_0\cdot 1+{\theta }_1{x}^{(i)}$. Forget that column and you force the line through the origin.

Build. $\mathbf{X}=\left[\begin{array}{cc}1& 0\\ 1& 1\\ 1& 2\end{array}\right],\mathbf{y}=\left[\begin{array}{c}1\\ 3\\ 3\end{array}\right]$.

Gram matrix. ${\mathbf{X}}^T\mathbf{X}=\left[\begin{array}{cc}3& 3\\ 3& 5\end{array}\right]$ (top-left $=$ count $=3$; the $3$'s are $\sum x=0+1+2$; bottom-right $=\sum x^2=0+1+4=5$).

Right side. ${\mathbf{X}}^T\mathbf{y}=\left[\begin{array}{c}\sum y\\ \sum xy\end{array}\right]=\left[\begin{array}{c}7\\ 0\cdot 1+1\cdot 3+2\cdot 3\end{array}\right]=\left[\begin{array}{c}7\\ 9\end{array}\right]$.

Invert. $\det =3\cdot 5-3\cdot 3=6$, so $({\mathbf{X}}^T\mathbf{X}{)}^{-1}=\frac{1}{6}\left[\begin{array}{cc}5& -3\\ -3& 3\end{array}\right]$.

Solve. $\mathbf{\theta }=\frac{1}{6}\left[\begin{array}{cc}5& -3\\ -3& 3\end{array}\right]\left[\begin{array}{c}7\\ 9\end{array}\right]=\frac{1}{6}\left[\begin{array}{c}35-27\\ -21+27\end{array}\right]=\frac{1}{6}\left[\begin{array}{c}8\\ 6\end{array}\right]=\left[\begin{array}{c}4/3\\ 1\end{array}\right]$.

Answer: ${\theta }_0=4/3\approx 1.33,{\theta }_1=1$. Line: $y=1.33+x$. See the fit and its residuals below.

$\mathbf{X}=\left[\begin{array}{ccc}1& 1& 0\\ 1& 0& 1\\ 1& 1& 1\end{array}\right],\mathbf{y}=\left[\begin{array}{c}3\\ 2\\ 6\end{array}\right]$

${\mathbf{X}}^T\mathbf{X}=\left[\begin{array}{ccc}3& 2& 2\\ 2& 2& 1\\ 2& 1& 2\end{array}\right]$, ${\mathbf{X}}^T\mathbf{y}=\left[\begin{array}{c}11\\ 9\\ 8\end{array}\right]$.

Solving the $3\times 3$ system ${\mathbf{X}}^T\mathbf{X}\mathbf{\theta }={\mathbf{X}}^T\mathbf{y}$ (columns are independent, so it is invertible) gives $\mathbf{\theta }=\left[\begin{array}{c}-1\\ 4\\ 3\end{array}\right]$.

Check: row 1 $=-1+4(1)+3(0)=3$ ✓, row 2 $=-1+0+3=2$ ✓, row 3 $=-1+4+3=6$ ✓. A perfect fit, because 3 points with 3 free parameters can be matched exactly.

$\mathbf{X}=\left[\begin{array}{ccc}1& 1& 3\\ 1& 2& 6\\ 1& 3& 9\end{array}\right]$ Column 3 $=3\times$ column 2, so the columns are linearly dependent. Thus $\mathrm{rank}(\mathbf{X})=2<3$, and ${\mathbf{X}}^T\mathbf{X}$ is singular — its determinant is $0$, so no inverse exists.

Consequence: infinitely many $\mathbf{\theta }$ give the same predictions (you can trade weight between the two proportional columns). The plain normal equation cannot be used. The correct tool is the Moore-Penrose Pseudoinverse ${\mathbf{X}}^{+}$ (computed via SVD), which returns the minimum-norm least-squares solution. Alternatively, drop the redundant feature, or add a Ridge Regression penalty to restore invertibility (see L4).

At the optimum the gradient is zero: ${\mathbf{X}}^T\mathbf{X}\mathbf{\theta }-{\mathbf{X}}^T\mathbf{y}=0$. Factor: ${\mathbf{X}}^T(\mathbf{X}\mathbf{\theta }-\mathbf{y})={\mathbf{X}}^T\mathbf{r}=0.$ Geometry: every column of $\mathbf{X}$ is perpendicular ("normal") to the residual. So $\hat{\mathbf{y}}$ is the orthogonal projection of $\mathbf{y}$ onto the column space of $\mathbf{X}$ — the closest reachable point, which is exactly why the error is minimised. This is the origin of the name normal equation. The projection picture is below.

${\mathbf{X}}^T\mathbf{X}$ is symmetric positive semidefinite, so its eigenvalues are $\ge 0$; being singular means at least one eigenvalue is exactly $0$. Adding $\lambda \mathbf{I}$ with $\lambda =1$ shifts every eigenvalue up by $1$: $\text{eig}({\mathbf{X}}^T\mathbf{X}+\lambda \mathbf{I})=\text{eig}({\mathbf{X}}^T\mathbf{X})+\lambda .$ Now all eigenvalues are $\ge 1>0$, so the matrix is positive definite and hence invertible. The formula is well-defined and yields a unique $\mathbf{\theta }$ — the ridge penalty picks the shrunken solution among the previously infinite ones. This is the invertibility bridge between the normal equation and Ridge Regression.

The normal equation inverts ${\mathbf{X}}^T\mathbf{X}$, an $(n+1)\times (n+1)=10001\times 10001$ matrix, costing $\approx O(n^3)\approx 10^{12}$ operations — brutal. With $n\gg m$ it is also rank-deficient ($\mathrm{rank}\le m=500$), so the plain inverse doesn't even exist.

Reveal: this is the underdetermined regime (more unknowns than equations). Prefer Gradient Descent (cost per step $O(mn)$) or a QR Decomposition / pseudoinverse route. Rule of thumb: normal equation shines when $n$ is small (roughly $\le 10^4$) and $m>n$; otherwise iterate.

The trick: treat $x^2$ as a new feature. Design matrix uses columns $1,x,x^2$: $\mathbf{X}=\left[\begin{array}{ccc}1& -1& 1\\ 1& 0& 0\\ 1& 1& 1\end{array}\right],\mathbf{y}=\left[\begin{array}{c}1\\ 0\\ 1\end{array}\right]$ ${\mathbf{X}}^T\mathbf{X}=\left[\begin{array}{ccc}3& 0& 2\\ 0& 2& 0\\ 2& 0& 2\end{array}\right]$, ${\mathbf{X}}^T\mathbf{y}=\left[\begin{array}{c}2\\ 0\\ 2\end{array}\right]$.

Solving gives $\mathbf{\theta }=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]$, i.e. $y=x^2$. Check: $(-1{)}^2=1$ ✓, $0^2=0$ ✓, $1^2=1$ ✓. The model is nonlinear in $x$ but linear in $\mathbf{\theta }$, which is all the normal equation needs.

Predictions are $\hat{\mathbf{y}}=({\theta }_0+{\theta }_1)\left[\begin{array}{c}1\\ 1\end{array}\right]$ — only the sum $s={\theta }_0+{\theta }_1$ matters. Least squares onto the direction $(1,1)$ projects $\mathbf{y}=(2,4)$: the best $s$ is the mean-matched value $s=\frac{2+4}{2}=3$.

Among all $({\theta }_0,{\theta }_1)$ with ${\theta }_0+{\theta }_1=3$, the minimum-norm one minimises ${\theta }_0^2+{\theta }_1^2$, which by symmetry is ${\theta }_0={\theta }_1=\frac{3}{2}$. $\mathbf{\theta }=\left[\begin{array}{c}3/2\\ 3/2\end{array}\right].$ Interpretation: the pseudoinverse spreads weight equally across redundant features rather than dumping it all on one — the least-committal, most stable choice.

$\mathbf{X}=\left[\begin{array}{cc}1& 2\end{array}\right]$ ($1\times 2$), so ${\mathbf{X}}^T\mathbf{X}=\left[\begin{array}{cc}1& 2\\ 2& 4\end{array}\right]$, $\det =1\cdot 4-2\cdot 2=0$ — singular. One equation, two unknowns: underdetermined, infinitely many lines pass through one point.

The pseudoinverse gives the minimum-norm fit. Since prediction is ${\theta }_0+2{\theta }_1=5$, minimise ${\theta }_0^2+{\theta }_1^2$ on that line: geometrically the closest point on the line ${\theta }_0+2{\theta }_1=5$ to the origin. That point is $\mathbf{\theta }=\frac{5}{1^2+2^2}\left[\begin{array}{c}1\\ 2\end{array}\right]=\left[\begin{array}{c}1\\ 2\end{array}\right]$. Check: $1+2(2)=5$ ✓.

The column space of $\mathbf{X}$ has dimension $n+1$. When $m=n+1$, that column space fills all of ${\mathbb{R}}^m$, so the projection of $\mathbf{y}$ is $\mathbf{y}$ itself — residual $\mathbf{r}=0$, an exact interpolation (this is why Exercise 2.2 and 4.3 fit perfectly). When $m>n+1$, the column space is a proper subspace of ${\mathbb{R}}^m$; $\mathbf{y}$ generally sticks out of it, so its projection $\hat{\mathbf{y}}\ne \mathbf{y}$ and a non-zero residual is unavoidable — the overdetermined regime. Least squares finds the closest reachable $\hat{\mathbf{y}}$, nothing more.

The matrix you invert in the normal equation has size ==$(n+1)\times (n+1)$==. A matrix ${\mathbf{X}}^T\mathbf{X}$ is invertible exactly when $\mathbf{X}$ has full column rank. At the optimum the residual is orthogonal (normal) to the column space of $\mathbf{X}$. Adding $\lambda \mathbf{I}$ shifts every eigenvalue up by ==$\lambda$==, restoring invertibility.

When should you NOT use the plain normal equation? ::: When ${\mathbf{X}}^T\mathbf{X}$ is singular/ill-conditioned, or when $n$ is very large (costly $O(n^3)$ inverse) — use pseudoinverse/ridge/gradient descent instead. What does the pseudoinverse return for a rank-deficient system? ::: The minimum-norm least-squares solution.