Calculus for Machine Learning

Studies rates of change. Given a function f(x), the derivative f'(x) tells us how fast f is changing at each point x. This is how we compute gradients to train models.

Studies accumulation. Given a rate of change, integration recovers the total quantity. In ML, integrals appear in probability distributions, expected values, and KL divergence.

Straight line. The simplest model. Linear regression is just finding the best m and b.

Parabola. MSE loss is quadratic — this is why gradient descent works so well on it (convex!).

S-curve. Squashes any input to (0, 1). Used in binary classification and gating.

The most popular activation. Zero for negatives, identity for positives. Simple derivative.

import numpy as np

# Define common functions
x = np.linspace(-5, 5, 100)

linear  = 2 * x + 1
quad    = x ** 2
sigmoid = 1 / (1 + np.exp(-x))
relu    = np.maximum(0, x)

# Evaluate at a point
sigmoid_at_0 = 1 / (1 + np.exp(0))
print(sigmoid_at_0)  # 0.5

import math

def linear(x, m=2, b=1):
    return m * x + b

def sigmoid(x):
    return 1 / (1 + math.exp(-x))

def relu(x):
    return max(0, x)

# Plot manually: generate (x, y) pairs
xs = [i * 0.1 for i in range(-50, 51)]
ys = [sigmoid(x) for x in xs]
print(sigmoid(0))  # 0.5

# lim x→0 sin(x)/x = 1
import numpy as np

# Approach from right
h_values = [0.1, 0.01, 0.001, 0.0001]
for h in h_values:
    print(f"h={h}: {np.sin(h)/h:.8f}")

# h=0.1:    0.99833417
# h=0.01:   0.99998333
# h=0.001:  0.99999983
# h=0.0001: 1.00000000
# → Approaches 1!

import math

def numerical_limit(f, a, h_values=None):
    """Estimate lim x→a f(x)"""
    if h_values is None:
        h_values = [10**-i
                    for i in range(1,9)]
    results = []
    for h in h_values:
        results.append(f(a + h))
    return results[-1]

# lim x→0 sin(x)/x
limit = numerical_limit(
    lambda x: math.sin(x)/x, 0
)
print(limit)  # ≈ 1.0

import numpy as np

# f(x) = x² → f'(x) = 2x
def f(x): return x**2
def f_prime(x): return 2*x

x = np.array([-2, -1, 0, 1, 2])
print(f"f(x)  = {f(x)}")
# [4, 1, 0, 1, 4]
print(f"f'(x) = {f_prime(x)}")
# [-4, -2, 0, 2, 4]

# At x=-2: slope is -4 (steeply downhill)
# At x=0:  slope is  0 (minimum!)
# At x=2:  slope is  4 (steeply uphill)

# Numerical derivative approximation
def numerical_derivative(f, x, h=1e-7):
    """Central difference formula"""
    return (f(x + h) - f(x - h)) / (2 * h)

def f(x): return x ** 2

# Test at several points
for x in [-2, -1, 0, 1, 2]:
    slope = numerical_derivative(f, x)
    print(f"x={x}: f'={slope:.4f}")

# x=-2: f'=-4.0000
# x=-1: f'=-2.0000
# x= 0: f'= 0.0000  ← minimum!
# x= 1: f'= 2.0000
# x= 2: f'= 4.0000

import numpy as np

def f(x): return x**2

x = 3.0; h = 1e-7

# Forward difference
fwd = (f(x+h) - f(x)) / h
# Backward difference
bwd = (f(x) - f(x-h)) / h
# Central difference (best!)
ctr = (f(x+h) - f(x-h)) / (2*h)

print(f"Forward:  {fwd:.10f}")
print(f"Backward: {bwd:.10f}")
print(f"Central:  {ctr:.10f}")
print(f"Exact:    {2*x:.10f}")
# Central is most accurate!

def numerical_gradient(f, params, h=1e-7):
    """Gradient for a list of params"""
    grad = []
    for i in range(len(params)):
        # Perturb param i
        params_plus = params[:]
        params_minus = params[:]
        params_plus[i] += h
        params_minus[i] -= h
        # Central difference
        g = (f(params_plus) - f(params_minus))
        g /= (2 * h)
        grad.append(g)
    return grad

# f(w0, w1) = w0² + 3*w1²
def loss(w):
    return w[0]**2 + 3*w[1]**2

g = numerical_gradient(loss, [2.0, 1.0])
print(g)  # [4.0, 6.0]

import numpy as np

# f(x) = 3x⁴ + 2x² − 7x + 5
# f'(x) = 12x³ + 4x − 7

def f(x):
    return 3*x**4 + 2*x**2 - 7*x + 5

def f_prime(x):
    return 12*x**3 + 4*x - 7

x = np.array([0, 1, 2])
print(f_prime(x))  # [-7, 9, 97]

def power_rule(coeff, exp):
    """Apply power rule: d/dx [c·x^n]"""
    if exp == 0:
        return (0, 0)  # constant
    new_coeff = coeff * exp
    new_exp = exp - 1
    return (new_coeff, new_exp)

# Differentiate 3x⁴
c, e = power_rule(3, 4)
print(f"{c}x^{e}")
# 12x^3

# Differentiate √x = x^0.5
c, e = power_rule(1, 0.5)
print(f"{c}x^{e}")
# 0.5x^-0.5

import numpy as np

# f(x) = x² · eˣ
# f'(x) = eˣ(2x + x²)

def f(x):
    return x**2 * np.exp(x)

def f_prime(x):
    return np.exp(x) * (2*x + x**2)

x = np.array([0, 1, 2])
print(f_prime(x))
# [0.0, 8.154, 43.510]

import math

def product_rule(f, f_p, g, g_p, x):
    """(fg)' = f'g + fg'"""
    return f_p(x)*g(x) + f(x)*g_p(x)

# f=x², f'=2x, g=eˣ, g'=eˣ
result = product_rule(
    lambda x: x**2,
    lambda x: 2*x,
    lambda x: math.exp(x),
    lambda x: math.exp(x),
    1.0
)
print(result)  # 8.1548

import numpy as np

# f(x) = sigmoid(3x + 2)
# outer: sigmoid(u), inner: 3x+2

def sigmoid(x):
    return 1 / (1 + np.exp(-x))

def sigmoid_prime(x):
    s = sigmoid(x)
    return s * (1 - s)

# Chain rule: σ'(3x+2) · 3
x = np.array([0, 1, 2])
inner = 3*x + 2
deriv = sigmoid_prime(inner) * 3
print(deriv)
# [0.315, 0.044, 0.005]

def chain_rule(outer_deriv, inner_deriv,
              inner_fn, x):
    """d/dx f(g(x)) = f'(g(x))·g'(x)"""
    u = inner_fn(x)        # g(x)
    df_du = outer_deriv(u)  # f'(g(x))
    du_dx = inner_deriv(x)  # g'(x)
    return df_du * du_dx

import math
# d/dx sin(x²) = cos(x²)·2x
result = chain_rule(
    lambda u: math.cos(u),  # f'
    lambda x: 2*x,          # g'
    lambda x: x**2,         # g
    1.0
)
print(result)  # ≈ 1.0806

import numpy as np

def sigmoid(x):
    return 1 / (1 + np.exp(-x))

def sigmoid_deriv(x):
    s = sigmoid(x)
    return s * (1 - s)

x = np.linspace(-6, 6, 5)
print("σ(x):  ", np.round(sigmoid(x), 4))
print("σ'(x): ", np.round(sigmoid_deriv(x), 4))
# σ(x):   [0.0025, 0.0474, 0.5, 0.9526, 0.9975]
# σ'(x):  [0.0025, 0.0452, 0.25, 0.0452, 0.0025]
# ↑ Derivative peaks at 0.25 (x=0) and vanishes at extremes!

import math

def sigmoid(x):
    # Numerically stable version
    if x >= 0:
        z = math.exp(-x)
        return 1 / (1 + z)
    else:
        z = math.exp(x)
        return z / (1 + z)

def sigmoid_deriv(x):
    s = sigmoid(x)
    return s * (1 - s)

# Verify max at x=0
print(sigmoid_deriv(0))   # 0.25
print(sigmoid_deriv(5))   # 0.0066
print(sigmoid_deriv(-5))  # 0.0066

import numpy as np

def tanh(x):
    return np.tanh(x)

def tanh_deriv(x):
    t = np.tanh(x)
    return 1 - t**2

x = np.array([-2, -1, 0, 1, 2])
print("tanh:  ", np.round(tanh(x), 4))
print("tanh': ", np.round(tanh_deriv(x), 4))
# tanh:   [-0.9640, -0.7616, 0.0, 0.7616, 0.9640]
# tanh':  [0.0707, 0.4200, 1.0, 0.4200, 0.0707]

import math

def tanh(x):
    ep = math.exp(x)
    em = math.exp(-x)
    return (ep - em) / (ep + em)

def tanh_deriv(x):
    t = tanh(x)
    return 1 - t**2

# Verify: max deriv at x=0 is 1.0
print(tanh_deriv(0))  # 1.0
print(tanh_deriv(3))  # 0.0099

import numpy as np

def relu(x): return np.maximum(0, x)
def relu_deriv(x): return (x > 0).astype(float)

def leaky_relu(x, a=0.01):
    return np.where(x > 0, x, a*x)
def leaky_relu_deriv(x, a=0.01):
    return np.where(x > 0, 1, a)

def gelu(x):
    return 0.5*x*(1+np.tanh(
        np.sqrt(2/np.pi)*(x+0.044715*x**3)))

x = np.array([-2, -1, 0, 1, 2])
print("ReLU:  ", relu(x))
# [0, 0, 0, 1, 2]

def relu(x):
    return max(0, x)

def relu_deriv(x):
    return 1 if x > 0 else 0

def leaky_relu(x, alpha=0.01):
    return x if x > 0 else alpha * x

def leaky_relu_deriv(x, alpha=0.01):
    return 1 if x > 0 else alpha

# The "dying ReLU" problem:
# If a neuron always receives
# negative input, its gradient
# is always 0 → it never learns!
# Leaky ReLU fixes this by
# allowing a small gradient (α)
# for negative inputs.

import numpy as np

def softmax(z):
    # Numerically stable
    e = np.exp(z - np.max(z))
    return e / e.sum()

def cross_entropy(y_true, y_pred):
    return -np.sum(y_true * np.log(y_pred + 1e-15))

# Example: 3-class classification
logits = np.array([2.0, 1.0, 0.1])
probs = softmax(logits)
print(probs)
# [0.659, 0.242, 0.099]
print(probs.sum())  # 1.0

# Gradient: simply ŷ − y!
y_true = np.array([1, 0, 0])  # one-hot
grad = probs - y_true
print(grad)
# [-0.341, 0.242, 0.099]

import math

def softmax(z):
    # Subtract max for stability
    m = max(z)
    exps = [math.exp(zi - m) for zi in z]
    total = sum(exps)
    return [e / total for e in exps]

def cross_entropy(y_true, y_pred):
    return -sum(
        yt * math.log(yp + 1e-15)
        for yt, yp in zip(y_true, y_pred)
    )

logits = [2.0, 1.0, 0.1]
probs = softmax(logits)
print([round(p, 3) for p in probs])
# [0.659, 0.242, 0.099]

import numpy as np

def mse(y_true, y_pred):
    return np.mean((y_pred - y_true)**2)

def mse_grad(y_true, y_pred):
    n = len(y_true)
    return 2 * (y_pred - y_true) / n

def binary_ce(y_true, y_pred):
    eps = 1e-15
    y_pred = np.clip(y_pred, eps, 1-eps)
    return -np.mean(
        y_true*np.log(y_pred) +
        (1-y_true)*np.log(1-y_pred)
    )

y = np.array([1, 0, 1])
yh = np.array([0.9, 0.1, 0.8])
print(f"MSE:  {mse(y, yh):.4f}")
print(f"BCE:  {binary_ce(y, yh):.4f}")

import math

def mse(y_true, y_pred):
    n = len(y_true)
    return sum((yp-yt)**2
        for yt,yp in zip(y_true,y_pred)) / n

def mse_grad(y_true, y_pred):
    n = len(y_true)
    return [2*(yp-yt)/n
        for yt,yp in zip(y_true,y_pred)]

def binary_ce(y_true, y_pred):
    eps = 1e-15
    n = len(y_true)
    total = 0
    for yt, yp in zip(y_true, y_pred):
        yp = max(eps, min(1-eps, yp))
        total -= (yt*math.log(yp) +
                  (1-yt)*math.log(1-yp))
    return total / n

import numpy as np
from scipy import integrate

# ∫₀² x² dx = [x³/3]₀² = 8/3 ≈ 2.667
result, error = integrate.quad(
    lambda x: x**2, 0, 2
)
print(f"{result:.4f}")  # 2.6667

# Numerical integration with trapezoid
x = np.linspace(0, 2, 1000)
y = x**2
area = np.trapz(y, x)
print(f"{area:.4f}")    # 2.6667

def riemann_sum(f, a, b, n=1000):
    """Left Riemann sum approximation"""
    dx = (b - a) / n
    total = 0
    for i in range(n):
        x = a + i * dx
        total += f(x) * dx
    return total

def trapezoidal(f, a, b, n=1000):
    """More accurate than Riemann"""
    dx = (b - a) / n
    total = (f(a) + f(b)) / 2
    for i in range(1, n):
        total += f(a + i * dx)
    return total * dx

# ∫₀² x² dx = 8/3
result = trapezoidal(
    lambda x: x**2, 0, 2
)
print(f"{result:.6f}")
# 2.666667

import numpy as np

# Verify FTC: ∫₀³ 2x dx = [x²]₀³ = 9
# Antiderivative F(x) = x²
F = lambda x: x**2
exact = F(3) - F(0)
print(exact)  # 9

# Numerical check
x = np.linspace(0, 3, 10000)
numerical = np.trapz(2*x, x)
print(f"{numerical:.4f}")  # 9.0000

def integrate_poly(coeffs, a, b):
    """Exact integral of polynomial
    coeffs: [c0, c1, c2, ...]
    p(x) = c0 + c1*x + c2*x² + ...
    """
    # Antiderivative coeffs
    F_coeffs = [0]  # constant term
    for i, c in enumerate(coeffs):
        F_coeffs.append(c / (i + 1))

    def F(x):
        return sum(c * x**(i)
            for i, c in enumerate(F_coeffs))

    return F(b) - F(a)

# ∫₀³ 2x dx (coeffs = [0, 2])
print(integrate_poly([0, 2], 0, 3))
# 9.0

import numpy as np

def simpsons(f, a, b, n=1000):
    """Simpson's 1/3 rule (most accurate)"""
    assert n % 2 == 0
    x = np.linspace(a, b, n+1)
    y = f(x)
    h = (b - a) / n
    return h/3 * (y[0] + y[-1]
        + 4*np.sum(y[1:-1:2])
        + 2*np.sum(y[2:-2:2]))

# ∫₀π sin(x) dx = 2
result = simpsons(np.sin, 0, np.pi)
print(f"{result:.10f}")
# 2.0000000000

import random, math

def monte_carlo_integrate(f, a, b, n=100000):
    """Approximate ∫ₐᵇ f(x) dx
    by sampling random points"""
    total = 0
    for _ in range(n):
        x = random.uniform(a, b)
        total += f(x)
    return (b - a) * total / n

# ∫₀π sin(x) dx = 2
result = monte_carlo_integrate(
    math.sin, 0, math.pi
)
print(f"{result:.4f}")
# ≈ 2.0000 (stochastic!)

import numpy as np

def kl_divergence(p, q):
    """KL(P||Q) for discrete dists"""
    p = np.asarray(p, dtype=float)
    q = np.asarray(q, dtype=float)
    # Avoid log(0)
    mask = p > 0
    return np.sum(p[mask] * np.log(p[mask] / q[mask]))

p = [0.4, 0.3, 0.3]
q = [0.33, 0.33, 0.34]
print(f"KL(P||Q) = {kl_divergence(p, q):.4f}")
# KL(P||Q) = 0.0133

import math

def expected_value(values, probs):
    """E[X] = Σ xᵢ · p(xᵢ)"""
    return sum(x*p for x,p
               in zip(values, probs))

# Fair die: E[X] = 3.5
vals = [1,2,3,4,5,6]
probs = [1/6]*6
print(expected_value(vals, probs))
# 3.5

# Variance: E[X²] − E[X]²
ex2 = expected_value([x**2 for x in vals], probs)
ex = expected_value(vals, probs)
print(f"Var = {ex2 - ex**2:.4f}")
# Var = 2.9167

import numpy as np

def f(x, y):
    return 3*x**2*y + 2*x*y**3 - 5

def partial_x(f, x, y, h=1e-7):
    return (f(x+h, y) - f(x-h, y)) / (2*h)

def partial_y(f, x, y, h=1e-7):
    return (f(x, y+h) - f(x, y-h)) / (2*h)

print(partial_x(f, 1, 2))  # 28.0 (6·1·2 + 2·8)
print(partial_y(f, 1, 2))  # 27.0 (3·1 + 6·1·4)

def partial_derivative(f, args, idx, h=1e-7):
    """Partial deriv w.r.t. args[idx]"""
    args_p = list(args)
    args_m = list(args)
    args_p[idx] += h
    args_m[idx] -= h
    return (f(*args_p) - f(*args_m)) / (2*h)

def f(x, y):
    return 3*x**2*y + 2*x*y**3 - 5

# ∂f/∂x at (1,2)
print(partial_derivative(f, [1,2], 0))
# 28.0
# ∂f/∂y at (1,2)
print(partial_derivative(f, [1,2], 1))
# 27.0

import numpy as np

def gradient(f, params, h=1e-7):
    """Compute ∇f numerically"""
    grad = np.zeros_like(params, dtype=float)
    for i in range(len(params)):
        p_plus = params.copy(); p_plus[i] += h
        p_minus = params.copy(); p_minus[i] -= h
        grad[i] = (f(p_plus) - f(p_minus)) / (2*h)
    return grad

# f(w) = w₀² + 3w₁² (bowl)
f = lambda w: w[0]**2 + 3*w[1]**2
w = np.array([4.0, 2.0])
print(gradient(f, w))  # [8.0, 12.0]
# Exact: [2w₀, 6w₁] = [8, 12] ✓

def gradient(f, params, h=1e-7):
    """Compute ∇f for list of params"""
    grad = []
    for i in range(len(params)):
        p_plus = params[:]
        p_minus = params[:]
        p_plus[i] += h
        p_minus[i] -= h
        g = (f(p_plus) - f(p_minus)) / (2*h)
        grad.append(g)
    return grad

def f(w):
    return w[0]**2 + 3*w[1]**2

print(gradient(f, [4.0, 2.0]))
# [8.0, 12.0]

import numpy as np

def gradient_descent(f, grad_f, w0,
                     lr=0.01, epochs=100):
    w = w0.copy()
    history = [f(w)]
    for _ in range(epochs):
        g = grad_f(w)
        w = w - lr * g
        history.append(f(w))
    return w, history

# Minimize f(w) = w₀² + 3w₁²
f = lambda w: w[0]**2 + 3*w[1]**2
grad = lambda w: np.array([2*w[0], 6*w[1]])

w_opt, hist = gradient_descent(
    f, grad, np.array([4.0, 2.0]),
    lr=0.1, epochs=50)
print(f"w = {w_opt}")
print(f"f(w) = {f(w_opt):.6f}")
# w ≈ [0.0, 0.0]  f(w) ≈ 0.0

def gradient_descent(f, grad_f, w0,
                     lr=0.01, epochs=100):
    w = w0[:]
    for _ in range(epochs):
        g = grad_f(w)
        w = [wi - lr*gi
             for wi, gi in zip(w, g)]
    return w

def f(w): return w[0]**2 + 3*w[1]**2
def grad_f(w): return [2*w[0], 6*w[1]]

w = gradient_descent(
    f, grad_f, [4.0, 2.0],
    lr=0.1, epochs=50)
print([round(x,6) for x in w])
# [0.0, 0.0]

import numpy as np

def sgd_momentum(grad_f, w0, lr=0.01,
                 beta=0.9, epochs=100):
    w = w0.copy()
    v = np.zeros_like(w)  # velocity
    for _ in range(epochs):
        g = grad_f(w)
        v = beta * v + g       # accumulate
        w = w - lr * v         # step
    return w

grad = lambda w: np.array([2*w[0], 6*w[1]])
w = sgd_momentum(grad, np.array([4.0, 2.0]),
                 lr=0.01, beta=0.9)
print(w)  # ≈ [0, 0]

def sgd_momentum(grad_f, w0, lr=0.01,
                 beta=0.9, epochs=100):
    w = w0[:]
    v = [0.0] * len(w)
    for _ in range(epochs):
        g = grad_f(w)
        v = [beta*vi + gi
             for vi,gi in zip(v,g)]
        w = [wi - lr*vi
             for wi,vi in zip(w,v)]
    return w

def grad_f(w):
    return [2*w[0], 6*w[1]]

w = sgd_momentum(grad_f, [4.0, 2.0])
print([round(x,4) for x in w])
# [0.0, 0.0]

import numpy as np

def adam(grad_f, w0, lr=0.001,
        b1=0.9, b2=0.999, eps=1e-8,
        epochs=200):
    w = w0.copy()
    m = np.zeros_like(w)
    v = np.zeros_like(w)
    for t in range(1, epochs+1):
        g = grad_f(w)
        m = b1*m + (1-b1)*g
        v = b2*v + (1-b2)*g**2
        m_hat = m / (1 - b1**t)
        v_hat = v / (1 - b2**t)
        w = w - lr * m_hat / (np.sqrt(v_hat) + eps)
    return w

grad = lambda w: np.array([2*w[0], 6*w[1]])
w = adam(grad, np.array([4.0, 2.0]))
print(np.round(w, 4))  # [0, 0]

import math

def adam(grad_f, w0, lr=0.001,
        b1=0.9, b2=0.999, eps=1e-8,
        epochs=200):
    w = w0[:]
    m = [0.0]*len(w)
    v = [0.0]*len(w)
    for t in range(1, epochs+1):
        g = grad_f(w)
        for i in range(len(w)):
            m[i] = b1*m[i]+(1-b1)*g[i]
            v[i] = b2*v[i]+(1-b2)*g[i]**2
            mh = m[i]/(1-b1**t)
            vh = v[i]/(1-b2**t)
            w[i] -= lr*mh/(math.sqrt(vh)+eps)
    return w

def grad_f(w):
    return [2*w[0], 6*w[1]]

w = adam(grad_f, [4.0, 2.0])
print([round(x,4) for x in w])
# [0.0, 0.0]

import numpy as np

def jacobian(f, x, h=1e-7):
    """Numerical Jacobian matrix"""
    x = np.asarray(x, dtype=float)
    f0 = np.asarray(f(x))
    n = len(x); m = len(f0)
    J = np.zeros((m, n))
    for j in range(n):
        xp = x.copy(); xp[j] += h
        xm = x.copy(); xm[j] -= h
        J[:, j] = (f(xp) - f(xm)) / (2*h)
    return J

# f(x,y) = [x²+y, xy]
f = lambda x: np.array([x[0]**2+x[1], x[0]*x[1]])
J = jacobian(f, [3.0, 2.0])
print(J)
# [[6. 1.]    ← [2x, 1]
#  [2. 3.]]  ← [y, x]

def jacobian(f, x, h=1e-7):
    n = len(x)
    f0 = f(x)
    m = len(f0)
    J = [[0.0]*n for _ in range(m)]
    for j in range(n):
        xp = x[:]; xm = x[:]
        xp[j] += h; xm[j] -= h
        fp = f(xp); fm = f(xm)
        for i in range(m):
            J[i][j] = (fp[i]-fm[i])/(2*h)
    return J

def f(x):
    return [x[0]**2+x[1], x[0]*x[1]]

J = jacobian(f, [3.0, 2.0])
for row in J: print(row)
# [6.0, 1.0]
# [2.0, 3.0]

import numpy as np

def hessian(f, x, h=1e-5):
    """Numerical Hessian matrix"""
    n = len(x)
    H = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            xpp=x.copy(); xpp[i]+=h; xpp[j]+=h
            xpm=x.copy(); xpm[i]+=h; xpm[j]-=h
            xmp=x.copy(); xmp[i]-=h; xmp[j]+=h
            xmm=x.copy(); xmm[i]-=h; xmm[j]-=h
            H[i,j] = (f(xpp)-f(xpm)-f(xmp)+f(xmm))
            H[i,j] /= (4*h*h)
    return H

# f(x,y) = x² + 3y² → H = [[2,0],[0,6]]
f = lambda x: x[0]**2 + 3*x[1]**2
print(hessian(f, np.array([1.0,1.0])))
# [[2. 0.]
#  [0. 6.]]

# Hessian eigenvalues tell us
# about the loss landscape shape:
#
# All positive → local minimum
# All negative → local maximum
# Mixed → saddle point
#
# Condition number (λ_max / λ_min)
# tells us how "elongated" the
# bowl is. High condition number
# → GD oscillates → slow.
# This is why Adam helps: it
# scales each dimension by 1/√v,
# effectively preconditioning.

import numpy as np
H = np.array([[2, 0], [0, 6]])
eigvals = np.linalg.eigvalsh(H)
print(f"Eigenvalues: {eigvals}")
# [2. 6.]
print(f"Condition: {max(eigvals)/min(eigvals)}")
# 3.0

Each node computes its output from its inputs and caches intermediate values (needed for backward pass).

Each node receives ∂L/∂output, computes ∂L/∂input using its local derivative, and passes it upstream.

import numpy as np

def taylor_exp(x, n_terms=10):
    """eˣ ≈ 1 + x + x²/2! + x³/3! +..."""
    result = np.zeros_like(x, dtype=float)
    for n in range(n_terms):
        result += x**n / np.math.factorial(n)
    return result

x = np.array([0, 1, 2])
print("Taylor:", taylor_exp(x))
print("Exact: ", np.exp(x))
# Taylor: [1.    2.718 7.389]
# Exact:  [1.    2.718 7.389]

def factorial(n):
    r = 1
    for i in range(2, n+1): r *= i
    return r

def taylor_exp(x, n_terms=15):
    return sum(x**n / factorial(n)
               for n in range(n_terms))

def taylor_sin(x, n_terms=10):
    return sum(
        (-1)**n * x**(2*n+1) / factorial(2*n+1)
        for n in range(n_terms))

print(taylor_exp(1))    # 2.71828...
print(taylor_sin(3.14159/2))
# ≈ 1.0

Propagates derivatives alongside values. Cost = O(n) passes for n inputs. Good when inputs << outputs (rare in ML).

One forward pass + one backward pass computes ALL gradients. Cost = O(1) backward passes regardless of # inputs. This is why backprop is used in ML.