Linear Regression

Linear regression is the workhorse of supervised learning: a model that predicts a continuous target as a linear combination of features, fit by minimising a squared-error loss. This lesson derives the cost function from first principles, shows the gradient that drives optimisation, walks through gradient descent with real numeric iterations, explains how to read the fitted coefficients, defines the coefficient of determination $R^2$, lists the four classical assumptions together with the symptom of each violation, and closes with an end-to-end worked example on a tiny dataset.

Learning Objectives

1. State the linear-regression hypothesis $\hat{y} = \theta^\top x$ and identify the role of each parameter vector entry.
2. Derive the mean-squared-error cost function, take its gradient with respect to the parameters, and set the gradient to zero to obtain the normal equation.
3. Run a few iterations of batch gradient descent on a one-feature problem and explain how the learning rate shapes the trajectory.
4. Interpret the fitted coefficients as feature effects, define $R^2$, and state why $R^2$ can mislead.
5. Enumerate the linearity, independence, homoscedasticity, and normality assumptions, and recognise the residual pattern that signals each violation.

1. The Linear Hypothesis

A linear regression model assumes the target is a linear combination of the features. For a single training example with feature vector $x = (1, x_1, x_2, \ldots, x_d)^\top$ — the leading $1$ is the constant term, sometimes called the bias unit — the prediction is

where $\theta = (\theta_0, \theta_1, \ldots, \theta_d)^\top$ is the parameter vector. "Linear" refers to linearity in the parameters, not in the features: the model $\hat{y} = \theta_0 + \theta_1 x + \theta_2 x^2$ is still linear regression, because $\theta_1$ and $\theta_2$ appear linearly. A model $\hat{y} = \theta_0 e^{\theta_1 x}$ is not, even though the relationship between $x$ and $\hat{y}$ is monotone.

For $n$ training examples stacked into a design matrix $X \in \mathbb{R}^{n \times (d+1)}$ whose $i$-th row is $x_i^\top$, the vector of predictions is

The columns of $X$ are the features; the rows are the observations. The convention that every row begins with a $1$ folds the intercept $\theta_0$ into the matrix product, so the same code works whether or not an intercept is included.

The hypothesis is intentionally simple — that is its strength. When the underlying relationship is approximately linear, no more expressive model is needed; the variance saved by not fitting extra parameters becomes prediction accuracy on new data.

2. The Mean-Squared-Error Cost Function

Pick $\theta$ by minimising a loss that penalises the gap between $\hat{y}_i$ and the observed $y_i$. The classical choice is the mean squared error

Three reasons make MSE the default regression loss.

1. Algebraic convenience. The square kills the sign so positive and negative residuals are penalised symmetrically; the gradient becomes linear in $\theta$; the optimum has a closed-form solution.
2. Probabilistic interpretation under Gaussian noise. If the targets are generated by $y = x^\top \theta + \varepsilon$ with $\varepsilon \sim \mathcal{N}(0, \sigma^2)$ independent across examples, minimising the negative log-likelihood is exactly minimising the MSE. The squared loss is therefore the maximum-likelihood estimator under additive Gaussian noise.
3. Convexity. The MSE is a positive-definite quadratic form in $\theta$; it has a single global minimum and no local minima to trap an optimiser.

The vector form, used in every line of production code that fits a linear model, is

where $\|v\|_2^2 = v^\top v$ is the squared Euclidean norm. The factor $1/n$ does not change the optimum; it is conventional so that $J$ is on the same scale as a per-example loss and so that gradients do not scale with the dataset size.

2.1 Why not mean absolute error?

Mean absolute error $\text{MAE} = \frac{1}{n}\sum |\hat{y}_i - y_i|$ is also a valid loss. It penalises outliers less harshly, which makes it preferable when the noise is heavy-tailed or when occasional large errors are acceptable. Its gradient is not defined at the residuals equal to zero, however, and the optimisation has no closed form. Use MAE when robustness to outliers matters; use MSE when outlier sensitivity is acceptable or desirable and when closed-form inference is convenient.

3. Gradient of the MSE

Differentiate $J(\theta)$ with respect to $\theta$ and set the result to zero to find the optimum. The scalar form expands as

Take the partial derivative with respect to one parameter $\theta_j$. The chain rule gives

Stacking the $d+1$ partials into a vector, using the design matrix, the gradient is

The two factors have clear interpretations. The residual vector $r = X\theta - y$ says which direction each prediction is wrong and by how much; the multiplication $X^\top r$ projects those residuals onto each parameter axis. A feature whose column is large in absolute value has a big effect on the gradient, which is why feature scaling matters for gradient descent: features measured in millions push the gradient in their direction much harder than features measured in fractions.

Setting the gradient to zero,

which is the normal equation. As long as $X^\top X$ is invertible — that is, no feature is a linear combination of the others — the unique solution is

4. Gradient Descent With Real Numbers

The normal equation is closed-form but cubic in the number of features. When $d$ is large, or when the data arrives in a stream and the model must be updated online, gradient descent is preferred. The update rule is

where $\alpha > 0$ is the learning rate. The negative gradient is the direction of steepest descent; $\alpha$ controls how far to step in that direction.

4.1 A one-feature dataset

Use the same five-point dataset introduced in lesson 1, restated here with simpler numbers for clarity.

With $n = 5$ the relevant sums are $\sum x_i y_i = 47{,}140$, $\sum x_i^2 = 2{,}125$, $\sum x_i = 97$, $\sum y_i = 2{,}200$. Means are $\bar{x} = 19.4$, $\bar{y} = 440$. The closed-form solution (rescaled back to thousands of dollars per 100 sq ft) is $\theta_0 \approx 84.21$, $\theta_1 \approx 18.34$, so the fitted line is

4.2 Iterating gradient descent on this surface

Start at $\theta_0 = 0$, $\theta_1 = 0$ and use learning rate $\alpha = 0.001$. The MSE surface for this dataset is a smooth, convex bowl with its minimum at the closed-form solution above. Each gradient step takes a small bite out of the loss.

Two observations. First, the loss falls rapidly at first — most of the error lives in the first ten iterations — and then trails off as the iterates approach the closed-form minimum near $\theta_0 \approx 17.74$, $\theta_1 \approx 25.36$. (The slight discrepancy with the rescaled figures above is the $1/n$ vs $1/(2n)$ convention; the minimum is the same.) Second, the iterates overshoot slightly before settling, a signature of a learning rate that is large enough to explore but small enough to converge.

4.3 What changes when $\alpha$ changes

The bowl shape is convex but anisotropic — long and narrow along one axis, short and steep along another — when features have very different scales. Standardising every feature to zero mean and unit variance before running gradient descent makes the bowl near-spherical and lets a single learning rate work for all parameters.

4.4 Gradient descent in Python

import numpy as np

x = np.array([10.0, 15.0, 18.0, 24.0, 30.0])
y = np.array([250.0, 380.0, 410.0, 540.0, 620.0])

X = np.column_stack([np.ones_like(x), x])

theta = np.zeros(2)
alpha = 1e-3
n = len(y)

for it in range(200):
    residual = X @ theta - y
    grad = (2.0 / n) * (X.T @ residual)
    theta -= alpha * grad

print(f"theta_0 = {theta[0]:.2f}, theta_1 = {theta[1]:.2f}")

The four lines inside the loop are the entire algorithm: form the residual vector, project it onto each parameter via $X^\top r$, scale by $2/n$, and subtract the result scaled by $\alpha$ from the current parameters. The same skeleton, with residual replaced by a stochastic estimate, becomes stochastic gradient descent.

5. The Closed-Form Solution (Normal Equation)

Setting the gradient to zero gives $\theta^\star = (X^\top X)^{-1} X^\top y$, also known as the ordinary least squares (OLS) solution. Three properties make it worth knowing.

1. Exact in one shot. No learning rate, no iteration count, no convergence worries. The output is the unique minimum of the convex MSE surface.
2. Cheap when $d$ is small. Cost is dominated by the matrix product $X^\top X$, which is $O(n d^2)$, plus the inversion or solve, which is $O(d^3)$. For $d$ up to a few thousand the closed form is faster than gradient descent to convergence.
3. Brittle when $X^\top X$ is ill-conditioned. Features that are near-collinear make $X^\top X$ nearly singular; numerical errors in the inversion dominate. The standard fix is to add a small ridge regulariser, giving the ridge regression solution $\theta^\star = (X^\top X + \lambda I)^{-1} X^\top y$, covered in a later lesson.

In production code the closed form is rarely written by hand. The standard tool is

import numpy as np
from sklearn.linear_model import LinearRegression

model = LinearRegression()
model.fit(x.reshape(-1, 1), y)
print(model.intercept_, model.coef_)   # ~17.74, [25.36]

Both np.linalg.lstsq and scikit-learn's LinearRegression solve the normal equation under the hood; the difference is that scikit-learn also stores metadata, supports multiple features, and plays nicely with the rest of the library's ecosystem.

6. Interpreting Coefficients

Each entry of $\theta^\star$ is read as the marginal effect of one feature on the predicted target, holding every other feature fixed.

For the worked example above, $\theta_1 = 25.36$ thousand dollars per 100 sq ft means that a 100-sq-ft increase in floor area is associated with a $25{,}360 increase in predicted sale price, all else equal. The intercept $\theta_0 = 17.74$ thousand dollars is the prediction when every feature is zero — for floor area, that is the extrapolation to a zero-square-foot "house", which is meaningless as a price but necessary for the algebra to work.

A few rules make coefficient interpretation more reliable.

A prediction interval for a new observation is wider than a confidence interval on the mean prediction: it includes both the uncertainty in $\theta^\star$ and the unavoidable noise $\sigma^2$. Reporting only the point estimate hides this distinction.

7. R-squared and Its Limits

The coefficient of determination is the fraction of variance in $y$ that the model explains:

A model that predicts the mean for every observation has $R^2 = 0$; a model with zero residual error has $R^2 = 1$. For simple linear regression, $R^2$ is also the square of the Pearson correlation between $x$ and $y$.

Three reasons $R^2$ can mislead.

1. It always rises when a feature is added. Adding a noisy feature to the design matrix increases $\text{SS}_{\text{res}}$ by a tiny amount or not at all, so $R^2$ creeps upward even for useless features. Adjusted $R^2$ penalises the model for each parameter added and should be reported alongside the raw value.
2. It says nothing about prediction quality on new data. A model with $R^2 = 0.95$ in-sample can have $R^2 = 0.30$ out-of-sample if it has overfit. Cross-validated $R^2$ or a held-out test set is the right diagnostic.
3. It is undefined outside the range of the data. $R^2$ compares to the training mean; if the test distribution has shifted, the right comparison is the test mean, which $R^2$ does not respect.

The right metric depends on the question. MSE and RMSE are on the same scale as $y$ and are easy to communicate to non-technical stakeholders. MAE is more robust to outliers. MAPE (mean absolute percentage error) is unit-free but blows up when any $y_i$ is near zero. Pick the metric that aligns with the cost structure of the application, then report it on a held-out set.

8. The Four Assumptions of Linear Regression

The OLS estimator has known statistical properties — unbiased, minimum variance among linear unbiased estimators — only under four assumptions on the data-generating process. Diagnose each one with a residual plot before trusting a fitted model.

When an assumption fails, two broad fixes exist. The first is to transform the data: take logs of a heavy-tailed target, take square roots of count data, or fit the model on a different scale. The second is to change the model: replace ordinary least squares with weighted least squares when heteroscedasticity is the problem, with generalised least squares for correlated errors, or with a tree-based model when the relationship is genuinely non-linear.

A diagnostic routine that catches most issues in five minutes:

1. Plot residuals vs fitted. Look for curvature (linearity) and funnel shape (homoscedasticity).
2. Plot residuals vs each feature. Look for structure that the model missed.
3. Plot the Q-Q plot of residuals. Look for heavy tails (normality).
4. For time-ordered data, plot residuals vs time and compute Durbin-Watson.
5. Compute Cook's distance per observation. Points with distance above $4/n$ are high-leverage outliers worth investigating.

8.1 When Linear Regression Is Not Enough

Linear regression is the wrong tool in three common situations.

• The relationship is non-linear. Use a basis expansion (polynomial features, splines) or switch to a tree-based model, kernel method, or neural network. Adding polynomial features keeps the linear-regression machinery and gains flexibility, at the cost of more parameters and more overfitting risk.
• The target is bounded. If $y$ is a probability, a count, or a non-negative quantity, the Gaussian noise assumption is wrong. Use logistic regression, Poisson regression, or a generalised linear model with the appropriate link.
• Features are highly correlated. Multicollinearity inflates the variance of the coefficients and makes them uninterpretable. Drop redundant features, combine them into a single index, or use ridge regression to regularise.

Reach for a more flexible model when the bias of linear regression is clearly limiting performance. Reach for a simpler model when the variance is the binding constraint. The art of applied ML is choosing where on the bias-variance curve the problem actually lives.

9. Worked Example: Sales Versus Advertising Spend

A small marketing team has tracked weekly sales (in thousands of dollars) against weekly spending on online ads (in thousands of dollars) for ten weeks. The data:

The relationship looks plausibly linear; the team wants a model to forecast sales given a planned ad budget.

9.1 Computing the coefficients

The required sums are

• $\sum x_i = 30.0$, $\sum y_i = 578$
• $\bar{x} = 3.0$, $\bar{y} = 57.8$
• $\sum x_i y_i = 1.2 \cdot 32 + 1.6 \cdot 38 + \cdots + 4.8 \cdot 81 = 1{,}912.0$
• $\sum x_i^2 = 1.44 + 2.56 + 4.00 + 5.76 + 7.84 + 10.24 + 12.96 + 16.00 + 19.36 + 23.04 = 103.20$

Plug into the closed-form slope

The intercept is

The fitted regression line is therefore

Interpretation: every extra $1{,}000 of weekly ad spend is associated with about $13{,}480 of additional weekly sales. The intercept says the baseline sales — with no ad spend — would be about $17{,}360, which captures organic demand and brand recognition.

9.2 Predictions, residuals, and $R^2$

For each of the ten weeks, the prediction and the residual:

Residuals oscillate gently around zero, with no obvious pattern — a good sign for the linearity and homoscedasticity assumptions. The residual sum of squares is

The total sum of squares, comparing each $y_i$ to the mean $\bar{y} = 57.8$, is

The coefficient of determination is therefore

The model explains about 99.7% of the variance in sales, an unusually tight fit. The RMSE, on the scale of $y$ (thousands of dollars), is

so the typical prediction error is about $830 in weekly sales.

9.3 A forecast for next week

If the marketing team plans to spend $5{,}000 on ads next week, the model's forecast is

or roughly $84{,}760 in sales. The right way to communicate this to the team is "around $85k, give or take $1k, assuming the relationship we have seen continues to hold for this budget level". The "give or take" is the residual standard error; the "assuming" is the linearity assumption; both belong in the conversation.

9.4 Cross-check with scikit-learn

The whole computation in twenty lines of Python:

import numpy as np
from sklearn.linear_model import LinearRegression

x = np.array([1.2, 1.6, 2.0, 2.4, 2.8, 3.2, 3.6, 4.0, 4.4, 4.8]).reshape(-1, 1)
y = np.array([32, 38, 45, 50, 56, 61, 67, 72, 76, 81], dtype=float)

model = LinearRegression().fit(x, y)
print(f"intercept = {model.intercept_:.2f}, slope = {model.coef_[0]:.2f}")

r_squared = model.score(x, y)
print(f"R^2 = {r_squared:.4f}")

y_pred = model.predict(x)
rmse = float(np.sqrt(np.mean((y - y_pred) ** 2)))
print(f"RMSE = {rmse:.2f} thousand dollars")

forecast = model.predict(np.array([[5.0]]))
print(f"forecast for $5k ad spend: ${float(forecast[0]):.2f}k")

The four printed lines match the hand calculation. In production code the script would split into train and test sets, report metrics on the test set rather than on the training set, and use np.linalg.lstsq or model.predict inside a pipeline that handles feature scaling, missing-value imputation, and cross-validation.

9.5 Sanity checks before shipping the model

The coefficients are credible. The residuals look like noise. $R^2$ is high. Before declaring victory, three further checks belong in the workflow.

1. Residual plot. Plot $y - \hat{y}$ against $\hat{y}$. A random cloud around zero confirms linearity and homoscedasticity; a curved shape means a non-linear feature is missing; a fan shape means the variance grows with the prediction.
2. Out-of-sample test. Hold out the last two weeks, refit on the first eight, and score on the held-out two. A test $R^2$ close to the training $R^2$ is the proof that the model generalises.
3. Domain plausibility. Is a 13.5:1 return on ad spend realistic for this business? If the team's historical return on ads is closer to 5:1, the model is probably picking up a confound (a new product launch that drove both sales and ad spend upward, for example). The numbers should pass a sniff test from a domain expert before they pass to a stakeholder.

A model that is statistically correct but operationally wrong is a worse outcome than no model at all: it gives false confidence to a decision that should have been made by other means.

Key Takeaways

• The linear-regression hypothesis $\hat{y} = \theta^\top x$ is a linear combination of features fit by minimising the mean-squared-error loss; "linear" means linear in the parameters, not necessarily in the features.
• The gradient of the MSE is $\nabla_\theta J = (2/n) X^\top (X\theta - y)$; setting it to zero gives the normal equation $\theta^\star = (X^\top X)^{-1} X^\top y$, the closed-form OLS solution.
• Gradient descent iteratively subtracts a learning-rate-scaled gradient from the parameters; feature scaling and a sensible learning rate are the keys to fast convergence.
• Each fitted coefficient is the marginal effect of one feature on the predicted target, holding the others fixed. Report confidence intervals alongside the point estimates.
• $R^2$ is the fraction of variance explained; it rises whenever a feature is added (use adjusted $R^2$ to compare models) and does not measure out-of-sample performance on its own.
• The four classical assumptions are linearity, independence of residuals, constant variance (homoscedasticity), and approximate normality of residuals; diagnose each with a residual plot before trusting the model.
• The worked sales-versus-ad-spend example fits $\hat{y} = 17.36 + 13.48\,x$ with $R^2 \approx 0.997$, RMSE of about $830, and a forecast of about $85k for the next week given a $5k ad budget.