22 - Calculus 2: Gradient, Integrals, Gradient Descent

Machine Learning Engineer course · Module 02 · Mathematics and statistics · Lesson 7

Goal. Extend the derivative to functions of many variables (partial derivatives and the gradient), prove that the gradient indicates the direction of steepest ascent, turn this idea into the algorithm that trains almost every model (gradient descent), see backpropagation as the chain rule applied to a neuron, and introduce integrals as areas. All with code that runs and that can be verified.

Myth to dispel. "A bigger gradient descent step always speeds up training." There is a threshold beyond which the steps diverge: we will derive it and measure it.

Introduction

In the previous lesson the function depended on a single variable. A neural network has millions of parameters: its "loss function" (the error to minimise) is a function of millions of variables. Three tools are needed: how to differentiate with respect to one variable at a time, how to combine these derivatives into a vector that indicates where to go, and how to use that vector to descend step by step.

Terms you will use (defined here at first use): loss function: a number that measures how much the model is wrong, to be made small; parameters: the numbers that the model adjusts during training; minimum: a point where the function does not go down any further.

The path starts from partial derivatives and the gradient, proves that the gradient points to the steepest ascent and derives gradient descent from it, measuring the maximum step and the effect of narrow valleys. It then applies the algorithm to linear regression and to backpropagation on a neuron, deals with convexity and local minima, and closes with integrals as areas. From the Calculus 1 lesson you need the derivative and the chain rule; the complete algorithm is taken up again in Gradient descent.

Functions of several variables and partial derivatives

The analogy: the landscape

A function f(x, y) assigns to every point of a plane an altitude: it is a landscape. If you are at a point and move east, the elevation rises or falls with a certain slope; if you move north you find another. The partial derivative with respect to x, ∂f/∂x, is the slope you find moving only along x, holding y still. It is computed like the derivative of one variable, treating the others as constants.

Example: f(x, y) = x² + 3y² + xy. Then ∂f/∂x = 2x + y (the term 3y² is constant with respect to x; xy has derivative y) and ∂f/∂y = 6y + x.

The gradient is the vector of the partial derivatives:

∇f(x, y) = (∂f/∂x, ∂f/∂y)

(∇ is read "nabla"). For a function of d variables it is a vector of d components.

import numpy as np

f = lambda p: p[0] ** 2 + 3 * p[1] ** 2 + p[0] * p[1]
grad_f = lambda p: np.array([2 * p[0] + p[1], 6 * p[1] + p[0]])        # partial derivatives computed by hand

def grad_numerico(funzione, p, h=1e-5):
    """Gradient by central differences: one partial derivative per coordinate."""
    p = np.asarray(p, dtype=float)
    g = np.zeros_like(p)
    for i in range(p.size):
        passo = np.zeros_like(p)
        passo[i] = h                                  # we move only along coordinate i
        g[i] = (funzione(p + passo) - funzione(p - passo)) / (2 * h)
    return g

p0 = np.array([1.5, -0.7])
print("analytic gradient :", grad_f(p0))
print("numerical gradient:", np.round(grad_numerico(f, p0), 6))

Output:

analytic gradient : [ 2.3 -2.7]
numerical gradient: [ 2.3 -2.7]

What the code does. grad_f implements the two formulas derived above. grad_numerico computes each component by central difference (previous lesson), moving by ±h only along the coordinate i: it is the definition of partial derivative, translated into code. The two vectors coincide: this comparison is the gradient check, which we will meet again.

The gradient indicates the steepest ascent

Why precisely the gradient? Let us move from p in the direction of a vector of length 1 (unit vector) u for a small stretch t. By the first-order Taylor expansion (the linearisation of the previous lesson, now in several variables):

f(p + t·u) ≈ f(p) + t · (∇f(p) · u)

where · is the dot product (the sum of the products of the components). The quantity D_u f = ∇f · u is the directional derivative: by how much f grows per unit of path in the direction u.

Which u makes it maximal? The dot product of two vectors is |∇f| · |u| · cos θ, with θ the angle between them. Since |u| = 1:

D_u f = |∇f| · cos θ ≤ |∇f|,

with equality only for cos θ = 1, that is θ = 0: u parallel to the gradient. Therefore:

• the gradient points in the direction of maximum growth, and its length is the maximum slope;
• −∇f points in the direction of steepest descent;
• in the directions perpendicular to the gradient (cos θ = 0) the function, to first order, does not change: these are the level curves (the lines of constant elevation on a map).

The three consequences are the same inequality read in three directions: maximum ascent along the gradient, maximum descent along its opposite, no first-order change along the level curves. For optimisation the second one matters, because it says how to choose the next step knowing only the local slope; the listing below checks it numerically by trying 3601 directions.

angoli = np.linspace(0, 2 * np.pi, 3601)
versori = np.stack([np.cos(angoli), np.sin(angoli)], axis=1)         # 3601 directions of length 1
g = grad_f(p0)
derivate_direzionali = versori @ g                                    # D_u f = ∇f · u for each u

migliore = versori[np.argmax(derivate_direzionali)]
print("direction of maximum growth:", np.round(migliore, 4))
print("normalised gradient        :", np.round(g / np.linalg.norm(g), 4))
print("maximum slope:", round(float(derivate_direzionali.max()), 4), "| |∇f| =", round(float(np.linalg.norm(g)), 4))

Output:

direction of maximum growth: [ 0.6481 -0.7615]
normalised gradient        : [ 0.6485 -0.7612]
maximum slope: 3.5468 | |∇f| = 3.5468

What the code does. It generates 3601 directions of length 1 (the rows of versori), computes for each the directional derivative with a single matrix-vector product (@), and looks for the best direction. It coincides with the normalised gradient and the maximum slope is |∇f|, as proved.

Gradient descent

The analogy: the hiker in the fog

You are on a mountain, in thick fog, and want to get down to the valley. You cannot see the landscape: you only feel the slope under your feet. Strategy: at each step look in which direction the ground falls most steeply and take a step in that direction. Repeat. This is gradient descent (Cauchy proposed it in 1847):

xₖ₊₁ = xₖ − η · ∇f(xₖ)

where η ("eta", the learning rate) is the length of the step in proportion to the slope. The minus sign because we descend.

How large must the step be? A case solved to the end

Take the simplest function with a minimum, f(x) = (x − 3)², with minimum at x = 3. The derivative is f'(x) = 2(x − 3). One step gives

xₖ₊₁ − 3 = xₖ − 3 − η·2(xₖ − 3) = (1 − 2η)(xₖ − 3).

The distance from the minimum is multiplied by (1 − 2η) at each step. Hence:

• if |1 − 2η| < 1, that is 0 < η < 1, the distance tends to zero: it converges;
• if η = 0.5 the factor is 0: it reaches the minimum in a single step;
• if 0.5 < η < 1 the factor is negative: it oscillates around the minimum, but gets closer;
• if η > 1 the factor has modulus > 1: the distance grows, the steps diverge.

The four regimes depend only on the modulus and the sign of the factor (1 − 2η): monotone if positive, oscillating if negative, divergent if the modulus exceeds 1. The general result is the same for any parabola, with 2 replaced by the curvature, and it is the reason why the maximum step depends on the curvature of the function and cannot be chosen once and for all.

def discesa_1d(eta, x0=-5.0, passi=8):
    x, traiettoria = x0, [x0]
    for _ in range(passi):
        x = x - eta * 2 * (x - 3)              # x_{k+1} = x_k - η f'(x_k)
        traiettoria.append(x)
    return traiettoria

for eta in (0.1, 0.5, 0.9, 1.05):
    t = discesa_1d(eta)
    print(f"η = {eta:<4}: factor (1-2η) = {1 - 2 * eta:+.2f} | trajectory: {[round(v, 2) for v in t[:7]]}")

Output:

η = 0.1 : factor (1-2η) = +0.80 | trajectory: [-5.0, -3.4, -2.12, -1.1, -0.28, 0.38, 0.9]
η = 0.5 : factor (1-2η) = +0.00 | trajectory: [-5.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0]
η = 0.9 : factor (1-2η) = -0.80 | trajectory: [-5.0, 9.4, -2.12, 7.1, -0.28, 5.62, 0.9]
η = 1.05: factor (1-2η) = -1.10 | trajectory: [-5.0, 11.8, -6.68, 13.65, -8.71, 15.88, -11.17]

What the code does. The function discesa_1d applies the update for 8 steps starting from x = −5. With η = 0.1 the factor is 0.8: slow, monotone approach. With η = 0.5 the factor is 0: after the first step you are already at 3. With η = 0.9 the factor is −0.8: it oscillates above and below 3 while approaching. With η = 1.05 the factor is −1.1: the oscillation widens, the values diverge. The number printed as "factor" predicts exactly what is seen.

Narrow valleys and conditioning

In the real world the variables have very different curvatures. Take f(x, y) = x² + 10y²: a narrow and long "valley", steep along y (curvature 20) and gentle along x (curvature 2). The gradient is (2x, 20y), and a step gives

x ← (1 − 2η)·x, y ← (1 − 20η)·y.

To converge both |1 − 2η| < 1 and |1 − 20η| < 1 are needed: the second is the more restrictive and imposes η < 0.1. In general, the maximum step is 2 / (largest curvature). But with η close to that limit the gentle direction converges very slowly: the factor (1 − 2η) stays close to 1. The ratio between maximum and minimum curvature (here 10) is the condition number: the larger it is, the slower or more zig-zagging the descent.

def discesa_2d(eta, passi=200, p0=(-2.6, 0.8), tol=1e-6):
    p = np.array(p0)
    for k in range(passi):
        p = p - eta * np.array([2 * p[0], 20 * p[1]])
        if np.linalg.norm(p) < tol:
            return f"converges in {k + 1} steps"
        if np.linalg.norm(p) > 1e6:
            return f"DIVERGES after {k + 1} steps"
    return f"after {passi} steps |p| = {np.linalg.norm(p):.2e}"

for eta in (0.01, 0.05, 0.09, 0.099, 0.11):
    print(f"η = {eta:<5}: {discesa_2d(eta)}")

Output:

η = 0.01 : after 200 steps |p| = 4.57e-02
η = 0.05 : converges in 141 steps
η = 0.09 : converges in 75 steps
η = 0.099: after 200 steps |p| = 1.41e-02
η = 0.11 : DIVERGES after 77 steps

What the code does. It applies the descent to the function x² + 10y² with various steps and reports whether it converges (norm below 10⁻⁶) or diverges (norm above 10⁶). The optimal step for this problem is a compromise: η = 0.01 is safe but after 200 steps has not yet reached the tolerance (slow gentle direction); with 0.05 and 0.09 it converges in 141 and 75 steps, with zig-zag along y; close to the limit (0.099) the factor 1 − 20η ≈ −0.98 makes the oscillation almost undamped and it does not converge in 200 steps; 0.11 exceeds the limit 0.1 and diverges. It is the same mathematics that governs all neural networks, only in millions of dimensions.

Practical consequence: scaling the variables. If two features have very different scales (one in units, the other in thousands), the loss function has a narrow valley like this one, and the descent is slow or unstable. Standardising the features (NumPy lesson) makes the curvatures similar. It is the deeper reason why in the scikit-learn lesson we scaled before using linear models trained with these methods.

Linear regression with gradient descent

We take up the least-squares line of the statistics lesson, but solve it without the closed formula. The loss is the mean squared error

L(a, b) = (1/n) Σ (yᵢ − a − b·xᵢ)².

With the chain rule (the derivative of a square is 2·(...) times the inner derivative), denoting the residual rᵢ = yᵢ − a − b·xᵢ:

∂L/∂a = −(2/n) Σ rᵢ, ∂L/∂b = −(2/n) Σ xᵢ·rᵢ.

rng = np.random.default_rng(0)
n = 200
x_grezza = rng.uniform(0, 100, n)                           # feature with a large scale
y = 5.0 + 0.3 * x_grezza + rng.normal(0, 2.0, n)            # true line: a = 5, b = 0.3

def gd_regressione(x, y, eta, passi):
    a = b = 0.0
    for _ in range(passi):
        r = y - (a + b * x)
        a -= eta * (-2 * r.mean())
        b -= eta * (-2 * (x * r).mean())
    return a, b

# closed formula from the statistics lesson
b_chiusa = ((x_grezza - x_grezza.mean()) * (y - y.mean())).sum() / ((x_grezza - x_grezza.mean()) ** 2).sum()
a_chiusa = y.mean() - b_chiusa * x_grezza.mean()
print(f"closed formula          : a = {a_chiusa:.4f}, b = {b_chiusa:.4f}")

x_std = (x_grezza - x_grezza.mean()) / x_grezza.std()
a_s, b_s = gd_regressione(x_std, y, eta=0.1, passi=200)
# convert back to the original scale: y = a_s + b_s*(x - m)/s  =>  b = b_s/s, a = a_s - b_s*m/s
print(f"GD on standardised x    : a = {a_s - b_s * x_grezza.mean() / x_grezza.std():.4f}, b = {b_s / x_grezza.std():.4f}  (200 steps, η = 0.1)")

for eta in (1e-3, 1e-4):
    a_g, b_g = gd_regressione(x_grezza, y, eta=eta, passi=200)
    esito = f"a = {a_g:.4f}, b = {b_g:.4f}" if abs(a_g) < 1e6 else f"DIVERGES (a = {a_g:.1e})"
    print(f"GD on raw x, η = {eta:g}: {esito}")

Output:

closed formula          : a = 5.0942, b = 0.2949
GD on standardised x    : a = 5.0942, b = 0.2949  (200 steps, η = 0.1)
GD on raw x, η = 0.001: DIVERGES (a = -1.5e+162)
GD on raw x, η = 0.0001: a = 0.0533, b = 0.3662

What the code does. gd_regressione computes at each step the residuals and updates a and b with the two partial derivatives (the means correspond to the 1/n Σ). It compares three results with the closed formula: the descent on the standardised x converges in 200 steps with η = 0.1 and, converted back to the original scale, coincides with the formula. On the raw variable (values up to 100) the step must be small: with η = 10⁻³ the descent diverges (the code flags it), with 10⁻⁴ it approaches but after 200 steps it is still far from the solution, in particular on the intercept. The reason is the conditioning from before.

The "in pieces" version: stochastic descent

With millions of examples computing the gradient on all the data at each step is too costly. Instead a mini-batch is used: a small random subset (say 32 examples) to estimate the gradient. It is a noisy but much cheaper estimate, and on average it points in the right direction (the mean of the gradients on random subsets is the full gradient). It is stochastic gradient descent (SGD), linked to the idea of stochastic approximation of Robbins and Monro (1951).

def sgd_regressione(x, y, eta, epoche, batch, seme=1):
    rng = np.random.default_rng(seme)
    a = b = 0.0
    for _ in range(epoche):
        ordine = rng.permutation(len(x))                  # one "epoch" = one pass over all the data, in random order
        for inizio in range(0, len(x), batch):
            idx = ordine[inizio:inizio + batch]
            r = y[idx] - (a + b * x[idx])
            a -= eta * (-2 * r.mean())
            b -= eta * (-2 * (x[idx] * r).mean())
    return a, b

a_sgd, b_sgd = sgd_regressione(x_std, y, eta=0.05, epoche=30, batch=16)
print("SGD (batch 16, 30 epochs):", round(a_sgd - b_sgd * x_grezza.mean() / x_grezza.std(), 3), round(b_sgd / x_grezza.std(), 4))

Output:

SGD (batch 16, 30 epochs): 4.931 0.2985

What the code does. At each epoch (a full pass over the data) it shuffles the order of the examples and takes them in groups of 16; for each group it computes the gradient on the group alone and takes a step. The result gets very close to the closed formula, but with a little noise: it is the price of the economy. The exact values depend on the seed and the batch size.

Backpropagation: the chain rule on a neuron

A classification neuron computes z = w·x + b, then p = σ(z) (the sigmoid, previous lesson), which is interpreted as the probability of class 1. The error for an example with label y ∈ {0, 1} is the log-loss (cross-entropy):

L = −[ y·ln p + (1 − y)·ln(1 − p) ].

To train, ∂L/∂w and ∂L/∂b are needed. The chain is w → z → p → L, and by the chain rule the derivatives multiply along the path:

∂L/∂w = (∂L/∂p) · (∂p/∂z) · (∂z/∂w).

The three factors: ∂L/∂p = −y/p + (1 − y)/(1 − p) = (p − y)/(p(1 − p)); ∂p/∂z = p(1 − p) (derivative of the sigmoid); ∂z/∂w = x. The product simplifies the factor p(1 − p):

∂L/∂z = p − y, ∂L/∂w = (p − y)·x, ∂L/∂b = p − y.

An elegant result: the gradient with respect to the input of the neuron is simply "predicted minus true". Backpropagation (Rumelhart, Hinton and Williams 1986) is this: the values are computed forwards (forward pass), then the derivatives are multiplied backwards along the graph, reusing the values already computed.

rng = np.random.default_rng(3)
# synthetic data: two groups of points in the plane
n = 300
X0 = rng.normal([-1.0, -1.0], 1.0, size=(n // 2, 2))
X1 = rng.normal([1.2, 1.0], 1.0, size=(n // 2, 2))
X = np.vstack([X0, X1])
y = np.concatenate([np.zeros(n // 2), np.ones(n // 2)])

sigmoide = lambda z: 1 / (1 + np.exp(-z))

def perdita(w, b):
    p = np.clip(sigmoide(X @ w + b), 1e-12, 1 - 1e-12)
    return float(-(y * np.log(p) + (1 - y) * np.log(1 - p)).mean())

def gradiente(w, b):
    p = sigmoide(X @ w + b)
    err = p - y                                   # ∂L/∂z for each example
    return X.T @ err / len(y), err.mean()         # means of (p - y)·x  and of (p - y)

# gradient check: comparison with finite differences
w0, b0 = np.array([0.3, -0.2]), 0.1
gw, gb = gradiente(w0, b0)
h = 1e-6
num_w = np.array([(perdita(w0 + h * e, b0) - perdita(w0 - h * e, b0)) / (2 * h) for e in np.eye(2)])
num_b = (perdita(w0, b0 + h) - perdita(w0, b0 - h)) / (2 * h)
print("backpropagated gradient :", np.round(gw, 6), round(float(gb), 6))
print("numerical gradient      :", np.round(num_w, 6), round(float(num_b), 6))

# training with gradient descent
w, b = np.zeros(2), 0.0
for passo in range(401):
    gw, gb = gradiente(w, b)
    w, b = w - 0.5 * gw, b - 0.5 * gb
    if passo in (0, 10, 100, 400):
        accuratezza = ((sigmoide(X @ w + b) > 0.5) == y).mean()
        print(f"step {passo:>3}: loss {perdita(w, b):.4f} | accuracy {accuratezza:.3f}")
print("learned weights:", np.round(w, 3), "| bias:", round(float(b), 3))

Output:

backpropagated gradient : [-0.428585 -0.549415] 0.025666
numerical gradient      : [-0.428585 -0.549415] 0.025666
step   0: loss 0.4656 | accuracy 0.917
step  10: loss 0.2291 | accuracy 0.923
step 100: loss 0.1845 | accuracy 0.937
step 400: loss 0.1834 | accuracy 0.937
learned weights: [2.099 1.77 ] | bias: -0.457

What the code does. perdita computes the mean log-loss on the whole set (with np.clip to avoid ln(0): floats lesson). gradiente implements exactly the formulas derived above, in vector form: X.T @ err / n is the mean of the vectors (p − y)·x. The gradient check compares the backpropagated gradient with the numerical one: the two rows coincide up to the sixth digit, so the derivative is right. Then the loop applies the descent (η = 0.5) and prints loss and accuracy at four moments: the loss falls and the accuracy rises towards the optimal value for two overlapping groups. The model is, in practice, a logistic regression: the same algorithm as scikit-learn's LogisticRegression, written by hand.

Convexity and local minima

A function is convex if the segment joining any two points of its graph lies above the graph (a cup). For a convex function every local minimum is a global minimum and gradient descent, with an adequate step, finds it. The squared error of linear regression and the log-loss of logistic regression are convex (you would study it with the matrix of second derivatives; here we take it as a fact: see Boyd and Vandenberghe).

Deep neural networks are not convex: there are several valleys. Depending on the starting point you end up in different minima. Let us show it with a one-variable function with two minima, f(x) = x⁴ − 3x² + x.

f4 = lambda x: x ** 4 - 3 * x ** 2 + x
df4 = lambda x: 4 * x ** 3 - 6 * x + 1

for partenza in (-2.0, -0.2, 0.2, 2.0):
    x = partenza
    for _ in range(2000):
        x -= 0.01 * df4(x)
    print(f"start {partenza:>5}: ends at x = {x:.4f}, f(x) = {f4(x):.4f}")

Output:

start  -2.0: ends at x = -1.3008, f(x) = -3.5139
start  -0.2: ends at x = -1.3008, f(x) = -3.5139
start   0.2: ends at x = 1.1309, f(x) = -1.0702
start   2.0: ends at x = 1.1309, f(x) = -1.0702

What the code does. It applies the descent from four starting points. The starts to the left of about 0.17 end in the deep minimum (on the left, f ≈ −3.5); those to the right end in the higher minimum (on the right, f ≈ −1.1). The minimum found depends on where you start. This is why the initialisation of the parameters and the use of several starts or of noise (like that of SGD) are important.

Integrals as areas

The analogy: distance from speed

If you know the speed of the car at every instant, how do you reconstruct how much road it has covered? You divide the time into tiny intervals Δt, in each the speed is almost constant: the distance covered is about speed × Δt. You add up all the pieces. This is the idea of the integral: the area under a curve f between a and b is approximated with n thin rectangles of base Δx = (b − a)/n and height f(xᵢ):

∫ₐᵇ f(x) dx = lim_{n→∞} Σ f(xᵢ)·Δx (Riemann sum).

The fundamental theorem

Derivative and integral are inverse operations. If F' = f, then

∫ₐᵇ f(x) dx = F(b) − F(a).

The intuition is in the sum: divide [a, b] into small intervals with endpoints x₀ = a, x₁, …, xₙ = b. For each interval, F(xᵢ₊₁) − F(xᵢ) ≈ F'(xᵢ)·Δx = f(xᵢ)·Δx (definition of derivative). Summing over all the intervals, the intermediate terms cancel (F(x₁) appears with a minus sign and then plus...): what remains is F(b) − F(a), the telescoping sum. So the area is the "total variation" of F.

from scipy import integrate
import math

f = lambda x: x ** 2
a, b = 0.0, 3.0                                  # ∫ x² dx from 0 to 3 = 3³/3 = 9
print("exact value F(b) - F(a) = 9")
for n in (4, 16, 64, 256, 1024):
    xs = np.linspace(a, b, n, endpoint=False)    # left endpoints of the subintervals
    riemann = (f(xs) * (b - a) / n).sum()
    trapezi = np.trapezoid(f(np.linspace(a, b, n + 1)), dx=(b - a) / n)
    print(f"n = {n:>4}: Riemann sum {riemann:.5f} (error {abs(riemann - 9):.2e}) | trapezoids {trapezi:.5f} (error {abs(trapezi - 9):.2e})")

valore, errore = integrate.quad(f, a, b)
print("scipy.integrate.quad:", valore, "(error estimate:", f"{errore:.1e})")

Output:

exact value F(b) - F(a) = 9
n =    4: Riemann sum 5.90625 (error 3.09e+00) | trapezoids 9.28125 (error 2.81e-01)
n =   16: Riemann sum 8.17383 (error 8.26e-01) | trapezoids 9.01758 (error 1.76e-02)
n =   64: Riemann sum 8.79016 (error 2.10e-01) | trapezoids 9.00110 (error 1.10e-03)
n =  256: Riemann sum 8.94733 (error 5.27e-02) | trapezoids 9.00007 (error 6.87e-05)
n = 1024: Riemann sum 8.98682 (error 1.32e-02) | trapezoids 9.00000 (error 4.29e-06)
scipy.integrate.quad: 9.000000000000002 (error estimate: 1.0e-13)

What the code does. It approximates the area under x² with rectangles (Riemann sum, using the left endpoint) and with trapezoids (which "cut" the curve better), for increasing n. The error of the Riemann sum falls as 1/n, that of the trapezoids as 1/n²: the trapezoids are much more efficient. scipy.integrate.quad uses an adaptive algorithm and also returns an error estimate.

Integrals and probability

Probability, for a continuous variable, is an area (random variables lesson): P(a ≤ X ≤ b) = ∫ₐᵇ f(x) dx, and the expected value is E[X] = ∫ x·f(x) dx.

densita_normale = lambda x: math.exp(-x * x / 2) / math.sqrt(2 * math.pi)
area, _ = integrate.quad(densita_normale, -1.96, 1.96)
print("∫ standard normal density between -1.96 and 1.96 =", round(area, 6))
print("total area over (-∞, +∞)                         =", round(integrate.quad(densita_normale, -np.inf, np.inf)[0], 6))

lam = 0.5
valore_atteso, _ = integrate.quad(lambda x: x * lam * math.exp(-lam * x), 0, np.inf)
print("E[X] for the exponential with λ = 0.5 (theory 1/λ = 2):", round(valore_atteso, 6))

Output:

∫ standard normal density between -1.96 and 1.96 = 0.950004
total area over (-∞, +∞)                         = 1.0
E[X] for the exponential with λ = 0.5 (theory 1/λ = 2): 2.0

What the code does. It computes as an integral what in lesson 4 we obtained with the function erf: the area of the normal density between ±1.96 is 0.95, and the total area 1. Then it verifies that the expected value of the exponential, defined as the integral ∫ x·λe^(−λx) dx, is 1/λ = 2. np.inf indicates an infinite endpoint.

Typical mistakes (symptom, cause, remedy)

• Loss that rises or nan. Cause: η too large, or scale of the variables. Remedy: reduce η (by a factor of 10), standardise the features.
• Loss that falls very slowly. Cause: η too small or narrow valley. Remedy: increase η almost to the limit, adaptive optimisers, scaling.
• Wrong gradient without noticing. Cause: a badly written derivative. Remedy: gradient check with finite differences on a small example.
• Forgetting the minus sign (you climb instead of descending). Remedy: check that the loss falls in the first steps.
• Confusing a non-convex function with a convex one. Remedy: several starts; look at the validation loss.
• Using ln(0) in the log-loss. Remedy: np.clip or stable functions.
• Integrating numerically with too few points or with singularities. Remedy: increase n, use quad and check the error estimate.
• Forgetting to divide by n in the mean of the gradients: it changes the scale of the optimal step.

The first two items are the two faces of the same parameter, η, and are told apart by looking at the loss curve: it rises or gives nan in the first case, it barely falls in the second. The third, the badly written gradient, is the most insidious because the algorithm runs without errors and merely produces worse results: this is why the gradient check should be done before every training.

How to work in practice

Start from the gradient check every time you write a gradient by hand, record the loss at every epoch and plot it (the shape says whether η is suitable), standardise the variables before training and fix the seed to make shuffling and initialisation reproducible. In production use the libraries (scikit-learn, PyTorch): writing the gradient by hand serves to understand, not to replace them.

There are few new terms. The partial derivative is the slope with respect to one variable with the others held still, the gradient ∇f is the vector of the partial derivatives and the directional derivative measures the variation per unit of path in a given direction. The learning rate η is the length of the step; an epoch is a full pass over the data and a mini-batch the subset that estimates the gradient. A convex function has a single minimum; backpropagation computes the gradients with the chain rule backwards along the graph; the Riemann sum approximates the area under a curve with rectangles.

Mini-project

1. Implement gradient descent for logistic regression on scikit-learn's breast-cancer dataset (standardised), compare accuracy and coefficients with LogisticRegression(penalty=None).
2. For the function f(x, y) = x² + 10y² verify for several values of η that the convergence limit is 0.1 and that with η = 0.1 exactly the behaviour is oscillatory.
3. Plot (with Matplotlib) the loss as a function of the step for η = 0.01, 0.1, 0.5, 1.0 in the synthetic logistic regression.

The first point uses a dataset included in scikit-learn, so it requires no download; the second repeats for two variables the check that the lesson did with one; the third makes visible the loss curves of the four regimes (slow, good, oscillating, divergent) that so far you have read as numbers.

Challenge

• Easy. Compute by hand the gradient of f(x, y) = x²y + sin(y) and verify it with grad_numerico.
• Medium. Prove that, for f(x) = (a/2)(x − m)², gradient descent converges if and only if 0 < η < 2/a, and find the η that converges in one step.
• Hard. Add momentum to the descent: v ← β·v + ∇f, x ← x − η·v, with β = 0.9. Show on x² + 10y² that it reduces the zig-zag and compare the number of steps.
• Expert. Extend backpropagation to a network with one hidden layer of 3 neurons (sigmoid) and one output; derive the gradients by hand with the chain rule, implement them in NumPy and verify them with the gradient check.
• Enterprise. Choose a real regression problem from your organisation. Train it with gradient descent written by hand, with and without standardisation, and document the learning rate chosen, the number of epochs and the final loss.

The Medium challenge is the algebraic generalisation of the case solved in the text, the Hard one shows the simplest remedy to the zig-zag of narrow valleys, the Expert one is the logical step from logistic regression to deep networks. If you pick only one, pick the Expert: it directly prepares the lesson on complete backpropagation.

What you can do now

You can compute a gradient and verify it with finite differences, prove why it points to the steepest ascent, predict whether gradient descent converges from the factor 1 − 2η (or from the maximum curvature in narrow valleys) and train a linear regression and a logistic neuron by hand. You can also read backpropagation as the chain rule and compute an area with Riemann sums and with scipy.integrate.quad. The next step is linear algebra, the language in which X @ w replaces the loops.

The same series also includes: 19 - Random Variables and Distributions, 25 - Matrix Calculus and Automatic Differentiation: From the Jacobian to an Autograd Written From Scratch, 26 - Constrained Optimisation: Lagrange Multipliers, Duality and the KKT Conditions (which refer back to this lesson).