Matrix Formulation

In simple linear regression, we modeled the relationship between a single predictor variable and a response. The mathematics was elegant: a slope, an intercept, and clear geometric intuition. But real-world problems rarely involve a single predictor.

Consider predicting house prices. You don't just have square footage—you have the number of bedrooms, bathrooms, age of the house, lot size, proximity to schools, crime rates, and dozens of other variables. Each adds predictive power, but also adds mathematical complexity.

The key insight: When we have multiple predictors, we transition from scalar algebra to matrix algebra. This isn't merely a notational convenience—it's a fundamental shift that unlocks computational efficiency, theoretical elegance, and connections to geometry that would be impossible to express otherwise.

Let's establish the model formally. In multiple linear regression, we assume that a response variable $y$ is a linear combination of $p$ predictor variables plus random noise:

$$y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_p x_{ip} + \epsilon_i$$

where:

• $y_i$ is the response for observation $i$
• $x_{i1}, x_{i2}, \ldots, x_{ip}$ are the $p$ predictor values for observation $i$
• $\beta_0$ is the intercept
• $\beta_1, \ldots, \beta_p$ are the regression coefficients
• $\epsilon_i$ is the random error term

For $n$ observations, we have $n$ such equations, one for each data point. Writing them all out individually becomes unwieldy—this is where matrix notation becomes indispensable.

The summation form:

We can write the model more compactly using summation notation:

$$y_i = \beta_0 + \sum_{j=1}^{p} \beta_j x_{ij} + \epsilon_i \quad \text{for } i = 1, 2, \ldots, n$$

This captures all $n$ equations but still requires index manipulation. To derive estimators and prove properties, we need something even more compact. Enter matrix formulation.

The design matrix (also called the model matrix or data matrix) is the central object that organizes all predictor information. We construct it by stacking observations as rows and features as columns.

Definition: Given $n$ observations and $p$ predictor variables, the design matrix $\mathbf{X}$ is an $n \times (p+1)$ matrix:

$$\mathbf{X} = \begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1p} \ 1 & x_{21} & x_{22} & \cdots & x_{2p} \ \vdots & \vdots & \vdots & \ddots & \vdots \ 1 & x_{n1} & x_{n2} & \cdots & x_{np} \end{bmatrix}$$

Critical observation: The first column is all 1s. This is the intercept column, and it's essential for the matrix formulation to work correctly.

Reading the design matrix:

• Row i contains all feature values for observation $i$ (preceded by a 1 for the intercept)
• Column j (for $j \geq 1$) contains all values of predictor $j$ across observations
• Column 0 is the intercept column of ones

The design matrix encodes the entire training dataset's feature information.

Now we can express the entire multiple regression model as a single matrix equation:

$$\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}$$

where:

• $\mathbf{y} = \begin{bmatrix} y_1 \ y_2 \ \vdots \ y_n \end{bmatrix}$ is the $n \times 1$ response vector

• $\boldsymbol{\beta} = \begin{bmatrix} \beta_0 \ \beta_1 \ \vdots \ \beta_p \end{bmatrix}$ is the $(p+1) \times 1$ coefficient vector

• $\boldsymbol{\epsilon} = \begin{bmatrix} \epsilon_1 \ \epsilon_2 \ \vdots \ \epsilon_n \end{bmatrix}$ is the $n \times 1$ error vector

This single equation replaces n equations. Every derivation, proof, and algorithm flows from this compact representation.

Unpacking the multiplication:

To see why this works, consider the product $\mathbf{X}\boldsymbol{\beta}$. The $i$-th element is:

$$(\mathbf{X}\boldsymbol{\beta})i = 1 \cdot \beta_0 + x{i1} \cdot \beta_1 + x_{i2} \cdot \beta_2 + \cdots + x_{ip} \cdot \beta_p$$

which is exactly $\beta_0 + \sum_{j=1}^{p} \beta_j x_{ij}$—the linear predictor for observation $i$. The matrix form is not an approximation or simplification; it's an exact encoding.

Once we have estimated coefficients $\hat{\boldsymbol{\beta}}$ (we'll derive these in the next page), the fitted values or predictions are:

$$\hat{\mathbf{y}} = \mathbf{X}\hat{\boldsymbol{\beta}}$$

The residual vector captures the difference between observed and predicted values:

$$\mathbf{e} = \mathbf{y} - \hat{\mathbf{y}} = \mathbf{y} - \mathbf{X}\hat{\boldsymbol{\beta}}$$

Important distinction:

• $\boldsymbol{\epsilon}$ represents the true (unobservable) errors in the population model
• $\mathbf{e}$ represents the residuals, which are the observable estimates of those errors

Predicting new data:

For a new observation with features $\mathbf{x}_{\text{new}} = [1, x_1, x_2, \ldots, x_p]^\top$, the predicted response is:

$$\hat{y}{\text{new}} = \mathbf{x}{\text{new}}^\top \hat{\boldsymbol{\beta}}$$

The same coefficient vector is used; only the feature values change. This is the core of the linear model's utility—parameters learned from training data apply directly to new predictions.

The goal of Ordinary Least Squares (OLS) is to find the coefficient vector $\boldsymbol{\beta}$ that minimizes the sum of squared residuals. In scalar form, we minimize:

$$\text{SSE} = \sum_{i=1}^{n} (y_i - \hat{y}i)^2 = \sum{i=1}^{n} (y_i - \beta_0 - \beta_1 x_{i1} - \cdots - \beta_p x_{ip})^2$$

In matrix form, this becomes remarkably elegant:

$$\text{SSE}(\boldsymbol{\beta}) = (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^\top(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) = |\mathbf{y} - \mathbf{X}\boldsymbol{\beta}|^2$$

Expanding the objective:

Let's expand the matrix expression to connect with optimization:

$$\begin{aligned} \text{SSE}(\boldsymbol{\beta}) &= (\mathbf{y} - \mathbf{X}\boldsymbol{\beta})^\top(\mathbf{y} - \mathbf{X}\boldsymbol{\beta}) &= \mathbf{y}^\top\mathbf{y} - \mathbf{y}^\top\mathbf{X}\boldsymbol{\beta} - \boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{y} + \boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta} &= \mathbf{y}^\top\mathbf{y} - 2\boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{y} + \boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta} \end{aligned}$$

Note that $\mathbf{y}^\top\mathbf{X}\boldsymbol{\beta}$ and $\boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{y}$ are both scalars and are transposes of each other, hence equal. This quadratic form in $\boldsymbol{\beta}$ is the starting point for deriving the normal equations.

The matrix formulation isn't just elegant notation—it provides fundamental advantages across theory, computation, and intuition.

A matrix that appears repeatedly in regression theory is the Gram matrix $\mathbf{X}^\top\mathbf{X}$. Understanding its structure is crucial.

Definition: The Gram matrix is a $(p+1) \times (p+1)$ symmetric matrix:

$$\mathbf{X}^\top\mathbf{X} = \begin{bmatrix} \mathbf{x}_0^\top\mathbf{x}_0 & \mathbf{x}_0^\top\mathbf{x}_1 & \cdots & \mathbf{x}_0^\top\mathbf{x}_p \ \mathbf{x}_1^\top\mathbf{x}_0 & \mathbf{x}_1^\top\mathbf{x}_1 & \cdots & \mathbf{x}_1^\top\mathbf{x}_p \ \vdots & \vdots & \ddots & \vdots \ \mathbf{x}_p^\top\mathbf{x}_0 & \mathbf{x}_p^\top\mathbf{x}_1 & \cdots & \mathbf{x}_p^\top\mathbf{x}_p \end{bmatrix}$$

where $\mathbf{x}_j$ denotes column $j$ of the design matrix (treating columns as vectors).

Properties of the Gram matrix:

1. Symmetric: $(\mathbf{X}^\top\mathbf{X})^\top = \mathbf{X}^\top\mathbf{X}$

2. Positive semi-definite: For any vector $\mathbf{v}$, $\mathbf{v}^\top(\mathbf{X}^\top\mathbf{X})\mathbf{v} = |\mathbf{X}\mathbf{v}|^2 \geq 0$

3. Captures inner products: Element $(j, k)$ is the inner product $\langle \mathbf{x}j, \mathbf{x}k \rangle = \sum{i=1}^{n} x{ij} x_{ik}$

4. Diagonal entries: $(\mathbf{X}^\top\mathbf{X})_{jj} = |\mathbf{x}_j|^2$ = sum of squares of column $j$