Linear regression from scratch with NumPy

Libraries hide how a model learns. Linear regression is the best place to lift the lid: it is the simplest model that actually trains, and every idea in it — a prediction, a loss, a gradient, a learning rate — shows up again in neural networks. We will build it in about twenty lines and then check it against scikit-learn.

The goal

Given inputs x and targets y, find a line:

prediction = w * x + b

w is the slope (the weight) and b is the intercept (the bias). Training means finding the w and b that make predictions closest to the real values.

Make some data

import numpy as np

rng = np.random.default_rng(42)

X = rng.uniform(0, 10, size=50)
y = 3 * X + 2 + rng.normal(0, 1.5, size=50)   # true line: slope 3, intercept 2, plus noise

We know the true answer is roughly w = 3, b = 2, so we can tell whether the algorithm works.

Step 1 — Measure how wrong we are

We need one number that says how bad a line is. The standard choice is mean squared error:

def loss(w, b, X, y):
    pred = w * X + b
    return ((pred - y) ** 2).mean()

Squaring does two jobs: it makes every error positive, and it punishes large mistakes more than small ones. A perfect line scores 0.

Picture loss as a bowl. Every (w, b) pair is a point on the bowl; the bottom is the best line. Training is rolling downhill.

Step 2 — Which way is downhill?

The gradientThe direction in which a function increases fastest; moving the opposite way lowers it tells us how the loss changes when we nudge each parameter. For mean squared error:

d(loss)/dw = (2/n) * sum( (pred - y) * x )
d(loss)/db = (2/n) * sum(  pred - y      )

Intuition: if predictions are too high on average (pred - y positive), the gradient is positive and we should decrease the parameter. The x in the w gradient means points far from zero have more influence over the slope.

Step 3 — Gradient descent

Repeat: compute predictions, compute gradients, take a small step downhill.

w, b = 0.0, 0.0
lr = 0.01          # learning rate: step size
n = len(X)

for epoch in range(1000):
    pred = w * X + b
    error = pred - y

    dw = (2 / n) * (error * X).sum()
    db = (2 / n) * error.sum()

    w -= lr * dw
    b -= lr * db

    if epoch % 200 == 0:
        print(epoch, round(loss(w, b, X, y), 3), round(w, 3), round(b, 3))

print(f"w = {w:.3f}, b = {b:.3f}")

Typical output with this seed (each line is printed after that epoch's update):

0 388.693 2.372 0.356
200 1.346 3.094 1.142
400 1.301 3.051 1.431
600 1.292 3.032 1.556
800 1.291 3.024 1.611
w = 3.021, b = 1.635

The loss falls from the hundreds to about 1.3 within the first couple of hundred steps, then crawls. We recovered a slope near 3 and an intercept near 2. The intercept is still creeping toward its final value after 1000 steps; more on that below.

The learning rate

Learning rateHow big a step gradient descent takes each update is the one setting you must get right:

• Too small: works, but needs thousands of extra steps.
• Too large: the steps overshoot the bottom and the loss grows. If you see the loss turning into nan or inf, lower it.

Try lr = 0.1 on this unscaled data and watch it blow up. That is a good five-minute experiment.

Why the intercept is slow: scale your inputs

X ranges from 0 to 10, so the gradient for w is much larger than for b, and one learning rate cannot suit both. Standardising the input fixes it:

X_scaled = (X - X.mean()) / X.std()

w2, b2 = 0.0, 0.0
for epoch in range(200):
    error = w2 * X_scaled + b2 - y
    w2 -= 0.1 * (2 / n) * (error * X_scaled).sum()
    b2 -= 0.1 * (2 / n) * error.sum()

# convert back to the original units
w_orig = w2 / X.std()
b_orig = b2 - w2 * X.mean() / X.std()
print(w_orig, b_orig)   # about 3.018 and 1.653

A bigger learning rate and fewer steps reach the exact answer. Feature scaling is not cosmetic; it changes how fast training converges.

Check against scikit-learn

from sklearn.linear_model import LinearRegression

model = LinearRegression().fit(X.reshape(-1, 1), y)
print(model.coef_[0], model.intercept_)   # about 3.018 and 1.653

Same answer. scikit-learn solves it directly with linear algebra rather than by stepping, which is why it is instant. For huge datasets and for neural networks the direct solution is not available, which is why gradient descent matters.

You can reproduce the direct solution yourself:

A = np.column_stack([X, np.ones(n)])
w_exact, b_exact = np.linalg.lstsq(A, y, rcond=None)[0]

How good is the fit?

R-squaredThe share of the variation in y that the model explains; 1 is perfect, 0 is no better than predicting the mean gives a quick summary:

ss_res = ((y - (w * X + b)) ** 2).sum()
ss_tot = ((y - y.mean()) ** 2).sum()
print(1 - ss_res / ss_tot)   # about 0.98

Remember that this score is measured on the data we trained with. To judge a model honestly you need a held-out set — see train/test split and data leakage.

From here to neural networks

Everything above is the whole training loop of deep learning in miniature:

1. Predict
2. Measure the loss
3. Compute gradients
4. Update the parameters

A neural network just has more parameters and needs backpropagation to compute the gradients through many layers.

FAQ

What is linear regression in simple terms?

Linear regression fits a straight line to data so you can predict a number from one or more inputs. Training finds the slope and intercept that minimise the average squared error.

What is gradient descent?

Gradient descent is an optimisation method that repeatedly moves the parameters in the direction that lowers the loss. The learning rate controls the step size.

Why does my loss become nan?

The learning rate is too high, so each step overshoots and the values grow until they overflow. Lower the learning rate or standardise your inputs.

Should I use gradient descent or the closed-form solution?

For small linear models the closed-form solution (np.linalg.lstsq or scikit-learn) is faster and exact. Gradient descent is the method that scales to huge data and to neural networks.

Next reading on Sythra Articles

• Backpropagation worked out with real numbers
• Your first machine learning model with scikit-learn
• How to detect overfitting