Multiple Linear Regression As A Matrix Problem

From Several Predictors to One Matrix Equation

In a multiple linear regression model, each observation has a response and measurements of several predictors. The model assigns one coefficient to each predictor, along with a separate intercept. Writing the observations together turns the problem into a single matrix equation: the design matrix records the predictor values, the coefficient vector records the unknowns, and the response vector records the data to be fitted.

Suppose there are \(n\) observations and \(p\) predictors. For observation \(i\), write its response as \(y_i\), and its predictor values as \(x_{i1},\dots,x_{ip}\). The coefficient vector has \(p+1\) entries: one intercept coefficient and one coefficient for each predictor. All vectors in this tutorial are columns in \(\mathbb{R}^n\) or \(\mathbb{R}^{p+1}\), as appropriate. The intercept is represented by adding a first column of ones to the design matrix.

The matrix equation \(\mathbf{X}\mathbf{b}=\mathbf{y}\) asks whether one choice of coefficients fits every observation exactly. Usually there are more observations than coefficients, and no exact solution exists. Least squares instead chooses \(\mathbf{b}\) to minimize the squared error \(\|\mathbf{y}-\mathbf{X}\mathbf{b}\|^2\). The vector \(\mathbf{r}=\mathbf{y}-\mathbf{X}\mathbf{b}\) is the residual vector; here it is used to express the fitting condition, rather than to interpret the individual coefficients.

The Normal Equations

For any candidate coefficient vector \(\mathbf{b}\), the residual is \(\mathbf{y}-\mathbf{X}\mathbf{b}\). A least-squares solution is a vector that minimizes the squared norm of this residual. The normal equations give a system of equations that every such minimizer satisfies—and, in fact, characterizes the minimizers.

Proof. Put \(\mathbf{r}=\mathbf{y}-\mathbf{X}\mathbf{b}\). First suppose \(\mathbf{X}^{\top}\mathbf{r}=\mathbf{0}\), which is equivalent to the stated normal equations. For any other candidate \(\mathbf{c}\), let \(\mathbf{d}=\mathbf{c}-\mathbf{b}\). Then \(\mathbf{y}-\mathbf{X}\mathbf{c}=\mathbf{r}-\mathbf{X}\mathbf{d}\), and expansion of the squared norm gives

The middle equality holds because \(\mathbf{r}^{\top}\mathbf{X}\mathbf{d}=(\mathbf{X}^{\top}\mathbf{r})^{\top}\mathbf{d}=0\). Thus \(\mathbf{b}\) minimizes the squared residual norm. Conversely, suppose \(\mathbf{b}\) is a minimizer. For any \(\mathbf{h}\in\mathbb{R}^q\), define \(a=\mathbf{r}^{\top}\mathbf{X}\mathbf{h}\). If \(\mathbf{X}\mathbf{h}=\mathbf{0}\), then \(a=0\). If \(\mathbf{X}\mathbf{h}\ne\mathbf{0}\) and \(a\ne0\), choose \(t=a/\|\mathbf{X}\mathbf{h}\|^2\). The candidate \(\mathbf{b}+t\mathbf{h}\) has squared residual norm

This contradicts minimality, so \(a=0\) for every \(\mathbf{h}\). Taking \(\mathbf{h}\) to be each standard basis vector of \(\mathbb{R}^q\) shows that every entry of \(\mathbf{X}^{\top}\mathbf{r}\) is zero. Hence \(\mathbf{X}^{\top}(\mathbf{y}-\mathbf{X}\mathbf{b})=\mathbf{0}\), which rearranges to the normal equations. \(\square\)

The matrix \(\mathbf{X}^{\top}\mathbf{X}\) is \(q\times q\), and \(\mathbf{X}^{\top}\mathbf{y}\) is a column with \(q\) entries, so the equation is dimensionally consistent. In a regression with an intercept, the first column of \(\mathbf{X}\) consists of ones. The first normal equation therefore involves the sum of the residual entries; the other equations involve their products with the corresponding predictor columns.

Solving a Multiple-Regression Example

When the Coefficients Are Not Unique

Solving the normal equations does not always produce one coefficient vector. If a predictor column is a linear combination of other columns, different coefficient vectors can produce the same fitted vector. The relevant condition is linear independence of all columns of the design matrix, including the intercept column.

Proof. A least-squares solution exists for every \(\mathbf{X}\) and \(\mathbf{y}\), by the existence result for least-squares problems. Suppose first that the columns of \(\mathbf{X}\) are linearly independent, and let \(\mathbf{b}\) and \(\mathbf{c}\) both be least-squares solutions. By the normal equations, \(\mathbf{X}^{\top}\mathbf{X}\mathbf{b}=\mathbf{X}^{\top}\mathbf{y}=\mathbf{X}^{\top}\mathbf{X}\mathbf{c}\). Set \(\mathbf{d}=\mathbf{c}-\mathbf{b}\); then \(\mathbf{X}^{\top}\mathbf{X}\mathbf{d}=\mathbf{0}\). Multiplying on the left by \(\mathbf{d}^{\top}\) gives \(\|\mathbf{X}\mathbf{d}\|^2=0\), so \(\mathbf{X}\mathbf{d}=\mathbf{0}\). Independence of the columns means the only solution of this homogeneous equation is \(\mathbf{d}=\mathbf{0}\). Therefore \(\mathbf{c}=\mathbf{b}\). Conversely, if the columns are dependent, there is a nonzero vector \(\mathbf{d}\) with \(\mathbf{X}\mathbf{d}=\mathbf{0}\). For any least-squares solution \(\mathbf{b}\), the vector \(\mathbf{b}+\mathbf{d}\) gives \(\mathbf{X}(\mathbf{b}+\mathbf{d})=\mathbf{X}\mathbf{b}\), so it gives the same residual and is also a minimizer. Since \(\mathbf{d}\ne\mathbf{0}\), it differs from \(\mathbf{b}\), and the solution is not unique. \(\square\)

Common Mistake: Leaving Out the Intercept

A model with an intercept must include a column of ones. Omitting it changes the set of vectors that the model can fit; the result is not simply a shorter way to write the same model. For the four-observation example above, suppose the intercept column is mistakenly removed. The resulting two-column matrix has columns \(x\) and \(z\), so its normal equations use \(\begin{bmatrix}6&0\\0&4\end{bmatrix}\) and right-hand side \(\begin{bmatrix}14\\-4\end{bmatrix}\). They give the coefficients \((7/3,-1)^{\top}\), not \((2,-1)^{\top}\).

This calculation is internally consistent for a model forced to have zero intercept, but it does not solve the stated model that includes an intercept. For example, the resulting fitted vector is \((-10/3,1,10/3,11/3)^{\top}\), whose residual is \((1/3,4,-7/3,4/3)^{\top}\). The residual entries sum to \(10/3\ne0\), so the residual is not orthogonal to the omitted column of ones. The correction is to restore the intercept column before forming \(\mathbf{X}^{\top}\mathbf{X}\) and \(\mathbf{X}^{\top}\mathbf{y}\).

In practice, the matrix setup determines the model being fitted: rows must keep the observations aligned, columns must represent the intended predictors, and the intercept must be included when required. Once the design matrix and response vector are assembled, the normal equations provide a systematic route to a least-squares solution. If the columns are independent, that solution is unique; if they are dependent, the fitted vector may still be determined even though its coefficient representation is not.