Linear Regression From Scratch in Python: How a Model Actually Learns

📎 This article includes 1 downloadable practice file ↓

Every neural network — including the biggest LLMs — learns the same way a straight line does: make a guess, measure the error, nudge the parameters to reduce it, repeat. Linear regression is the smallest version of that loop, so it is the best place to see it clearly.

The data

We create data where we know the true answer: y = 3x + 7 plus some noise. If our code works, it should rediscover 3 and 7.

import numpy as np
rng = np.random.default_rng(1)
x = rng.uniform(0, 10, 100)
y = 3 * x + 7 + rng.normal(0, 2, 100)

The model and the error

The model has two parameters: slope w and intercept b. Its prediction is w*x + b. We measure how wrong it is with the mean squared error: the average of (prediction − actual)².

Gradient descent

The gradient tells us which way to move each parameter to reduce the error. For mean squared error the formulas are short:

w, b = 0.0, 0.0          # start with a bad guess
lr = 0.01                # learning rate: size of each nudge

for epoch in range(3000):
    pred = w * x + b
    err = pred - y
    loss = (err ** 2).mean()
    dw = 2 * (err * x).mean()       # slope of the loss with respect to w
    db = 2 * err.mean()             # ... and with respect to b
    w -= lr * dw
    b -= lr * db
    if epoch % 500 == 0:
        print(f"epoch {epoch:4d}  loss {loss:8.2f}  w {w:.2f}  b {b:.2f}")

print("learned:", round(w, 2), round(b, 2))

The loss drops quickly at first and then slowly. w lands near 3 and b near 7 — not exactly, because of the noise we added.

The learning rate matters

Try all three. Tuning the learning rate is still one of the first things people do when training large models.

Check against the libraries

print(np.polyfit(x, y, 1))                 # [slope, intercept], exact solution

from sklearn.linear_model import LinearRegression
m = LinearRegression().fit(x.reshape(-1, 1), y)
print(m.coef_[0], m.intercept_)

For linear regression there is an exact formula, so libraries do not need gradient descent. We used it anyway because it is the method that scales: a neural network has no exact formula, only the loop.

Several inputs

With more features (ad spend, discount, month) the model becomes y = X @ w + b, where w is a vector. The loop is identical:

X = np.column_stack([x, rng.uniform(0, 5, 100)])
y2 = 3 * X[:, 0] - 2 * X[:, 1] + 7 + rng.normal(0, 2, 100)
w, b = np.zeros(2), 0.0
for _ in range(5000):
    err = X @ w + b - y2
    w -= 0.01 * 2 * X.T @ err / len(y2)
    b -= 0.01 * 2 * err.mean()
print(w.round(2), round(b, 2))             # close to [3, -2] and 7

From here to neural networks

A neural network is many of these linear steps stacked, with a simple non-linear function between them. Same gradient descent, more parameters. See the next step: a neural network from scratch in NumPy.