Gradient Descent, Line by Line: Twelve Lines of Python and the Calculus Behind Them

This function draws the best straight line through a cloud of points. It starts with a flat line at zero, and 900 times over it asks “which way is less wrong?” and nudges the line that way.

Nothing in it is more complicated than multiply, add and average. The rest of this post takes it apart one line at a time.

def gradient_descent_training(X, y, alpha=0.0001, epochs=900):
    m, b = 0.0, 0.0
    for epoch in range(epochs):
        n = len(X)
        y_pred = m * X + b
        error = y_pred - y
        dm = (2/n) * np.dot(error, X)
        db = (2/n) * np.sum(error)
        m -= alpha * dm
        b -= alpha * db
    return m, b

What goes in, what comes out

It returns two numbers, m and b, the slope and intercept of the line it settled on. Together they are the trained model — to predict a new value, you compute m * x_new + b.

Name	What it is	Example
X	The inputs, as a NumPy array — one number per data point	[0, 0.17, 0.34, ...]
y	The true outputs you observed, same length as X	[-0.39, -2.66, 1.35, ...]
alpha	The learning rate: how big a step to take each time	0.0001
epochs	How many times to repeat the whole process	900

That is worth sitting with for a second. A trained machine learning model is not a mysterious object. Here it is literally two numbers, and the function’s entire job is to find good values for them.

Line by line

m, b = 0.0, 0.0

The starting guess. A slope of 0 and an intercept of 0 is the horizontal line through the origin — almost certainly a terrible fit, and that is fine. The whole point is that the function improves it from wherever it starts.

for epoch in range(epochs):

Do everything below 900 times. One pass is called an epoch. Notice that epoch is never used inside the loop; it is just a counter.

y_pred = m * X + b

What the current line predicts, for every data point at once.

This is the line that surprises people coming from plain Python. X is a NumPy array, so m * X + b does the arithmetic to every element, and y_pred comes back as a whole array of predictions. No loop needed. If X has 60 points, this makes 60 predictions in one line.

error = y_pred - y

How wrong each prediction is, one number per data point. Positive means the line sits above that point; negative means below.

The order matters. It must be prediction minus truth. Write it the other way round and every step later in the function goes the wrong direction — the line will run away from the data instead of toward it. This is the single most common bug in code like this.

dm and db

The two lines that do the actual thinking. They get a section to themselves below.

m -= alpha * dm and b -= alpha * db

The nudge. Move each parameter a little way in the direction that reduces error. Also gets its own section.

return m, b

Hand back the two numbers the loop arrived at.

The two gradient lines

dm = (2/n) * np.dot(error, X)
db = (2/n) * np.sum(error)

These answer one question: if I nudged m (or b) up a tiny bit, would the total error get better or worse, and by how much?

That is all a gradient is. A number saying which way is uphill, and how steep the hill is.

Take db first — it is the easy one

np.sum(error) adds up all the errors. Divided by n that would be the average error, and the 2 is a leftover from the calculus. So:

• Predictions too high on average → db is positive → the update pushes b down.
• Predictions too low on average → db is negative → the update pushes b up.

That is exactly what common sense says to do with the height of a line. The calculus agrees with you.

Now dm

np.dot(error, X) is the dot product: multiply the two arrays element by element, then add up the results. It is identical to np.sum(error * X), just faster and more idiomatic.

np.dot([2, -1, 3], [10, 20, 30])   # = 2*10 + (-1)*20 + 3*30 = 90

So dm is the average error weighted by x. Each point’s error is scaled by how far out along the x-axis it sits.

Why weight by x? Because changing the slope pivots the line around the y-axis. It barely moves the line near x = 0, but swings it a long way out at x = 100. Points far from the origin therefore have far more say in whether the slope is right, and multiplying by X is precisely how that gets accounted for.

The calculus, in full

This section derives dm and db from scratch. You need one rule from calculus (the chain rule) and the idea of a partial derivative. If you skip it, take the two plain-English descriptions above as definitions — nothing later depends on the derivation.

Step 1: write down what is being minimised

The loss is the mean squared error:

L(m, b) = (1/n) * sum over i of ( m*x_i + b - y_i )**2

Read the notation carefully, because this is where most confusion starts.

L is a function of m and b, not of x. The data x_i and y_i are fixed — they were measured, they are not going to change. The only things free to vary are the two parameters. So when we differentiate, we differentiate with respect to m and b, and every x_i and y_i is treated as a constant.

This is the mental flip that makes gradient descent make sense. You are used to functions where x is the variable and m is a fixed coefficient. Here it is the other way around.

Step 2: name the inside

Each term in that sum is something squared. Give the something a name:

u_i = m*x_i + b - y_i          # the error on point i
L   = (1/n) * sum of u_i**2

So u_i is exactly what the code calls error. The loss is the average of u_i squared.

Step 3: differentiate with respect to m

The sum rule lets us handle one term at a time, so focus on u_i**2. The chain rule says: to differentiate a function of a function, differentiate the outside, keep the inside, then multiply by the derivative of the inside.

d/dm of (u_i**2)  =  2 * u_i * (du_i/dm)

Now what is du_i/dm? Differentiate u_i = m*x_i + b - y_i with respect to m, holding everything else constant:

d/dm of (m*x_i)  =  x_i        # x_i is a constant multiplier
d/dm of (b)      =  0          # b does not depend on m
d/dm of (-y_i)   =  0          # data is constant

so du_i/dm = x_i

Substituting back:

d/dm of (u_i**2)  =  2 * u_i * x_i

Put the sum and the 1/n back on:

dL/dm  =  (2/n) * sum over i of ( u_i * x_i )
       =  (2/n) * sum over i of ( error_i * x_i )

That sum of products, element by element, is a dot product — which is the line in the code:

dm = (2/n) * np.dot(error, X)

Step 4: differentiate with respect to b

Identical, except the inside derivative changes:

d/db of (m*x_i)  =  0
d/db of (b)      =  1
d/db of (-y_i)   =  0

so du_i/db = 1

Which gives:

dL/db  =  (2/n) * sum over i of ( u_i * 1 )
       =  (2/n) * sum over i of error_i

db = (2/n) * np.sum(error)

The entire difference between the two formulas is that du/dm = x_i while du/db = 1. That one factor of x_i is why the slope gradient is weighted by x and the intercept gradient is not — and, further down, why the two parameters converge at such different speeds.

Why the gradient points uphill

A derivative is a rate of change: dL/dm is how much L goes up per unit increase in m.

So if dL/dm is positive, increasing m increases the loss. You want the opposite, so you decrease m. If dL/dm is negative, increasing m decreases the loss, so you increase m. In both cases you move opposite to the sign of the derivative — which is precisely what m -= alpha * dm does. That single line is the whole reason gradient descent works, and this is the argument behind it.

The pair (dL/dm, dL/db) taken together is the gradient, written as an upside-down triangle in textbooks. Read as a vector on the (m, b) plane, it points in the direction of steepest increase in loss. Negate it and you have the direction of steepest decrease. That vector is what the two update lines step along.

Verify it numerically — the trick worth knowing

You never have to trust a derivation. A derivative is the limit of a slope, so you can approximate it by just nudging the parameter and seeing what the loss does:

dL/dm  ≈  ( L(m+h, b) - L(m-h, b) ) / (2h)     for small h, e.g. 1e-5

Run this against the formula on the three-point dataset below, at m = 0.5, b = 0.5:

Quantity	From the formula	From nudging
dL/dm	-16.000000	-16.000000
dL/db	-7.000000	-7.000000

They agree to six decimals. In real work this is called gradient checking, and it is how people debug hand-written backpropagation. If your analytic gradient and your numerical gradient disagree, the analytic one has a bug — and it is usually a dropped factor or a sign.

The 2/n constant, and what happens without it

The 2 comes from differentiating a square. The 1/n comes from the loss being a mean rather than a sum.

Neither is load-bearing. Both are positive constants, so dropping them scales every gradient by the same factor and does not change which direction is downhill — it only changes how far each step goes. Halving the 2 and doubling alpha gives identical behaviour. Some textbooks define the loss with a 1/(2n) out front precisely so the 2 cancels and the derivative comes out clean. That is cosmetic, and if you see it elsewhere, that is why.

The 1/n does do one useful thing though: it keeps the gradient the same size regardless of how many data points you have. Without it, doubling your dataset would double your gradients, and a learning rate tuned on 100 points would explode on 10,000.

Why the minimum exists at all

Differentiate a second time:

d2L/dm2 = (2/n) * sum of x_i**2

A sum of squares is always positive, so the second derivative is always positive, so the loss curve is bowl-shaped. In two parameters the same argument makes it a genuine bowl: convex, with exactly one lowest point and no false valleys to get trapped in.

This matters more than it might look. It means gradient descent on this problem cannot fail for structural reasons. If it fails, it is your alpha or your code.

The shortcut it makes possible

At the bottom of a bowl the ground is flat, so both partial derivatives are zero. Setting them to zero gives two equations in two unknowns:

(2/n) * sum( (m*x_i + b - y_i) * x_i ) = 0
(2/n) * sum(  m*x_i + b - y_i        ) = 0

Solve them and you get the closed-form least squares answer directly, with no loop at all:

m = sum( (x_i - x_mean)*(y_i - y_mean) ) / sum( (x_i - x_mean)**2 )
b = y_mean - m * x_mean

So why iterate? Because this shortcut exists only for a handful of models. The moment the loss is anything harder — a logistic regression, a neural network — those equations have no closed-form solution. Gradient descent needs nothing except a loss you can differentiate, which is why it scaled all the way to modern machine learning and the formula did not.

It is still worth computing both on a small problem. Watching a loop grind its way to the number algebra produces instantly is the most convincing demonstration that the loop is doing real mathematics.

Why minus, and what alpha does

m -= alpha * dm
b -= alpha * db

The minus sign

The gradient points uphill — toward more error. We want less error. So we go the other way, and subtracting is how.

Change -= to += and the function does the opposite of its job: it climbs, the error grows every epoch, and the numbers eventually overflow. Try it once, deliberately. It is a good way to convince yourself the sign is doing real work.

Alpha, the learning rate

The gradient says which way. Alpha says how far. This matters more than almost anything else in the function:

Alpha	What happens
Too small	Correct direction, but the line creeps. You run out of epochs before you arrive.
About right	Drops fast, then settles near the best answer.
Too large	Each step overshoots the target and lands further past it than the last. The error explodes and you end up with inf, then nan.

The blindfold picture is the usual one, and it is a good one. You are standing on a hillside in fog, trying to reach the bottom. You cannot see the valley, but you can feel which way the ground slopes. So you feel the slope, take a step downhill, and repeat. The gradient is the feel of the ground; alpha is the length of your stride.

There is no universally correct alpha. It depends on the data, and finding a good one is part of the job.

One epoch, with real numbers

Three data points that sit exactly on the line y = 2x + 1, so you know what the answer should be. Use alpha = 0.1.

X = np.array([1., 2., 3.])
y = np.array([3., 5., 7.])

Epoch 1, starting from m = 0, b = 0:

Step	Value
y_pred = 0*X + 0	[0, 0, 0]
error = y_pred - y	[-3, -5, -7] — every prediction is too low
np.dot(error, X)	-3(1) + -5(2) + -7(3) = -34
dm = (2/3)(-34)	-22.67
db = (2/3)(-15)	-10.0
m = 0 - 0.1(-22.67)	2.267
b = 0 - 0.1(-10.0)	1.0

Both gradients came out negative, because every prediction was too low. Subtracting a negative pushes both parameters up — which is right. After one epoch the line is already y = 2.27x + 1.0, against a true y = 2x + 1.

What follows:

Epoch	m	b	Mean squared error
start	0.000	0.000	27.67
1	2.267	1.000	0.33
2	2.018	0.893	0.0053
3	2.044	0.908	0.0013
900	2.000	1.000	~0

Notice the overshoot at epoch 1 and the correction back at epoch 2. That wobble is normal — it is the line rocking past the answer and settling back, and a smaller alpha would make it gentler.

Also notice how much of the work happens in the first step: the error falls from 27.67 to 0.33 immediately, then spends 900 epochs polishing. Nearly every training curve you will ever see has this shape.

The whole thing as a picture

        ┌─────────────────────────┐
        │  Start: m = 0, b = 0    │
        └────────────┬────────────┘
                     ↓
        ┌─────────────────────────┐
   ┌──→ │  Predict                │
   │    │  y_pred = m*X + b       │
   │    └────────────┬────────────┘
   │                 ↓
   │    ┌─────────────────────────┐
   │    │  Measure                │
   │    │  error = y_pred - y     │
   │    └────────────┬────────────┘
   │                 ↓
   │    ┌─────────────────────────┐
   │    │  Ask which way is down  │
   │    │  dm, db                 │
   │    └────────────┬────────────┘
   │                 ↓
   │    ┌─────────────────────────┐
   │    │  Step that way          │
   │    │  m -= alpha*dm          │
   │    │  b -= alpha*db          │
   │    └────────────┬────────────┘
   │                 ↓
   │          ⟨ done 900 epochs? ⟩
   └───── no ────────┤
                    yes
                     ↓
              Return m, b

Four steps in a circle: predict, measure, ask which way, step. Nine hundred laps and you have a trained model.

This loop is the reason the function is worth understanding in detail. It is not specific to straight lines. Swap the model for a neural network with a billion parameters and the four steps are unchanged — only the gradient calculation in step three gets harder, which is the problem backpropagation solves. Every deep learning framework you will ever use is an industrial-strength version of what is written here.

Four things to watch for in this code

1. The default alpha is very small. On data where x runs from 0 to 10, alpha=0.0001 with 900 epochs is not enough to finish. In a test run the slope reached 2.33 against a best-fit 2.54, and the intercept only got to 0.29 against a best-fit of -1.07 — it had barely left its starting value. Reaching the answer took about 50,000 epochs. At alpha=0.01, 900 epochs was plenty. This is not a bug in the code, but it does mean that running it as-is may look like it “doesn’t work”. Try a few values of alpha and watch what changes.

2. The intercept lags behind the slope. dm is multiplied by X, so when x ranges over 0 to 10 it is roughly ten times larger than db. Both get the same alpha, so whatever step size suits the slope is far too small for the intercept. The standard fix is to rescale X to have mean 0 and standard deviation 1 before training.

3. n = len(X) is inside the loop. It recomputes the same number 900 times. Harmless, but it belongs above the loop.

4. Nothing records the error. The function returns m and b but keeps no history, so you cannot plot a loss curve or tell whether it converged, wobbled or exploded. Appending the mean squared error to a list each epoch costs one line and makes the function much easier to debug.

Check your understanding

1. Why is error = y_pred - y rather than y - y_pred?
2. np.dot(error, X) uses X but np.sum(error) does not. Why the difference?
3. Every prediction is too high. What are the signs of dm and db, and which way do m and b move?
4. You change -= to +=. Describe what the printed error does over 900 epochs.
5. The function returns m = 0.02, b = 0.01 after 900 epochs on data that clearly slopes upward steeply. What is the most likely cause?