Ordinary Least Squares

We’ve now spent a number of lectures exploring how to build effective models – we introduced the SLR and constant models, selected cost functions to suit our modeling task, and applied transformations to improve the linear fit.

Throughout all of this, we considered models of one feature ($\hat{y}_i = \theta_0 + \theta_1 x_i$) or zero features ($\hat{y}_i = \theta_0$). As data scientists, we usually have access to datasets containing many features. To make the best models we can, it will be beneficial to consider all of the variables available to us as inputs to a model, rather than just one. In today’s lecture, we’ll introduce multiple linear regression as a framework to incorporate multiple features into a model. We will also learn how to accelerate the modeling process – specifically, we’ll see how linear algebra offers us a powerful set of tools for understanding model performance.

OLS Problem Formulation¶

Multiple Linear Regression¶

Multiple linear regression is an extension of simple linear regression that adds additional features to the model. The multiple linear regression model takes the form:

Our predicted value of $y$, $\hat{y}$, is a linear combination of the single observations (features), $x_i$, and the parameters, $\theta_i$.

As an example, you could predict cats’ ages $\hat{y}$ as a linear model of two features: length $x_1$ and weight $x_2$. This would give us the equation

where $\theta_0$ is the intercept, $\theta_1$ is the parameter for length, and $\theta_2$ is the parameter for weight.

We can explore this idea further by looking at a dataset containing aggregate per-player data from the 2018-19 NBA season, downloaded from Kaggle.

Let’s say we are interested in predicting the number of points (PTS) an athlete will score in a basketball game this season.

Suppose we want to fit a linear model by using some characteristics, or features of a player. Specifically, we’ll focus on field goals, assists, and 3-point attempts.

• FG, the average number of (2-point) field goals per game

• AST, the average number of assists per game

• 3PA, the average number of 3-point field goals attempted per game

Because we are now dealing with many parameter values, we’ve collected them all into a parameter vector with dimensions $(p+1) \times 1$ to keep things tidy. Remember that $p$ represents the number of features we have (in this case, 3).

We are working with two vectors here: a row vector representing the observed data, and a column vector containing the model parameters. The multiple linear regression model is equivalent to the dot (scalar) product of the observation vector and parameter vector.

Notice that we have inserted 1 as the first value in the observation vector. When the dot product is computed, this 1 will be multiplied with $\theta_0$ to give the intercept of the regression model. We call this 1 entry the intercept or bias term.

Given that we have three features here, we can express this model as:

Our features are represented by $x_1$ (FG), $x_2$ (AST), and $x_3$ (3PA) with each having corresponding parameters, $\theta_1$, $\theta_2$, and $\theta_3$.

In statistics, this model + loss is called Ordinary Least Squares (OLS). The solution to OLS is the minimizing loss for parameters $\hat{\theta}$, also called the least squares estimate.

From One Feature To Many Features¶

In an earlier chapter, we worked with the Simple Linear Regression (SLR) model $\hat{y} = \theta_0 + \theta_1x$, along the actual observations $y$. We can represent the dataset that corresponds to this model as follows:

We have also worked with the Constant Model $\hat{y} = \theta$, along with the actual observations $y$. We can represent the dataset that corresponds to this model even more simply, as there is no feature data to consider:

Accordingly, we can use a similar representation for our new Multiple Regression Model:

From here, we can pull out just the first row, in order to get the first observation $\{x_{11}, x_{12}, \cdots, x_{1p}, y_1\}$. Similarly, we can pull out just the second column $\{x_{12}, x_{22}, \cdots, x_{n2}\}$, in order to get the second feature.

Recall that this is our model for Multiple Linear Regression:

We can then write all of our predictions $\hat{y}_i$ like this:

Linear Algebra Approach¶

We now know how to generate a single prediction from multiple observed features. Data scientists usually work at scale – that is, they want to build models that can produce many predictions, all at once. The vector notation we introduced above gives us a hint on how we can expedite multiple linear regression. We want to use the tools of linear algebra.

Let’s think about how we can apply what we did above. To accommodate for the fact that we’re considering several feature variables, we’ll adjust our notation slightly. Each observation can now be thought of as a row vector with an entry for each of $p$ features.

To make a prediction from the first observation in the data, we take the dot product of the parameter vector and first observation vector. To make a prediction from the second observation, we would repeat this process to find the dot product of the parameter vector and the second observation vector. If we wanted to find the model predictions for each observation in the dataset, we’d repeat this process for all $n$ observations in the data.

Expanding our transpose vectors, we get the following:

We can then vectorize our predictions $\hat{y}$ and our parameters $\theta$. We will keep the $n$ row vectors, each with dimension $p + 1$.

We can collect all of our $x_i^T$ row vectors into a single matrix.

Briefly expanding our parameter vector $\theta$, we can see that our predictions take the following form:

Let’s collapse our parameter vector back into $\theta$ and call our matrix of observations $\mathbb{X}$.

Finally, we can call our prediction vector $\hat{\mathbb{Y}}$. We can now represent our Multiple Linear Regression model for many observations fully in terms of matrices.

The matrix $\mathbb{X}$ is known as the design matrix. It contains all observed data for each of our $p$ features, where each row corresponds to one observation, and each column corresponds to a feature. It often (but not always) contains an additional column of all ones to represent the intercept or bias column.

Our true output can also be represented as a vector, which we will call $\mathbb{R}^{p + 1}$. As an example, our true output in the basketball example was the actual recorded PTS (5.1, 1.7, etc.). This is distinct from the predictions our model would make for $PTS$ based on our features (FG, AST, 3PA)

To review what is happening in the design matrix: each row represents a single observation. For example, a student in Data 100. Each column represents a feature. For example, the ages of students in Data 100. This convention allows us to easily transfer our previous work in DataFrames over to this new linear algebra perspective.

To summarize where we have ended up, our model now takes the form $\mathbb{\hat{Y}} = \mathbb{X} \theta$. $\mathbb{\hat{Y}}$ is the prediction vector with $n$ elements ($\mathbb{\hat{Y}} \in \mathbb{R}^{n}$); it contains the prediction made by the model for each of the $n$ input observations. $\mathbb{X}$ is the design matrix with dimensions $\mathbb{X} \in \mathbb{R}^{n \times (p + 1)}$, and $\theta$ is the parameter vector with dimensions $\theta \in \mathbb{R}^{(p + 1)}$. Our true output $\mathbb{Y}$ contains $n$ elements ($\mathbb{Y} \in \mathbb{R}^{n}$).

Mean Squared Error¶

We now have a new approach to understanding models in terms of vectors and matrices. To accompany this new convention, we should update our understanding of risk functions and model fitting.

Recall our definition of MSE:

At its heart, the MSE is a measure of distance – it gives an indication of how “far away” the predictions are from the true values, on average.

We can express the MSE as a squared L2 norm if we rewrite it in terms of the prediction vector, $\hat{\mathbb{Y}}$, and true target vector, $\mathbb{Y}$:

Here, the superscript 2 outside of the parentheses means that we are squaring the norm. If we plug in our linear model $\hat{\mathbb{Y}} = \mathbb{X} \theta$, we find the MSE cost function in vector notation:

Under the linear algebra perspective, our new task is to fit the optimal parameter vector $\theta$ such that the cost function is minimized. Equivalently, we wish to minimize the norm

We can restate this goal in three ways:

• Minimize the mean squared error for the linear model $\mathbb{\hat{Y}} = \mathbb{X} \theta$

• Minimize the distance between the vector of true values, $\mathbb{Y}$, and the vector of predicted values, $\mathbb{\hat{Y}}$

• Minimize the length of the residual vector, defined as:

  $$e = \mathbb{Y} - \mathbb{\hat{Y}} = \begin{bmatrix} y_1 - \hat{y}_1 \\ y_2 - \hat{y}_2 \\ \vdots \\ y_n - \hat{y}_n \end{bmatrix}$$

A Note on Terminology for Multiple Linear Regression¶

There are several equivalent terms in the context of regression. The ones we use most often for this course are bolded.

• $x$ can be called a

  • Feature(s)

  • Covariate(s)

  • Independent variable(s)

  • Explanatory variable(s)

  • Predictor(s)

  • Input(s)

  • Regressor(s)

• $y$ can be called an

  • Output

  • Outcome

  • Response

  • Dependent variable

• $\hat{y}$ can be called a

  • Prediction

  • Predicted response

  • Estimated value

• $\theta$ can be called a

  • Weight(s)

  • Parameter(s)

  • Coefficient(s)

• $\hat{\theta}$ can be called a

  • Estimator(s)

  • Optimal parameter(s)

• A datapoint $(x, y)$ is also called an observation.

Geometric Derivation¶

Rethinking Matrix-Vector Multiplication¶

There are 2 ways we can think about matrix-vector multiplication

1. So far, we’ve thought of our model as horizontally stacked predictions per datapoint:

2. Notice that we can express our predictions $\mathbb{\hat{Y}}$ in the following way:

Notice that the $j$th element (zero-indexed) of every sum in the vector on the right side contains $\theta_j$. So, the zeroth element of every sum contains $\theta_0$, the second element of every sum (again, zero-indexed) contains $\theta_2$, etc. It would be nice to group these elements with shared $\theta$ values together!

Because adding vectors of the same length produces a new vector of element-wise sums, we can split this already-summed vector into the sum of vectors:

Note that any $j$th summed column vector contains exclusively elements of the form $\theta_j x_{ij}$ (except for the first one, which contains only elements of the form $\theta_0 = \theta_0 \times 1$). Thus, you can factor out $\theta_j$ from every $j$-th summed column vector:

3. It is sometimes very helpful to think of matrix-vector multiplication as performed by columns. We can also think of $\mathbb{Y}$ as a linear combination of feature vectors, scaled by parameters.

Up until now, we’ve mostly thought of our model as a scalar product between horizontally stacked observations and the parameter vector. We can also think of $\hat{\mathbb{Y}}$ as a linear combination of feature vectors, scaled by the parameters. We use the notation $\mathbb{X}_{:, i}$ to denote the $i$th column of the design matrix.

This new approach is useful because it allows us to take advantage of the properties of linear combinations.

Because the prediction vector, $\hat{\mathbb{Y}} = \mathbb{X} \theta$, is a linear combination of the columns of $\mathbb{X}$, we know that the predictions are contained in the span of $\mathbb{X}$. That is, we know that $\mathbb{\hat{Y}} \in \text{Span}(\mathbb{X})$.

Column Space And The Normal Equation¶

The diagram below is a simplified view of $\text{Span}(\mathbb{X})$, assuming that each column of $\mathbb{X}$ has length $n$. Notice that the columns of $\mathbb{X}$ define a subspace of $\mathbb{R}^n$, where each point in the subspace can be reached by a linear combination of $\mathbb{X}$'s columns. The prediction vector $\mathbb{\hat{Y}}$ lies somewhere in this subspace.

Examining this diagram, we find a problem. The vector of true values, $\mathbb{Y}$, could theoretically lie anywhere in $\mathbb{R}^n$ space – its exact location depends on the data we collect out in the real world. However, our multiple linear regression model can only make predictions in the subspace of $\mathbb{R}^n$ spanned by $\mathbb{X}$. Remember the model fitting goal we established in the previous section: we want to generate predictions such that the distance between the vector of true values, $\mathbb{Y}$, and the vector of predicted values, $\mathbb{\hat{Y}}$, is minimized. This means that we want $\mathbb{\hat{Y}}$ to be the vector in $\text{Span}(\mathbb{X})$ that is closest to $\mathbb{Y}$.

Another way of rephrasing this goal is to say that we wish to minimize the length of the residual vector $e$, as measured by its $L_2$ norm.

The vector in $\text{Span}(\mathbb{X})$ that is closest to $\mathbb{Y}$ is always the orthogonal projection of $\mathbb{Y}$ onto $\text{Span}(\mathbb{X}).$ Thus, we should choose the parameter vector $\theta$ that makes the residual vector orthogonal to any vector in $\text{Span}(\mathbb{X})$. You can visualize this as the vector created by dropping a perpendicular line from $\mathbb{Y}$ onto the span of $\mathbb{X}$.

Remember our goal is to find $\hat{\theta}$ such that we minimize the objective function $R(\theta)$. Equivalently, this is the $\hat{\theta}$ such that the residual vector $e = \mathbb{Y} - \mathbb{X} \hat{\theta}$ is orthogonal to $\text{Span}(\mathbb{X})$.

Looking at the definition of orthogonality of $\mathbb{Y} - \mathbb{X}\hat{\theta}$ to $span(\mathbb{X})$, we can write:

Let’s then rearrange the terms:

And finally, we end up with the normal equation:

Any vector $\theta$ that minimizes MSE on a dataset must satisfy this equation.

If $\mathbb{X}^T \mathbb{X}$ is invertible, we can conclude:

This is called the least squares estimate of $\theta$: it is the value of $\theta$ that minimizes the squared loss.

Note that the least squares estimate was derived under the assumption that $\mathbb{X}^T \mathbb{X}$ is invertible. This condition holds true when $\mathbb{X}^T \mathbb{X}$ is full column rank, which, in turn, happens when $\mathbb{X}$ is full column rank. We will get to the proof for why $\mathbb{X}$ needs to be full column rank in the OLS Properties section.

OLS Properties¶

Residuals¶

When using the optimal parameter vector, our residuals $e = \mathbb{Y} - \hat{\mathbb{Y}}$ are orthogonal to $span(\mathbb{X})$.

The Bias/Intercept Term¶

For all linear models with an intercept term, the sum of residuals is zero.

Existence of a Unique Solution¶

To summarize:

Uniqueness of the OLS Solution¶

In most settings, the number of observations (n) is much greater than the number of features (p).

In practice, instead of directly inverting matrices, we can use more efficient numerical solvers to directly solve a system of linear equations using the normal equation shown below. Note that at least one solution always exists because intuitively, we can always draw a line of best fit for a given set of data, but there may be multiple lines that are “equally good”. (Formal proof is beyond this course.)

The Least Squares estimate $\hat{\theta}$ is unique if and only if $\mathbb{X}$ is full column rank.

Therefore, if $\mathbb{X}$ is not full column rank, we will not have unique estimates. This can happen for two major reasons.

1. If our design matrix $\mathbb{X}$ is “wide”:

   • If n < p, then we have way more features (columns) than observations (rows).

   • Then $rank(\mathbb{X}) = min(n, p+1) < p+1$, so $\hat{\theta}$ is not unique.

   • Intuition

     • With one data point ($n=1$), there are $\infty$ SLR models ($p = 1$) that minimize MSE. Try drawing a picture!

   • More intuition

     • Given two non-overlapping vectors ($p = 1$) in 2D space ($n = 2$) there is exactly one linear combination that can reach each point in 2D space.

     • If you have three of these vectors ($p = 2$), there are $\infty$ possible combinations.

   • Typically we have $n$ >> $p$, so this is less of an issue.

2. If our design matrix $\mathbb{X}$ has features that are linear combinations of other features:

   • By definition, rank of $\mathbb{X}$ is number of linearly independent columns in $\mathbb{X}$.

   • Example: If “Width”, “Height”, and “Perimeter” are all columns,

     • Perimeter = 2 * Width + 2 * Height $\rightarrow$ $\mathbb{X}$ is not full rank.

   • Important with one-hot encoding (to discuss later).

Evaluating Model Performance¶

Our geometric view of multiple linear regression has taken us far! We have identified the optimal set of parameter values to minimize MSE in a model of multiple features. Now, we want to understand how well our fitted model performs.

RMSE¶

One measure of model performance is the Root Mean Squared Error, or RMSE. The RMSE is simply the square root of MSE. Taking the square root converts the value back into the original, non-squared units of $y_i$, which is useful for understanding the model’s performance. A low RMSE indicates more “accurate” predictions – that there is a lower average loss across the dataset.

Multiple $R^2$¶

For SLR, we used the correlation coefficient to capture the association between the target variable and a single feature variable. In a multiple linear model setting, we will need a performance metric that can account for multiple features at once. Multiple $R^2$, also called the coefficient of determination, is the proportion of variance of our fitted values (predictions) $\hat{y}_i$ to our true values $y_i$. It ranges from 0 to 1 and is effectively the proportion of variance in the observations that the model explains.

Note that for OLS with an intercept term, for example $\hat{y} = \theta_0 + \theta_1x_1 + \theta_2x_2 + \cdots + \theta_px_p$, $R^2$ is equal to the square of the correlation between $y$ and $\hat{y}$. On the other hand for SLR, $R^2$ is equal to $r^2$, the correlation between $x$ and $y$. The proof of these last two properties is out of scope for this course.

As we add more features, our fitted values tend to become closer and closer to our actual values. Thus, $R^2$ increases.

Returning to our basketball example:

• The SLR model (AST only) explains 45.7% of the variance in the true $y$.

• The AST & 3PA model explains 60.9% of the variance in the true $y$.

Adding more features doesn’t always mean our model is better, though! We’ll see why later in the course.

Demo: Predicting Penguin Bill Depth¶

Can we predict bill depth from flipper length and body mass? We can now apply what we have learned about Multiple Linear Regression models and Ordinary Least Squares to do just that. We will conclude by introducing the scikit-learn library for simple machine learning tasks, and using it to recreate the same prediction of bill depth.

Consider the penguins dataset.

Suppose we could measure flippers and weight easily, but not bill dimensions. How can we predict bill depth from flipper length and/or body mass?

For demo purposes, we’ll drop all columns except the variables whose relationships we’re interested in modeling.

Suppose we want to create a linear regression model that predicts a penguin’s bill depth $y$ using both their flipper length $x_1$ and body mass $x_2$, plus an intercept term.

Another way to write this:

Using The Solution To The Normal Equation¶

We saw in last lecture that we can model the above multiple linear regression equation using matrix multiplication:

The optimal $\hat{\theta}$ that minimizes MSE also solves the normal equation:

If $\mathbb{X}$ is full column rank, then there is a unique solution to $\hat{\theta}$:

Note: This derivation is one of the most challenging in the course, especially for those of you who are learning linear algebra during the same semester. Be kind to yourself!

Let’s try this in code. We have to add a bias column to the design matrix if we want to include an intercept in the model:

Recall the solution to the normal equation if $\mathbb{X}$ is full column rank:

So, our fitted model is as follows:

Note: The code above is inefficient. We won’t go into this in Data 100, but it’s better to use np.linalg.solve rather than computing the explicit matrix inverse:

Make Predictions¶

We’ll save the predictions in a column so we can compare them against predictions from sklearn in the next section.

Machine Learning Frameworks¶

We have introduced the Multiple Linear Regression model $\hat{\mathbb{Y}} = \mathbb{X}\theta$. We have expressed the Mean Squared Error loss function in vector notation:

We then used geometry to minimize average loss by fitting our model with optimal parameter values $\hat{\theta}$. We touched upon the use of matrix calculus to fit the model (out of scope for this course), and we used calculus to find the optimal parameter values for SLR in the previous chapters.

In research and industry, it is common to use libraries for creating and training models. In Data 100, we will use scikit-learn, commonly called sklearn.

The sklearn Workflow¶

At a high level, there are three steps to creating an sklearn model:

1. Initialize a new model instance.

• Make a “copy” of the model template

my_model = lm.LinearRegression()

2. Fit the model to the training data.

• Save the optimal model parameters

my_model.fit(X, y)

3. Use fitted models to make predictions.

my_model.predict(X_new)

To extract the fitted model’s parameters, use my_model.coef_ (to extract parameters for all features without the bias $\theta_1, \theta_2, \cdots, \theta_p$) and my_model.intercept_ (to extract the parameter for the bias $\theta_0$).

Let’s try this out.

Using sklearn to fit our Multiple Linear Regression Model¶

Do we obtain the same fitted OLS parameters provided by the normal equation using sklearn? We can use the sklearn workflow we just described to find out.

An alternate approach to optimize our model is to use the sklearn.linear_model.LinearRegression class.(Documentation)

1. Create an sklearn object.

   First we create a model. At this point the model has not been fit so it has no parameters.

2. fit the object to data.

Then we “fit” the model, which means computing the parameters that minimize the loss function.

The first argument of the .fit method should be a matrix, or DataFrame, with dimension $n \times p$ containing your feature data (the dimension is not $n \times p + 1$ because LinearRegression objects add an intercept by default). The second argument should be a vector, or Series, of dimension $n$ containing the observations we’re trying to predict.

Note: The LinearRegression class is hard-coded to use the MSE as its loss function.

3. Analyze fit or call predict.

Now that our model is trained, we can ask it questions.

The code below asks the model to estimate the bill depth (in mm) for three different penguins with different flipper lengths and body masses.

Here are the predictions for all of the penguins in our dataset:

Let’s save the scikit-learn model predictions as a column in our DataFrame:

Looking at the predictions above, we see that they are exactly the same as what we calculated with our analytic formula!

Analyze parameters.

We can ask the model for its intercept and slope with .intercept_ and .coef_, respectively.

So, our fitted model is as follows:

The parameters are the same as with our analytic formula.

They look the same!

Fitting General Models¶

The solution we found for the normal equation can only help us fit our model, or find our optimal parameters $\hat{\theta}$ when we use Multiple Linear Regression as our model and Mean Squared Error as our loss function. How do these software packages, like scikit-learn, fit general models?

We made some big assumptions with the calculus-based and geometric techniques for minimizing loss:

• Calculus:

  • We can’t always algebraically solve for the zero points of the derivative.

• Geometric:

  • We don’t always use a linear model with MSE loss.

To design more complex models with different loss functions, we need a new optimization technique: gradient descent.

In the next chapter, we will use this iterative algorithm to numerically compute the parameters $\hat{\theta}$ that minimize the loss.

[Bonus] Time Complexity of the Normal Equation¶

It can be easy to assume that a closed-form solution is always the fastest approach, but evaluating the normal equation on large datasets can be extremely computationally taxing and can take a lot of resources to calculate. Here we’ll explore why this is true.

We know that the solution to the normal equation is $\hat{\theta} = (\mathbb{X}^T\mathbb{X})^{-1}\mathbb{X}^T\mathbb{Y}$. Let’s break this down and determine the computational complexity for this solution.

Let $n$ be the number of samples (rows) and $d$ be the number of features (columns).

• Computing $(\mathbb{X}^{\top}\mathbb{X})$ takes $O(nd^2)$ time, and it’s inverse takes another $O(d^3)$ time to calculate; overall, $(\mathbb{X}^{\top}\mathbb{X})^{-1}$ takes $O(nd^2) + O(d^3)$ time.

• $\mathbb{X}^{\top}\mathbb{Y}$ takes $O(nd)$ time.

• Multiplying $(\mathbb{X}^{\top}\mathbb{X})^{-1}$ and $\mathbb{X}^{\top}\mathbb{Y}$ takes $O(d^2)$ time.

In total, calculating the solution to the normal equation takes $O(nd^2) + O(d^3) + O(nd) + O(d^2)$ time. We can see that $O(nd^2) + O(d^3)$ dominates the complexity — this can be problematic for high-dimensional models and very large datasets.