Multivariable Calculus & Gradients

In single-variable calculus, the derivative of f(x) at a point x₀ is defined as the limit:

$$f'(x_0) = \lim_{h \to 0} \frac{f(x_0 + h) - f(x_0)}{h}$$

The key idea: perturb the one input by a tiny amount h and measure how the output changes, normalised by h.

For a function of two variables, say f(x, y), there is no single "direction" along which to perturb anymore — perturbations can go in any direction in the xy-plane. The simplest thing we can do is hold one variable fixed and differentiate with respect to the other. This gives us a partial derivative.

Definition. The partial derivative of f with respect to x at (x₀, y₀) is:

$$\frac{\partial f}{\partial x}(x_0, y_0) = \lim_{h \to 0} \frac{f(x_0 + h, y_0) - f(x_0, y_0)}{h}$$

Notation: ∂f/∂x, fₓ, or D₁f. The ∂ symbol (pronounced "del" or "partial") signals that other variables are held constant.

Geometric picture. The graph of f(x, y) is a surface in 3D. If you slice that surface with the vertical plane y = y₀, you get a curve. The partial derivative ∂f/∂x at (x₀, y₀) is simply the slope of that curve at x₀. Similarly, slice with x = x₀ to get the curve whose slope is ∂f/∂y.

Example. Let f(x, y) = x²y + sin(y).

• ∂f/∂x = 2xy (treat y as a constant; sin(y) is a constant, so its x-derivative is 0)
• ∂f/∂y = x² + cos(y) (treat x as a constant; d/dy of x²y is x²; d/dy of sin(y) is cos(y))

Higher-order partials. We can differentiate again:

• ∂²f/∂x² = 2y (differentiate ∂f/∂x = 2xy with respect to x)
• ∂²f/∂x∂y = 2x (differentiate ∂f/∂x = 2xy with respect to y)
• ∂²f/∂y∂x = 2x (differentiate ∂f/∂y = x² + cos(y) with respect to x)

Clairaut's theorem (Schwarz's theorem): if the mixed partial derivatives are continuous, the order of differentiation doesn't matter: ∂²f/∂x∂y = ∂²f/∂y∂x. This is almost always true in practice for the functions we encounter.

import sympy as sp

x, y = sp.symbols('x y')
f = x**2 * y + sp.sin(y)

df_dx = sp.diff(f, x)      # 2*x*y
df_dy = sp.diff(f, y)      # x**2 + cos(y)
d2f_dxdy = sp.diff(f, x, y)  # 2*x  (Clairaut: same as diff(f,y,x))

print(df_dx, df_dy, d2f_dxdy)

Partial derivatives tell us about change along the coordinate axes. But what about change in an arbitrary direction? And is there a single object that encodes all the first-order directional information about f at a point?

Yes. That object is the gradient.

Definition. For f: ℝⁿ → ℝ, the gradient at a point x ∈ ℝⁿ is the vector of all partial derivatives:

$$\nabla f(\mathbf{x}) = \begin{pmatrix} \partial f/\partial x_1 \ \partial f/\partial x_2 \ \vdots \ \partial f/\partial x_n \end{pmatrix}$$

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

The gradient is a vector field. At each point (x, y) in the domain, ∇f(x, y) gives a vector. Plot all these vectors together and you get a picture of how f changes throughout its domain — a map showing the slope-landscape.

The key geometric facts:

1. ∇f(x) points in the direction of steepest ascent of f at x. Move in the direction of ∇f and f increases as fast as possible.
2. ∇f(x) is perpendicular to the level set of f through x. The level set is {y : f(y) = f(x)}. Strikingly, the direction of maximum change is orthogonal to the surface of no change. (This is why level curves on a topographic map are perpendicular to the paths of steepest ascent.)
3. |∇f(x)| is the rate of steepest ascent.

Example. f(x, y) = x² + y². The level sets are circles. The gradient is ∇f = (2x, 2y) = 2(x, y), which is the radial outward vector — perpendicular to the circles. Moving away from the origin is steepest ascent. Moving along a circle changes f not at all.

Analogy. Imagine f is the elevation of a hilly landscape. The gradient at your location points in the compass direction you should walk to climb as steeply as possible. The contour lines on a topo map (level sets) are exactly perpendicular to this direction.

import numpy as np
import matplotlib.pyplot as plt

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

# Numerical gradient at a grid of points
X, Y = np.meshgrid(np.linspace(-2, 2, 20), np.linspace(-2, 2, 20))
GX = 2 * X   # analytical: partial of (x^2+y^2) w.r.t. x
GY = 2 * Y   # analytical: partial of (x^2+y^2) w.r.t. y

plt.contour(X, Y, f(X, Y), levels=10, cmap='Blues')
plt.quiver(X, Y, GX, GY, alpha=0.7)
plt.title("Gradient field of f(x,y)=x²+y² — arrows perpendicular to contours")
plt.axis('equal'); plt.show()

The partial derivatives ∂f/∂x and ∂f/∂y measure change in the x and y directions respectively. But what is the rate of change in an arbitrary direction u?

Definition. The directional derivative of f at x in the direction of a unit vector u (|u| = 1) is:

$$D_\mathbf{u} f(\mathbf{x}) = \lim_{h \to 0} \frac{f(\mathbf{x} + h\mathbf{u}) - f(\mathbf{x})}{h}$$

This is just the ordinary derivative of the single-variable function g(t) = f(x + tu) evaluated at t = 0.

The theorem that connects everything:

$$D_\mathbf{u} f(\mathbf{x}) = \nabla f(\mathbf{x}) \cdot \mathbf{u} = |\nabla f(\mathbf{x})| \cos\theta$$

where θ is the angle between ∇f and u. This follows from differentiating g(t) = f(x + tu) using the chain rule: g'(0) = ∇f(x)·u.

Why the gradient points uphill. The expression |∇f|cos(θ) is maximised when θ = 0 (i.e., u is in the same direction as ∇f), giving a maximum rate of change of |∇f|. It is minimised (steepest descent) when θ = π. It is zero when θ = π/2, i.e., u is perpendicular to the gradient — which is exactly the direction along the level set.

Example. f(x, y) = 3x + y², so ∇f = (3, 2y). At the point (1, 2), ∇f = (3, 4). The directional derivative in the direction u = (1/√2, 1/√2) is: D_u f = (3)(1/√2) + (4)(1/√2) = 7/√2 ≈ 4.95.

Important edge case. When f is not differentiable at a point (e.g., f(x, y) = |x| + |y| at the origin), directional derivatives in all directions may still exist, but they won't be given by the gradient formula — the gradient doesn't exist there.

import numpy as np

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

# Directional derivative numerically: finite difference along u
def directional_deriv(f, x0, y0, u, h=1e-7):
    u = np.array(u, dtype=float)
    u = u / np.linalg.norm(u)  # ensure unit vector
    return (f(x0 + h*u[0], y0 + h*u[1]) - f(x0, y0)) / h

# Analytical gradient at (1, 2): (3, 4)
grad = np.array([3.0, 4.0])
u = np.array([1.0, 1.0]) / np.sqrt(2)

analytical = np.dot(grad, u)               # 7/sqrt(2) ≈ 4.950
numerical  = directional_deriv(f, 1, 2, u) # same, up to float error

print(f"Analytical: {analytical:.6f}, Numerical: {numerical:.6f}")

Partial derivatives measure change along coordinate axes. But what if you perturb all inputs simultaneously — perhaps you move in an arbitrary direction in ℝⁿ? The right object is the total derivative (also called the differential or Fréchet derivative).

Linear approximation. For a function f: ℝⁿ → ℝ that is differentiable at x₀, the best linear approximation near x₀ is:

$$f(\mathbf{x}_0 + \boldsymbol{\delta}) \approx f(\mathbf{x}_0) + \nabla f(\mathbf{x}_0)^\top \boldsymbol{\delta}$$

for small perturbations δ. This is the multivariable version of f(x₀+h) ≈ f(x₀) + f'(x₀)h.

The gradient ∇f(x₀)ᵀ acts as a linear map from ℝⁿ to ℝ: it takes a perturbation vector and returns a scalar (the first-order change in f). In linear algebra terms, this linear map is a row vector (a 1×n matrix).

What differentiability really means. f is differentiable at x₀ if there exists a linear map L (i.e., a row vector) such that:

$$\lim_{\boldsymbol{\delta} \to 0} \frac{|f(\mathbf{x}_0 + \boldsymbol{\delta}) - f(\mathbf{x}_0) - L \boldsymbol{\delta}|}{|\boldsymbol{\delta}|} = 0$$

When this holds, L = ∇f(x₀)ᵀ. The condition says: the linear approximation is not just close to f near x₀, it is close relative to the size of the perturbation, meaning the error is o(|δ|).

Caution: partial derivatives existing does not imply differentiability. The classic counterexample is f(x, y) = xy/(x² + y²) for (x,y) ≠ (0,0) and f(0,0) = 0. Both partial derivatives at (0,0) are 0, yet f is not even continuous there (its value along y = x is 1/2). Differentiability requires partial derivatives to exist and to correctly predict change in all directions.

Chain rule (the most important theorem in machine learning). If g: ℝⁿ → ℝᵐ and f: ℝᵐ → ℝ, and we compose h = f ∘ g, then:

$$\nabla h(\mathbf{x}) = J_g(\mathbf{x})^\top \nabla f(g(\mathbf{x}))$$

where Jg is the Jacobian matrix of g (defined in the next concept). This is backpropagation.

# Linear approximation demonstration in Python
import numpy as np

def f(x, y): return np.sin(x) * np.exp(y)

# At (pi/4, 0): f = sin(pi/4)*exp(0) = 1/sqrt(2) ≈ 0.7071
# Gradient: df/dx = cos(x)*exp(y), df/dy = sin(x)*exp(y)
# At (pi/4, 0): grad = (cos(pi/4), sin(pi/4)) = (1/sqrt(2), 1/sqrt(2))

x0, y0 = np.pi/4, 0.0
grad = np.array([np.cos(x0)*np.exp(y0), np.sin(x0)*np.exp(y0)])

delta = np.array([0.01, 0.02])  # small perturbation

f_approx = f(x0, y0) + grad @ delta       # linear approximation
f_exact  = f(x0 + delta[0], y0 + delta[1]) # true value
error    = abs(f_exact - f_approx)

print(f"Approx: {f_approx:.8f}")
print(f"Exact:  {f_exact:.8f}")
print(f"Error:  {error:.2e}  (should be O(delta^2) ≈ {np.linalg.norm(delta)**2:.2e})")

So far we have considered f: ℝⁿ → ℝ (scalar-valued). But many functions map ℝⁿ → ℝᵐ — they have multiple outputs. The Jacobian is the matrix of all partial derivatives that generalises the gradient to this vector-valued case.

Definition. For f: ℝⁿ → ℝᵐ where f = (f₁, …, fₘ)ᵀ, the Jacobian at x is the m×n matrix:

$$J_f(\mathbf{x}) = \frac{\partial \mathbf{f}}{\partial \mathbf{x}} = \begin{pmatrix} \partial f_1/\partial x_1 & \cdots & \partial f_1/\partial x_n \ \vdots & \ddots & \vdots \ \partial f_m/\partial x_1 & \cdots & \partial f_m/\partial x_n \end{pmatrix}$$

Row i is the gradient (transposed) of fᵢ. Column j tells you how all outputs respond to a unit change in input xⱼ.

Geometric meaning. The Jacobian is the best linear approximation of f near x:

$$\mathbf{f}(\mathbf{x} + \boldsymbol{\delta}) \approx \mathbf{f}(\mathbf{x}) + J_f(\mathbf{x}) \boldsymbol{\delta}$$

This is the multivariable generalization of the tangent line. J_f(x) maps tangent vectors at x (in ℝⁿ) to tangent vectors at f(x) (in ℝᵐ). This is the beginning of differential geometry.

Special cases:

• n=1, m=1: J_f = f'(x), the ordinary derivative
• n=n, m=1: J_f = ∇f(x)ᵀ, the gradient as a row vector
• n=m: |det(J_f)| is the local volume scaling factor (how the map stretches/squishes infinitesimal volumes)

The multivariable chain rule in matrix form. For h = f ∘ g:

$$J_h(\mathbf{x}) = J_f(g(\mathbf{x})) \cdot J_g(\mathbf{x})$$

This is matrix multiplication. In neural networks, each layer performs a linear map (whose Jacobian is the weight matrix) followed by an elementwise nonlinearity (whose Jacobian is a diagonal matrix of derivative values). The product of all these Jacobians — computed right-to-left in the backward pass — is backpropagation.

Example. f(x, y) = (x² + y, sin(x) + y²).

$$J_f = \begin{pmatrix} 2x & 1 \ \cos x & 2y \end{pmatrix}$$

import numpy as np

def f(v):
    x, y = v
    return np.array([x**2 + y, np.sin(x) + y**2])

def jacobian_numerical(f, v, h=1e-6):
    n = len(v)
    fv = f(v)
    m = len(fv)
    J = np.zeros((m, n))
    for j in range(n):
        vp = v.copy(); vp[j] += h
        J[:, j] = (f(vp) - fv) / h
    return J

v0 = np.array([1.0, 2.0])
J_num = jacobian_numerical(f, v0)
# Analytical: [[2x, 1], [cos(x), 2y]] at (1,2) = [[2, 1], [cos(1), 4]]
J_ana = np.array([[2*1, 1], [np.cos(1), 2*2]])
print("Numerical Jacobian:\n", J_num)
print("Analytical Jacobian:\n", J_ana)

Just as the second derivative f''(x) tells us about the curvature of a single-variable function (is the graph concave up or concave down?), the Hessian captures the curvature of a multivariable function f: ℝⁿ → ℝ.

Definition. The Hessian of f at x is the n×n matrix of all second partial derivatives:

$$H_f(\mathbf{x}) = \nabla^2 f(\mathbf{x}) = \begin{pmatrix} \partial^2 f/\partial x_1^2 & \partial^2 f/\partial x_1 \partial x_2 & \cdots \ \partial^2 f/\partial x_2 \partial x_1 & \partial^2 f/\partial x_2^2 & \cdots \ \vdots & & \ddots \end{pmatrix}$$

By Clairaut's theorem, the Hessian is symmetric: (Hf)ᵢⱼ = (Hf)ⱼᵢ whenever the mixed partials are continuous.

Second-order Taylor expansion. The Hessian allows a quadratic approximation:

$$f(\mathbf{x} + \boldsymbol{\delta}) \approx f(\mathbf{x}) + \nabla f(\mathbf{x})^\top \boldsymbol{\delta} + \frac{1}{2} \boldsymbol{\delta}^\top H_f(\mathbf{x}) \boldsymbol{\delta}$$

This is the multivariable analogue of f(x+h) ≈ f(x) + f'(x)h + ½f''(x)h².

Definiteness and critical points. At a critical point where ∇f = 0, the Hessian determines the nature of the point:

• H is positive definite (all eigenvalues > 0): local minimum
• H is negative definite (all eigenvalues < 0): local maximum
• H is indefinite (mixed-sign eigenvalues): saddle point
• H is positive/negative semi-definite (some eigenvalues = 0): inconclusive (higher-order terms needed)

Why eigenvalues? The quadratic form δᵀ H δ measures the curvature of f in direction δ. Its sign is determined by the eigenvalues of H. If all eigenvalues are positive, the curvature is positive in every direction — like the bottom of a bowl — so the point is a local minimum.

Condition number and optimization. The ratio λ_max/λ_min of the Hessian's eigenvalues is the condition number. When this is large (ill-conditioned), gradient descent oscillates wildly: the gradient in the steep direction overcorrects while barely moving in the shallow direction. This is why training neural networks is hard — the loss landscape is highly ill-conditioned.

Newton's method uses the Hessian to take a smarter step: rather than moving against the gradient by a fixed step size, it solves for the minimum of the local quadratic approximation:

$$\mathbf{x}_{k+1} = \mathbf{x}_k - H_f(\mathbf{x}_k)^{-1} \nabla f(\mathbf{x}_k)$$

This converges quadratically (number of correct digits doubles each step) near a minimum, compared to linear convergence for gradient descent.

import numpy as np

def f(v):
    x, y = v
    return x**4 - 2*x**2 + y**2 + x*y   # has interesting critical points

def gradient(v):
    x, y = v
    return np.array([4*x**3 - 4*x + y, 2*y + x])

def hessian_numerical(f, v, h=1e-4):
    n = len(v)
    H = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            vpp = v.copy(); vpp[i] += h; vpp[j] += h
            vpm = v.copy(); vpm[i] += h; vpm[j] -= h
            vmp = v.copy(); vmp[i] -= h; vmp[j] += h
            vmm = v.copy(); vmm[i] -= h; vmm[j] -= h
            H[i,j] = (f(vpp)-f(vpm)-f(vmp)+f(vmm)) / (4*h*h)
    return H

v0 = np.array([1.0, -0.5])  # near a local minimum
H = hessian_numerical(f, v0)
eigenvalues = np.linalg.eigvalsh(H)  # symmetric → real eigenvalues
print("Hessian:\n", H)
print("Eigenvalues:", eigenvalues)   # both positive → local min

Now we arrive at the payoff. Given f: ℝⁿ → ℝ (think: the loss function of a neural network), we want to find a point x* where f is small. Gradient descent exploits the geometric meaning of the gradient to do this iteratively.

The algorithm.

1. Start at some initial point x₀.
2. Compute the gradient ∇f(xₖ).
3. Move opposite to the gradient (since ∇f points uphill, −∇f points downhill): xₖ₊₁ = xₖ − α ∇f(xₖ)
4. Repeat until ‖∇f(xₖ)‖ is small.

The scalar α > 0 is the learning rate (or step size). It controls how far we move each step.

Why does this work? From the directional derivative formula, moving in direction u changes f by approximately |∇f| cos(θ) per unit distance. The direction u = −∇f/|∇f| (i.e., opposite the gradient) maximises the rate of decrease. So each step decreases f as fast as possible locally.

The learning rate dilemma. If α is too small, convergence is slow. If α is too large, we overshoot the minimum and may diverge. The right α depends on the curvature (Hessian). A safe choice is α < 2/λ_max, where λ_max is the largest eigenvalue of the Hessian. This is why understanding the Hessian matters even for first-order methods.

Convergence for convex functions. For a convex function with L-smooth gradient (‖∇f(x) − ∇f(y)‖ ≤ L‖x − y‖), gradient descent with α = 1/L converges: f(xₖ) − f(x*) ≤ O(1/k).

What does 'convex' mean geometrically? A function is convex if its graph lies above every tangent hyperplane. Formally, f(λx + (1-λ)y) ≤ λf(x) + (1-λ)f(y) for all x, y and λ ∈ [0,1]. For smooth functions, convexity is equivalent to the Hessian being positive semidefinite everywhere.

Stochastic gradient descent (SGD). In ML, f is a sum over training examples: f(θ) = (1/N) Σᵢ ℓᵢ(θ). Computing the full gradient costs O(N). SGD estimates the gradient using a single random example (or a minibatch), drastically reducing per-step cost. The gradient estimate is noisy, but this actually helps escape sharp local minima.

Beyond SGD: Adam. Modern optimizers like Adam maintain running estimates of the gradient and its elementwise square root (second moment), effectively adapting the learning rate per parameter. The gradient is still the fundamental object; the adaptations just make it easier to follow.

import numpy as np

def f(v):    return v[0]**2 + 5*v[1]**2  # ill-conditioned bowl
def grad(v): return np.array([2*v[0], 10*v[1]])

def gradient_descent(f, grad, x0, lr=0.1, steps=50):
    x = x0.copy()
    history = [x.copy()]
    for _ in range(steps):
        g = grad(x)
        x = x - lr * g
        history.append(x.copy())
        if np.linalg.norm(g) < 1e-8:
            break
    return np.array(history)

x0 = np.array([3.0, 2.0])
history = gradient_descent(f, grad, x0, lr=0.09)
print(f"Start:  f={f(x0):.4f}  at {x0}")
print(f"Finish: f={f(history[-1]):.6f} at {history[-1]}")
# Note: with lr=0.2 (>2/lambda_max=0.2), it diverges on x[1]!
# This shows why learning rate < 2/lambda_max matters.

Computing gradients by hand for a neural network with millions of parameters is impossible. But computing them exactly (not approximately like finite differences) is what makes deep learning tractable. Automatic differentiation (autodiff) does this by applying the chain rule to a computational graph.

The idea. Every function computation can be broken into elementary operations (add, multiply, sin, exp, …). Record these in a graph where nodes are values and edges are operations. Then apply the chain rule systematically through this graph.

Two modes:

1. Forward mode. Propagate derivatives forward alongside values. For each input xᵢ, compute ∂y/∂xᵢ for all outputs y in one pass. Cost: n forward passes for n inputs. Efficient when n << m.

2. Reverse mode (backpropagation). First do a forward pass computing all values. Then propagate gradients backward: for each output yⱼ, compute ∂yⱼ/∂xᵢ for all inputs in one backward pass. Cost: one forward + one backward for each output. Efficient when m << n.

Why reverse mode dominates ML. Neural networks have n = millions (parameters) and m = 1 (scalar loss). One backward pass gives ∇L(θ) — the gradient of the loss with respect to all parameters — at the cost of ≈2× the forward pass. This is the reason deep learning is computationally feasible.

PyTorch example. PyTorch builds a dynamic computational graph during the forward pass (each operation records its inputs and the rule for differentiating through it). Calling .backward() traverses this graph in reverse, accumulating gradients in .grad attributes.

Mathematical foundation. Each node computes a function g(u, v) → z. During the backward pass, if we know ∂L/∂z (the incoming gradient), we propagate:

• ∂L/∂u += (∂L/∂z)(∂z/∂u) [chain rule]
• ∂L/∂v += (∂L/∂z)(∂z/∂v)

The local Jacobians ∂z/∂u and ∂z/∂v are cheap to compute (they're just elementwise derivatives for operations like ReLU). The backpropagation algorithm is nothing more than the chain rule, applied left-to-right through the graph from outputs to inputs.

import torch

# Build a tiny neural net forward pass and let PyTorch differentiate it
x = torch.tensor([1.0, 2.0, 3.0], requires_grad=False)
W = torch.tensor([[0.1, 0.2, 0.3],
                   [0.4, 0.5, 0.6]], requires_grad=True)   # 2x3 weight matrix
b = torch.tensor([0.1, -0.1], requires_grad=True)          # bias

# Forward pass: z = W @ x + b, loss = sum(z^2)
z = W @ x + b
loss = (z ** 2).sum()

# Backward pass: computes dL/dW and dL/db via reverse-mode autodiff
loss.backward()

print("dL/dW:\n", W.grad)  # shape (2,3) — the Jacobian contracted with upstream grad
print("dL/db:\n", b.grad)  # shape (2,)