Understanding the Mechanisms beneath Gradient Descent

How does gradient descent find the best parameters for a model? This notebook explores how a linear regression model recovers unknown slope and intercept by minimizing a loss function, what the error surface looks like in 2D and 3D, and how partial derivatives and gradient vectors guide each optimization step.

pytorch
gradient descent
linear regression
SciML
neural-networks
Author

Mei-Chin Pang

Published

April 7, 2026

import torch
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

A neural network, whether shallow or deep, implements a function that maps inputs to expected outputs. This function is learned through a training process that iteratively searches for a set of weights enabling the network to model the variations in the training data.

The simplest such function is a linear mapping from a single input to a single output, represented by the equation of a line with slope m and y-intercept c:

y = mx + c

Varying each parameter produces a different linear model that defines a different input-output mapping. Learning the mapping therefore involves finding the parameter values that result in the minimum error between predicted and target outputs. This error is calculated by a loss function (also called a cost function or error function), and the process of minimizing it is referred to as function optimization.

Gradient descent is the optimization algorithm at the heart of this process. At each step it:

  1. Computes the loss for the current weights.
  2. Calculates the gradient: the direction and magnitude of steepest increase of the loss with respect to each weight.
  3. Updates the weights by a small step opposite to the gradient, pushing them toward lower loss.

Repeating these steps drives the weights toward the minimum of the loss surface, recovering the parameters that best fit the data.

Linear Regression with Gradient Descent

Step 1: Generating the Training Data

The code below creates a synthetic dataset of 20 (x, y) pairs from a known linear relationship, simulating a scenario where we have observed data and want to recover the underlying parameters.

  1. true_m = 2.5, true_c = 4.0: The ground-truth slope and intercept that we will try to recover.

  2. torch.manual_seed(42): Fixes the random seed so the results are reproducible.

  3. X = torch.linspace(0, 10, N).unsqueeze(1): Creates 20 evenly spaced x-values in [0, 10], reshaped to a column vector of shape (20, 1).

  4. noise = torch.randn(N, 1) * 1.5: Draws 20 random values from a standard normal distribution, scaled by 1.5, to simulate measurement noise.

  5. Y = true_m * X + true_c + noise: Computes the noisy target values y = 2.5x + 4 + \epsilon. These are the simulated “observed” outputs our model will try to fit.

  6. w_true = torch.tensor([true_m, true_c]): Stores the true coefficients [m, c] as a tensor for later comparison with the learned values.

# Simulate ground-truth parameters
true_m = 2.5
true_c = 4.0

# Create sample dataset for 20 points
torch.manual_seed(42)

# (20, 1)
N = 20
X = torch.linspace(0, 10, N).unsqueeze(1)           

# Add Gaussian noise
noise = torch.randn(N, 1) * 1.5

# Simulate observed Y values
Y = true_m * X + true_c + noise

for i in range(N):
    print(f"X[{i}] = {X[i].item():6.2f},  Y[{i}] = {Y[i].item():7.2f}")

w_true = torch.tensor(
    [true_m, true_c],
    dtype=torch.float32)

print("True coefficients w (slope m and intercept c):")
print(w_true)
X[0] =   0.00,  Y[0] =    6.89
X[1] =   0.53,  Y[1] =    7.55
X[2] =   1.05,  Y[2] =    7.98
X[3] =   1.58,  Y[3] =    4.79
X[4] =   2.11,  Y[4] =    8.13
X[5] =   2.63,  Y[5] =   12.20
X[6] =   3.16,  Y[6] =   13.10
X[7] =   3.68,  Y[7] =   15.73
X[8] =   4.21,  Y[8] =   15.06
X[9] =   4.74,  Y[9] =   14.81
X[10] =   5.26,  Y[10] =   16.42
X[11] =   5.79,  Y[11] =   18.84
X[12] =   6.32,  Y[12] =   19.44
X[13] =   6.84,  Y[13] =   21.17
X[14] =   7.37,  Y[14] =   22.04
X[15] =   7.89,  Y[15] =   25.03
X[16] =   8.42,  Y[16] =   24.59
X[17] =   8.95,  Y[17] =   25.77
X[18] =   9.47,  Y[18] =   28.89
X[19] =  10.00,  Y[19] =   28.07
True coefficients w (slope m and intercept c):
tensor([2.5000, 4.0000])

Step 2: Preparing Inputs, Weights, and Optimizer

This code sets up everything needed before the training loop begins.

  1. XX = torch.hstack([X, torch.ones_like(X)]): Builds the design matrix of shape (N, 2). Each row for a sample x_i becomes [x_i \;\; 1]. This lets us express the linear model mx + c as a single matrix multiplication \hat{y} = \mathbf{X}\mathbf{w}.
Column index Multiplies Parameter
0 X (the feature) slope m (w_0)
1 ones (bias term) intercept c (w_1)

  1. w_pred = torch.randn(2, 1, requires_grad=True): Initialises the two unknown coefficients with random values. The flag requires_grad=True tells PyTorch to track every operation on w_pred so gradients can be computed later.

  2. optimizer = torch.optim.NAdam([w_pred], lr=0.01): Creates an NAdam optimizer (a variant of Adam with Nesterov momentum) that will update w_pred using the gradients computed by autograd. The learning rate lr=0.01 controls the step size of each update.

# Prepare input as a design matrix of shape (N, 2)
# Each row: [x_i, 1]  so that  y_hat = XX @ w  =  m*x + c
XX = torch.hstack([X, torch.ones_like(X)])

# Initialise the unknown coefficients randomly
# requires_grad=True tells PyTorch to track operations for autograd
w_pred = torch.randn(2, 1, requires_grad=True)

print("Design matrix XX (first 5 rows):")
print(XX[:5])
print("-" * 79)

print("Target tensor Y (first 5 rows):")
print(Y[:5])
print("-" * 79)

# NAdam optimizer
optimizer = torch.optim.NAdam([w_pred], lr=0.01)

print("Coefficients w before training:")
print(w_pred)
print("-" * 79)
Design matrix XX (first 5 rows):
tensor([[0.0000, 1.0000],
        [0.5263, 1.0000],
        [1.0526, 1.0000],
        [1.5789, 1.0000],
        [2.1053, 1.0000]])
-------------------------------------------------------------------------------
Target tensor Y (first 5 rows):
tensor([[6.8904],
        [7.5467],
        [7.9827],
        [4.7891],
        [8.1260]])
-------------------------------------------------------------------------------
Coefficients w before training:
tensor([[-0.7658],
        [-0.7506]], requires_grad=True)
-------------------------------------------------------------------------------

Step 3: Run the Optimization for 5000 Iterations

The code below runs the training loop for 5000 iterations. Each iteration performs four actions:

  1. optimizer.zero_grad(): Clears the gradients from the previous iteration. PyTorch accumulates gradients by default, so without this reset the gradients would keep growing and the optimizer would take incorrect steps.

  2. y_pred = XX @ w_pred: Computes the predicted \hat{y} values via matrix multiplication \mathbf{X}\mathbf{w}. This is the forward pass.

  3. mse = torch.mean(torch.square(Y - y_pred)): Calculates the Mean Squared Error loss between the true Y and the predictions. This single number measures how far off the current weights are.

  4. mse.backward(): Triggers backpropagation: PyTorch walks the computation graph in reverse and fills w_pred.grad with \frac{\partial\,\text{MSE}}{\partial\mathbf{w}}.

  5. optimizer.step(): The NAdam optimizer uses those gradients to nudge w_pred in the direction that reduces the loss.

After 5000 such updates, w_pred should converge close to the true coefficients [m, c] = [2.5, 4.0], and the final lines print both sets of values plus the remaining error.

# Run optimizer
for i in range(5000):
    optimizer.zero_grad()
    y_pred = XX @ w_pred
    mse = torch.mean(torch.square(Y - y_pred))
    mse.backward()
    optimizer.step()

print("Coefficients w after training:")
print(w_pred)
print("-" * 79)

print("True coefficients w:")
print(w_true)
print("-" * 79)

print("Errors in coefficients:")
print(w_true - w_pred.detach().squeeze())
Coefficients w after training:
tensor([[2.3612],
        [5.0163]], requires_grad=True)
-------------------------------------------------------------------------------
True coefficients w:
tensor([2.5000, 4.0000])
-------------------------------------------------------------------------------
Errors in coefficients:
tensor([ 0.1388, -1.0163])

Step 4: Plotting the Results

The plot below shows the 20 training samples alongside the true line and the line recovered by gradient descent. A close match confirms that the optimizer successfully recovered the coefficients.

  1. x_plot = np.linspace(..., 200): Creates 200 evenly spaced points spanning the data range (with a 1-unit margin on each side) for drawing smooth curves.

  2. y_true_line = true_m * x_plot + true_c: Evaluates the ground-truth line y = 2.5x + 4.

  3. w_np = w_pred.detach().numpy().flatten(): Extracts the learned weights from the PyTorch tensor into a NumPy array. .detach() removes the tensor from the computation graph so NumPy can read it.

  4. y_pred_line = w_np[0] * x_plot + w_np[1]: Evaluates the learned line using the recovered slope (w_0) and intercept (w_1).

  5. The plot overlays:

    • Green scatter points: the 20 noisy training samples.
    • Black solid line: the true relationship y = 2.5x + 4.
    • Red dashed line: the learned line from gradient descent, with its coefficients shown in the legend.
# Smooth x-values for plotting the curves
x_plot = np.linspace(
    X.min().item() - 1,
    X.max().item() + 1,
    200).reshape(-1, 1)

# True line
y_true_line = true_m * x_plot + true_c

# Learned line using recovered coefficients
w_np = w_pred.detach().numpy().flatten()
y_pred_line = w_np[0] * x_plot + w_np[1]

fig, ax = plt.subplots(figsize=(8, 5))

# Plot training samples
ax.scatter(
    X.numpy(),
    Y.numpy(),
    color="#2d6a4f",
    edgecolors="#1b4332",
    s=60,
    zorder=3,
    label="Training samples")

# Plot the true line with known coefficients
ax.plot(
    x_plot,
    y_true_line,
    color="black",
    linewidth=2,
    label=f"True: $y = {true_m}x + {true_c}$")

# Plot the learned line with recovered coefficients
ax.plot(
    x_plot,
    y_pred_line,
    color="red",
    linewidth=2,
    linestyle="--",
    label=f"Learned: $y = {w_np[0]:.2f}x + {w_np[1]:.2f}$")

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("Linear Regression via Gradient Descent")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Visualizing the Error Surface

To build intuition for why gradient descent works, we can visualize how the Sum of Squared Errors (SSE) changes as we vary each parameter individually.

  1. m_range / c_range: Create 200 candidate values centered on the true slope and intercept, respectively.

  2. SSE sweep for slope (sse_m): For each candidate m, the intercept is held fixed at the true value c = 4.0 and the SSE \sum(y_i - (m x_i + c))^2 is computed. This produces a U-shaped curve whose minimum sits at the true slope m = 2.5.

  3. SSE sweep for intercept (sse_c): Likewise, the slope is fixed at m = 2.5 and the intercept is varied.

  4. Normalization: Both SSE curves are divided by their maximum so they share a common [0, 1] scale for easy comparison.

  5. Side-by-side plots: The left panel shows SSE vs w_0 (slope) and the right panel shows SSE vs w_1 (intercept). Both curves are convex with a single minimum, which is what makes gradient descent guaranteed to converge for linear regression.

# SSE error profiles: 
# sweep one parameter while holding the other at its true value

# Range of values to sweep
m_range = np.linspace(true_m - 4, true_m + 4, 200)
c_range = np.linspace(true_c - 8, true_c + 8, 200)

X_np = X.numpy()
Y_np = Y.numpy()

# SSE when sweeping m (slope), with c fixed at true_c
sse_m = np.array(
    [np.sum(
        (Y_np - (m_val * X_np + true_c)) ** 2) 
        for m_val in m_range])

# SSE when sweeping c (intercept), with m fixed at true_m
sse_c = np.array(
    [np.sum(
        (Y_np - (true_m * X_np + c_val)) ** 2) 
        for c_val in c_range])

# Normalize to [0, 1] for comparison for two parameters
sse_m_norm = sse_m / sse_m.max()
sse_c_norm = sse_c / sse_c.max()

# Find the minimum positions
m_min_idx = np.argmin(sse_m_norm)
c_min_idx = np.argmin(sse_c_norm)
m_at_min = m_range[m_min_idx]
c_at_min = c_range[c_min_idx]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))

# Plot SSE vs slope (m)
ax1.plot(m_range, sse_m_norm, color="royalblue", linewidth=2)

# Mark the minimum point on the curve
ax1.axvline(
    m_at_min, color="red", linewidth=1.5,
    linestyle="--", label="Global min")
ax1.text(
    m_at_min + 0.3, 0.05,
    f"$m$ = {m_at_min:.2f}", color="red", fontsize=10)

ax1.set_xlabel(r"$w_0$ (slope, $m$)")
ax1.set_ylabel("SSE")
ax1.set_ylim(-0.05, 1.05)
ax1.grid(True, alpha=0.3)
ax1.set_title("SSE vs Slope (m)")

# Plot SSE vs intercept (c)
ax2.plot(c_range, sse_c_norm, color="royalblue", linewidth=2)

# Mark the minimum point on the curve
ax2.axvline(
    c_at_min, color="red", linewidth=1.5,
    linestyle="--", label="Global min")
ax2.text(
    c_at_min + 0.5, 0.05,
    f"$c$ = {c_at_min:.2f}", color="red", fontsize=10)

ax2.set_xlabel(r"$w_1$ ($y$-intercept, $c$)")
ax2.set_ylabel("SSE")
ax2.set_ylim(-0.05, 1.05)
ax2.grid(True, alpha=0.3)
ax2.set_title("SSE vs Intercept (c)")


fig.suptitle(
    "Error (SSE) profiles when sweeping slope and intercept",
    y=1.02)
plt.tight_layout()
plt.show()

NoteWhy do the SSE minima not occur exactly at the true parameters (2.5, 4.0)?

The global minimum of each SSE curve is the parameter value that best fits the noisy samples, not the true underlying line. Because we added Gaussian noise (noise = torch.randn(N, 1) * 1.5), the 20 observed points are randomly shifted away from the true line y = 2.5x + 4. The least-squares solution on those noisy points will be slightly different from (2.5, 4.0).

With only 20 samples and noise scale 1.5, the finite sample “pulls” the best-fit slope and intercept away from the true values. If we increase N to thousands or reduce the noise, the SSE minima would converge toward the true parameters.

This is the standard distinction between:

  • True (population) parameters: m = 2.5, c = 4.0 (i.e., the data-generating values).
  • Estimated (sample) parameters: the values that minimize the error on the finite noisy dataset (i.e., what the SSE plots show). The gradient descent process is trying to find the best fit to the observed data, which is why it converges to the SSE minima rather than the true parameters.

3D Error Surface and Gradient Descent Trajectory

The previous step swept one parameter at a time. Here we vary both the slope and intercept simultaneously to produce a full 3D error surface, then overlay the path that gradient descent actually follows.

  1. MM, CC = torch.meshgrid(m_grid, c_grid, indexing="xy"): Creates an 80×80 grid of (m, c) pairs covering the search space.

  2. SSE computation: For every grid point, the residuals y_i - (m x_i + c) are computed and squared, giving the SSE at that (m, c). The result is normalized to [0, 1].

  3. Global minimum: torch.argmin locates the grid point with the lowest SSE. Its coordinates (m_\text{min}, c_\text{min}) are stored for plotting as a black dot.

  4. ax.plot_surface(..., cmap="RdBu"): Renders the SSE surface as a colored 3D plot, where blue denotes the low-error regions, whereas red represents the high-error regions.

  5. Gradient descent trajectory: Starting from the corner of the grid (m_0, c_0), we manually run 600 steps of gradient descent with learning rate 0.001. At each step, the analytical gradients \frac{\partial\,\text{SSE}}{\partial m} = -2\sum r_i x_i and \frac{\partial\,\text{SSE}}{\partial c} = -2\sum r_i are used to update the parameters. The resulting path is drawn as a grey dashed line sliding down the surface toward the minimum.

  6. Black marker: The global minimum is marked with a scatter point, confirming that the gradient descent path converges to the bowl’s lowest point (i.e., the best-fit parameters).

NoteWhy do we compute the gradients manually here instead of using autograd?

The manual gradient descent implemented here serves a different purpose from Step 3. In Step 3, we used PyTorch’s autograd + NAdam to learn the weights. Here, we need to record the full (m, c, SSE) trajectory at every iteration so we can draw the path on the 3D surface (i.e., a visualization that autograd doesn’t produce automatically).

The analytical gradients come from differentiating the SSE directly:

\text{SSE} = \sum_{i=1}^{N} r_i^2 \quad \text{where } r_i = y_i - (mx_i + c)

Applying the chain rule:

\frac{\partial\,\text{SSE}}{\partial m} = \sum_{i} 2r_i \cdot \frac{\partial r_i}{\partial m} = \sum_{i} 2r_i \cdot (-x_i) = -2\sum r_i x_i

\frac{\partial\,\text{SSE}}{\partial c} = \sum_{i} 2r_i \cdot \frac{\partial r_i}{\partial c} = \sum_{i} 2r_i \cdot (-1) = -2\sum r_i

Each gradient tells us: “which direction and how much should I move this parameter to reduce SSE the fastest?” The update rule m -= lr * grad_m then takes a small step in the opposite direction of steepest ascent (i.e., steepest descent). We use vanilla gradient descent (not NAdam) here with a small learning rate (0.001) so the trajectory has many visible steps that trace a smooth curve down the bowl, making the visualization clearer.

# 3D SSE surface: vary both slope (w0) and intercept (w1)

m_grid = torch.linspace(true_m - 4, true_m + 4, 80)
c_grid = torch.linspace(true_c - 8, true_c + 8, 80)
MM, CC = torch.meshgrid(m_grid, c_grid, indexing="xy")

# Compute SSE at every (m, c) combination
SSE = torch.zeros_like(MM)
for i in range(MM.shape[0]):
    for j in range(MM.shape[1]):
        residuals = Y.squeeze() - (MM[i, j] * X.squeeze() + CC[i, j])
        SSE[i, j] = torch.sum(residuals ** 2)

# Normalize to [0, 1]
SSE_norm = SSE / SSE.max()

# Find global minimum
min_idx = torch.argmin(SSE_norm)
min_i, min_j = divmod(min_idx.item(), SSE_norm.shape[1])
m_min = MM[min_i, min_j].item()
c_min = CC[min_i, min_j].item()
sse_min = SSE_norm[min_i, min_j].item()

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection="3d")

# Surface plot of SSE over the (m, c) grid
# RdBu_r: reversed so blue = low error, red = high error
ax.plot_surface(
    MM.numpy(), CC.numpy(), SSE_norm.numpy(),
    cmap="RdBu_r", alpha=0.8)

# Gradient descent trajectory on the surface
m_path, c_path, sse_path = [m_grid[0].item()], [c_grid[0].item()], []
m_gd = m_grid[0].item()
c_gd = c_grid[0].item()
lr_gd = 0.001

for _ in range(600):
    res = Y_np - (m_gd * X_np + c_gd)
    sse_val = np.sum(res ** 2)
    sse_path.append(sse_val / SSE.max().item())
    grad_m = -2 * np.sum(res * X_np)
    grad_c = -2 * np.sum(res)
    m_gd -= lr_gd * grad_m
    c_gd -= lr_gd * grad_c
    m_path.append(m_gd)
    c_path.append(c_gd)

# Final SSE for last point
res = Y_np - (m_gd * X_np + c_gd)
sse_path.append(np.sum(res ** 2) / SSE.max().item())

ax.plot(
    m_path, c_path, sse_path,
    color="grey",
    linewidth=3,
    linestyle="--",
    label="Gradient descent path")

# Mark the global minimum
ax.scatter(
    [m_min], [c_min], [sse_min],
    color="black", s=100, zorder=5,
    label=f"Global min ($m$={m_min:.2f}, $c$={c_min:.2f})")

ax.set_xlabel(
    r"$w_0$ (slope, $m$)",
    fontweight="bold",
    fontsize=12,
    labelpad=1)

ax.set_ylabel(
    r"$w_1$ (intercept, $c$)",
    fontweight="bold",
    fontsize=12,
    labelpad=1)

ax.set_zlabel(
    "SSE",
    fontweight="bold",
    fontsize=12,
    labelpad=0.1,
    rotation=90)

ax.set_title(
    "3D SSE Error Surface for Parameters Slope and Intercept",
    pad=15)

ax.legend(fontsize=10)
# ax.view_init(elev=30, azim=225)
plt.tight_layout()
plt.show()

Gradient Vectors and Partial Derivatives

Partial Derivatives

When a function depends on more than one variable, we can ask: how fast does the output change if we vary just one of the inputs while keeping the others fixed? The answer is a partial derivative.

The partial derivative of a function f with respect to the variable x is written \frac{\partial f}{\partial x}, and we find it by differentiating f with respect to x while treating every other variable as a constant.

To illustrate, consider two different functions with two input variables x and y:

f_1(x, y) = x + y \qquad\qquad f_2(x, y) = x^2 + y^2

The code below plots each function as a 3D surface (left) alongside its contour map (right).

  • f_1 (top row): A flat plane tilted equally in the x- and y-directions. Its contours are evenly spaced parallel lines, reflecting the fact that both partial derivatives are constants (\frac{\partial f_1}{\partial x} = 1, \frac{\partial f_1}{\partial y} = 1), which means that the slope is the same everywhere.

  • f_2 (bottom row): An elliptic paraboloid centered at the origin. Its contours are concentric circles because \frac{\partial f_2}{\partial x} = 2x and \frac{\partial f_2}{\partial y} = 2y. The slope grows as we move away from the center, so the contours become more tightly packed further out.

In general, \frac{\partial f}{\partial x} gives the rate of change of f in the x-direction, whereas \frac{\partial f}{\partial y} gives the rate of change in the y-direction.

# Grid over [-10, 10] x [-10, 10]
xg = torch.linspace(-10, 10, 200)
yg = torch.linspace(-10, 10, 200)
XG, YG = torch.meshgrid(xg, yg, indexing="xy")

# f1(x, y) = x + y   
# (a plane)
F1 = XG + YG

# f2(x, y) = x^2 + y^2  
# (a paraboloid)
F2 = XG ** 2 + YG ** 2

# Convert to NumPy for plotting
xg_np, yg_np = XG.numpy(), YG.numpy()
f1_np, f2_np = F1.numpy(), F2.numpy()

fig = plt.figure(figsize=(8, 10))

# --- Top-left: 3D surface of f1 ---
ax1 = fig.add_subplot(2, 2, 1, projection="3d")
ax1.plot_surface(xg_np, yg_np, f1_np, cmap="RdBu_r", alpha=0.9)
ax1.set_xlabel("x")
ax1.set_ylabel("y")
ax1.set_title(r"$f_1(x, y) = x + y$")

# --- Top-right: Contours of f1 ---
ax2 = fig.add_subplot(2, 2, 2)
cs1 = ax2.contour(xg_np, yg_np, f1_np, levels=15, cmap="RdBu_r")
ax2.set_xlabel("x")
ax2.set_ylabel("y")
ax2.set_title(r"Contours of $f_1(x, y) = x + y$")
ax2.grid(True, alpha=0.3)
ax2.set_aspect("equal")

# --- Bottom-left: 3D surface of f2 ---
ax3 = fig.add_subplot(2, 2, 3, projection="3d")
ax3.plot_surface(xg_np, yg_np, f2_np, cmap="RdBu_r", alpha=0.9)
ax3.set_xlabel("x")
ax3.set_ylabel("y")
ax3.set_title(r"$f_2(x, y) = x^2 + y^2$")

# --- Bottom-right: Contours of f2 ---
ax4 = fig.add_subplot(2, 2, 4)
cs2 = ax4.contour(xg_np, yg_np, f2_np, levels=15, cmap="RdBu_r")
ax4.set_xlabel("x")
ax4.set_ylabel("y")
ax4.set_title(r"Contours of $f_2(x, y) = x^2 + y^2$")
ax4.grid(True, alpha=0.3)
ax4.set_aspect("equal")

fig.suptitle(
    r"The functions $f_1$ and $f_2$ and their corresponding contours",
    fontsize=13, y=0.02)

plt.tight_layout()
plt.show()

The Gradient Vector

When we collect all the partial derivatives of a function into a single vector, we get the gradient, written \nabla f:

\nabla f(x, y) = \begin{bmatrix} \frac{\partial f}{\partial x} \\[6pt] \frac{\partial f}{\partial y} \end{bmatrix}

The gradient points in the direction of steepest increase of f, and its magnitude tells us how steep that increase is; the negative gradient -\nabla f therefore points toward steepest decrease, which is exactly the direction gradient descent follows.

Example: Gradient of f_1(x, y) = x + y

The partial derivatives are \frac{\partial f_1}{\partial x} = 1 and \frac{\partial f_1}{\partial y} = 1, so the gradient is:

\nabla f_1 \begin{aligned} = \begin{bmatrix} \frac{\partial f_1}{\partial x} \\ \frac{\partial f_1}{\partial y} \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \end{bmatrix} \end{aligned}

This is a constant vector, which is the same everywhere in the (x, y) plane. It tells us the function increases at the same rate regardless of where we stand, which is consistent with f_1 being a flat tilted plane with the parallel and evenly spaced contours.

Example: Gradient of f_2(x, y) = x^2 + y^2

The partial derivatives are \frac{\partial f_2}{\partial x} = 2x and \frac{\partial f_2}{\partial y} = 2y, so the gradient is:

\nabla f_2 \begin{aligned} = \begin{bmatrix} \frac{\partial f_2}{\partial x} \\ \frac{\partial f_2}{\partial y} \end{bmatrix} = \begin{bmatrix} 2x \\ 2y \end{bmatrix} \end{aligned}

Unlike f_1, this gradient depends on the position (x, y). At the origin (0, 0) the gradient is the zero vector (i.e. at the bottom of the bowl where there is no slope). As we move away from the origin the gradient grows in magnitude, pointing radially outward, which matches the tightening concentric contours we saw in the plot above. For instance, at (1, 1) the gradient is (2, 2), while at (2, 1) it is (4, 2), which is steeper in the x-direction because we are further from the centre along that axis.

Connecting it back to our SSE

In the linear-regression example above, the loss function \text{SSE}(m, c) depends on two parameters. Its gradient is:

\nabla \text{SSE} = \begin{bmatrix} \frac{\partial\,\text{SSE}}{\partial m} \\[6pt] \frac{\partial\,\text{SSE}}{\partial c} \end{bmatrix} = \begin{bmatrix} -2\sum r_i\, x_i \\[6pt] -2\sum r_i \end{bmatrix}

At each training step, gradient descent evaluates this vector and moves (m, c) in the opposite direction, sliding down the 3D error bowl until it reaches the minimum.

Computing Partial Derivatives with Symbolic Differentiation

Of course, we can also verify our partial derivatives symbolically using SymPy, a Python library for symbolic mathematics.

The code below takes f_2(x, y) = x^2 + y^2 and uses diff(f2, x) and diff(f2, y) to compute the partial derivatives automatically:

\frac{\partial f_2}{\partial x} = 2x \qquad\qquad \frac{\partial f_2}{\partial y} = 2y

It then evaluates the gradient vector \nabla f_2 = (2x,\; 2y) at two points using .subs():

  • At (1, 1): \nabla f_2 = (2, 2), a moderate slope in both directions.
  • At (2, 1): \nabla f_2 = (4, 2), steeper in the x-direction because we are further from the centre of the bowl.

This confirms that the gradient grows in magnitude as we move away from the minimum at the origin, consistent with the tightening contour spacing we saw in the plot above.

from sympy.abc import x, y
from sympy import diff, pprint

f2 = x**2 + y**2
df2dx = diff(f2, x)
df2dy = diff(f2, y)

print("Partial derivative of")
pprint(f2)
print("-" * 79)

print("with respect to x is")
pprint(df2dx)
print("-" * 79)

print("and with respect to y is")
pprint(df2dy)
print("-" * 79)

print("gradient at (1,1) is ({},{})".format(df2dx.subs([(x,1),(y,1)]),
df2dy.subs([(x,1),(y,1)])))
print("-" * 79)

print("gradient at (2,1) is ({},{})".format(df2dx.subs([(x,2),(y,1)]),
df2dy.subs([(x,2),(y,1)])))
Partial derivative of
 2    2
x  + y 
-------------------------------------------------------------------------------
with respect to x is
2⋅x
-------------------------------------------------------------------------------
and with respect to y is
2⋅y
-------------------------------------------------------------------------------
gradient at (1,1) is (2,2)
-------------------------------------------------------------------------------
gradient at (2,1) is (4,2)

References

  1. Brownlee, Jason, Stefania Cristina, and Mehreen Saeed. Calculus for machine learning. Machine Learning Mastery, 2022.