Linear Regression

Learning Objectives

After reading this page, you should be able to:

1. Derive the vectorized gradient for linear regression variants.
2. Derive the direct solution for linear regression variants.
3. Explain the advantages and limitations of direct solution and gradient descent.

1 Introduction

So far, we have been discussing supervised learning. The two models we have discussed are k-nearest neighbours and decision trees.

As with any supervised learning task, we are given a training set consisting of inputs \(\mathbf{x} \in \mathbb{R}^D\) and the corresponding target labels \(\mathbf{t}\). We have set aside a separate validation set for tuning hyperparameters. We also talked about setting aside a separate labeled test set with the goal of determining how well our model generalizes to unseen data.

Today, we will start our discussion on linear models. In this section, we will discuss regression. In the next section, we will discuss classification.

Question: What do you recall is the main difference between regression and classification?

• For regression, the target is a real number.
• For classification, the target is a categorical variable.

Starting with this section, we will take a modular approach to machine learning. Initially, it may not be obvious what the “modules” are, but we will revisit this idea at the end of the section.

2 Regression Setup

Let’s recall the regression setup and become familiar with our notation. We have a training data set with \(N\) data points.

\[(\mathbf{x}^{(1)}, t^{(1)}), (\mathbf{x}^{(2)}, t^{(2)}), \dots, (\mathbf{x}^{(N)}, t^{(N)})\]

For the \(i\)-th data point \((\mathbf{x}^{(i)}, t^{(i)})\),

• \(\mathbf{x}^{(i)} \in \mathbb{R}^D\) is a \(D\)-dimensional feature vector, and
• \(t^{(i)}\) is a scalar target value for each feature vector.

We assume that there is an underlying function \(f: \mathbb{R}^D \rightarrow \mathbb{R}\) that maps from each \(D\)-dimensional feature vector to the corresponding scalar target value. We will learn a model for this function. Given a feature vector \(\mathbf{x}\), our model will produce a prediction \(y\). We hope that this prediction \(y\) is close to the corresponding target value \(t\).

There are many examples of regression problems. The table below shows some examples. The main distinction between regression and classification problems is the target value. The target value for a regression problem is a continuous value (or a real number). The target value for a classification problem is a discrete value or a categorical variable.

problem	feature vector \(\mathbf{x}^{(i)}\)	target scalar \(t^{(i)}\)
Predicting housing price	Square footage, location, # of bedrooms & bathrooms	Price of the house
Predicting rainfall	Temperature, humidity, wind, etc.	Amount of rainfall
Predicting company revenue	Previous sales data	Company’s revenue

3 Motivating Example

As a motivating example, consider the UTM pond. We want to predict the air temperature at the UTM pond given information from the UTM forest. To achieve this goal, we have collected various information regarding the air temperature, soil temperature, humidity and soil water content at the forest. Our goal is to predict the air temperature at the pond. Note that this problem has multiple real-valued features and a real-valued target.

4 Introducing the Linear Model

As we discussed in supervised learning, a hypothesis is a specific mapping from input to output. What makes learning possible is not considering all possible mappings, but restricting the set of hypotheses we are willing to consider. By choosing a particular model family, we impose assumptions about how inputs relate to outputs, which in turn defines and constrains the hypothesis space. These restrictions reduce an otherwise vast set of possible functions to a manageable space that we can search during learning.

You might have heard of the no free lunch theorem, which states that there is no universally best algorithm that works optimally across every problem. As a result, the effectiveness of any learning algorithm depends entirely on how well its assumptions align with the specific problem at hand. This is where inductive bias comes in. Your model’s inductive bias is precisely those assumptions it makes about the world. The theorem suggests that it is not only helpful but essential to make the right assumptions. Without them, you can’t learn anything useful from data.

4.1 Linear Model (Scalar Form)

For linear models, we restrict the mapping function \(f\) to have the specific form below.

\[\begin{align} y = f(\mathbf{x}) = w_1 x_1 + w_2 x_2 + \dots + w_D x_D + b = \sum_{j=1}^D w_j x_j + b \end{align}\]

Let’s break down the notation.

• \(y\) denotes the prediction output generated by our model.

• The weight vector \(\mathbf{w} = \begin{bmatrix} w_1 & w_2 & \dots & w_D \end{bmatrix}^\top\) controls how much each input feature contributes to the prediction.

• The bias \(b\) is an offset term that shifts the prediction.

• \(\mathbf{w}\) and \(b\) are the parameters of the model, which we will learn from data.

Why is this model “linear”? The function \(f(\mathbf{x})\) is linear in the features \(x_1, x_2, \dots, x_D\). Geometrically, this model describes very specific shapes: a straight line in 2D, a plane in 3D, and a hyperplane in higher dimensions.

Each weight \(w_j\) controls the magnitude and direction of how each feature \(x_j\) influences the prediction \(y\). A weight with a large magnitude means the model considers that feature particularly important for making predictions. As a result, small changes in this feature will have a big impact on the model’s output — the prediction is sensitive to that feature. We’ll revisit this idea when we discuss overfitting and see how such sensitivity can sometimes lead to problems.

Why do we include a bias term in the model? The bias \(b\) (also called the intercept) in a linear model allows us to shift the hyperplane. To see why this matters, consider a one-feature model: \(y = wx + b\). Here \(w\) controls the slope and \(b\) is the \(y\)-intercept. If we remove \(b\), the line is forced to pass through the origin \((0,0)\). This is a strong restriction and implies that the model cannot represent data that doesn’t go through the origin. By including \(b\), we allow the model to shift and align with the actual data, making it much more expressive. In other words, the bias term is essential for a flexible model and we always keep it.

4.2 Vectorizing the Linear Model

Before optimizing the model parameters, let’s rewrite the linear model in a more powerful form. So far, we have expressed the model in scalar form, with individual weights multiplying individual features. In practice, however, we typically represent these models using matrices and vectors. We use the term vectorization to refer to this process of converting an equation in scalar form to one in vectorized form.

Why Vectorization Matters

There are several advantages of expressing the model in a vectorized form. At the surface level, the vectorized form is mathematically elegant. It provides us with compact notation to handle thousands or millions of data points simultaneously.

However, the real motivation of vectorization is computational efficiency. Machine learning, especially deep learning, is computationally intensive. Training models often requires processing large datasets and performing millions of operations. GPUs (and even CPUs) are designed to perform many simple operations in parallel. Writing code using vectors and matrices allows us to take advantage of the GPU’s and CPU’s parallel computation capabilities.

Vectorization also allows us to use libraries such as NumPy and PyTorch, which are highly optimized for matrix and vector operations. In contrast, explicit Python loops are executed one step at a time and incur significant overhead. By avoiding loops, we allow these optimized routines to do the heavy lifting for us.

Computing the Prediction for One Example

Our ultimate goal is to compute predictions for all training examples using one vectorized equation with no summations. As a first step, let’s start by vectorizing the equation to compute the prediction for one training example.

Recall that we can write down the feature vector \(\mathbf{x} \in \mathbb{R}^{D \times 1}\) and the weight vector \(\mathbf{w} \in \mathbb{R}^{D \times 1}\) for one example as follows. We adopt the convention that every vector is defined as a column vector by default.

\[\begin{align} \mathbf{x} = \begin{bmatrix} x_1\\ x_2 \\ \vdots \\ x_D \\ \end{bmatrix}, \qquad \mathbf{w} = \begin{bmatrix} w_1 \\ w_2 \\ \vdots \\ w_D \\ \end{bmatrix} \end{align}\]

Having these vectors allows us to rewrite the linear model using a dot product (or an inner product) as shown below.

\[\begin{align} y = f(\mathbf{x}) & = \mathbf{w}^\top \mathbf{x} + b = \begin{bmatrix} w_1 & w_2 & \cdots & w_D \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_D \end{bmatrix} + b \\ & = {\mathbf{x}}^\top \mathbf{w} + b = \begin{bmatrix} x_1 & x_2 & \cdots & x_D \end{bmatrix} \begin{bmatrix} w_1 \\ w_2 \\ \vdots \\ w_D \\ \end{bmatrix} + b \end{align}\]

The two ways of writing the dot product (\(\mathbf{w}^\top \mathbf{x}\) and \({\mathbf{x}}^\top \mathbf{w}\)) are equivalent. The first expression (\(\mathbf{w}^\top \mathbf{x}\)) is more commonly used when we compute the prediction for one training example. However, the second expression (\({\mathbf{x}}^\top \mathbf{w}\)) will be useful in the later sections when we derive a vectorized equation for computing predictions for all training examples.

Treat the Bias as a Weight

In our linear model so far, notice that the bias \(b\) is treated separately from the weights \(\mathbf{w}\). While this reflects how we think about the model conceptually, it is mathematically inconvenient. The bias becomes a special case we have to handle separately in implementation. As a result, we often adopt a unified approach by treating the bias \(b\) as one of the weights.

We can do this by slightly changing our notation. Let’s introduce a dummy feature \(x_0 = 1\) (which is always \(1\)) and let its corresponding weight be the bias \(b\). (An alternative approach is to add a dummy feature at the end \(x_{D+1}=1\).) Then, we can incorporate the bias and the dummy feature into our weight and feature vectors respectively. Formally, we can rewrite the model as follows.

\[\begin{align} y & = \mathbf{w}^\top \mathbf{x} + b \cdot 1 = w_0 \cdot x_0 + \mathbf{w}^\top \mathbf{x} \\ & = \begin{bmatrix} w_0 & w_1 & \dots & w_D \end{bmatrix} \begin{bmatrix} x_0 \\ x_1\\ x_2 \\ \vdots \\ x_D \\ \end{bmatrix} = \mathbf{w'}^\top \mathbf{x'} \end{align}\]

In practice, we usually drop the primes for simplicity and continue writing the model as shown below. It is understood from context that \(\mathbf{x}\) includes the dummy feature \(x_0 = 1\) and \(\mathbf{w}\) includes the bias term \(w_0 = b\).

\[\begin{align} y = \mathbf{w}^\top \mathbf{x}. \end{align}\]

In this formulation, the bias \(b\) is no longer a parameter that requires special treatment — it is simply one of the weights that we want to optimize.

Computing Predictions for the Entire Training Set

So far, we can calculate the prediction for one example with the vectorized equation below. Our next step is to derive a vectorized equation for computing the predictions for all training examples.

\[\begin{align} y = \mathbf{w}^\top \mathbf{x} = \mathbf{x}^\top \mathbf{w}. \end{align}\]

To start, let’s define a prediction vector \(\mathbf{y} \in \mathbb{R}^{N \times 1}\) to represent the predictions for all training examples. As a sanity check, every example should have a prediction, so \(\mathbf{y}\) should be an \(N\)-dimensional vector. We will stack the individual predictions into a column vector, using the second expression above (\(\mathbf{x}^\top \mathbf{w}\)).

\[\begin{align} \mathbf{y} = \begin{bmatrix} y^{(1)} \\ y^{(2)} \\ \vdots \\ y^{(N)} \\ \end{bmatrix} = \begin{bmatrix} (\mathbf{x}^{(1)})^\top \mathbf{w} \\ (\mathbf{x}^{(2)})^\top \mathbf{w} \\ \vdots \\ (\mathbf{x}^{(N)})^\top \mathbf{w} \\ \end{bmatrix} \end{align}\]

Note that the \(i\)-th row of \(\mathbf{y}\) is a product of the two vectors below:

• \((\mathbf{x}^{(i)})^\top \in \mathbb{R}^{1 \times (D+1)}\) is a row vector corresponding to the features for the \(i\)-th training example, and

• \(\mathbf{w} \in \mathbb{R}^{(D+1) \times 1}\) is the column vector containing the weights.

We can factor out \(\mathbf{w}\) as a common term, as shown below. (Aside: To factor out a common term from a stack of products, that term must be positioned so it can be multiplied by each row. Since we’re stacking rows vertically, we need \(\mathbf{w}\) on the right as a column vector, so each row can multiply it.)

\[\begin{align} \mathbf{y} = \begin{bmatrix} y^{(1)} \\ y^{(2)} \\ \vdots \\ y^{(N)} \\ \end{bmatrix} = \begin{bmatrix} (\mathbf{x}^{(1)})^\top \mathbf{w} \\ (\mathbf{x}^{(2)})^\top \mathbf{w} \\ \vdots \\ (\mathbf{x}^{(N)})^\top \mathbf{w} \\ \end{bmatrix} = \begin{bmatrix} (\mathbf{x}^{(1)})^\top \\ (\mathbf{x}^{(2)})^\top \\ \vdots \\ (\mathbf{x}^{(N)})^\top \\ \end{bmatrix} \mathbf{w} \end{align}\]

What remains is to define a matrix containing the stacked feature vectors, which is the first term in the product above. We define the data matrix or design matrix \(\mathbf{X} \in \mathbb{R}^{N \times (D+1)}\). Each row represents the features for one example, and each column represents the values of one feature for all the examples.

\[\begin{align} \mathbf{X} = \begin{bmatrix} (\mathbf{x}^{(1)})^\top \\ (\mathbf{x}^{(2)})^\top \\ \vdots \\ (\mathbf{x}^{(N)})^\top \end{bmatrix} = \begin{bmatrix} x_0^{(1)} & x_1^{(1)} & x_2^{(1)} & \ldots & x_{D}^{(1)} \\ x_0^{(2)} & x_1^{(2)} & x_2^{(2)} & \ldots & x_{D}^{(2)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ x_0^{(N)} & x_1^{(N)} & x_2^{(N)} & \ldots & x_{D}^{(N)} \end{bmatrix} \end{align}\]

With this matrix, it is crucial to keep track of the indices: the superscripts represent the index of each example, and the subscripts represent the index of each feature. Notice that the data matrix already contains the dummy feature for each example. The final vectorized equation to compute the predictions for all the training examples is shown below, and we reiterate the definitions of \(\mathbf{y} \in \mathbb{R}^{N \times 1}, \mathbf{X} \in \mathbb{R}^{N \times (D+1)}, \mathbf{w} \in \mathbb{R}^{(D + 1)\times 1}\) for emphasis.

\[\begin{align} \mathbf{y} = \mathbf{X} \mathbf{w} \\ \end{align}\] \[\begin{align} \mathbf{X} = \begin{bmatrix} x_0^{(1)} & x_1^{(1)} & x_2^{(1)} & \ldots & x_{D}^{(1)} \\ x_0^{(2)} & x_1^{(2)} & x_2^{(2)} & \ldots & x_{D}^{(2)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ x_0^{(N)} & x_1^{(N)} & x_2^{(N)} & \ldots & x_{D}^{(N)} \end{bmatrix}, \quad \mathbf{w} = \begin{bmatrix} w_0 \\ w_1 \\ w_2 \\ \vdots \\ w_D \\ \end{bmatrix}, \quad \mathbf{y} = \begin{bmatrix} y^{(1)} \\ y^{(2)} \\ \vdots \\ y^{(N)} \\ \end{bmatrix} \end{align}\]

As a sanity check, note that it is valid to multiply \(\mathbf{X}\) and \(\mathbf{w}\) together since the second dimension of \(\mathbf{X}\) matches the first dimension of \(\mathbf{w}\), producing a \(\mathbb{R}^{N \times 1}\) column vector.

4.3 The Loss Function

In the linear model, the weights \(\mathbf{w}\) (including biases) are model parameters. You have already seen examples of hyperparameters, such as \(k\) in k-Nearest-Neighbours and maximum tree depth in decision trees. We define model parameter and hyperparameter formally below.

Definition: A model parameter is a value that the learning algorithm automatically learns from the training data.

Definition: A hyperparameter is a setting chosen before training and is typically tuned using the validation data.

Recall that the hypothesis space consists of all possible functions of our chosen form. For our linear regression model, each choice of \(\mathbf{w}\) will give us a different hyperplane. There are infinitely many possible parameter combinations, so infinitely many hypotheses to choose from.

How do we find the best parameters from infinite possibilities? We need to mathematically quantify how badly a model’s prediction matches the target, using a loss function. You have seen an example of a loss function before. Entropy/uncertainty is a loss function we optimize when building a decision tree.

A loss function quantifies how bad the hypothesis is for an example if our model predicts \(y\) but the target is \(t\). In other words, the loss function evaluates how badly our model is doing at predicting the target. If the prediction gets far away from the target, the loss increases.

Definition: The loss function measures the discrepancy between a model’s prediction and the ground-truth target value for one example.

For regression, we typically use the squared error loss. This loss function takes the difference between the prediction and the target, squares the result, and multiplies by \(1/2\).

\[\begin{align} \mathcal{L}(y,t) = \frac{1}{2} (y - t)^2 \end{align}\]

We also call the difference \(y - t\) the residual, and we want to minimize the residual for each example. Having the \(1/2\) factor simplifies derivative calculations since the power of 2 cancels out with \(1/2\).

We often choose squared error as a loss function because it corresponds to a natural statistical assumption: the prediction errors are normally distributed around the true values. Under this assumption, minimizing squared error arises directly from maximum likelihood estimation. While we do not emphasize statistical theory in this course, this connection provides a principled justification for squared error as a sensible objective.

Squared error also has an intuitive geometric advantage over alternatives such as absolute error. With absolute error, many different prediction functions can achieve the same total loss. For example, a line can be vertically shifted as long as there are equal numbers of points above and below it. Squared error, however, penalizes larger deviations much more heavily, which discourages such ambiguity and pulls the solution toward a single line that runs through the “middle” of the data. This behavior aligns well with our intuitive notion of what a good fit should look like.

4.4 The Cost Function

The loss function measures how bad a single prediction is. But we have multiple training examples. To aggregate the losses for all examples, we can compute the average loss.

We can define a cost function that is the average loss across all training data points. The cost function calculates the loss for each training point, sums them together, and takes the average over the \(N\) examples, as shown below.

\[\begin{align} \mathcal{E}(\mathbf{w}) = \frac{1}{N} \sum_{i = 1}^N \mathcal{L}(y^{(i)}, t^{(i)}) \end{align}\]

Definition: The cost function measures the average loss across all training examples, quantifying the overall performance of the model on the entire dataset.

The expression above emphasizes that cost is a function of the model parameters \(\mathbf{w}\). This might look different from how we defined the loss in terms of the prediction and the target. Ultimately, both loss and cost are functions of the model parameters, and the expression above makes this connection explicit. This reminder will be useful once we start taking derivatives of the loss/cost function.

In machine learning literature, “loss” and “cost” are often used interchangeably. However, by their definitions, loss refers to the error for a single training example, while cost refers to the average loss across the entire training dataset.

If we plug in the squared error loss into the cost function, we will get the mean squared error (MSE) cost function.

\[\begin{align} \mathcal{E}(\mathbf{w}) &= \frac{1}{N} \sum_{i = 1}^N \mathcal{L}(y^{(i)}, t^{(i)}) \\ & = \frac{1}{2N} \sum_{i = 1}^N (y^{(i)} - t^{(i)})^2 \\ & = \frac{1}{2N} \sum_{i = 1}^N (\mathbf{w}^\top \mathbf{x}^{(i)} - t^{(i)})^2 \label{mse_plugin} \end{align}\]

In equation \(\eqref{mse_plugin}\), we plug in \(y^{(i)} = \mathbf{w}^\top \mathbf{x}^{(i)}\) to highlight the dependence of the cost function on the model parameters \(\mathbf{w}\).

4.5 Vectorizing MSE Cost

To find optimal model parameters that minimize the MSE cost function, it would be useful to express the function in a vectorized form. We will show you a few equivalent vectorized expressions of the MSE cost function.

To start, we can write the MSE cost as an inner product of the residuals, as shown below.

\[\begin{align} \mathcal{E}(\mathbf{w}) = \frac{1}{2N} \sum_{i = 1}^N (y^{(i)} - t^{(i)})^2 = \frac{1}{2N} (\mathbf{y} - \mathbf{t})^\top (\mathbf{y} - \mathbf{t}) \end{align}\]

Note that the prediction vector \(\mathbf{y} \in \mathbb{R}^{N \times 1}\) and the target vector \(\mathbf{t} \in \mathbb{R}^{N \times 1}\) have the same dimensions.

Question: Explain why we can convert the summation on the left to the dot product on the right.

\[\frac{1}{2N} \sum_{i = 1}^N (y^{(i)} - t^{(i)})^2 = \frac{1}{2N} (\mathbf{y} - \mathbf{t})^\top (\mathbf{y} - \mathbf{t})\]

Answer: Let’s ignore the constant at the front and focus on transforming the summation. We start with the expression below.

\[\begin{equation} \sum_{i = 1}^N (y^{(i)} - t^{(i)}) (y^{(i)} - t^{(i)}) \end{equation}\]

To make the expression more compact, let’s define a vector \(\mathbf{r} \in \mathbb{R}^{N \times 1}\) where its \(i\)-th component is equal to the \(i\)-th residual.

\[\begin{equation} r^{(i)} = y^{(i)} - t^{(i)} \end{equation}\]

Then, the original expression becomes:

\[\begin{equation} \sum_{i = 1}^N r^{(i)} r^{(i)} = \begin{bmatrix} r^{(1)} & \dots & r^{(N)} \end{bmatrix} \begin{bmatrix} r^{(1)} \\ \vdots \\ r^{(N)} \end{bmatrix} = \mathbf{r}^\top \mathbf{r} \end{equation}\]

From the previous section, we can write the prediction vector as a product of the data matrix and the weight vector, leading to the following equivalent expression.

\[\begin{align} \mathcal{E}(\mathbf{w}) = \frac{1}{2N} (\mathbf{y} - \mathbf{t})^\top (\mathbf{y} - \mathbf{t}) = \frac{1}{2N} (\mathbf{X} \mathbf{w} - \mathbf{t})^\top (\mathbf{X} \mathbf{w} - \mathbf{t}) \end{align}\]

Finally, we can write the MSE cost function using the \(L^2\) norm notation below.

\[\begin{align} \mathcal{E}(\mathbf{w}) = \frac{1}{2N} \sum_{i = 1}^N (y^{(i)} - t^{(i)})^2 = \frac{1}{2N} \lVert \mathbf{X} \mathbf{w} - \mathbf{t} \rVert_2 ^2 \\ \end{align}\]

In the equation above, the superscript \(2\) denotes exponentiation and the \(\lVert \rVert_2\) notation with a subscript of \(2\) denotes the \(L^2\) norm. \(L^2\) norm is also called the Euclidean norm since it measures the length or magnitude of a vector using Euclidean distance. For a vector \(\mathbf{v} \in \mathbb{R}^{n \times 1}\), its \(L^2\) norm is defined as the square root of the sum of the squares of its components.

\[\begin{align} \lVert \mathbf{v} \rVert_2 = \sqrt{v_1^2 + v_2^2 + \dots + v_n^2} = \sqrt{\sum_{i=1}^n v_i^2} \end{align}\]

Therefore, we can write the squared \(L^2\) norm of \(\mathbf{v}\) in two equivalent ways as follows.

\[\begin{align} \lVert \mathbf{v} \rVert_2 ^2 = \sum_{i=1}^n v_i^2 = \begin{bmatrix} v_1 & \dots & v_n \end{bmatrix} \begin{bmatrix} v_1 \\ \vdots \\ v_n \end{bmatrix} = \mathbf{v}^\top \mathbf{v} \end{align}\]

Once you plug in \(\mathbf{v} = \mathbf{X} \mathbf{w} - \mathbf{t}\), you can see that the dot product \((\mathbf{X} \mathbf{w} - \mathbf{t})^\top (\mathbf{X} \mathbf{w} - \mathbf{t})\) is equivalent to \(\lVert \mathbf{X} \mathbf{w} - \mathbf{t} \rVert_2 ^2\).

5 Finding Optimal Parameters

Our goal is to find parameters \(\mathbf{w}\) that minimize the cost function \(\mathcal{E}(\mathbf{w})\).

\[\begin{align} \mathbf{w}^* = \arg \min_{\mathbf{w}} \mathcal{E}(\mathbf{w}) \end{align}\]

There are two main approaches to solve this optimization problem, solving for a direct solution and performing gradient descent.

We can use calculus to find the minimum of the cost function. The process involves three steps.

1. Compute the gradient (the partial derivative of the cost with respect to each parameter).

2. Set the gradient equal to zero.

3. Solve the resulting system of equations for the parameters.

This approach gives a closed-form solution that contains the optimal parameters in one shot.

5.1 Deriving Gradient for the Cost

We are going to take a series of steps to derive the gradient vector for linear regression. We start with the cost function below, defined as the average loss over the training set.

\[\begin{align} \mathcal{E}(\mathbf{w}) = \frac{1}{N} \sum_{i = 1}^N \mathcal{L}^{(i)}(\mathbf{w}) \end{align}\]

Here \(\mathcal{L}^{(i)}(\mathbf{w}) = \mathcal{L}(y^{(i)}, t^{(i)})\) denotes the loss for the \(i\)-th training example. Given this, to derive the gradient for the cost function, it is sufficient to derive the gradient for the loss function.

\[\begin{align} \nabla_{\mathbf{w}} \mathcal{E}(\mathbf{w}) = \begin{bmatrix} \displaystyle \frac{\partial \mathcal{E}}{\partial w_0} \\ \vdots \\ \displaystyle \frac{\partial \mathcal{E}}{\partial w_D} \end{bmatrix} = \frac{1}{N} \sum_{i = 1}^N \nabla_{\mathbf{w}} \mathcal{L}^{(i)}(\mathbf{w}) \end{align}\]

Next, let’s derive the derivative of the loss with respect to each weight, \(\frac{\partial \mathcal{L}^{(i)}(\mathbf{w})}{\partial w_j}, j = 0, 1, \dots, D\). Recall that the loss function is defined below.

\[\begin{align} \mathcal{L}^{(i)}(\mathbf{w}) = \mathcal{L}(y^{(i)}, t^{(i)}) = \frac{1}{2} (y^{(i)} - t^{(i)})^2 = \frac{1}{2} (\mathbf{w}^\top \mathbf{x}^{(i)} - t^{(i)})^2 \end{align}\]

We will compute the derivative of the loss function with respect to each weight using the chain rule.

\[\begin{align} & \frac{\partial \mathcal{L}^{(i)}}{\partial y^{(i)}} = \frac{\partial }{\partial y^{(i)}} \frac{1}{2} (y^{(i)} - t^{(i)})^2 = (y^{(i)} - t^{(i)}) = \left( \mathbf{w}^\top \mathbf{x}^{(i)} - t^{(i)} \right)\\ & \frac{\partial y^{(i)}}{\partial w_j} = \frac{\partial }{\partial w_j} (w_0 x_0^{(i)} + \dots + w_j x_j^{(i)} + \dots + w_D x_D^{(i)}) = x_j^{(i)} \end{align}\]

Therefore,

\[\begin{align} & \frac{\partial \mathcal{L}^{(i)} (\mathbf{w})}{\partial w_j} = \frac{\partial \mathcal{L}^{(i)}}{\partial y^{(i)}} \frac{\partial y^{(i)}}{\partial w_j} = \left( \mathbf{w}^\top \mathbf{x}^{(i)} - t^{(i)} \right) x_j^{(i)} \end{align}\]

Given the derivative of the loss with respect to each weight, we can write down the gradient vector for the loss function below.

\[\begin{align} \nabla_{\mathbf{w}} \mathcal{L}^{(i)}(\mathbf{w}) = \begin{bmatrix} \displaystyle \frac{\partial \mathcal{L}^{(i)}}{\partial w_0} \\ \vdots \\ \displaystyle \frac{\partial \mathcal{L}^{(i)}}{\partial w_D} \end{bmatrix} & = \begin{bmatrix} \left( \mathbf{w}^\top \mathbf{x}^{(i)} - t^{(i)} \right) x_0^{(i)} \\ \vdots \\ \left( \mathbf{w}^\top \mathbf{x}^{(i)} - t^{(i)} \right) x_D^{(i)} \end{bmatrix} \\ &= \left( \mathbf{w}^\top \mathbf{x}^{(i)} - t^{(i)} \right) \begin{bmatrix} x_0^{(i)} \\ \vdots \\ x_D^{(i)} \end{bmatrix} \\ & = \left( \mathbf{w}^\top \mathbf{x}^{(i)} - t^{(i)} \right) \mathbf{x}^{(i)} \end{align}\]

Next, we can write down the gradient vector for the cost function by adding back the average to the gradient vector for the loss function.

\[\begin{align} \nabla_{\mathbf{w}} \mathcal{E} (\mathbf{w}) = \frac{1}{N} \sum_{i = 1}^N \nabla_{\mathbf{w}} \mathcal{L}^{(i)}(\mathbf{w}) = \frac{1}{N} \sum_{i = 1}^N \left( \mathbf{w}^\top \mathbf{x}^{(i)} - t^{(i)} \right) \mathbf{x}^{(i)} \end{align}\]

In the next section, we will fully vectorize the gradient vector by eliminating the summation. The resulting vectorized gradient vector is shown below.

\[\begin{align} \nabla_{\mathbf{w}} \mathcal{E}(\mathbf{w}) = \frac{1}{N} \sum_{i = 1}^N \left( \mathbf{w}^\top \mathbf{x}^{(i)} - t^{(i)} \right) \mathbf{x}^{(i)} = \frac{1}{N} \mathbf{X}^\top \left( \mathbf{X} \mathbf{w} - \mathbf{t} \right) \end{align}\]

5.2 Vectorizing Gradient for the Cost

We will show that the following equation is valid.

\[ \begin{equation} \label{eq:test} \begin{aligned} \nabla_{\mathbf{w}} \mathcal{E} (\mathbf{w}) = \frac{1}{N} \sum_{i = 1}^N \left( \mathbf{w}^\top \mathbf{x}^{(i)} - t^{(i)} \right) \mathbf{x}^{(i)} \label{step_3} = \frac{1}{N} \mathbf{X}^\top \left( \mathbf{X} \mathbf{w} - \mathbf{t} \right) \end{aligned} \end{equation} \]

Recall the following notation: \(N\) denotes the number of training examples. \(D\) denotes the number of features. \(\mathbf{X} \in \mathbb{R}^{N \times (D+1)}\) denotes the data matrix. \(\mathbf{t} \in \mathbb{R}^{N \times 1}\) denotes the target vector. \(\mathbf{w} \in \mathbb{R}^{(D+1) \times 1}\) denotes the weight vector (including the bias).

Equation \(\eqref{eq:test}\) has three operations that we need to vectorize:

1. The operation between the residuals \(\left( \mathbf{w}^\top \mathbf{x}^{(i)} - t^{(i)} \right)\) and the features \(\mathbf{x}^{(i)}\).

2. The operation between the predictions \(\mathbf{w}^\top \mathbf{x}^{(i)}\) and the targets \(t^{(i)}\).

3. The operation between the weights \(\mathbf{w}\) and the features \(\mathbf{x}^{(i)}\).

Step 1: Vectorizing operation between residuals and features

To simplify this expression temporarily, let’s define the vector \(\mathbf{r} \in \mathbb{R}^{N \times 1}\) with its \(i\)-th component equal to the \(i\)-th residual.

\[\begin{align} \mathbf{r} = \begin{bmatrix} r^{(1)} \\ r^{(2)} \\ \vdots \\r^{(N)} \end{bmatrix}, \quad r^{(i)} = \mathbf{w}^\top \mathbf{x}^{(i)} - t^{(i)} \end{align}\]

With these definitions, the gradient vector \(\nabla_{\mathbf{w}} \mathcal{E}(\mathbf{w})\) becomes the following. \[\begin{align} \nabla_{\mathbf{w}} \mathcal{E}(\mathbf{w}) = \frac{1}{N} \sum_{i=1}^N \left(\mathbf{w}^\top \mathbf{x}^{(i)} - t^{(i)} \right) \mathbf{x}^{(i)} = \frac{1}{N} \sum_{i=1}^N r^{(i)} \mathbf{x}^{(i)}. \end{align}\]

Note that \(r^{(i)} \in \mathbb{R}\) is a scalar and \(\mathbf{x}^{(i)} \in \mathbb{R}^{(D+1) \times 1}\) is a vector. The summation over \(i = 1, \dots, N\) adds a dimension to both components. \(r^{(i)} \in \mathbb{R}\) becomes a vector \(\mathbf{r} \in \mathbb{R}^{N \times 1}\) (defined as a column vector by our convention). \(\mathbf{x}^{(i)} \in \mathbb{R}^{(D+1) \times 1}\) becomes a matrix with one dimension over the \((D+1)\) features and one dimension over the \(N\) examples. We will define \(\mathbf{X} \in \mathbb{R}^{N \times (D+1)}\) to be consistent with our previous definition of the data matrix.

Our goal is to multiply the residual vector \(\mathbf{r} \in \mathbb{R}^{N \times 1}\) and the data matrix \(\mathbf{X} \in \mathbb{R}^{N \times (D+1)}\) to produce the gradient vector \(\nabla_{\mathbf{w}} \mathcal{E}(\mathbf{w}) \in \mathbb{R}^{(D+1) \times 1}\). Recall that the gradient vector should have one derivative for each of the \((D+1)\) weights.

The matrix/vector multiplication requires the inner dimensions to match. Observing that the dimension over \(N\) should disappear after the product, one of these two products should work: \(\mathbf{r}^\top \mathbf{X}\) or \(\mathbf{X}^\top \mathbf{r}\).

By our convention, the resulting gradient vector should be a column vector. \(\mathbf{r}^\top \mathbf{X}\) produces a row vector, whereas \(\mathbf{X}^\top \mathbf{r}\) produces a column vector. Therefore, the vectorized equation is as follows.

\[\begin{align} \nabla_{\mathbf{w}} \mathcal{E}(\mathbf{w}) = \frac{1}{N} \sum_{i=1}^N r^{(i)} \mathbf{x}^{(i)} = \frac{1}{N} \mathbf{X}^\top \mathbf{r} \end{align}\]

Step 2: Vectorizing operation between predictions and targets

Next, we are going to vectorize the operation between the predictions and the targets inside the residual vector. Recall that each component of the residual vector is defined below. For this step, we will write the first term as the prediction \(y^{(i)}\) so that we don’t yet worry about how to vectorize the operations for computing the prediction.

\[\begin{align} r^{(i)} = \mathbf{w}^\top \mathbf{x}^{(i)} - t^{(i)} = y^{(i)} - t^{(i)} \end{align}\]

We can write the residual vector \(\mathbf{r}\) as follows, where \(\mathbf{y} \in \mathbb{R}^{N \times 1}\) is the prediction vector and \(\mathbf{t} \in \mathbb{R}^{N \times 1}\) is the target vector.

\[\begin{align} \mathbf{r} = \begin{bmatrix} r^{(1)} \\ r^{(2)} \\ \vdots \\ r^{(N)} \end{bmatrix} = \begin{bmatrix} y^{(1)} - t^{(1)} \\ y^{(2)} - t^{(2)} \\ \vdots \\ y^{(N)} - t^{(N)} \end{bmatrix} = \begin{bmatrix} y^{(1)} \\ y^{(2)} \\ \vdots \\ y^{(N)} \end{bmatrix} - \begin{bmatrix} t^{(1)} \\ t^{(2)} \\ \vdots \\ t^{(N)} \end{bmatrix} = \mathbf{y} - \mathbf{t} \end{align}\]

The gradient vector computation becomes the following. \[\begin{align} \nabla_{\mathbf{w}} \mathcal{E}(\mathbf{w}) = \frac{1}{N} \mathbf{X}^\top \mathbf{r} = \frac{1}{N} \mathbf{X}^\top (\mathbf{y} - \mathbf{t}) \end{align}\]

Step 3: Vectorizing operation between weights and features

What remains is to write the prediction vector \(\mathbf{y}\) in terms of the data matrix \(\mathbf{X}\) and the weight vector \(\mathbf{w}\). We will use the result of the detailed derivation in Section 4.2, shown below.

\[\begin{align} \mathbf{y} = \mathbf{X} \mathbf{w} \end{align}\]

Plugging this into the gradient vector, we have the final vectorized expression below.

\[\begin{align} \nabla_{\mathbf{w}} \mathcal{E}(\mathbf{w}) = \frac{1}{N} \mathbf{X}^\top (\mathbf{y} - \mathbf{t}) = \frac{1}{N} \mathbf{X}^\top (\mathbf{X} \mathbf{w} - \mathbf{t}) \end{align}\]

5.3 Solving for Direct Solution

Now that we have the gradient of the cost function, we will set it to be zero (which means every component is zero) and solve for the optimal weights.

\[\begin{align} \nabla_{\mathbf{w}} \mathcal{E} (\mathbf{w}) = \frac{1}{N} \mathbf{X}^\top (\mathbf{X} \mathbf{w} - \mathbf{t}) = \mathbf{0} \\ \Rightarrow \mathbf{X}^\top \mathbf{X} \mathbf{w} - \mathbf{X}^\top \mathbf{t} = \mathbf{0} \\ \Rightarrow \mathbf{X}^\top \mathbf{X} \mathbf{w} = \mathbf{X}^\top \mathbf{t} \\ \Rightarrow \mathbf{w} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{t} \end{align}\]

Note that in the last step, since \(\mathbf{X}^\top \mathbf{X}\) is a matrix, we cannot “divide” by \(\mathbf{X}^\top \mathbf{X}\). Instead, we need to multiply both sides with its inverse.

6 Direct Solution Properties

Now that we know how to solve for the optimal weights of linear regression using the direct solution approach, let’s discuss the advantages and limitations of this approach.

The direct (closed-form) solution approach works well in simple linear regression for several reasons. A closed-form solution exists for the optimal weights of linear regression. Moreover, the optimal weights are unique because the squared-loss objective is convex with a single global minimum. Therefore, the direct solution approach allows us to compute the unique optimal weights directly without relying on iterative algorithms.

However, the direct solution approach also has important limitations. First, it applies only to a small set of problems. Basic linear regression with squared loss is one of the few cases where a closed-form solution exists. Second, it becomes computationally expensive with many features. The direct solution requires computing a matrix inversion \((\mathbf{X}^\top \mathbf{X})^{-1}\), which is an \(O(D^3)\) operation. Third, it does not generalize to most of the complex machine learning models we will study later.

For these reasons, we will discuss a more general and scalable optimization approach on the next page, called gradient descent.

7 A Modular Approach to Machine Learning

On this page, we saw that each machine learning model can be broken down into three modular components:

1. A model architecture specifies the form of the relationship between inputs and outputs. An example is representing the prediction as a linear combination of features.

2. A loss (or cost) function measures how well the model fits the data by assigning a numerical penalty to prediction errors. An example is the squared error for regression.

3. An optimization algorithm determines how the model’s parameters are adjusted in order to reduce the loss. The direct solution approach is an example.

This modular approach is useful because it clearly distinguishes three parts: (1) what we assume about the data, (2) how we measure errors, and (3) how we learn from the data. By keeping some components fixed and changing others, we can construct many different machine learning algorithms within a single, unified framework. This modular way of thinking will carry forward as we move from simple models like linear regression to more complex learning methods.

8 Summary

This page is devoted to introducing the regression setup, the linear regression model for solving a regression problem, and the derivation of the closed-form solution for the optimal weights. Along the way, we introduced several important concepts: how regression differs from classification, what a model family and an inductive bias are, how model parameters differ from hyperparameters, and how a loss function differs from a cost function.

Much of our effort went into vectorization. We rewrote the model, the cost function, and finally the gradient, until each one became a matrix or vector expression with no summations left in it. Practically speaking, vectorized code lets optimized libraries and parallel hardware do the heavy lifting for us. Conceptually, vectorization connects to Fundamental Idea #4: ML Describes Geometric Processes. Once the entire training set is one data matrix, our model becomes a transformation acting on the whole data set rather than a rule applying to one data point at a time.

This page presents the first example of Fundamental Idea #1: Learning is Optimization. We chose a model family, defined an objective that states what a good fit means, and solved for the parameters that optimize it. This is also where the modular view introduced above pays off. A model, a loss, and an optimizer are separate choices. On the next page, we will keep the first two components and replace the direct solution with gradient descent.