What you’ll learn

• the update rule $\boldsymbol{\theta} \leftarrow \boldsymbol{\theta} - \alpha\nabla J$, derived rather than asserted
• two full iterations computed by hand, matching NumPy to four decimals
• the exact threshold at which a learning rate stops converging and starts exploding
• batch, stochastic and mini-batch descent — the real trade-off between them
• why feature scaling turns a slow ravine into a fast bowl
• learning schedules, and when to stop

Intuition

You are standing on a hillside in fog. You cannot see the valley floor, but you can feel which way the ground slopes under your feet. Take a step downhill. Feel again. Step again.

That is the entire algorithm. The gradient is the direction of steepest ascent, so you move against it. The learning rate is how long a stride you take. Everything else — momentum, Adam, schedules — is refinement of those two sentences.

  flowchart TD
  A["Start from random parameters"] --> B["Compute predictions"]
  B --> C["Compute cost J"]
  C --> D["Compute gradient of J"]
  D --> E["Step against the gradient"]
  E --> F{"Converged?"}
  F -->|no| B
  F -->|yes| G["Return theta"]

The Normal Equation from Multiple Linear Regression already solves linear regression exactly, so why iterate at all? Because inverting an $n \times n$ matrix costs roughly $O(n^3)$. At 100,000 features that is hopeless, and for models with no closed form at all — logistic regression, neural networks, gradient boosting — iteration is the only option. Gradient descent is the workhorse of the entire field; linear regression is just the clearest place to learn it.

The math

The update rule

For a cost $J(\boldsymbol{\theta})$, the gradient collects the partial derivatives:

and each iteration steps against it:

$\alpha$ is the learning rate. It is the only hyperparameter, and it decides everything.

The gradient of MSE

With $J(\boldsymbol{\theta}) = \frac{1}{m}\lVert\mathbf{X}\boldsymbol{\theta} - \mathbf{y}\rVert^2$, differentiating gives a strikingly simple result:

Read it as English: the residual vector, projected back onto the features, scaled by two over the sample count. Every feature gets pushed in proportion to how much it contributed to the current error.

When does it actually converge?

This is the part most tutorials skip. For a quadratic cost with Hessian $\mathbf{H} = \frac{2}{m}\mathbf{X}^\top\mathbf{X}$, gradient descent converges if and only if

where $\lambda_{\max}$ is the largest eigenvalue. Above that threshold each step overshoots by more than it corrects and the parameters diverge geometrically. The condition number $\kappa = \lambda_{\max}/\lambda_{\min}$ then decides how fast convergence happens: a $\kappa$ of 1 is a round bowl and descent goes straight in, while a $\kappa$ of 1,000 is a narrow ravine that forces a zig-zag.

Both facts have the same practical consequence, and it is the reason for the whole preprocessing phase: scale your features.

Worked example by hand

The five-point dataset with the feature centred ($x_c = x - 3$), so the Hessian is diagonal and the arithmetic stays exact. Parameters are $[b, w]$, starting from $[0, 0]$, with $\alpha = 0.2$.

The Hessian, and the stability limit.

Eigenvalues 2 and 4, so descent diverges for any $\alpha \geq 2/4 = 0.5$. Our 0.2 is comfortably inside. The optimum is $[b, w] = [4, 0.6]$.

Iteration 0. With $\boldsymbol{\theta} = [0, 0]$ every prediction is 0, so the residual is $-\mathbf{y}$:

Iteration 1. Predictions are now $1.6 + 0.48x_c$, giving residuals $[-1.36, -2.88, -3.4, -1.92, -2.44]$:

The cost falling.

The slope reaches its final value almost immediately while the intercept takes far longer. That is the condition number in action: the eigenvalue along $w$ is twice the eigenvalue along $b$, so the same $\alpha$ makes twice the progress in that direction.

The learning rate

figureSame start, three learning rates, 70 steps eachmatplotlib

figureThe learning curve tells you immediately which rate is wrongmatplotlib

Reading the plots

A practical recipe: try $\alpha \in \{0.001, 0.01, 0.1, 1\}$, plot all four cost curves on one log axis, and take the largest rate that still descends smoothly.

Notice the steps shrinking on their own as the curve flattens. The learning rate never changes — the gradient shrinks near the minimum, so descent automatically slows down as it arrives.

Why scaling matters

figureThe same optimisation problem, before and after scalingmatplotlib

Unscaled features produce an elongated cost surface, and the negative gradient points across the valley rather than down it. Standardising makes the contours near-circular and lets descent walk straight in — often an order of magnitude fewer iterations, from the same code.

from sklearn.datasets import load_diabetes
from sklearn.linear_model import SGDRegressor
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
 
X, y = load_diabetes(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=42
)
 
scaled = make_pipeline(
    StandardScaler(),
    SGDRegressor(max_iter=2000, tol=1e-4, eta0=0.01, random_state=0),
).fit(X_train, y_train)
 
raw = SGDRegressor(
    max_iter=2000, tol=1e-4, eta0=0.01, random_state=0
).fit(X_train, y_train)
 
print(f"scaled   R^2 = {scaled.score(X_test, y_test):.4f}")  # 0.4767
print(f"unscaled R^2 = {raw.score(X_test, y_test):.4f}")     # 0.4554
print(f"epochs used  = {scaled[-1].n_iter_}")                # 27

from sklearn.datasets import load_diabetes
from sklearn.linear_model import SGDRegressor
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
 
X, y = load_diabetes(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=42
)
 
scaled = make_pipeline(
    StandardScaler(),
    SGDRegressor(max_iter=2000, tol=1e-4, eta0=0.01, random_state=0),
).fit(X_train, y_train)
 
raw = SGDRegressor(
    max_iter=2000, tol=1e-4, eta0=0.01, random_state=0
).fit(X_train, y_train)
 
print(f"scaled   R^2 = {scaled.score(X_test, y_test):.4f}")  # 0.4767
print(f"unscaled R^2 = {raw.score(X_test, y_test):.4f}")     # 0.4554
print(f"epochs used  = {scaled[-1].n_iter_}")                # 27

The scaled pipeline converges in 27 epochs and lands on 0.4767 — indistinguishable from the exact Normal Equation answer of 0.4773. The unscaled version burns all 2,000 epochs, emits a convergence warning, and still ends up worse.

The three variants

How much data does one gradient use?

  flowchart TD
  A["Batch GD
all m rows per step"] --> A2["Smooth, expensive steps
settles exactly on the minimum"]
  B["Stochastic GD
1 random row per step"] --> B2["Cheap, noisy steps
bounces around the minimum forever"]
  C["Mini-batch GD
32 to 256 rows per step"] --> C2["Balanced; vectorises well on GPUs
the default everywhere"]

The stochastic gradient is an unbiased estimate of the batch gradient — right on average, wrong on any individual step. That noise is not purely a cost: it lets SGD escape shallow local minima that would trap batch descent on non-convex problems, which is exactly why deep learning uses it.

Learning schedules

Because SGD never settles, the standard fix is to shrink the learning rate over time:

Early on the steps are long enough to cover ground; later they are short enough to sit still. Decay too fast and it freezes before arriving; too slowly and it keeps bouncing. SGDRegressorSGDRegressor exposes this as learning_rate="invscaling"learning_rate="invscaling" with eta0eta0 and power_tpower_t.

From scratch

import numpy as np
 
 
def batch_gradient_descent(X, y, lr=0.2, n_iter=50):
    """Plain batch gradient descent. X must already include the bias column."""
    m = len(y)
    theta = np.zeros(X.shape[1])
    history = []
 
    for _ in range(n_iter):
        residual = X @ theta - y
        history.append(float((residual**2).mean()))
        gradient = (2 / m) * X.T @ residual
        theta = theta - lr * gradient
 
    return theta, history
 
 
x_centred = np.array([1.0, 2.0, 3.0, 4.0, 5.0]) - 3.0
X = np.c_[np.ones(5), x_centred]
y = np.array([2.0, 4.0, 5.0, 4.0, 5.0])
 
theta, history = batch_gradient_descent(X, y, lr=0.2, n_iter=50)
 
print(theta.round(4))            # [4.  0.6]
print([round(h, 4) for h in history[:4]])
# [17.2, 6.2688, 2.5548, 1.2265]
 
# The stability limit, computed rather than guessed
hessian = (2 / len(y)) * X.T @ X
print(f"max stable lr = {2 / np.linalg.eigvalsh(hessian).max():.3f}")   # 0.500

import numpy as np
 
 
def batch_gradient_descent(X, y, lr=0.2, n_iter=50):
    """Plain batch gradient descent. X must already include the bias column."""
    m = len(y)
    theta = np.zeros(X.shape[1])
    history = []
 
    for _ in range(n_iter):
        residual = X @ theta - y
        history.append(float((residual**2).mean()))
        gradient = (2 / m) * X.T @ residual
        theta = theta - lr * gradient
 
    return theta, history
 
 
x_centred = np.array([1.0, 2.0, 3.0, 4.0, 5.0]) - 3.0
X = np.c_[np.ones(5), x_centred]
y = np.array([2.0, 4.0, 5.0, 4.0, 5.0])
 
theta, history = batch_gradient_descent(X, y, lr=0.2, n_iter=50)
 
print(theta.round(4))            # [4.  0.6]
print([round(h, 4) for h in history[:4]])
# [17.2, 6.2688, 2.5548, 1.2265]
 
# The stability limit, computed rather than guessed
hessian = (2 / len(y)) * X.T @ X
print(f"max stable lr = {2 / np.linalg.eigvalsh(hessian).max():.3f}")   # 0.500

The first four cost values reproduce the hand-computed table exactly, and after 50 iterations the parameters have reached the Normal Equation’s answer of $[4, 0.6]$.

Pitfalls

🧪 Try It Yourself

Exercise 1 – One batch gradient step

Exercise 2 – Watch the cost fall

Exercise 3 – Compute the stability limit

Exercise 4 – Push it past the limit

Exercise 5 – Scaling changes everything for SGD

Recap

• The update rule is $\boldsymbol{\theta} \leftarrow \boldsymbol{\theta} - \alpha\nabla J$, and for MSE the gradient is $\frac{2}{m}\mathbf{X}^\top(\mathbf{X}\boldsymbol{\theta} - \mathbf{y})$.
• Convergence requires $\alpha < 2/\lambda_{\max}(\mathbf{H})$ — on the worked example, 0.5.
• Two hand iterations at $\alpha = 0.2$ take the cost from 17.20 to 6.27 to 2.55, matching NumPy exactly.
• Scaling controls the condition number, which controls how many iterations you pay for.
• Batch is smooth and expensive, stochastic is noisy and cheap, mini-batch is the default.
• Plot cost against iteration on every run; the shape of that curve diagnoses the learning rate immediately.

Exercise 6 – Find the learning rate that breaks it

Next

Continue to Regularization - Ridge and Lasso Regression — add a penalty term to the same cost and watch overfitting recede.