Linear Regression from Scratch
Fit a straight line to data three ways: the normal equation, gradient descent, and scikit-learn. All three give the same answer, shown with real numbers, plus MSE, R-squared and the traps to know.
Linear regression fits a line y = w x + b that minimises the mean squared error between predictions and actual values, known as ordinary least squares. The normal equation solves it exactly and gradient descent iteratively, reaching the same coefficients. Interviewers expect R-squared explained: the share of variance explained relative to always predicting the mean.
Linear regression fits a line y = w x + b that minimises the squared error between predictions and the real values. You can solve it exactly with the normal equation, or iteratively with gradient descent, and both reach the same coefficients. R-squared reports how much of the variation in y the line explains.
- Ordinary least squares linear regression chooses the slope and intercept that minimise the mean squared error of the residuals.
- The normal equation solves (X^T X) beta = X^T y exactly; using solve is faster and more accurate than computing a matrix inverse.
- On 50 points from y = 3x + 2 with noise, the normal equation, 200 steps of gradient descent and scikit-learn all gave intercept 1.537 and slope 3.024.
- The example fit had an MSE of 2.294 and an R-squared of 0.969, meaning the line explains about 97% of the variation in y.
- R-squared on training data always rises as features are added, even useless ones, so linear regression fits should be checked on held-out data.
- Scaling features before gradient descent turns a long narrow loss valley into a round bowl, so a single learning rate works.
01The model and what "best fit" means
Given inputs x and a numeric target y, the simplest model is a straight line: prediction = w times x plus b. The slope w says how much y changes when x goes up by one. The intercept b is the prediction when x is zero.
To choose w and b, we need a score for "how wrong is this line". The standard choice is the mean squared error (MSE): take each residual (actual minus predicted), square it, and average. Squaring makes every error positive and punishes big misses harder than small ones. The best line is the one with the lowest MSE. This is called ordinary least squares.
02The data used below
import numpy as np
rng = np.random.default_rng(42)
x = rng.uniform(0, 10, 50)
y = 3 * x + 2 + rng.normal(0, 2, 50) # true slope 3, intercept 2, noise sd 2Because we generated the data, we know the right answer is slope 3 and intercept 2. Real data never tells you that, but it lets us check the methods. We should get close to those values but not exactly, because the noise pulls the fit around.
03Method 1: the normal equation
Stack a column of ones next to x so the intercept is handled like any other weight. Then the best weights satisfy (X^T X) beta = X^T y.
X = np.column_stack([np.ones_like(x), x])
beta = np.linalg.solve(X.T @ X, X.T @ y)
print(beta.round(3)) # [1.537 3.024]The intercept is 1.537 and the slope is 3.024, close to the true 2 and 3. Use solve rather than computing the matrix inverse, which is slower and less accurate.
This is exact, but the cost grows quickly with the number of features (roughly cubic), and X^T X can be unsolvable when features are perfectly correlated. For millions of rows or thousands of features, people use iterative methods.
04Method 2: gradient descent
The MSE surface for linear regression is a smooth bowl. At any point the gradient tells you which way is uphill, so you step the opposite way by a small amount (the learning rate), and repeat.
xs = (x - x.mean()) / x.std() # scale the feature first
Xs = np.column_stack([np.ones_like(xs), xs])
w = np.zeros(2)
lr = 0.1
for _ in range(200):
grad = 2 / len(y) * Xs.T @ (Xs @ w - y)
w -= lr * grad
slope = w[1] / x.std() # convert back to original units
intercept = w[0] - slope * x.mean()
print(round(intercept, 3), round(slope, 3)) # 1.537 3.024After 200 steps it lands on the same 1.537 and 3.024. Two practical points. Scaling the feature first makes the bowl round instead of a long narrow valley, so a single learning rate works. And the learning rate matters: too large and the steps overshoot and diverge, too small and it takes forever. Gradient descent is the engine behind neural networks too, which is why it is worth knowing in this simple setting.
05Checking the fit: MSE and R-squared
pred = X @ beta
mse = ((y - pred) ** 2).mean()
r2 = 1 - ((y - pred) ** 2).sum() / ((y - y.mean()) ** 2).sum()
print(round(mse, 3), round(r2, 3)) # 2.294 0.969
from sklearn.linear_model import LinearRegression
m = LinearRegression().fit(x.reshape(-1, 1), y)
print(round(m.intercept_, 3), round(m.coef_[0], 3)) # 1.537 3.024MSE is 2.294, so a typical miss is about 1.5 units (the square root, RMSE). R-squared is 0.969: the line explains about 97% of the variation in y. R-squared compares your model with a model that always predicts the mean, so 0 means no better than that and 1 is a perfect fit.
scikit-learn gives identical coefficients, as expected, since it solves the same problem. Library code is what you use in practice. Writing it once by hand is what makes the library's behaviour make sense.
06Assumptions and traps
Linear regression assumes the relationship is roughly linear, the errors are independent with similar spread across the range, and the features are not perfectly correlated with each other. Plot residuals against predictions. A curve or funnel shape means one of those is broken.
Outliers have a big pull because errors are squared. One extreme point can tilt the whole line.
Extrapolation is risky. The line says nothing reliable outside the range of x you trained on.
With many features and few rows, the model can fit noise. R-squared on the training data always rises as you add features, even useless ones. Check performance on held-out data (see the topics on overfitting and cross-validation). Regularised variants, ridge and lasso, add a penalty on large weights to control this.
Correlation is not causation either: a coefficient describes association in your data, not what would happen if you changed x.
07Common interview questions
Q: What does the loss function in linear regression measure?
A: The mean squared error between predictions and actual values.
Q: Normal equation or gradient descent?
A: The normal equation is exact and fine for a modest number of features. Gradient descent scales to huge data and is the basis for models that have no closed form.
Q: Why scale features for gradient descent?
A: Features on different scales make the loss surface a narrow valley, which forces a tiny learning rate and slow convergence.
Q: What does R-squared of 0.97 mean?
A: The model explains about 97% of the variance in the target, relative to predicting the mean. It does not prove the model is correct or will generalise.
Q: Effect of outliers?
A: They distort the fit because squared error weights them heavily. Remove them if they are errors, use robust regression, or transform the target.
Q: What is multicollinearity?
A: Features that are strongly correlated with each other. Predictions may still be fine, but individual coefficients become unstable and hard to interpret.