Part 3: Calculus and Optimization

Part of the Mathematics for Programming 101 Series

The Debugging Session That Taught Me Calculus

I was training a neural network for image classification. Loss was decreasing, then suddenly spiking to infinity.

Epoch 1: loss = 2.305
Epoch 2: loss = 1.823
Epoch 3: loss = 1.456
Epoch 4: loss = 1.203
Epoch 5: loss = inf

Stack Overflow said: "Lower your learning rate."

I tried 0.01, 0.001, 0.0001. Still exploded. Sometimes earlier, sometimes later, but always the same infinity.

Then I actually looked at the gradients. They were growing exponentially through the layers. The problem wasn't the learning rate—it was gradient explosion in deep networks.

The solution: Gradient clipping, better initialization, and residual connections.

But here's the key: I only understood the problem when I understood the calculus. Derivatives weren't just theory—they were the diagnostic tool I needed.

What is Calculus in Programming?

Calculus is the mathematics of change. In programming:

• Derivatives = Rate of change = Slope = Gradient

• Integrals = Accumulation = Total change over time

You use calculus when:

• Training machine learning models

• Optimizing any function

• Understanding algorithm behavior

• Simulating physics or motion

• Analyzing trends and patterns

Derivatives: Measuring Change

The Intuition

A derivative measures: "If I change the input slightly, how much does the output change?"

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    """Simple function: f(x) = x²"""
    return x**2

def numerical_derivative(f, x, h=1e-7):
    """
    Approximate derivative using finite differences
    f'(x) ≈ [f(x+h) - f(x)] / h
    """
    return (f(x + h) - f(x)) / h

# Calculate derivative at x = 3
x = 3
derivative_at_3 = numerical_derivative(f, x)
analytical_derivative = 2 * x  # f'(x) = 2x for f(x) = x²

print(f"Numerical derivative at x={x}: {derivative_at_3:.6f}")
print(f"Analytical derivative at x={x}: {analytical_derivative}")
print(f"Difference: {abs(derivative_at_3 - analytical_derivative):.10f}")

# Visualize
x_vals = np.linspace(0, 5, 100)
y_vals = f(x_vals)

plt.figure(figsize=(10, 6))
plt.plot(x_vals, y_vals, label='f(x) = x²')

# Tangent line at x=3
tangent_slope = 2 * x
tangent_y = f(x) + tangent_slope * (x_vals - x)
plt.plot(x_vals, tangent_y, 'r--', label=f'Tangent at x={x} (slope={tangent_slope})')

plt.scatter([x], [f(x)], color='red', s=100, zorder=5)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()
plt.grid(True)
plt.title('Derivative = Slope of Tangent Line')
plt.savefig('derivative_visualization.png', dpi=100, bbox_inches='tight')
print("Visualization saved!")

Key insight: The derivative at a point is the slope of the tangent line at that point.

Common Derivatives

# Power rule: d/dx(x^n) = n*x^(n-1)
def power_rule(x, n):
    return n * (x ** (n - 1))

# Examples
print(f"d/dx(x²) at x=3: {power_rule(3, 2)}")  # 2*3 = 6
print(f"d/dx(x³) at x=2: {power_rule(2, 3)}")  # 3*4 = 12

# Exponential: d/dx(e^x) = e^x
def exp_derivative(x):
    return np.exp(x)

# Logarithm: d/dx(ln(x)) = 1/x
def log_derivative(x):
    return 1 / x

# Sine: d/dx(sin(x)) = cos(x)
def sin_derivative(x):
    return np.cos(x)

# Cosine: d/dx(cos(x)) = -sin(x)
def cos_derivative(x):
    return -np.sin(x)

Chain Rule: The Key to Deep Learning

The chain rule is how backpropagation works.

If you have composed functions: h(x) = f(g(x)), then:

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

def chain_rule_example():
    """
    Example: h(x) = (2x + 1)³
    Let g(x) = 2x + 1 and f(u) = u³
    Then h(x) = f(g(x))
    
    h'(x) = f'(g(x)) * g'(x) = 3*(2x+1)² * 2
    """
    def h(x):
        return (2*x + 1)**3
    
    def h_derivative(x):
        # Outer derivative: 3u²
        outer = 3 * (2*x + 1)**2
        # Inner derivative: 2
        inner = 2
        return outer * inner
    
    x = 2
    numerical = numerical_derivative(h, x)
    analytical = h_derivative(x)
    
    print(f"Numerical: {numerical:.6f}")
    print(f"Analytical: {analytical}")
    print(f"Match? {np.isclose(numerical, analytical)}")

chain_rule_example()

Why this matters for neural networks:

# Simple neural network with one hidden layer
def neural_net(x, w1, w2):
    """
    x → (linear) → (relu) → (linear) → output
    f(x) = w2 * max(0, w1 * x)
    """
    hidden = np.maximum(0, w1 * x)  # ReLU activation
    output = w2 * hidden
    return output

# To train this, we need gradients via chain rule:
# ∂output/∂w1 = ∂output/∂hidden * ∂hidden/∂w1

Partial Derivatives and Gradients

When you have multiple inputs, you need partial derivatives.

def f(x, y):
    """Function of two variables: f(x,y) = x²y + y³"""
    return x**2 * y + y**3

def partial_x(x, y, h=1e-7):
    """∂f/∂x: derivative with respect to x, holding y constant"""
    return (f(x + h, y) - f(x, y)) / h

def partial_y(x, y, h=1e-7):
    """∂f/∂y: derivative with respect to y, holding x constant"""
    return (f(x, y + h) - f(x, y)) / h

# Gradient: vector of all partial derivatives
def gradient(x, y):
    """∇f = [∂f/∂x, ∂f/∂y]"""
    return np.array([partial_x(x, y), partial_y(x, y)])

# Example
x, y = 2, 3
grad = gradient(x, y)

print(f"∂f/∂x at ({x},{y}): {grad[0]:.3f}")
print(f"∂f/∂y at ({x},{y}): {grad[1]:.3f}")
print(f"Gradient: {grad}")

# Analytical gradient for verification
analytical_grad = np.array([2*x*y, x**2 + 3*y**2])
print(f"Analytical gradient: {analytical_grad}")

Gradient = direction of steepest ascent Negative gradient = direction of steepest descent (for minimization)

Gradient Descent: Optimizing Functions

The algorithm that powers ML training:

def gradient_descent_1d(f, df, x_init, learning_rate=0.1, iterations=100):
    """
    Minimize function f using its derivative df
    
    Args:
        f: function to minimize
        df: derivative of f
        x_init: starting point
        learning_rate: step size
        iterations: number of steps
    """
    x = x_init
    history = [x]
    
    for i in range(iterations):
        # Compute gradient
        grad = df(x)
        
        # Update: move in opposite direction of gradient
        x = x - learning_rate * grad
        
        history.append(x)
        
        if i % 20 == 0:
            print(f"Iteration {i}: x={x:.4f}, f(x)={f(x):.4f}")
    
    return x, history

# Example: Minimize f(x) = (x - 3)²
def f(x):
    return (x - 3)**2

def df(x):
    return 2 * (x - 3)

# Start from x=0, should converge to x=3
final_x, history = gradient_descent_1d(f, df, x_init=0, learning_rate=0.1)

print(f"\nFinal x: {final_x:.6f} (should be ~3.0)")
print(f"Final f(x): {f(final_x):.6f} (should be ~0.0)")

Multidimensional Gradient Descent

def gradient_descent_nd(f, grad_f, x_init, learning_rate=0.01, iterations=1000):
    """
    Minimize function f in multiple dimensions
    
    Args:
        f: function to minimize
        grad_f: function that computes gradient
        x_init: starting point (array)
        learning_rate: step size
        iterations: number of steps
    """
    x = x_init.copy()
    history = [x.copy()]
    
    for i in range(iterations):
        # Compute gradient
        grad = grad_f(x)
        
        # Update
        x = x - learning_rate * grad
        
        history.append(x.copy())
        
        if i % 200 == 0:
            print(f"Iteration {i}: x={x}, f(x)={f(x):.4f}")
    
    return x, np.array(history)

# Example: Minimize f(x,y) = (x-1)² + (y-2)²
def f_2d(x):
    return (x[0] - 1)**2 + (x[1] - 2)**2

def grad_f_2d(x):
    return np.array([
        2 * (x[0] - 1),  # ∂f/∂x
        2 * (x[1] - 2)   # ∂f/∂y
    ])

# Start from origin, should converge to (1, 2)
x_init = np.array([0.0, 0.0])
final_x, history = gradient_descent_nd(f_2d, grad_f_2d, x_init, learning_rate=0.1, iterations=100)

print(f"\nFinal x: {final_x} (should be ~[1, 2])")
print(f"Final f(x): {f_2d(final_x):.6f} (should be ~0.0)")

Real Example: Linear Regression from Scratch

class LinearRegression:
    """Linear regression using gradient descent"""
    
    def __init__(self, learning_rate=0.01):
        self.lr = learning_rate
        self.weights = None
        self.bias = None
    
    def fit(self, X, y, iterations=1000):
        """Train using gradient descent"""
        n_samples, n_features = X.shape
        
        # Initialize parameters
        self.weights = np.zeros(n_features)
        self.bias = 0
        
        # Gradient descent
        for i in range(iterations):
            # Forward pass: predictions
            y_pred = X @ self.weights + self.bias
            
            # Compute loss (Mean Squared Error)
            loss = np.mean((y_pred - y)**2)
            
            # Compute gradients
            dw = (2/n_samples) * X.T @ (y_pred - y)
            db = (2/n_samples) * np.sum(y_pred - y)
            
            # Update parameters
            self.weights -= self.lr * dw
            self.bias -= self.lr * db
            
            if i % 100 == 0:
                print(f"Iteration {i}: Loss = {loss:.4f}")
        
        return self
    
    def predict(self, X):
        """Make predictions"""
        return X @ self.weights + self.bias

# Generate sample data
np.random.seed(42)
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X.flatten() + np.random.randn(100) * 0.5

# Train model
model = LinearRegression(learning_rate=0.1)
model.fit(X, y, iterations=500)

print(f"\nLearned weights: {model.weights}")
print(f"Learned bias: {model.bias:.3f}")
print(f"True values: weights=[3.0], bias=4.0")

# Predictions
y_pred = model.predict(X)
final_loss = np.mean((y_pred - y)**2)
print(f"Final MSE: {final_loss:.4f}")

Backpropagation: Chain Rule in Action

Building a simple neural network with backpropagation:

class SimpleNeuralNet:
    """Two-layer neural network with backpropagation"""
    
    def __init__(self, input_dim, hidden_dim, output_dim, learning_rate=0.01):
        # Initialize weights
        self.W1 = np.random.randn(input_dim, hidden_dim) * 0.01
        self.b1 = np.zeros(hidden_dim)
        self.W2 = np.random.randn(hidden_dim, output_dim) * 0.01
        self.b2 = np.zeros(output_dim)
        self.lr = learning_rate
    
    def relu(self, x):
        """ReLU activation"""
        return np.maximum(0, x)
    
    def relu_derivative(self, x):
        """Derivative of ReLU"""
        return (x > 0).astype(float)
    
    def forward(self, X):
        """Forward pass - save intermediate values for backprop"""
        # Layer 1
        self.Z1 = X @ self.W1 + self.b1
        self.A1 = self.relu(self.Z1)
        
        # Layer 2
        self.Z2 = self.A1 @ self.W2 + self.b2
        
        return self.Z2
    
    def backward(self, X, y, output):
        """Backpropagation - compute gradients using chain rule"""
        m = X.shape[0]
        
        # Output layer gradient
        dZ2 = output - y  # Simplified for MSE loss
        dW2 = (1/m) * self.A1.T @ dZ2
        db2 = (1/m) * np.sum(dZ2, axis=0)
        
        # Hidden layer gradient (chain rule!)
        dA1 = dZ2 @ self.W2.T
        dZ1 = dA1 * self.relu_derivative(self.Z1)
        dW1 = (1/m) * X.T @ dZ1
        db1 = (1/m) * np.sum(dZ1, axis=0)
        
        return dW1, db1, dW2, db2
    
    def train_step(self, X, y):
        """Single training step"""
        # Forward pass
        output = self.forward(X)
        
        # Compute loss
        loss = np.mean((output - y)**2)
        
        # Backward pass
        dW1, db1, dW2, db2 = self.backward(X, y, output)
        
        # Update weights
        self.W1 -= self.lr * dW1
        self.b1 -= self.lr * db1
        self.W2 -= self.lr * dW2
        self.b2 -= self.lr * db2
        
        return loss
    
    def fit(self, X, y, epochs=1000):
        """Train the network"""
        for epoch in range(epochs):
            loss = self.train_step(X, y)
            
            if epoch % 100 == 0:
                print(f"Epoch {epoch}: Loss = {loss:.4f}")
        
        return self

# Example: XOR problem (non-linearly separable)
X = np.array([[0, 0], [0, 1], [1, 0], [1, 1]])
y = np.array([[0], [1], [1], [0]])

# Train network
nn = SimpleNeuralNet(input_dim=2, hidden_dim=4, output_dim=1, learning_rate=0.5)
nn.fit(X, y, epochs=5000)

# Test predictions
predictions = nn.forward(X)
print("\nPredictions:")
for i in range(len(X)):
    print(f"Input: {X[i]}, Target: {y[i][0]}, Prediction: {predictions[i][0]:.3f}")

Optimization Algorithms

Momentum

Problem: Gradient descent oscillates in ravines. Solution: Accumulate velocity from past gradients.

def gradient_descent_momentum(grad_f, x_init, learning_rate=0.01, 
                               momentum=0.9, iterations=1000):
    """Gradient descent with momentum"""
    x = x_init.copy()
    velocity = np.zeros_like(x)
    
    for i in range(iterations):
        grad = grad_f(x)
        
        # Update velocity (exponentially weighted average of gradients)
        velocity = momentum * velocity - learning_rate * grad
        
        # Update position
        x = x + velocity
        
        if i % 200 == 0:
            print(f"Iteration {i}: x={x}")
    
    return x

# Test with same function as before
x_momentum = gradient_descent_momentum(grad_f_2d, np.array([0.0, 0.0]), 
                                       learning_rate=0.1, iterations=100)
print(f"Final x with momentum: {x_momentum}")

Adam Optimizer

Combines momentum and adaptive learning rates:

def adam_optimizer(grad_f, x_init, learning_rate=0.001, 
                   beta1=0.9, beta2=0.999, epsilon=1e-8, iterations=1000):
    """Adam optimizer implementation"""
    x = x_init.copy()
    m = np.zeros_like(x)  # First moment (mean)
    v = np.zeros_like(x)  # Second moment (variance)
    
    for t in range(1, iterations + 1):
        grad = grad_f(x)
        
        # Update biased first moment estimate
        m = beta1 * m + (1 - beta1) * grad
        
        # Update biased second moment estimate
        v = beta2 * v + (1 - beta2) * (grad**2)
        
        # Bias correction
        m_hat = m / (1 - beta1**t)
        v_hat = v / (1 - beta2**t)
        
        # Update parameters
        x = x - learning_rate * m_hat / (np.sqrt(v_hat) + epsilon)
        
        if t % 200 == 0:
            print(f"Iteration {t}: x={x}")
    
    return x

# Test Adam
x_adam = adam_optimizer(grad_f_2d, np.array([0.0, 0.0]), iterations=100)
print(f"Final x with Adam: {x_adam}")

Debugging with Calculus

Gradient Checking

Verify backpropagation implementation:

def gradient_check(f, grad_f, x, epsilon=1e-7):
    """
    Check if analytical gradient matches numerical gradient
    
    Args:
        f: function
        grad_f: analytical gradient function
        x: point to check
        epsilon: small value for numerical approximation
    """
    analytical_grad = grad_f(x)
    numerical_grad = np.zeros_like(x)
    
    # Compute numerical gradient for each dimension
    for i in range(len(x)):
        x_plus = x.copy()
        x_plus[i] += epsilon
        
        x_minus = x.copy()
        x_minus[i] -= epsilon
        
        numerical_grad[i] = (f(x_plus) - f(x_minus)) / (2 * epsilon)
    
    # Compare
    difference = np.linalg.norm(analytical_grad - numerical_grad)
    relative_difference = difference / (np.linalg.norm(analytical_grad) + 
                                       np.linalg.norm(numerical_grad))
    
    print(f"Analytical gradient: {analytical_grad}")
    print(f"Numerical gradient: {numerical_grad}")
    print(f"Relative difference: {relative_difference:.10f}")
    
    if relative_difference < 1e-7:
        print("✓ Gradient check passed!")
    else:
        print("✗ Gradient check failed!")
    
    return relative_difference

# Test with our 2D function
x_test = np.array([1.5, 2.5])
gradient_check(f_2d, grad_f_2d, x_test)

Diagnosing Training Issues

def diagnose_training(losses, gradients):
    """Diagnose common training problems"""
    
    # Check for exploding gradients
    max_grad = np.max(np.abs(gradients))
    if max_grad > 10:
        print(f"⚠ Exploding gradients detected! Max: {max_grad:.2f}")
        print("Solutions: Gradient clipping, lower learning rate, batch norm")
    
    # Check for vanishing gradients
    min_grad = np.min(np.abs(gradients))
    if min_grad < 1e-7:
        print(f"⚠ Vanishing gradients detected! Min: {min_grad:.2e}")
        print("Solutions: ReLU instead of sigmoid, residual connections, better initialization")
    
    # Check for stagnant loss
    recent_losses = losses[-10:]
    if len(recent_losses) > 1:
        loss_change = np.std(recent_losses)
        if loss_change < 1e-6:
            print(f"⚠ Loss plateaued! Std: {loss_change:.2e}")
            print("Solutions: Increase learning rate, different optimizer, check data")
    
    # Check for oscillating loss
    if len(losses) > 3:
        oscillation = np.mean(np.abs(np.diff(losses[-20:])))
        avg_loss = np.mean(losses[-20:])
        if oscillation / avg_loss > 0.1:
            print(f"⚠ Loss oscillating! Oscillation/Loss ratio: {oscillation/avg_loss:.2%}")
            print("Solutions: Lower learning rate, momentum, better batch size")

Key Takeaways

• Derivatives measure rates of change—essential for optimization

• Chain rule enables backpropagation in neural networks

• Gradients point in direction of steepest ascent (negative for descent)

• Gradient descent iteratively minimizes functions

• Advanced optimizers (momentum, Adam) converge faster and more reliably

• Gradient checking verifies backpropagation implementations

• Understanding calculus helps debug training issues

What's Next

In the next article, we'll explore probability and statistics—the mathematics behind uncertainty, decision-making, A/B testing, and recommendation systems.

You'll learn:

• Probability distributions in practice

• Bayesian thinking for better decisions

• Statistical testing that actually works

• Building anomaly detection systems

Continue to Part 4: Probability and Statistics →

---

Navigation

• ← Part 2: Linear Algebra

• Back to Series Overview

• Part 4: Probability and Statistics →