Linear Regression

Table of Contents

1. Learning Objectives
2. Prerequisites
3. Notation
4. Core Intuition
5. The Linear Model and Notation
6. Least Squares from Maximum Likelihood
7. The Normal Equations: Derivation
8. Geometric Interpretation: Orthogonal Projection
9. Existence, Uniqueness, and Rank Conditions
10. Gradient Descent for Linear Regression
11. Regularization: Ridge and Lasso
12. The Bias–Variance Decomposition
13. Statistical Properties of the OLS Estimator
14. Worked Examples
15. Connection to the Broader Curriculum
16. Common Pitfalls and Misconceptions
17. Research Perspective
18. Summary of Takeaways
19. Exercises

---

Learning Objectives

After reading this chapter, you should be able to:

1. State the linear regression model in scalar, vector, and matrix form, and explain the role of the design matrix $\mathbf{X}$.
2. Derive the least-squares objective from a Gaussian noise model via maximum likelihood.
3. Prove the normal equations $\mathbf{X}^T\mathbf{X}\hat{\boldsymbol{\theta}} = \mathbf{X}^T\mathbf{y}$ and solve for $\hat{\boldsymbol{\theta}}$ when the Gram matrix is invertible.
4. Interpret ordinary least squares (OLS) as orthogonal projection of $\mathbf{y}$ onto the column space of $\mathbf{X}$.
5. Characterize when the OLS solution is unique in terms of rank and linear independence of features.
6. Derive the gradient descent update for linear regression and connect it to the condition number of $\mathbf{X}^T\mathbf{X}$.
7. Derive closed-form Ridge regression and explain why $\lambda \mathbf{I}$ restores invertibility.
8. State the bias–variance decomposition for squared error and identify the sources of each term.

---

Prerequisites

This chapter assumes familiarity with:

• Vectors, Spans & Linear Independence — span, rank, linear independence, inner products
• Matrix Operations & Linear Transformations — column space, rank-nullity, matrix multiplication
• Gradient Descent — first-order optimization, convexity, convergence

Optional but helpful: Positive Definite Matrices for regularization theory, and Maximum Likelihood Estimation for the probabilistic derivation.

---

Notation

• $\mathbf{X} \in \mathbb{R}^{N \times (d+1)}$ — Design matrix
• $\mathbf{y} \in \mathbb{R}^N$ — Target vector
• $\boldsymbol{\theta} \in \mathbb{R}^{d+1}$ — Parameter vector (including bias)
• $\tilde{\mathbf{x}}_i$ — Augmented feature vector with leading 1
• $\hat{\boldsymbol{\theta}}_{\mathrm{OLS}}$ — Ordinary least squares estimator
• $\lambda$ — Regularization strength (Ridge/Lasso)

---

Core Intuition

Linear regression is the simplest supervised learning model — and the most important. Every neural network layer begins as a linear map; attention computes linear combinations; diffusion models predict noise through linear projections in latent space. Understanding linear regression deeply means understanding the template that all of supervised learning follows: specify a model, define a loss, optimize.

The question is deceptively simple. Given pairs $(\mathbf{x}_i, y_i)$ for $i = 1, \ldots, N$, find parameters $\boldsymbol{\theta}$ such that predictions $\hat{y}_i = f(\mathbf{x}_i; \boldsymbol{\theta})$ are close to observed targets $y_i$. When $f$ is linear in $\boldsymbol{\theta}$, the problem admits:

• A closed-form solution (the normal equations)
• A geometric interpretation (projection onto a subspace)
• A probabilistic foundation (Gaussian noise $\implies$ least squares)
• A convex optimization landscape (single global minimum)

These properties make linear regression the canonical entry point to machine learning theory. Logistic regression replaces the Gaussian likelihood with Bernoulli; neural networks replace linear $f$ with compositions of linear maps and nonlinearities; but the loss–gradient–update pattern established here persists throughout.

Series context. This opens Chapter 4 (Supervised Learning) in Volume II: Classical Machine Learning. Subsequent articles cover Logistic Regression, Bayesian Linear Regression, and the Bias–Variance Tradeoff.

---

The Linear Model and Notation

Definition 1 (Linear Regression Model). Given feature vectors $\mathbf{x}_i \in \mathbb{R}^d$ and scalar targets $y_i \in \mathbb{R}$ for $i = 1, \ldots, N$, the linear regression model predicts

$$\hat{y}_i = \boldsymbol{\theta}^T \tilde{\mathbf{x}}_i = \theta_0 + \sum_{j=1}^{d} \theta_j x_{ij} \tag{1}$$

where $\tilde{\mathbf{x}}_i = (1, x_{i1}, \ldots, x_{id})^T \in \mathbb{R}^{d+1}$ is the augmented feature vector (bias absorbed via a leading coordinate of 1), and $\boldsymbol{\theta} = (\theta_0, \theta_1, \ldots, \theta_d)^T \in \mathbb{R}^{d+1}$ is the parameter vector.

Definition 2 (Design Matrix). Stack augmented features into the design matrix

$$\mathbf{X} = \begin{bmatrix} \tilde{\mathbf{x}}_1^T \\ \vdots \\ \tilde{\mathbf{x}}_N^T \end{bmatrix} \in \mathbb{R}^{N \times (d+1)}, \quad \mathbf{y} = \begin{bmatrix} y_1 \\ \vdots \\ y_N \end{bmatrix} \in \mathbb{R}^N. \tag{2}$$

Then predictions for all samples simultaneously are

$$\hat{\mathbf{y}} = \mathbf{X}\boldsymbol{\theta}. \tag{3}$$

Important equation. Equation (3) is the matrix form of every linear model in machine learning. The design matrix $\mathbf{X}$ maps parameters to predictions; its column space $\mathcal{C}(\mathbf{X})$ is the set of all vectors reachable by some choice of $\boldsymbol{\theta}$.

Definition 3 (Hypothesis Class). The hypothesis class for linear regression is

$$\mathcal{H} = \{ f_{\boldsymbol{\theta}} : \mathbf{x} \mapsto \boldsymbol{\theta}^T \tilde{\mathbf{x}} \mid \boldsymbol{\theta} \in \mathbb{R}^{d+1} \}. \tag{4}$$

This is the set of all affine functions on $\mathbb{R}^d$. It is a $(d+1)$-dimensional vector space when features are non-degenerate.

---

Least Squares from Maximum Likelihood

We now derive why squared error is the natural loss function — not by fiat, but from a probabilistic model.

Assumption 1 (Gaussian Noise Model). Observations arise as

$$y_i = \boldsymbol{\theta}^T \tilde{\mathbf{x}}_i + \varepsilon_i, \quad \varepsilon_i \stackrel{\text{i.i.d.}}{\sim} \mathcal{N}(0, \sigma^2). \tag{5}$$

Equivalently, $y_i \mid \mathbf{x}_i, \boldsymbol{\theta} \sim \mathcal{N}(\boldsymbol{\theta}^T \tilde{\mathbf{x}}_i, \sigma^2)$.

Proposition 1. Under Assumption 1, the maximum likelihood estimator of $\boldsymbol{\theta}$ is identical to the minimizer of the sum of squared residuals:

$$\hat{\boldsymbol{\theta}}_{\text{MLE}} = \arg\min_{\boldsymbol{\theta}} \|\mathbf{X}\boldsymbol{\theta} - \mathbf{y}\|_2^2. \tag{6}$$

Proof. The log-likelihood for independent Gaussian observations is

$$\ell(\boldsymbol{\theta}) = -\frac{N}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^{N}(y_i - \boldsymbol{\theta}^T \tilde{\mathbf{x}}_i)^2. \tag{7}$$

Since $\sigma^2$ is constant with respect to $\boldsymbol{\theta}$, maximizing $\ell(\boldsymbol{\theta})$ is equivalent to minimizing $\sum_i (y_i - \boldsymbol{\theta}^T \tilde{\mathbf{x}}_i)^2 = \|\mathbf{X}\boldsymbol{\theta} - \mathbf{y}\|_2^2$. $\blacksquare$

Definition 4 (Empirical Risk / MSE Loss). The mean squared error objective is

$$J(\boldsymbol{\theta}) = \frac{1}{2N}\|\mathbf{X}\boldsymbol{\theta} - \mathbf{y}\|_2^2 = \frac{1}{2N}\sum_{i=1}^{N}(\hat{y}_i - y_i)^2. \tag{8}$$

The factor $\frac{1}{2N}$ is a convention: the $\frac{1}{2}$ cancels upon differentiation; the $\frac{1}{N}$ averages over samples. Minimizers are identical regardless of these constants.

Definition 5 (Ordinary Least Squares). The OLS estimator is

$$\hat{\boldsymbol{\theta}}_{\text{OLS}} = \arg\min_{\boldsymbol{\theta} \in \mathbb{R}^{d+1}} J(\boldsymbol{\theta}). \tag{9}$$

---

The Normal Equations: Derivation

We derive the closed-form solution step by step.

Theorem 1 (Normal Equations). A necessary condition for $\hat{\boldsymbol{\theta}}$ to minimize $J(\boldsymbol{\theta})$ is

$$\mathbf{X}^T\mathbf{X}\,\hat{\boldsymbol{\theta}} = \mathbf{X}^T\mathbf{y}. \tag{10}$$

If $\mathbf{X}^T\mathbf{X}$ is invertible, the unique minimizer is

$$\hat{\boldsymbol{\theta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}. \tag{11}$$

Proof. Expand the objective:

$$J(\boldsymbol{\theta}) = \frac{1}{2N}(\mathbf{X}\boldsymbol{\theta} - \mathbf{y})^T(\mathbf{X}\boldsymbol{\theta} - \mathbf{y}) = \frac{1}{2N}\left(\boldsymbol{\theta}^T\mathbf{X}^T\mathbf{X}\boldsymbol{\theta} - 2\boldsymbol{\theta}^T\mathbf{X}^T\mathbf{y} + \mathbf{y}^T\mathbf{y}\right). \tag{12}$$

Differentiating with respect to $\boldsymbol{\theta}$ (using $\nabla_{\boldsymbol{\theta}}(\boldsymbol{\theta}^T\mathbf{A}\boldsymbol{\theta}) = 2\mathbf{A}\boldsymbol{\theta}$ for symmetric $\mathbf{A}$):

$$\nabla_{\boldsymbol{\theta}} J = \frac{1}{N}\left(\mathbf{X}^T\mathbf{X}\boldsymbol{\theta} - \mathbf{X}^T\mathbf{y}\right). \tag{13}$$

Setting $\nabla_{\boldsymbol{\theta}} J = \mathbf{0}$ yields (10). Since $J$ is convex (as $\mathbf{X}^T\mathbf{X}$ is positive semidefinite), any critical point is a global minimum. If $\mathbf{X}^T\mathbf{X}$ is positive definite (hence invertible), the solution is unique and given by (11). $\blacksquare$

Definition 6 (Moore–Penrose Solution). When $\mathbf{X}^T\mathbf{X}$ is singular, (10) has infinitely many solutions. The minimum-norm least-squares solution is

$$\hat{\boldsymbol{\theta}} = \mathbf{X}^+\mathbf{y}, \tag{14}$$

where $\mathbf{X}^+ = (\mathbf{X}^T\mathbf{X})^+\mathbf{X}^T$ is the pseudoinverse. This is developed fully in Singular Value Decomposition.

Important equation. Equation (11) is the normal equation (or closed-form solution). It requires inverting a $(d+1) \times (d+1)$ matrix — cost $O(d^3)$ — making it impractical when $d$ is in the millions, even though forming $\mathbf{X}^T\mathbf{X}$ costs only $O(Nd^2)$.

---

Geometric Interpretation: Orthogonal Projection

The least-squares problem has an elegant geometric interpretation that clarifies many properties of OLS.

Theorem 2 (Projection Characterization). Let $\hat{\mathbf{y}} = \mathbf{X}\hat{\boldsymbol{\theta}}$ be the OLS fitted values. Then $\hat{\mathbf{y}}$ is the orthogonal projection of $\mathbf{y}$ onto $\mathcal{C}(\mathbf{X})$, the column space of $\mathbf{X}$. The residual $\mathbf{r} = \mathbf{y} - \hat{\mathbf{y}}$ is orthogonal to every column of $\mathbf{X}$:

$$\mathbf{X}^T\mathbf{r} = \mathbf{0}. \tag{15}$$

Proof. By definition, $\hat{\mathbf{y}} \in \mathcal{C}(\mathbf{X})$. The residual is $\mathbf{r} = \mathbf{y} - \mathbf{X}\hat{\boldsymbol{\theta}}$. From the normal equations (10): $\mathbf{X}^T(\mathbf{y} - \mathbf{X}\hat{\boldsymbol{\theta}}) = \mathbf{0}$, which is precisely (15). This is the orthogonality condition characterizing projection onto a subspace. $\blacksquare$

Corollary 1 (Pythagorean Decomposition). The total squared norm decomposes as

$$\|\mathbf{y}\|^2 = \|\hat{\mathbf{y}}\|^2 + \|\mathbf{r}\|^2, \tag{16}$$

since $\hat{\mathbf{y}} \perp \mathbf{r}$.

Interpretation. OLS finds the point in the column space of $\mathbf{X}$ closest to $\mathbf{y}$ in Euclidean distance. When $N > d+1$ (overdetermined system), we cannot fit $\mathbf{y}$ exactly; we find the best approximation. When $N < d+1$ (underdetermined), infinitely many $\boldsymbol{\theta}$ yield zero training error; the pseudoinverse selects the minimum-norm parameter vector.

---

Existence, Uniqueness, and Rank Conditions

Theorem 3 (Uniqueness Condition). The OLS solution $\hat{\boldsymbol{\theta}}$ is unique if and only if $\mathbf{X}$ has full column rank, i.e., $\text{rank}(\mathbf{X}) = d+1$.

Proof. $\mathbf{X}^T\mathbf{X}$ is invertible if and only if $\mathbf{X}$ has full column rank (see Positive Definite Matrices). Full column rank means the columns of $\mathbf{X}$ are linearly independent — no feature (including the bias column) is a linear combination of the others. $\blacksquare$

Proposition 2 (Training Error and Rank). If $\text{rank}(\mathbf{X}) = N$ (full row rank, requires $N \leq d+1$), there exists $\boldsymbol{\theta}$ such that $\mathbf{X}\boldsymbol{\theta} = \mathbf{y}$ exactly — zero training error.

When features are collinear (multicollinearity), $\mathbf{X}^T\mathbf{X}$ is near-singular: small perturbations in $\mathbf{y}$ cause large changes in $\hat{\boldsymbol{\theta}}$. The condition number $\kappa(\mathbf{X}^T\mathbf{X}) = \lambda_{\max}/\lambda_{\min}$ quantifies this sensitivity via the eigenvalue spread of the Gram matrix.

---

Gradient Descent for Linear Regression

When $d$ is large, forming and inverting $\mathbf{X}^T\mathbf{X}$ is prohibitive. Gradient descent provides an iterative alternative.

From (13), the gradient of the unscaled loss $\tilde{J}(\boldsymbol{\theta}) = \frac{1}{2}\|\mathbf{X}\boldsymbol{\theta} - \mathbf{y}\|^2$ is

$$\nabla_{\boldsymbol{\theta}} \tilde{J} = \mathbf{X}^T(\mathbf{X}\boldsymbol{\theta} - \mathbf{y}). \tag{17}$$

The update rule with learning rate $\alpha > 0$:

$$\boldsymbol{\theta}_{t+1} = \boldsymbol{\theta}_t - \alpha \,\mathbf{X}^T(\mathbf{X}\boldsymbol{\theta}_t - \mathbf{y}). \tag{18}$$

Proposition 3 (Convergence in the Eigenbasis). Let $\mathbf{X}^T\mathbf{X} = \mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^T$ be the eigendecomposition with eigenvalues $\lambda_1 \geq \cdots \geq \lambda_{d+1} > 0$. Gradient descent with step size $\alpha < 2/\lambda_1$ converges to $\hat{\boldsymbol{\theta}}$. The convergence rate per step in direction $\mathbf{q}_j$ (eigenvector) is governed by $|1 - \alpha\lambda_j|$.

When $\kappa(\mathbf{X}^T\mathbf{X})$ is large (ill-conditioned features), convergence is slow — a theme developed in Gradient Descent and addressed by preconditioning and adaptive methods (Adam, etc.).

Proposition 4 (Stochastic Gradient Descent). Replacing the full gradient (17) with an unbiased estimate based on a mini-batch $\mathcal{B} \subset \{1,\ldots,N\}$:

$$\tilde{\nabla}_{\boldsymbol{\theta}} = \frac{1}{|\mathcal{B}|}\sum_{i \in \mathcal{B}} \tilde{\mathbf{x}}_i(\boldsymbol{\theta}^T\tilde{\mathbf{x}}_i - y_i) \tag{19}$$

yields an unbiased gradient estimator with reduced per-step cost $O(|\mathcal{B}| \cdot d)$ instead of $O(Nd)$.

---

Regularization: Ridge and Lasso

When $\mathbf{X}^T\mathbf{X}$ is singular or ill-conditioned, regularization stabilizes estimation.

Definition 7 (Ridge Regression / L2 Penalty). The Ridge objective adds an $\ell_2$ penalty:

$$J_{\text{Ridge}}(\boldsymbol{\theta}) = \frac{1}{2N}\|\mathbf{X}\boldsymbol{\theta} - \mathbf{y}\|^2 + \lambda \sum_{j=1}^{d} \theta_j^2 = \frac{1}{2N}\|\mathbf{X}\boldsymbol{\theta} - \mathbf{y}\|^2 + \lambda \|\boldsymbol{\theta}_{1:d}\|_2^2, \tag{20}$$

where $\boldsymbol{\theta}_{1:d} = (\theta_1, \ldots, \theta_d)^T$ excludes the bias (typically not penalized).

Theorem 4 (Ridge Closed Form). The Ridge estimator is

$$\hat{\boldsymbol{\theta}}_{\text{Ridge}} = \left(\mathbf{X}^T\mathbf{X} + 2N\lambda \,\mathbf{D}\right)^{-1}\mathbf{X}^T\mathbf{y}, \tag{21}$$

where $\mathbf{D} = \text{diag}(0, 1, \ldots, 1) \in \mathbb{R}^{(d+1)\times(d+1)}$.

Proof. Setting the gradient of (20) to zero:

$$\frac{1}{N}\mathbf{X}^T(\mathbf{X}\boldsymbol{\theta} - \mathbf{y}) + 2\lambda \mathbf{D}\boldsymbol{\theta} = \mathbf{0}. \tag{22}$$

Rearranging: $(\mathbf{X}^T\mathbf{X} + 2N\lambda \mathbf{D})\boldsymbol{\theta} = \mathbf{X}^T\mathbf{y}$. For any $\lambda > 0$, the matrix $\mathbf{X}^T\mathbf{X} + 2N\lambda \mathbf{D}$ is positive definite and hence invertible. $\blacksquare$

Definition 8 (Lasso / L1 Penalty). The Lasso objective uses an $\ell_1$ penalty:

$$J_{\text{Lasso}}(\boldsymbol{\theta}) = \frac{1}{2N}\|\mathbf{X}\boldsymbol{\theta} - \mathbf{y}\|^2 + \lambda \|\boldsymbol{\theta}_{1:d}\|_1. \tag{23}$$

Unlike Ridge, Lasso has no closed-form solution. The $\ell_1$ penalty induces sparsity: many $\theta_j$ are driven exactly to zero, performing implicit feature selection. The geometry of the $\ell_1$ ball (corners on coordinate axes) explains this — the level sets of $\|\boldsymbol{\theta}\|_1$ are non-smooth at axes, so optima often lie on sparse faces.

Proposition 5 (Bayesian Interpretation of Ridge). Ridge regression corresponds to MAP estimation under a Gaussian prior $\boldsymbol{\theta}_{1:d} \sim \mathcal{N}(\mathbf{0}, \frac{1}{2\lambda}\mathbf{I})$ with the Gaussian noise model (5). See Bayesian Linear Regression.

---

The Bias–Variance Decomposition

Understanding generalization requires decomposing expected prediction error.

Theorem 5 (Bias–Variance Decomposition). Consider the model $y = f(\mathbf{x}) + \varepsilon$ with $\mathbb{E}[\varepsilon] = 0$, $\text{Var}(\varepsilon) = \sigma^2$. Let $\hat{f}(\mathbf{x})$ be an estimator trained on a random dataset $\mathcal{D}$. For a fixed test point $\mathbf{x}$, the expected squared error decomposes as

$$\mathbb{E}_{\mathcal{D}, \varepsilon}\left[(y - \hat{f}(\mathbf{x}))^2\right] = \underbrace{\left(\mathbb{E}_{\mathcal{D}}[\hat{f}(\mathbf{x})] - f(\mathbf{x})\right)^2}_{\text{Bias}^2} + \underbrace{\mathbb{E}_{\mathcal{D}}\left[(\hat{f}(\mathbf{x}) - \mathbb{E}_{\mathcal{D}}[\hat{f}(\mathbf{x})])^2\right]}_{\text{Variance}} + \underbrace{\sigma^2}_{\text{Irreducible noise}}. \tag{24}$$

Proof. Add and subtract $\mathbb{E}[\hat{f}]$ and expand; cross terms vanish by independence of $\varepsilon$ from $\mathcal{D}$. See Bias–Variance Tradeoff for the full derivation. $\blacksquare$

Interpretation for linear regression:

• Low complexity (few features, large $\lambda$): high bias, low variance — underfitting
• High complexity (many features, $\lambda = 0$ with collinearity): low bias, high variance — overfitting
• Ridge increases bias slightly but reduces variance substantially when $\mathbf{X}^T\mathbf{X}$ is ill-conditioned

The effective degrees of freedom of Ridge regression is $\text{tr}(\mathbf{X}(\mathbf{X}^T\mathbf{X} + \lambda \mathbf{I})^{-1}\mathbf{X}^T)$, a smooth measure of model complexity.

---

Statistical Properties of the OLS Estimator

Under standard assumptions, OLS enjoys strong statistical guarantees.

Assumption 2 (Gauss–Markov Conditions). (i) Linearity: $\mathbb{E}[y_i \mid \mathbf{x}_i] = \boldsymbol{\theta}^T\tilde{\mathbf{x}}_i$; (ii) Homoscedasticity: $\text{Var}(\varepsilon_i \mid \mathbf{x}_i) = \sigma^2$; (iii) Uncorrelated errors: $\text{Cov}(\varepsilon_i, \varepsilon_j) = 0$ for $i \neq j$; (iv) Full rank: $\text{rank}(\mathbf{X}) = d+1$.

Theorem 6 (Gauss–Markov). Under Assumption 2, the OLS estimator $\hat{\boldsymbol{\theta}}$ is the Best Linear Unbiased Estimator (BLUE): among all linear unbiased estimators, it has minimum variance.

Proposition 6 (Variance of OLS). Under homoscedastic Gaussian errors:

$$\text{Cov}(\hat{\boldsymbol{\theta}}) = \sigma^2 (\mathbf{X}^T\mathbf{X})^{-1}. \tag{25}$$

Large entries in $(\mathbf{X}^T\mathbf{X})^{-1}$ — caused by collinear features — inflate the variance of individual coefficient estimates.

---

Worked Examples

Example 1: Simple Linear Regression

Given $N = 3$ points with scalar feature: $(1, 2), (2, 3), (3, 5)$. Augmented design matrix and target:

$$\mathbf{X} = \begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \end{bmatrix}, \quad \mathbf{y} = \begin{bmatrix} 2 \\ 3 \\ 5 \end{bmatrix}. \tag{26}$$

Compute:

$$\mathbf{X}^T\mathbf{X} = \begin{bmatrix} 3 & 6 \\ 6 & 14 \end{bmatrix}, \quad \mathbf{X}^T\mathbf{y} = \begin{bmatrix} 10 \\ 23 \end{bmatrix}. \tag{27}$$ $$\hat{\boldsymbol{\theta}} = \begin{bmatrix} 3 & 6 \\ 6 & 14 \end{bmatrix}^{-1} \begin{bmatrix} 10 \\ 23 \end{bmatrix} = \frac{1}{6}\begin{bmatrix} 14 & -6 \\ -6 & 3 \end{bmatrix}\begin{bmatrix} 10 \\ 23 \end{bmatrix} = \begin{bmatrix} 0.33 \\ 1.5 \end{bmatrix}. \tag{28}$$

The fitted line is $\hat{y} = 0.33 + 1.5x$.

Example 2: Orthogonality of Residuals

For the solution above, $\hat{\mathbf{y}} = \mathbf{X}\hat{\boldsymbol{\theta}}$ and $\mathbf{r} = \mathbf{y} - \hat{\mathbf{y}}$. Direct computation confirms $\mathbf{X}^T\mathbf{r} = \mathbf{0}$ — the residual is orthogonal to both the bias column and the feature column.

Example 3: Ridge Stabilization

Suppose two features are nearly collinear: $\mathbf{x}_1 \approx \mathbf{x}_2$. Then $\mathbf{X}^T\mathbf{X}$ has a very small eigenvalue $\lambda_{\min}$, and $(\mathbf{X}^T\mathbf{X})^{-1}$ has a very large eigenvalue $1/\lambda_{\min}$. Adding $\lambda \mathbf{I}$ shifts all eigenvalues by $2N\lambda$, bounding $\kappa(\mathbf{X}^T\mathbf{X} + 2N\lambda \mathbf{I})$ and stabilizing $\hat{\boldsymbol{\theta}}$.

---

Connection to the Broader Curriculum

Linear regression is the Rosetta Stone of machine learning:

• Linear model $\mathbf{X}\boldsymbol{\theta}$ — Neural network layers
• Squared loss — Cross-entropy, Huber, contrastive losses
• Normal equations — No closed form; backpropagation instead
• Ridge/Lasso — Weight decay, dropout, all regularization
• Bias–variance — Universal generalization framework
• Design matrix rank — Representation rank, LoRA low-rank structure

Logistic Regression replaces the Gaussian likelihood with Bernoulli and the identity link with the logistic function. Backpropagation generalizes the gradient computation to arbitrary computation graphs. Self-Attention computes similarity-weighted linear combinations — regression over values indexed by queries.

---

Common Pitfalls and Misconceptions

Pitfall 1: Applying the normal equation when $N < d$. When there are more parameters than samples, $\mathbf{X}^T\mathbf{X}$ is singular. The normal equation as written in (11) does not apply; use the pseudoinverse or regularization.

Pitfall 2: Regularizing the bias term. Penalizing $\theta_0$ shifts predictions systematically and is usually undesirable. Always exclude the intercept from $\ell_1$ and $\ell_2$ penalties unless features are centered and targets are demeaned.

Pitfall 3: Confusing $R^2$ with model quality. High $R^2$ on training data does not imply good generalization. $R^2 = 1 - \|\mathbf{r}\|^2/\|\mathbf{y} - \bar{y}\mathbf{1}\|^2$ measures fit on the training set only.

Pitfall 4: Ignoring feature scaling before gradient descent. Features on different scales yield an ill-conditioned $\mathbf{X}^T\mathbf{X}$, slowing convergence. Standardization ($\tilde{x}_j = (x_j - \mu_j)/\sigma_j$) improves conditioning without changing the OLS solution (if the bias is unpenalized).

Pitfall 5: Treating linear regression as "always the right first model." Linearity in parameters does not mean linearity in features. Polynomial, interaction, and basis-expansion features can capture nonlinearity while preserving the closed-form solution.

Pitfall 6: Assuming homoscedasticity. When $\text{Var}(\varepsilon_i)$ depends on $\mathbf{x}_i$, OLS is still unbiased but no longer efficient. Weighted least squares or robust methods are appropriate.

---

Research Perspective

Linear regression dates to Gauss and Legendre's work on least squares in the early nineteenth century for astronomical orbit prediction. The probabilistic formulation through maximum likelihood connects to Fisher's foundations of statistical inference in the 1920s.

In modern machine learning, linear models remain indispensable:

• Generalized Linear Models (GLMs) extend the linear predictor to exponential-family distributions (Nelder & Wedderburn, 1972).
• Kernel Ridge Regression implicitly maps features to infinite-dimensional spaces while retaining a closed-form solution via the representer theorem.
• Random Features (Rahimi & Recht, 2007) approximate kernel methods with explicit finite-dimensional linear models.
• Scaling laws for language models show that linear trends in log-loss vs. compute persist across orders of magnitude — the simplest model class still governs the largest systems.

The resurgence of interest in interpretability has renewed focus on linear probes and linear concept directions in embedding spaces, demonstrating that linear algebra remains central even in the era of billion-parameter models.

---

Summary of Takeaways

• Linear model — $\hat{y} = \boldsymbol{\theta}^T\tilde{\mathbf{x}}$ — Simplest supervised hypothesis class
• MSE loss — $\frac{1}{2N}\ — \mathbf{X}\boldsymbol{\theta} - \mathbf{y}\ — ^2$ — MLE under Gaussian noise
• Normal equations — $\mathbf{X}^T\mathbf{X}\hat{\boldsymbol{\theta}} = \mathbf{X}^T\mathbf{y}$ — Closed-form OLS
• Projection — $\hat{\mathbf{y}} = \text{proj}_{\mathcal{C}(\mathbf{X})}(\mathbf{y})$ — Geometric interpretation
• Uniqueness — $\text{rank}(\mathbf{X}) = d+1$ — Full column rank required
• Ridge — $(\mathbf{X}^T\mathbf{X} + 2N\lambda \mathbf{D})^{-1}\mathbf{X}^T\mathbf{y}$ — Stabilizes ill-conditioned problems
• Bias–variance — Error = Bias$^2$ + Variance + $\sigma^2$ — Generalization decomposition

Next article: Logistic Regression → — where the linear predictor is composed with the logistic function and Bernoulli likelihood replaces the Gaussian.

---

Exercises

Exercise 1 (Derivation). Starting from the Gaussian log-likelihood (7), verify that $\hat{\sigma}^2_{\text{MLE}} = \frac{1}{N}\|\mathbf{X}\hat{\boldsymbol{\theta}} - \mathbf{y}\|^2$ after substituting the MLE for $\boldsymbol{\theta}$.

Exercise 2 (Projection). Prove that among all $\mathbf{v} \in \mathcal{C}(\mathbf{X})$, the vector $\hat{\mathbf{y}} = \mathbf{X}\hat{\boldsymbol{\theta}}$ uniquely minimizes $\|\mathbf{y} - \mathbf{v}\|_2$. Use the Pythagorean theorem with the orthogonal decomposition $\mathbf{y} = \hat{\mathbf{y}} + \mathbf{r}$.

Exercise 3 (Rank). Construct a design matrix $\mathbf{X} \in \mathbb{R}^{4 \times 3}$ with rank 2. Show that the normal equations have infinitely many solutions and characterize the solution set as an affine subspace of $\mathbb{R}^3$.

Exercise 4 (Ridge). Prove that the eigenvalues of $\mathbf{X}^T\mathbf{X} + \lambda \mathbf{I}$ are $\lambda_j + \lambda$ where $\lambda_j$ are eigenvalues of $\mathbf{X}^T\mathbf{X}$. Conclude that Ridge always produces a unique solution for $\lambda > 0$.

Exercise 5 (Bias–variance). For the simple model $y = \theta x + \varepsilon$ with a single observation $(x_1, y_1)$ and known $x_1 \neq 0$, compute $\mathbb{E}[\hat{\theta}]$ and $\text{Var}(\hat{\theta})$ where $\hat{\theta} = y_1/x_1$. Verify the bias–variance decomposition at a test point $x_*$.

Exercise 6 (Condition number). Let $\mathbf{X}$ have columns $\mathbf{x}_1, \mathbf{x}_2$ with $\|\mathbf{x}_1\| = \|\mathbf{x}_2\| = 1$ and $\mathbf{x}_1^T\mathbf{x}_2 = \rho$. Compute $\mathbf{X}^T\mathbf{X}$ and its condition number as a function of $\rho$. What happens as $\rho \to 1$?

Exercise 7 (Conceptual). Explain why centering features ($\tilde{x}_{ij} = x_{ij} - \bar{x}_j$) does not change the fitted values $\hat{\mathbf{y}}$, but does change the interpretation of $\theta_0$. How does centering affect the Ridge penalty if the bias is not excluded?

Exercise 8 (Connection). Show that the gradient descent update (18) can be written as $\boldsymbol{\theta}_{t+1} = (\mathbf{I} - \alpha \mathbf{X}^T\mathbf{X})\boldsymbol{\theta}_t + \alpha \mathbf{X}^T\mathbf{y}$. Analyze convergence via the spectral radius of $(\mathbf{I} - \alpha \mathbf{X}^T\mathbf{X})$.