Calculus and Optimization for Machine Learning

Phase 0 · Week 3–4 · ~20–25 hours

Calculus is the engine of learning. Not metaphorically — literally. Every weight update in every neural network, every convergence proof, every adaptive optimizer is a calculation of how much the loss function changes when you nudge a parameter. If linear algebra gives you the objects ML operates on (vectors and matrices), calculus gives you the dynamics — the mechanism by which a model gets better over time. This file treats calculus not as abstract mathematics but as the derivation of the update rule that runs every training loop you’ve ever written.

---

Why Calculus Is Not Optional (And Why You Should Derive, Not Memorize)

Many practitioners memorize: “gradient descent update is w = w - lr * gradient.” That’s not understanding — that’s transcription. The question is: where does that formula come from? The answer requires calculus, and understanding the answer changes how you debug training failures. When your loss isn’t decreasing, it’s not a magic hyperparameter problem — it’s a geometry problem about the loss surface, and calculus is the tool for reasoning about it.

---

Core Concepts

1. Derivatives as Local Linear Approximation

A derivative f'(x) at a point x₀ tells you the slope of the best linear approximation to f at that point. The fundamental equation:

f(x₀ + h) ≈ f(x₀) + f'(x₀) · h (for small h)

This is not just a definition — it’s the reason gradient descent works. If you move a small step h in the direction that decreases f, you get a proportional decrease. Gradient descent is the repeated application of this local linear approximation.

import numpy as np
import matplotlib.pyplot as plt

def numerical_derivative(f, x, h=1e-5):
    """Numerical derivative via central difference (more accurate than forward diff)."""
    return (f(x + h) - f(x - h)) / (2 * h)

# Test on a known function
f = lambda x: x**3 - 2*x**2 + x
f_prime_analytical = lambda x: 3*x**2 - 4*x + 1  # by hand

x = 2.0
analytical = f_prime_analytical(x)
numerical = numerical_derivative(f, x)

print(f"Analytical derivative at x=2: {analytical}")   # 5.0
print(f"Numerical derivative at x=2:  {numerical:.6f}")  # 5.000000
print(f"Error: {abs(analytical - numerical):.2e}")       # ~1e-10

Key takeaway: Numerical derivatives are how you verify your analytical gradient implementations (a technique called “gradient checking” — critical for debugging custom loss functions).

---

2. The Chain Rule — Backpropagation’s Foundation

The chain rule for h(x) = f(g(x)):

h'(x) = f'(g(x)) · g'(x)

This is it. Backpropagation in a neural network is the chain rule applied to a composed function. Nothing more.

# Chain rule: derivative of sin(x²) at x = π/4
# Let f(u) = sin(u), g(x) = x²
# h(x) = f(g(x)) = sin(x²)
# h'(x) = f'(g(x)) * g'(x) = cos(x²) * 2x

x = np.pi / 4
h_prime_analytical = np.cos(x**2) * 2 * x

h = lambda x: np.sin(x**2)
h_prime_numerical = numerical_derivative(h, x)

print(f"Analytical (chain rule): {h_prime_analytical:.6f}")
print(f"Numerical:               {h_prime_numerical:.6f}")

A 3-layer network (the actual chain rule application):

Loss = L(a₃)
a₃ = f₃(z₃)
z₃ = W₃ @ a₂ + b₃
a₂ = f₂(z₂)
...

∂L/∂W₁ = ∂L/∂a₃ · ∂a₃/∂z₃ · ∂z₃/∂a₂ · ∂a₂/∂z₂ · ∂z₂/∂W₁

Every backprop algorithm is this chain — computed right-to-left (hence “back-propagation”).

---

3. Partial Derivatives and Gradients

For a function of multiple variables f(x₁, x₂, ..., xₙ), the gradient is the vector of all partial derivatives:

∇f = [∂f/∂x₁, ∂f/∂x₂, ..., ∂f/∂xₙ]

The gradient points in the direction of steepest ascent. Gradient descent moves in the negative gradient direction.

def gradient_of_mse(X, y, w):
    """
    Manually compute gradient of MSE loss w.r.t. weights w.
    
    MSE = (1/n) * ||Xw - y||²
    ∂MSE/∂w = (2/n) * Xᵀ(Xw - y)
    
    Derivation:
    L = (1/n) * (Xw - y)ᵀ(Xw - y)
      = (1/n) * (wᵀXᵀXw - 2yᵀXw + yᵀy)
    ∂L/∂w = (2/n) * (XᵀXw - Xᵀy)
           = (2/n) * Xᵀ(Xw - y)
    """
    n = len(y)
    predictions = X @ w
    residuals = predictions - y
    gradient = (2 / n) * X.T @ residuals
    return gradient

# Verify with numerical gradient
np.random.seed(42)
n, d = 100, 3
X = np.random.randn(n, d)
w_true = np.array([1.0, -2.0, 0.5])
y = X @ w_true + 0.1 * np.random.randn(n)
w_test = np.random.randn(d)

mse = lambda w: np.mean((X @ w - y)**2)
analytical_grad = gradient_of_mse(X, y, w_test)

# Numerical gradient check
numerical_grad = np.zeros(d)
h = 1e-5
for i in range(d):
    w_plus = w_test.copy(); w_plus[i] += h
    w_minus = w_test.copy(); w_minus[i] -= h
    numerical_grad[i] = (mse(w_plus) - mse(w_minus)) / (2 * h)

print(f"Analytical gradient: {analytical_grad}")
print(f"Numerical gradient:  {numerical_grad}")
print(f"Max error: {np.max(np.abs(analytical_grad - numerical_grad)):.2e}")
# Should be < 1e-6

---

4. Gradient Descent — Full Implementation

def gradient_descent(X, y, lr=0.01, n_iter=1000, tolerance=1e-6):
    """
    Gradient descent for linear regression.
    
    Returns: weight history, loss history
    """
    n, d = X.shape
    w = np.zeros(d)  # initialization
    
    loss_history = []
    w_history = [w.copy()]
    
    for iteration in range(n_iter):
        # Forward pass: compute predictions
        predictions = X @ w
        
        # Compute loss
        loss = np.mean((predictions - y)**2)
        loss_history.append(loss)
        
        # Backward pass: compute gradient
        gradient = gradient_of_mse(X, y, w)
        
        # Update step
        w = w - lr * gradient
        w_history.append(w.copy())
        
        # Convergence check
        if iteration > 0 and abs(loss_history[-2] - loss_history[-1]) < tolerance:
            print(f"Converged at iteration {iteration}")
            break
    
    return np.array(w_history), np.array(loss_history)

w_history, loss_history = gradient_descent(X, y, lr=0.01, n_iter=2000)

print(f"Learned weights: {w_history[-1]}")
print(f"True weights:    {w_true}")
print(f"Final MSE: {loss_history[-1]:.6f}")

# Visualize convergence
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.semilogy(loss_history)
plt.xlabel('Iteration')
plt.ylabel('MSE Loss (log scale)')
plt.title('Loss Convergence')
plt.grid(True)

plt.subplot(1, 2, 2)
# Plot weight trajectory (first two weights)
plt.plot(w_history[:, 0], w_history[:, 1], 'b.-', alpha=0.5)
plt.scatter([w_true[0]], [w_true[1]], c='red', s=100, zorder=5, label='True weights')
plt.xlabel('w₀')
plt.ylabel('w₁')
plt.title('Weight Trajectory')
plt.legend()
plt.tight_layout()
plt.savefig('gradient_descent_convergence.png', dpi=150)

---

5. Jacobian and Hessian — Why They Matter

Jacobian: For a vector-valued function f: Rⁿ → Rᵐ, the Jacobian J is the matrix of all partial derivatives. Shape: (m, n).

• In backprop: the Jacobian of a layer’s output w.r.t. its input is what gets multiplied during gradient propagation.

• In batch normalization: the Jacobian of normalized outputs is notoriously complex to derive (this is why BN is tricky to implement from scratch).

Hessian: For a scalar function f: Rⁿ → R, the Hessian H is the matrix of second derivatives.

Hᵢⱼ = ∂²f / (∂xᵢ ∂xⱼ)

def hessian_numerical(f, x, h=1e-4):
    """Numerical Hessian via finite differences."""
    n = len(x)
    H = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            x_pp = x.copy(); x_pp[i] += h; x_pp[j] += h
            x_pm = x.copy(); x_pm[i] += h; x_pm[j] -= h
            x_mp = x.copy(); x_mp[i] -= h; x_mp[j] += h
            x_mm = x.copy(); x_mm[i] -= h; x_mm[j] -= h
            H[i, j] = (f(x_pp) - f(x_pm) - f(x_mp) + f(x_mm)) / (4 * h**2)
    return H

# Quadratic: f(x, y) = x² + 2y² + xy
f = lambda x: x[0]**2 + 2*x[1]**2 + x[0]*x[1]
x0 = np.array([1.0, 2.0])
H = hessian_numerical(f, x0)
print(f"Hessian:\n{H}")
# Expected: [[2, 1], [1, 4]] (constants for quadratic)

Why the Hessian matters in ML:

• Eigenvalues of the Hessian at a loss minimum describe the loss landscape curvature

• Large max eigenvalue / small min eigenvalue (high condition number) → poorly conditioned optimization → slow convergence → need for adaptive learning rates

• This is the actual mathematical reason Adam outperforms SGD on some tasks

---

6. Visualizing Loss Landscapes

def visualize_loss_surface(X, y):
    """Visualize the MSE loss landscape for 2-parameter linear regression."""
    w0_range = np.linspace(-1, 3, 100)
    w1_range = np.linspace(-4, 0, 100)
    W0, W1 = np.meshgrid(w0_range, w1_range)
    
    loss_surface = np.zeros_like(W0)
    for i in range(len(w0_range)):
        for j in range(len(w1_range)):
            w = np.array([W0[j, i], W1[j, i]])
            # Use only first 2 features for visualization
            loss_surface[j, i] = np.mean((X[:, :2] @ w - y)**2)
    
    plt.figure(figsize=(10, 4))
    plt.subplot(1, 2, 1)
    plt.contourf(W0, W1, loss_surface, levels=50, cmap='viridis')
    plt.colorbar(label='MSE Loss')
    plt.xlabel('w₀'); plt.ylabel('w₁')
    plt.title('Loss Surface (contour)')
    
    plt.subplot(1, 2, 2)
    ax = plt.subplot(1, 2, 2, projection='3d') if False else plt.subplot(1, 2, 2)
    plt.contour(W0, W1, np.log(loss_surface + 1e-10), levels=30, cmap='plasma')
    plt.xlabel('w₀'); plt.ylabel('w₁')
    plt.title('Log-Loss Surface')
    plt.tight_layout()
    plt.savefig('loss_landscape.png', dpi=150)

---

What Most Learners Get Wrong

They treat calculus as abstract symbols rather than as the update rule.

The derivative ∂L/∂w is not a number you compute once. It is a function that you evaluate at the current weights to determine which direction to move. Understanding this distinction — between a derivative function and a derivative value — is why some people can derive new loss functions and others can only use predefined ones.

Second misconception: believing that gradient descent finds a global minimum. For non-convex loss surfaces (every neural network), gradient descent finds a local minimum, not a global one. The practical reason deep networks work is that most local minima in high-dimensional spaces are approximately equally good (a result from random matrix theory). Understanding this requires calculus + linear algebra together.

Third gap: not knowing how to gradient-check. Every time you implement a new loss function or layer, you should run a gradient check (numerical vs analytical). The code for this is above. If you skip this step, you will spend hours debugging training that silently fails because of a sign error in your backprop.

---

Practice Problems with Acceptance Criteria

Problem 1 — Manual Gradient Derivation Derive the gradient of binary cross-entropy loss L = -[y log(p) + (1-y) log(1-p)] with respect to the pre-sigmoid logit z (where p = sigmoid(z)).

• Acceptance criteria: Arrive at the clean result ∂L/∂z = p - y. Verify numerically with gradient checking (np.allclose within 1e-6).

Problem 2 — Learning Rate Ablation Implement gradient descent on the 2-parameter MSE surface above. Run with lr ∈ {0.0001, 0.001, 0.01, 0.1, 1.0}. Plot convergence curves for all 5.

• Acceptance criteria: Correctly identify which learning rates converge, which diverge, and which oscillate. Explain the geometry of why lr=1.0 diverges using the Hessian eigenvalue.

Problem 3 — Gradient of L2-Regularized Loss Derive and implement the gradient of L_reg = MSE + λ||w||² w.r.t. w. Run gradient descent with λ ∈ {0, 0.01, 0.1} and plot how the learned weights change.

• Acceptance criteria: Gradient implementation verified numerically. Plot shows weight shrinkage increasing with λ. Explain what happens as λ → ∞.

---

Resources

Resource	Focus	Time	Free?
3Blue1Brown “Essence of Calculus”	Geometric intuition for derivatives	~3 hrs	✅
Khan Academy Multivariable Calculus	Partial derivatives, gradients	~5 hrs	✅
Stewart “Early Transcendentals” Ch. 14–15	Rigorous multivariable treatment	~8 hrs	❌
Andrej Karpathy “micrograd” repo	Backprop from scratch in 100 lines	~3 hrs	✅
Practice problems above	Cement with code	~5 hrs	—

Highest-ROI single resource: Karpathy’s micrograd (github.com/karpathy/micrograd). It’s 100 lines of Python that implement full automatic differentiation and backprop. Reading and re-implementing it is worth 20 hours of textbook calculus for building backprop intuition.

---

Return to README.md · Next: 03_probability_and_statistics.md