Linear Regression

In this lecture, we learn the linear regression model: how it predicts a continuous target such as a house price, how we measure its error with RMSE, and how we find its best parameters with the normal equation and with scikit-learn. The next lecture finds the same parameters with gradient descent.

This lecture has three more pages:

• Deriving the Normal Equation (optional) — where the normal equation comes from, and what happens when it fails.
• Interpreting a Linear Regression Model — linear regression on the housing data, and what the weights, the bias, and the RMSE mean in dollars, plus which features matter most.
• OLS Regression and R-squared (optional) — the same model from a statistics perspective with statsmodels: p-values, confidence intervals, and R-squared.

Dataset

We use the California housing dataset from the Matplotlib lecture: housing.csv.zip (also on Kaggle). Unzip it into your ml4biz folder and run the code below in a notebook there — all the packages are already installed.

Classification vs. Regression

• Classification: the target value is discrete/categorical, e.g., we try to predict the class label for a given data point, such as Survived/Died, Churn/Loyal.
• Regression: the target value is continuous, e.g., we aim to predict a value for the target variable, such as house price and sales.

Functional Relation vs. Statistical Relation

• For functional relations, variables can be expressed by a mathematical formula $Y = f(X)$, such as the linear relationship between the gross sales of a product ($s$) and its quantity sold ($q$) at a fixed price of $10: $s = 10q$

  

• For statistical relations, the data points (observations) do not fall directly on the curve of the relationship, such as income and housing price.

  

The goal of machine learning is to learn (estimate) the functional relationships from the statistical relationships to make predictions:

The code for the figures on this page is in the Appendix.

Linear Regression Definition

Linear regression is one of the simplest machine learning models that aims to learn the linear relationship between the features and target.

A linear regression model can be informally defined as follows:

• Target is also called the dependent variable
• Features are also called the independent variables
• $w_1, \ldots, w_d$ and $b$ are the model parameters, including:
  • the feature weights $w_1, w_2, \ldots, w_d$
  • the bias term (also called intercept term) $b$, which captures all other factors that influence the dependent/target variable other than the features.

Univariate and Multivariate Regression

A linear regression model assumes that the relationship between the target variable and the features is linear.

• Univariate: one feature. When there is only one feature (independent variable), it's called Univariate Linear Regression, such as $y = 5x + 3$, which is a line as shown in the example above.
• Multivariate: multiple features. When there are two or more features (independent variables), it's called Multivariate Linear Regression.
  • when there are two features, such as $y = 3x_1 + 5x_2 + 4$, it is a plane in 3D space
  • when there are three or more, such as $y = 3x_1 + 5x_2 + 7x_3 + 9$, it is a hyperplane in higher dimensional space (hard to visualize)

The Housing Example

Here we are going to predict the median house value. For simplicity, we only use three features — population, median_income, and housing_median_age:

import numpy as np
import pandas as pd

df = pd.read_csv('housing.csv')
df = df[['population', 'median_income', 'housing_median_age', 'median_house_value']]
df.head()

Output:
   population  median_income  housing_median_age  median_house_value
0       322.0         8.3252                41.0            452600.0
1      2401.0         8.3014                21.0            358500.0
2       496.0         7.2574                52.0            352100.0
3       558.0         5.6431                52.0            341300.0
4       565.0         3.8462                52.0            342200.0

Then, the linear regression model is simply:

Our goal is to find the best values for $w_1, w_2, w_3$ and $b$ via training. Each row of the training data gives one equation:

• first row: $452600 = w_1 \cdot 322 + w_2 \cdot 8.3252 + w_3 \cdot 41.0 + b$
• second row: $358500 = w_1 \cdot 2401 + w_2 \cdot 8.3014 + w_3 \cdot 21.0 + b$
• ...

Vector Notation

In machine learning, vector notation is used to write definitions and algorithms more concisely.

A scalar is a single number, such as $3$ or $-5.2$. One house's median income, or its price, is a scalar. A vector is an ordered list of scalars, such as all the feature values of one house.

Row Vectors

A row vector is just a 1D NumPy array. We write vectors as bold lowercase letters:

import numpy as np

u = np.array([1, 3, -5])
v = np.array([4, -2, -1])

This course follows the notation of Dive into Deep Learning: a plain lowercase letter such as $x$ is a scalar, a bold lowercase letter such as $\mathbf{x}$ is a vector, and a bold capital such as $\mathbf{X}$ is a matrix. In LaTeX, write \mathbf{x}. Some books write vectors with an arrow instead, $\vec{u}$, so you may see both.

Column Vectors and the Dot Product

You need some basic linear algebra to understand vector/matrix calculations, such as the dot product and the transpose.

In machine learning, vectors are usually column vectors instead of row vectors. This follows the linear algebra convention of writing a matrix acting on a vector from the left, $\mathbf{y} = \mathbf{A}\mathbf{x}$: for the shapes to match, $\mathbf{x}$ has to be a column of shape $(n, 1)$. The convention is much older than machine learning and was inherited from math and physics. It also keeps chains of operations readable: applying $\mathbf{B}$ and then $\mathbf{A}$ is written $\mathbf{A}\mathbf{B}\mathbf{x}$, the same order as function notation $f(g(x))$.

The transpose ($^\top$) of a column vector is a row vector and vice versa. For a matrix, the transpose flips it over its diagonal, so rows become columns: a matrix of shape $(n, d)$ becomes $(d, n)$.

The dot product can be written as a matrix multiplication:

# a row vector/1D array has shape (3,)
# a column vector/2D array has shape (3, 1)

# u and v are 1D arrays
u = np.array([1, 3, -5])
v = np.array([4, -2, -1])
print(u, v, u.shape, v.shape)

print('*' * 50)
# reshape u and v to be column vectors
u = u.reshape(3, 1)
v = v.reshape(3, 1)
print(u)
print(v)
print(u.shape, v.shape)

print('*' * 50)
# dot product
print(np.dot(u.T, v))

Output:
[ 1  3 -5] [ 4 -2 -1] (3,) (3,)
**************************************************
[[ 1]
 [ 3]
 [-5]]
[[ 4]
 [-2]
 [-1]]
(3, 1) (3, 1)
**************************************************
[[3]]

Linear Regression in Vector Notation

We use the following notation:

• $n$ is the number of instances (examples) in the training dataset, e.g., $n = 16512$ for the housing dataset after an 80/20 train/test split
• $d$ is the number of features, e.g., $d = 3$ for the housing example
• $\mathbf{x}^{(i)}$ is the vector of all the feature values of the $i^{th}$ instance, and $y^{(i)}$ is its target value
• $\mathbf{w}$ is the vector of feature weights $w_1, \ldots, w_d$, and $b$ is the bias, a scalar; together they are the model parameters

Recall the first two rows of the housing data:

• first row: $452600 = w_1 \cdot 322 + w_2 \cdot 8.3252 + w_3 \cdot 41.0 + b$
• second row: $358500 = w_1 \cdot 2401 + w_2 \cdot 8.3014 + w_3 \cdot 21.0 + b$

The weights form the column vector $\mathbf{w}$, and the first row of feature values forms the column vector $\mathbf{x}^{(1)}$:

The first row's equation is then a dot product plus the bias, $y^{(1)} = \mathbf{w}^\top \mathbf{x}^{(1)} + b$, with $y^{(1)} = 452600$. In general, for an instance with features $x_1, \ldots, x_d$ the model predicts $\hat{y} = w_1 x_1 + \cdots + w_d x_d + b = \mathbf{w}^\top \mathbf{x} + b$, so the prediction for the $i^{th}$ instance is:

To handle all rows at once, stack the transposed feature vectors into the matrix $\mathbf{X}$ — one row per instance — and the targets into the column vector $\mathbf{y}$:

Now, we can put all rows together in vectorized form — $\mathbf{X} \in \mathbb{R}^{n \times d}$ has shape $(n, 3)$ and $\mathbf{w}$ has shape $(3, 1)$, so $\mathbf{X}\mathbf{w}$ gives one value per row, and the scalar $b$ is added to every row:

Row $i$ of $\mathbf{X}\mathbf{w}$ is $(\mathbf{x}^{(i)})^\top \mathbf{w}$, which equals $\mathbf{w}^\top \mathbf{x}^{(i)}$ because a dot product does not depend on order. So this computes $\hat{y}^{(i)} = \mathbf{w}^\top \mathbf{x}^{(i)} + b$ for every row at once. In the matrix form, $\mathbf{X}$ must come first so the shapes match: $(n, 3) \times (3, 1)$ gives $(n, 1)$.

The bias $b$ is one number shared by all instances. Adding a scalar to every entry of a vector is called broadcasting, and NumPy does the same when you write X @ w + b.

Evaluation Metric for Linear Regression

In linear regression we try to learn/estimate a linear relationship. For example, given the following two linear relationships, how do we decide which one is better?

We need to define the prediction error as follows:

• The ground truth for the $i^{th}$ instance is $y^{(i)}$
• The prediction is $\hat{y}^{(i)} = \mathbf{w}^\top \mathbf{x}^{(i)} + b$
• The prediction error, or residual, is $y^{(i)} - \hat{y}^{(i)}$

For $n$ instances, we square the errors, average them, and take the square root to get one of the most commonly used evaluation metrics of a regression model, the Root Mean Square Error (RMSE):

Our goal is to minimize RMSE.

NOTE: the RMSE is measured on the same scale with the same units as $y$.

Back to the two lines above: the green line has an RMSE of 0.898 and the red line 1.009, so the green line is better. In fact, no line has a lower RMSE on this data — the green line is the one the normal equation finds below.

In practice, it is often simpler to minimize the Mean Square Error (MSE) instead — the $\mathbf{w}$ and $b$ that minimize MSE also minimize RMSE:

Another common evaluation metric is the Mean Absolute Error (MAE), the average size of the errors, ignoring their sign:

For example, if the predictions for three houses are off by $10k, -$20k, and $30k:

• MAE $= (10 + 20 + 30) / 3 = 20$, so on average we are off by $20k
• RMSE $= \sqrt{(10^2 + 20^2 + 30^2) / 3} = \sqrt{467} \approx 21.6$, a bit higher, because squaring makes the $30k error count more

A big gap between MAE and RMSE means a few predictions are very wrong. Use RMSE when big mistakes are especially costly; use MAE when every dollar of error costs the same, or when you need a number that is easy to explain. Linear regression is trained to minimize MSE, so RMSE is its usual metric.

Now that we know the goal, how can we find the parameters that minimize RMSE/MSE?

• The normal equation (this lecture)
• Gradient descent (next lecture)

The Normal Equation

There exists a "closed-form" solution to find the optimal values of $\mathbf{w}$ and $b$ that minimize MSE, which is called the Normal Equation (if you are interested in where it comes from, see Deriving the Normal Equation (optional)).

To write it in one formula, we absorb the bias into the weights by appending a column of 1s to $\mathbf{X}$ (as its first column) and making $b$ the first entry of $\mathbf{w}$. Then $\mathbf{X}\mathbf{w}$ already includes $b$, because each row computes $1 \cdot b + w_1 x_1 + \cdots + w_d x_d$. With this extended $\mathbf{X}$ and $\mathbf{w}$:

• $\mathbf{w}^*$ is the value of $\mathbf{w}$ that minimizes the cost function, i.e., MSE; its first entry is the bias $b$ and the rest are the weights $w_1, \ldots, w_d$
• $\mathbf{X}$ is the matrix of features with a column of 1s added for the bias, so it has shape $(n, d + 1)$, and $\mathbf{y}$ is the vector of targets
• $^\top$ is the transpose (see Column Vectors and the Dot Product) and $^{-1}$ is the matrix inverse: for a square matrix $\mathbf{A}$, $\mathbf{A}^{-1}$ is the matrix with $\mathbf{A}^{-1}\mathbf{A} = \mathbf{I}$ (the identity matrix), the matrix version of $1/a$. Not every matrix has one.

The normal equation is one way to find the parameters of a linear regression model. Gradient Descent shows another, and explains the difference between a model and the algorithm used to train it.

Let's try it on simulated data, where we know the ideal model is $y = 5 + 3x$:

# fit linear regression using the normal equation
import numpy as np

np.random.seed(42)  # fix the random seed so that each run generates the same set of random numbers
X = 2 * np.random.rand(100, 1)  # generate 100 random numbers between 0 and 2 with shape (100, 1)
y = 5 + 3 * X + np.random.randn(100, 1)  # add noise from a standard normal distribution

X_with_bias = np.c_[np.ones((100, 1)), X]  # add a column of 1s for the bias b
print(X_with_bias[:2])
print('*' * 50)

# np.linalg.inv() calculates the inverse matrix and .dot() does matrix multiplication
# w_best[0] is the bias b, w_best[1] is the weight w1
# so our best model is: y = 5.215 + 2.770 * x
w_best = np.linalg.inv(X_with_bias.T.dot(X_with_bias)).dot(X_with_bias.T).dot(y)
print(w_best)
print('*' * 50)

# given a new set of X_new, we predict y_pred using the dot product of X and w_best
X_new = np.array([[0.6], [0.9], [1.3]])
X_new_with_bias = np.c_[np.ones((3, 1)), X_new]  # add a column of 1s for the bias b
y_pred = X_new_with_bias.dot(w_best)
print(y_pred)

Output:
[[1.         0.74908024]
 [1.         1.90142861]]
**************************************************
[[5.21509616]
 [2.77011339]]
**************************************************
[[6.87716419]
 [7.70819821]
 [8.81624356]]

The estimates (5.215 and 2.770) are close to, but not exactly, the ideal 5 and 3, because of the noise in the data.

With scikit-learn

For the above example, we can use scikit-learn to get the same result. scikit-learn solves the same least-squares problem, but instead of inverting $\mathbf{X}^\top \mathbf{X}$, it uses a more numerically stable method called Singular Value Decomposition (SVD). SVD also works when $\mathbf{X}^\top \mathbf{X}$ has no inverse.

from sklearn.linear_model import LinearRegression

lin_reg = LinearRegression()
lin_reg.fit(X, y)  # no need to add the bias column: LinearRegression fits the intercept for you

print(lin_reg.predict(X_new))
print('*' * 50)
print(lin_reg.intercept_, lin_reg.coef_)  # intercept_ is b, coef_ is w: exactly the same as our calculation above!

Output:
[[6.87716419]
 [7.70819821]
 [8.81624356]]
**************************************************
[5.21509616] [[2.77011339]]

Appendix: Plotting Code

The code below draws the figures on this page. It generates its own data, so you can run it on its own.

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8')

# figure 1: a functional relation - sales = 10 * quantity
quantity = np.arange(0, 10)
sales = 10 * quantity

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(quantity, sales, marker='o')
ax.set_xlabel('Quantity sold (q)')
ax.set_ylabel('Gross sales (s)')
ax.set_title('Functional relation: s = 10q')

# the simulated data from the normal equation example: y = 5 + 3x + noise
np.random.seed(42)
X = 2 * np.random.rand(100, 1)
y = 5 + 3 * X + np.random.randn(100, 1)

# figure 2: a statistical relation - the points scatter around a line
fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(X, y, s=15)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Statistical relation: the points do not fall on a line')

# the best line from the normal equation: y = 5.215 + 2.770x
X_with_bias = np.c_[np.ones((100, 1)), X]
w_best = np.linalg.inv(X_with_bias.T.dot(X_with_bias)).dot(X_with_bias.T).dot(y)  # w_best[0] is b, w_best[1] is w1
x_line = np.array([0, 2])
y_best = w_best[0, 0] + w_best[1, 0] * x_line

# figure 3: the line learned from the data
fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(X, y, s=15)
ax.plot(x_line, y_best, color='tab:green', linewidth=2,
        label=f'learned: y = {w_best[0, 0]:.3f} + {w_best[1, 0]:.3f}x')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.legend(loc='upper left')

# figure 4: which line is better? compare their RMSE
def rmse(intercept, slope):
    return np.sqrt(np.mean((intercept + slope * X - y) ** 2))

fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(X, y, s=15)
ax.plot(x_line, y_best, color='tab:green', linewidth=2,
        label=f'green: y = {w_best[0, 0]:.3f} + {w_best[1, 0]:.3f}x (RMSE {rmse(w_best[0, 0], w_best[1, 0]):.3f})')
ax.plot(x_line, 6 + 2 * x_line, color='tab:red', linewidth=2,
        label=f'red: y = 6 + 2x (RMSE {rmse(6, 2):.3f})')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.legend(loc='upper left')

# figure 5: residuals - the vertical distance from each point to the line
order = np.argsort(X[:, 0])[::8]  # 13 points spread across the x range, to keep the figure readable
X_few, y_few = X[order, 0], y[order, 0]
y_pred_few = w_best[0, 0] + w_best[1, 0] * X_few
i = np.argmax(y_few - y_pred_few)  # label the point with the largest positive residual

fig, ax = plt.subplots(figsize=(8, 5))
ax.vlines(X_few, y_pred_few, y_few, color='tab:red', linewidth=1.5, label='residuals')
ax.plot(x_line, y_best, color='tab:green', linewidth=2, label='model')
ax.scatter(X_few, y_few, s=40, color='tab:blue', zorder=3, label='actual values $y^{(i)}$')
ax.scatter(X_few, y_pred_few, s=40, color='tab:orange', zorder=3, label=r'predicted values $\hat y^{(i)}$')
ax.annotate(r'residual $= y^{(i)} - \hat y^{(i)}$',
            xy=(X_few[i], (y_few[i] + y_pred_few[i]) / 2),
            xytext=(X_few[i] + 0.1, y_few[i] + 1.0),
            arrowprops=dict(arrowstyle='->', color='black'))
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.legend(loc='lower right')