Multivariable Calculus

calculus in \(\mathbb{R}^n\): functions of several variables, the surfaces they define, and the fields that flow over them. the programme is the same as one variable — linearise locally, integrate globally — but the derivative becomes a matrix, the chain rule becomes matrix multiplication, and the fundamental theorem splits into three named theorems (green, stokes, gauss) that are secretly one (Courant, Richard, 1996). this page is also the mathematical spine of machine learning: gradients, hessians, jacobians and constrained optima are chapter 5 of (Deisenroth, Marc Peter and Faisal, A. Aldo and Ong, Cheng Soon, 2020).

partial derivatives and the gradient#

definitions#

for \(f : \mathbb{R}^n \to \mathbb{R}\), the partial derivative freezes all coordinates but one:

\begin{equation} \frac{\partial f}{\partial x_i}(\mathbf{x}) = \lim_{h \to 0} \frac{f(\mathbf{x} + h\,\mathbf{e}_i) - f(\mathbf{x})}{h}, \end{equation}

and the gradient collects them, \(\nabla f = \big(\tfrac{\partial f}{\partial x_1}, \dots, \tfrac{\partial f}{\partial x_n}\big)^{\top}\). but the honest definition of differentiability is the one-variable idea verbatim — existence of a best linear approximation:

\begin{equation} f(\mathbf{x} + \mathbf{h}) = f(\mathbf{x}) + \nabla f(\mathbf{x})^{\top} \mathbf{h} + o(\lVert \mathbf{h} \rVert). \end{equation}

partials existing does not imply differentiability (partials only probe axis directions); continuous partials do. 𐃏

for differentiable \(f\), the directional derivative along a unit vector \(\mathbf{v}\) is

\begin{equation} D_{\mathbf{v}} f(\mathbf{x}) = \lim_{t \to 0} \frac{f(\mathbf{x} + t\mathbf{v}) - f(\mathbf{x})}{t} = \nabla f(\mathbf{x})^{\top} \mathbf{v}. \end{equation}

steepest ascent, with proof#

claim. among unit vectors, \(D_{\mathbf{v}} f\) is maximised by \(\mathbf{v} = \nabla f / \lVert \nabla f \rVert\), with maximum value \(\lVert \nabla f \rVert\); and \(\nabla f\) is orthogonal to the level set through \(\mathbf{x}\).

proof. cauchy–schwarz gives \(D_{\mathbf{v}} f = \nabla f^{\top} \mathbf{v} \leq \lVert \nabla f \rVert \, \lVert \mathbf{v} \rVert = \lVert \nabla f \rVert\), with equality iff \(\mathbf{v}\) is a positive multiple of \(\nabla f\) — so the gradient direction wins, and the opposite direction is steepest descent. for orthogonality: let \(\gamma(t)\) be any differentiable curve inside the level set \(\{f = c\}\) with \(\gamma(0) = \mathbf{x}\). then \(f(\gamma(t)) \equiv c\), and the chain rule gives \(\nabla f^{\top} \gamma’(0) = 0\) — the gradient kills every tangent direction. \(\blacksquare\)

this pair of facts is the whole geometry of gradient descent: move perpendicular to the contour lines, at a rate set by their crowding.

jacobian, hessian, chain rule#

the derivative is a matrix#

for \(\mathbf{f} : \mathbb{R}^n \to \mathbb{R}^m\), the derivative at \(\mathbf{x}\) is the \(m \times n\) jacobian

\begin{equation} J_{\mathbf{f}}(\mathbf{x}) = \begin{pmatrix} \frac{\partial f_1}{\partial x_1} & \cdots & \frac{\partial f_1}{\partial x_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial f_m}{\partial x_1} & \cdots & \frac{\partial f_m}{\partial x_n} \end{pmatrix}, \qquad \mathbf{f}(\mathbf{x} + \mathbf{h}) = \mathbf{f}(\mathbf{x}) + J_{\mathbf{f}}(\mathbf{x})\,\mathbf{h} + o(\lVert \mathbf{h} \rVert). \end{equation}

a scalar \(f\) has \(J_f = \nabla f^{\top}\) (a row); a curve \(\gamma : \mathbb{R} \to \mathbb{R}^m\) has \(J_{\gamma} = \gamma’\) (a column). the chain rule is now literally matrix multiplication:

\begin{equation} J_{\mathbf{g} \circ \mathbf{f}}(\mathbf{x}) = J_{\mathbf{g}}(\mathbf{f}(\mathbf{x})) \; J_{\mathbf{f}}(\mathbf{x}), \end{equation}

composition of linear approximations. every backpropagation pass in a neural network is this equation evaluated right-to-left. 𐃏

second derivatives#

for scalar \(f\), the hessian collects second partials, \(H(\mathbf{x})_{ij} = \frac{\partial^2 f}{\partial x_i \partial x_j}\). by schwarz’s theorem, if the second partials are continuous then \(\partial_i \partial_j f = \partial_j \partial_i f\) — the hessian is symmetric, which is what lets the spectral theorem act on it in the second-derivative test below.

taylor to second order#

for \(f\) twice continuously differentiable,

\begin{equation} f(\mathbf{x} + \mathbf{h}) = f(\mathbf{x}) + \nabla f(\mathbf{x})^{\top} \mathbf{h} + \tfrac{1}{2}\, \mathbf{h}^{\top} H(\mathbf{x})\, \mathbf{h} + o(\lVert \mathbf{h} \rVert^{2}). \end{equation}

the local model of any smooth function is a quadratic: a plane tilted by \(\nabla f\), curved by \(H\). newton’s method in optimisation is nothing but minimising this quadratic exactly at each step.

optimisation#

unconstrained: critical points and the hessian test#

at an interior extremum of differentiable \(f\), \(\nabla f = \mathbf{0}\) (fermat, coordinate by coordinate). classification at a critical point \(\mathbf{x}^*\) reads off the taylor quadratic \(\tfrac12 \mathbf{h}^{\top} H \mathbf{h}\):

• \(H\) positive definite (all eigenvalues \(> 0\)): strict local minimum.
• \(H\) negative definite: strict local maximum.
• \(H\) indefinite (eigenvalues of both signs): saddle point.
• \(H\) singular: the test is silent — higher-order terms decide (e.g. \(f(x,y) = x^4 + y^4\) vs \(x^4 - y^4\) at the origin).

in two variables the eigenvalue conditions compress to the familiar discriminant: with \(D = f_{xx} f_{yy} - f_{xy}^2\), \(D > 0\) with \(f_{xx} > 0\) is a min, \(D > 0\) with \(f_{xx} < 0\) a max, \(D < 0\) a saddle.

constrained: lagrange multipliers#

problem. extremise \(f(\mathbf{x})\) subject to \(g(\mathbf{x}) = 0\), with \(\nabla g \neq \mathbf{0}\) on the constraint set.

derivation. let \(\mathbf{x}^*\) be a constrained extremum and \(\gamma(t)\) any differentiable curve lying in the constraint surface with \(\gamma(0) = \mathbf{x}^*\). two observations:

• \(t \mapsto f(\gamma(t))\) has an extremum at \(0\), so \(\nabla f(\mathbf{x}^*)^{\top} \gamma’(0) = 0\): \(\nabla f\) is orthogonal to every tangent direction of the surface.
• the tangent space of \(\{g = 0\}\) at \(\mathbf{x}^*\) is exactly the orthogonal complement of \(\nabla g(\mathbf{x}^*)\) (differentiate \(g(\gamma(t)) \equiv 0\); the implicit function theorem guarantees enough curves to fill the whole complement).

so \(\nabla f(\mathbf{x}^*)\) lies in the one-dimensional space spanned by \(\nabla g(\mathbf{x}^*)\):

\begin{equation} \boxed{\;\nabla f(\mathbf{x}^*) = \lambda\, \nabla g(\mathbf{x}^*)\;} \end{equation}

for some multiplier \(\lambda \in \mathbb{R}\). geometrically: at the optimum, the level curve of \(f\) is tangent to the constraint — if they crossed transversally, sliding along the constraint would still improve \(f\). 𐃏

worked example#

maximise \(f(x,y) = xy\) on the circle \(g(x,y) = x^2 + y^2 - 1 = 0\).

• stationarity: \(\nabla f = \lambda \nabla g\) gives \(y = 2\lambda x\) and \(x = 2\lambda y\).
• eliminate: substituting, \(y = 4\lambda^2 y\). if \(y = 0\) then \(x = 0\), off the constraint; so \(\lambda = \pm\tfrac12\) and \(y = \pm x\).
• constraint: \(2x^2 = 1\), so \(x = \pm \tfrac{1}{\sqrt{2}}\).
• classify: \(f = \tfrac12\) at \(\pm(\tfrac{1}{\sqrt2}, \tfrac{1}{\sqrt2})\) (maxima, \(\lambda = \tfrac12\)), \(f = -\tfrac12\) at \(\pm(\tfrac{1}{\sqrt2}, -\tfrac{1}{\sqrt2})\) (minima). the constraint set is compact, so evt guarantees these are global.

the numeric check at the bottom of the page sweeps the circle and lands on the same point, value and multiplier.

multiple integrals#

fubini#

the double integral over a rectangle \(R = [a,b] \times [c,d]\) is defined by riemann sums over grid cells, exactly as in one dimension. fubini’s theorem: for \(f\) continuous on \(R\),

\begin{equation} \iint_R f\,dA = \int_a^b \!\! \left( \int_c^d f(x,y)\,dy \right) dx = \int_c^d \!\! \left( \int_a^b f(x,y)\,dx \right) dy \end{equation}

— iterate in either order. over non-rectangular regions, describe the region by bounding curves and let the inner limits depend on the outer variable.

change of variables#

if \(T : U \to \mathbb{R}^n\) is injective and continuously differentiable with \(\det J_T \neq 0\), then

\begin{equation} \int_{T(U)} f(\mathbf{x})\,d\mathbf{x} = \int_U f(T(\mathbf{u}))\,\bigl|\det J_T(\mathbf{u})\bigr|\,d\mathbf{u}. \end{equation}

the jacobian determinant is the local volume-scaling factor of the map — an infinitesimal \(d\mathbf{u}\)-box lands as a parallelepiped of volume \(|\det J_T|\,d\mathbf{u}\). in one dimension this is substitution; the absolute value appears because multiple integrals are unoriented.

worked example: the gaussian via polar coordinates#

polar coordinates \(T(r, \theta) = (r\cos\theta, r\sin\theta)\) have

\begin{equation} J_T = \begin{pmatrix} \cos\theta & -r\sin\theta \\ \sin\theta & r\cos\theta \end{pmatrix}, \qquad \det J_T = r. \end{equation}

let \(I = \int_{-\infty}^{\infty} e^{-x^2} dx\). square it and switch to polar:

\begin{align*} I^2 &= \int_{-\infty}^{\infty} \!\! \int_{-\infty}^{\infty} e^{-(x^2 + y^2)}\,dx\,dy = \int_0^{2\pi} \!\! \int_0^{\infty} e^{-r^2}\, r\,dr\,d\theta \\ &= 2\pi \left[ -\tfrac{1}{2} e^{-r^2} \right]_0^{\infty} = 2\pi \cdot \tfrac{1}{2} = \pi, \end{align*}

so \(I = \sqrt{\pi}\). the factor \(r\) from the jacobian is precisely what makes the integrand elementary — the trick that has no one-variable counterpart, and the normalising constant of the entire theory of probability. 𐃏

vector calculus#

div, curl, and the two kinds of derivative of a field#

a vector field \(\mathbf{F} = (P, Q, R) : \mathbb{R}^3 \to \mathbb{R}^3\) has two first-order derivatives that matter:

\begin{equation} \operatorname{div} \mathbf{F} = \nabla \cdot \mathbf{F} = \frac{\partial P}{\partial x} + \frac{\partial Q}{\partial y} + \frac{\partial R}{\partial z}, \qquad \operatorname{curl} \mathbf{F} = \nabla \times \mathbf{F} = \begin{pmatrix} \frac{\partial R}{\partial y} - \frac{\partial Q}{\partial z} \\[0.4ex] \frac{\partial P}{\partial z} - \frac{\partial R}{\partial x} \\[0.4ex] \frac{\partial Q}{\partial x} - \frac{\partial P}{\partial y} \end{pmatrix}. \end{equation}

divergence measures net outflow per unit volume (source density); curl measures circulation per unit area (local rotation). two identities anchor the algebra: \(\operatorname{curl}(\nabla f) = \mathbf{0}\) and \(\operatorname{div}(\operatorname{curl} \mathbf{F}) = 0\).

line and surface integrals#

• line integral of a field along an oriented curve \(C\) parametrised by \(\gamma(t)\), \(t \in [a,b]\):

  \begin{equation} \int_C \mathbf{F} \cdot d\mathbf{r} = \int_a^b \mathbf{F}(\gamma(t)) \cdot \gamma’(t)\,dt \end{equation}

  — work done by the field along the path. for gradient fields, \(\int_C \nabla f \cdot d\mathbf{r} = f(\gamma(b)) - f(\gamma(a))\): path-independent, the ftc for curves.

• surface integral (flux) through an oriented surface \(S\) with unit normal \(\mathbf{n}\):

  \begin{equation} \iint_S \mathbf{F} \cdot \mathbf{n}\,dS \end{equation}

  — net flow of the field through the membrane.

the three big theorems, and why they are one#

(orientations: boundaries traversed so the region stays on the left; outward normals for gauss.) each says the same thing the one-variable ftc says: accumulating a derivative over a region equals evaluating the original object on the boundary. green is stokes for flat surfaces; the modern formulation absorbs all three (and the ftc itself) into the generalised stokes theorem

\begin{equation} \int_{\partial \Omega} \omega = \int_{\Omega} d\omega \end{equation}

for differential forms \(\omega\) — one theorem, every dimension (Courant, Richard, 1996). physics runs on the special cases: maxwell’s equations pair gauss with stokes, and every conservation law is a divergence theorem in disguise.

numerics: verifying the lagrange example#

sweep the constraint circle densely and compare against the analytic answer \(x = y = \tfrac{1}{\sqrt2}\), \(f = \tfrac12\), \(\lambda = \tfrac12\):

import numpy as np

# maximise f(x,y) = xy on the unit circle x^2 + y^2 = 1.
# lagrange says: grad f = lambda grad g  =>  (y, x) = lambda (2x, 2y)
# solving gives x = y = 1/sqrt(2), f = 1/2 (and the symmetric variants).

theta = np.linspace(0, 2 * np.pi, 2_000_001)     # dense sweep of the constraint
f = np.cos(theta) * np.sin(theta)

i = f.argmax()
x, y = np.cos(theta[i]), np.sin(theta[i])
print(f"numeric argmax : x = {x:.6f}, y = {y:.6f}")
print(f"analytic       : x = y = 1/sqrt(2) = {1/np.sqrt(2):.6f}")
print(f"numeric max f  : {f[i]:.8f}   (analytic: 0.5)")

# check the gradients are parallel at the optimum: (y, x) vs (2x, 2y)
grad_f = np.array([y, x])
grad_g = np.array([2 * x, 2 * y])
lam = grad_f @ grad_g / (grad_g @ grad_g)
print(f"lambda         : {lam:.6f}   (analytic: 0.5)")
print(f"residual |grad f - lambda grad g| = {np.linalg.norm(grad_f - lam * grad_g):.2e}")

numeric argmax : x = 0.707107, y = 0.707107
analytic       : x = y = 1/sqrt(2) = 0.707107
numeric max f  : 0.50000000   (analytic: 0.5)
lambda         : 0.500000   (analytic: 0.5)
residual |grad f - lambda grad g| = 1.57e-16

the argmax lands on the analytic point to six decimals, the least-squares multiplier equals the analytic \(\lambda\), and the stationarity residual is at machine precision — the gradients really are parallel there.

see also#

• calculus — the parent map
• single-variable calculus — where every theorem here comes from
• linear algebra — jacobians, hessians and definiteness live here
• optimisation — gradient descent, newton, and constraints at industrial scale
• differential equations — fields as dynamics

References#

Courant, Richard (1996). Differential and Integral Calculus, Springer.

Deisenroth, Marc Peter and Faisal, A. Aldo and Ong, Cheng Soon (2020). Mathematics for Machine Learning, Cambridge University Press.