Linear Regression from Scratch in NumPy: A Complete, Correct Walkthrough

Why build it by hand?

You can fit a line in one line with scikit-learn. So why implement linear regression yourself? Because the moment you write the gradient by hand, three things stop being magic: how a model learns, why the learning rate matters, and what "minimising a loss" actually does step by step. Every deep learning optimiser you will ever use is a descendant of the tiny loop we are about to write.

This post is a complete walkthrough. Every code block runs top to bottom with just NumPy (plus Matplotlib for the plots). By the end you will have two working implementations, gradient descent and the closed-form solution, and you will know when to reach for each.

The model in one equation

Linear regression assumes the target y is a straight-line function of the input x, plus some noise we cannot explain:

y = w · x + b

• w (the weight, or slope) says how much y changes when x goes up by one.
• b (the bias, or intercept) is the value of y when x is zero.

Training means finding the w and b that make the line pass as close as possible to the data. "As close as possible" needs a precise definition, and that is the cost function.

Step 0: make some data

We generate data that truly follows y = 4 + 3x and then add noise, so we know the answer we are hoping to recover.

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
X = 2 * np.random.rand(100, 1)          # 100 points in [0, 2)
y = 4 + 3 * X + np.random.randn(100, 1) # true line: b=4, w=3, plus noise

plt.scatter(X, y, s=15)
plt.xlabel("X"); plt.ylabel("y"); plt.title("Synthetic data")
plt.show()

The true parameters are w = 3 and b = 4. A correct implementation should land near those values.

Step 1: the cost function (MSE)

We measure error with Mean Squared Error: average the squared gap between what the model predicts and the truth.

MSE = (1 / n) · Σ (ŷ - y)²

Squaring does two useful things: it makes every error positive (over- and under-shooting both count), and it punishes big misses much harder than small ones. We add a factor of 1/2 out front so the derivative comes out clean later. It changes the value of the cost but not the location of its minimum, so the best w and b are unaffected.

To handle the bias b with the same matrix math as the weights, we prepend a column of ones to X. Then the bias is just "the weight on a feature that is always 1", and prediction becomes a single dot product.

# Prepend a column of ones so theta[0] is the bias, theta[1] the slope
X_b = np.c_[np.ones((X.shape[0], 1)), X]   # shape (100, 2)

def compute_cost(X_b, y, theta):
    n = len(y)
    predictions = X_b.dot(theta)           # shape (n, 1)
    return (1 / (2 * n)) * np.sum((predictions - y) ** 2)

Step 2: the gradient, derived

Gradient descent needs the slope of the cost with respect to each parameter. Take the derivative of the MSE above with respect to theta:

∂J/∂θ = (1 / n) · Xᵀ · (Xθ - y)

That is the whole thing. Xθ - y is the vector of errors, and multiplying by Xᵀ spreads each error back onto the feature that caused it. A positive gradient means "the cost goes up if you increase this parameter", so we step in the opposite direction.

def gradient_descent(X_b, y, theta, alpha, n_iters):
    n = len(y)
    history = []
    for _ in range(n_iters):
        errors = X_b.dot(theta) - y            # (n, 1)
        gradient = (1 / n) * X_b.T.dot(errors) # (2, 1)
        theta = theta - alpha * gradient       # step downhill
        history.append(compute_cost(X_b, y, theta))
    return theta, history

Notice the fix that trips up most from-scratch tutorials: theta is returned and reassigned, and the bias column is baked into X_b, so both parameters learn.

Step 3: train it

theta = np.random.randn(2, 1)   # random start: [bias, slope]
alpha = 0.1                     # learning rate
n_iters = 1000

theta, history = gradient_descent(X_b, y, theta, alpha, n_iters)
print(f"Learned bias  b = {theta[0, 0]:.3f}")   # ~4.2
print(f"Learned slope w = {theta[1, 0]:.3f}")   # ~2.8

You should see values close to b = 4 and w = 3. They are not exact, and they should not be: the data has noise, so the best-fit line is the one that balances all the scatter, not the invisible "true" line.

Watch it converge

Plotting the cost history is the single most useful debugging habit in all of machine learning. A healthy run drops fast and then flattens.

plt.plot(history)
plt.xlabel("Iteration"); plt.ylabel("Cost (MSE/2)")
plt.title("Gradient descent convergence")
plt.show()

If the curve shoots up to infinity, your learning rate is too high. If it crawls down and never flattens, it is too low. This one plot tells you which.

Step 4: the shortcut, the normal equation

Linear regression is special: you can solve for the exact best parameters in closed form, no iteration needed. This is the normal equation:

θ = (Xᵀ X)⁻¹ Xᵀ y

theta_exact = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
print(theta_exact.ravel())   # ~[4.2, 2.8], matches gradient descent

Both methods land in the same place, which is a good sanity check. So why bother with gradient descent at all?

| | Normal equation | Gradient descent | |---|---|---| | Iteration | None, one formula | Many small steps | | Cost | Inverts a matrix, ~O(features³) | O(iterations · features) | | Scales to millions of features? | No, the inverse blows up | Yes | | Needs a learning rate? | No | Yes, must be tuned | | Works for neural nets? | No | Yes, this is the whole idea |

Use the normal equation for small, tidy problems. Use gradient descent (and its descendants like Adam) for everything large or non-linear, which is most of modern ML.

In practice, prefer np.linalg.lstsq(X_b, y, rcond=None) over inverting XᵀX directly. It is more numerically stable when features are correlated, which makes XᵀX nearly singular.

Step 5: evaluate the fit

A learned theta is not "good" just because training finished. Two numbers tell you whether the line is actually useful.

predictions = X_b.dot(theta)

mse = np.mean((predictions - y) ** 2)
ss_res = np.sum((y - predictions) ** 2)
ss_tot = np.sum((y - y.mean()) ** 2)
r2 = 1 - ss_res / ss_tot

print(f"MSE       : {mse:.3f}")
print(f"R squared : {r2:.3f}")   # ~0.7 on this noisy data

MSE is the average squared error in the units of y squared. R squared is friendlier: it is the fraction of the variation in y your line explains, from 0 (no better than guessing the mean) to 1 (perfect). Around 0.7 here is honest for data this noisy.

One trap: feature scaling

Our X lived in [0, 2), so a learning rate of 0.1 worked fine. Give gradient descent features on wildly different scales (say, "age" from 0 to 100 next to "income" from 0 to 100000) and it will zig-zag painfully or diverge. The fix is to standardise every feature before training:

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

The normal equation does not care about scaling, one more reason it is convenient for small problems. Gradient descent almost always wants it.

Putting it together

You now have the full arc of a supervised learning algorithm in about 30 lines: define a model, measure its error with a cost function, derive the gradient, step downhill until the cost flattens, then check the fit with real metrics. Swap the straight line for a network of them and the cost for cross-entropy, and this exact loop trains a neural network.

Try it on your own data next. Load a CSV, pick one numeric column to predict from another, prepend the ones column, and run the same loop. When the cost curve flattens and R squared climbs, you will have built, not imported, a model that learns.