12. Multivariable Differential Calculus
Most quantities in the real world depend on more than one input: the temperature in a room depends on three coordinates, the cost of a product on many prices, and the error of a machine-learning model on millions of weights. This chapter extends limits, derivatives, linear approximation, and optimization from functions of one variable (Chapters 2 to 4) to functions of several variables, using the vector language of Chapter 11. Its central object, the gradient, is the tool behind steepest-ascent reasoning, tangent planes, and gradient-based optimization in science and engineering.
<a id="121-functions-of-several-variables"></a>
12.1 Functions of Several Variables
This subchapter defines functions with several inputs and develops the two ways to picture them: as a surface in space, and as a contour map of level curves.
Domain and Range in Several Variables
Theory
A function of one variable takes a number and returns a number. Many quantities need more than one number as input. The volume of a cylinder needs both a radius and a height, $V = \pi r^2 h$. The wind chill needs both the air temperature and the wind speed. So we allow the input to be a point in the plane or in space.
A function of two variables $f$ is a rule that assigns to each ordered pair $(x, y)$ in a set $D \subseteq \mathbb{R}^2$ exactly one real number $f(x, y)$. We write $f : D \to \mathbb{R}$ and often $z = f(x, y)$.
| Term | Meaning | One-variable analogue |
|---|---|---|
| Domain $D$ | The set of input points $(x, y)$ where $f$ is defined | An interval or union of intervals |
| Range | The set of output values ${ f(x, y) : (x, y) \in D }$ | Same idea |
| Independent variables | $x$ and $y$ | $x$ |
| Dependent variable | $z = f(x, y)$ | $y = f(x)$ |
As in Section 1.2, when a function is given by a formula and no domain is stated, its natural domain is the set of all points where the formula produces a real number. The restrictions are the same as before, now applied to expressions in two variables:
| Expression | Restriction | Typical domain shape |
|---|---|---|
| $\sqrt{g(x, y)}$ | $g(x, y) \ge 0$ | A closed region bounded by the curve $g = 0$ |
| $\ln g(x, y)$ | $g(x, y) > 0$ | An open region (boundary excluded) |
| $\frac{p(x, y)}{q(x, y)}$ | $q(x, y) \ne 0$ | The plane with a curve removed |
| $\arcsin g(x, y)$ | $-1 \le g(x, y) \le 1$ | A band between two curves |
The important change is geometric. In one variable a domain is a piece of a line. In two variables it is a region of the plane, and describing it means drawing it: a disk, a half-plane, the inside of a parabola, the plane minus a line. The boundary curve is found by turning the inequality into an equation, and then a test point tells you which side is included.
The range is usually harder. A good strategy is to think about which outputs the formula can produce as the input moves around the domain, often by spotting a single combination such as $x^2 + y^2$ that controls the output.
Examples
- $f(x, y) = \sqrt{9 - x^2 - y^2}$. We need $x^2 + y^2 \le 9$, so the domain is the closed disk of radius $3$ centered at the origin. The quantity $9 - x^2 - y^2$ ranges from $0$ (on the boundary circle) to $9$ (at the center), so the range is $[0, 3]$.
- $g(x, y) = \ln(y - x^2)$. We need $y > x^2$: the region strictly above the parabola $y = x^2$. Since $y - x^2$ takes every positive value, the range is all of $\mathbb{R}$.
- $h(x, y) = \frac{x + y}{x - y}$. The domain is every point except the line $y = x$. For the range, try to hit an arbitrary value $c$: $h = c$ means $(1 + c),y = (c - 1),x$. For $c \ne -1$ take $x = 1 + c$, $y = c - 1$ (then $x - y = 2 \ne 0$), and for $c = -1$ take $(0, 1)$. Every $c$ is reached, so the range is $\mathbb{R}$.
- $V(r, h) = \pi r^2 h$ as a formula has domain $\mathbb{R}^2$, but as a physical model its domain is $r \ge 0$, $h \ge 0$. As in one variable, the context can shrink the domain.
Domain of g(x, y) = ln(y - x^2): the open region above the parabola y = x^2
y
|
\ | /
\####|####/ # = in the domain (y > x^2)
\###|###/
\##|##/ the parabola itself is NOT included (dashed boundary)
------`-+-'-------- x
| test point (0, 1): 1 - 0 > 0, so the inside is in
The code below detects a natural domain the way a computer sees it: evaluate the formula on a grid, and the points that return NaN are outside. The fraction of valid grid points also estimates the area of the domain.
import numpy as np
# Natural domains as boolean masks on a grid over the square [-4, 4] x [-4, 4]
x, y = np.meshgrid(np.linspace(-4, 4, 801), np.linspace(-4, 4, 801))
with np.errstate(invalid='ignore', divide='ignore'):
f = np.sqrt(9 - x**2 - y**2) # NaN outside the disk x^2 + y^2 <= 9
g = np.log(y - x**2) # NaN or -inf on and below the parabola
print(round(np.isfinite(f).mean() * 64, 2), round(9 * np.pi, 2)) # domain area vs exact pi * 3^2
# 28.2 28.27
print(np.nanmin(f), np.nanmax(f)) # range of f is [0, 3]
# 0.0 3.0
# ASCII picture of the domain of g(x, y) = ln(y - x^2): '#' where y > x^2
for yv in np.linspace(3, -1, 9):
row = ''.join('#' if yv - xv**2 > 0 else '.' for xv in np.linspace(-2.5, 2.5, 41))
print(f'{yv:5.1f} {row}')
# 3.0 .......###########################.......
# 2.5 ........#########################........
# 2.0 .........#######################.........
# 1.5 ...........###################...........
# 1.0 .............###############.............
# 0.5 ...............###########...............
# 0.0 .........................................
# -0.5 .........................................
# -1.0 .........................................
Key Takeaways
- A function of two variables maps points $(x, y)$ of a region $D \subseteq \mathbb{R}^2$ to real numbers.
- The natural domain comes from the same restrictions as in one variable (roots, logarithms, denominators), but it is now a region: sketch its boundary curve and use a test point to pick the side.
- Track whether the boundary is included ($\ge$ versus $>$). This matters later for limits, continuity, and absolute extrema.
- For the range, look for a single combination of the variables (like $x^2 + y^2$) that controls the output.
🧪 Practice
- Find and sketch the domain of $f(x, y) = \sqrt{x + y} + \sqrt{1 - x}$. Is the boundary included?
- Find the domain and range of $f(x, y) = \frac{1}{\sqrt{16 - 4x^2 - y^2}}$.
- Find the domain of $f(x, y) = \arcsin(x^2 + y^2 - 2)$ and describe it in words. What is its range?
- Interview-style: A data pipeline computes $\log(\text{income} - \text{expenses})$ for every customer row and silently produces
NaNfor a few percent of them. Describe how you would detect, report, and handle the rows that fall outside the domain. Hint: the domain is the half-plane where income exceeds expenses. Decide explicitly what the model should do on the other side (clip, transform with something like $\operatorname{sign}(u)\ln(1 + \lvert u \rvert)$, or flag), rather than lettingNaNpropagate.
Graphs and Surfaces
Theory
The graph of a function of one variable is a curve: all points $(x, f(x))$ in the plane. For a function of two variables, the inputs already use up a plane, so the output needs a third axis. The graph of $f$ is the set
$${ (x, y, z) \in \mathbb{R}^3 : z = f(x, y),\ (x, y) \in D },$$
a surface in space. Think of the domain as a floor plan and $f(x, y)$ as the height of a roof (or the ground) above each point of the floor. Mountains, valleys, and ridges become features of this surface.
Each input point has exactly one output, so a vertical line through any point of $D$ meets the graph exactly once. This is the vertical line test for surfaces. A sphere $x^2 + y^2 + z^2 = 1$ fails it and is not the graph of one function, but its upper half $z = \sqrt{1 - x^2 - y^2}$ passes.
Sketching with traces. A surface is hard to draw directly. Instead we slice it with planes and draw the curves we get, called traces or cross-sections (the same method used for quadric surfaces in Section 11.1):
| Slice | Equation | What you get |
|---|---|---|
| Vertical plane $x = k$ | $z = f(k, y)$ | A curve in the $yz$-direction: the profile walking north-south |
| Vertical plane $y = k$ | $z = f(x, k)$ | A curve in the $xz$-direction: the profile walking east-west |
| Horizontal plane $z = k$ | $f(x, y) = k$ | A level curve (next topic) |
A catalog of common graphs.
| Function | Surface | Key feature |
|---|---|---|
| $f = ax + by + c$ | Plane | Straight-line traces in every vertical slice |
| $f = x^2 + y^2$ | Circular paraboloid | Bowl; parabola traces, circle level curves |
| $f = \sqrt{x^2 + y^2}$ | Cone | Straight-line traces through the origin; sharp tip |
| $f = \sqrt{r^2 - x^2 - y^2}$ | Upper hemisphere | Domain is a disk |
| $f = y^2 - x^2$ | Hyperbolic paraboloid (saddle) | Up in one direction, down in the other |
| $f = g(x)$ (no $y$) | Cylinder | The curve $z = g(x)$ extruded along the $y$-axis |
| $f = g(x^2 + y^2)$ | Surface of revolution | Rotate the profile $z = g(r)$ around the $z$-axis |
The last row is a useful shortcut. Whenever $f$ depends only on $r = \sqrt{x^2 + y^2}$, the surface is symmetric around the $z$-axis, and you only need to sketch one profile.
Examples
Sketch $z = 4 - x^2 - 2y^2$.
- Trace $x = k$: $z = (4 - k^2) - 2y^2$, a downward parabola in $y$ whose top drops as $\lvert k \rvert$ grows.
- Trace $y = k$: $z = (4 - 2k^2) - x^2$, a downward parabola in $x$.
- Trace $z = k$: $x^2 + 2y^2 = 4 - k$, an ellipse for $k < 4$, a point at $k = 4$, and empty for $k > 4$.
So the surface is an upside-down elliptic paraboloid with its peak at $(0, 0, 4)$, steeper in the $y$-direction.
z
| peak (0, 0, 4)
| _.-*-._
| .-' | '-. trace y = 0: z = 4 - x^2
| / | \ trace x = 0: z = 4 - 2y^2 (narrower)
| / | \
|/ | \
-----+---------+---------+------ x
-2 0 2
A table of heights is a "numerical picture" of a surface: reading along a row is walking along a trace. The saddle $z = x^2 - y^2$ goes up along the $x$-axis and down along the $y$-axis.
import numpy as np
import sympy as sp
x, y, k = sp.symbols('x y k', real=True)
f = 4 - x**2 - 2*y**2 # elliptic paraboloid opening downward
print(sp.expand(f.subs(x, k))) # trace in the plane x = k: parabola in y
# -k**2 - 2*y**2 + 4
print(sp.expand(f.subs(y, k))) # trace in the plane y = k: parabola in x
# -2*k**2 - x**2 + 4
print(sp.solve(sp.Eq(f, k), y)) # horizontal trace z = k: an ellipse, solved for y
# [-sqrt(-2*k - 2*x**2 + 8)/2, sqrt(-2*k - 2*x**2 + 8)/2]
xs = np.arange(-2, 3)
F = sp.lambdify((x, y), f)
S = sp.lambdify((x, y), x**2 - y**2) # saddle for comparison
print('paraboloid z = 4 - x^2 - 2y^2 saddle z = x^2 - y^2')
for yv in xs[::-1]:
print(f'y={yv:2d}', ' '.join(f'{F(xv, yv):4d}' for xv in xs), ' ', ' '.join(f'{S(xv, yv):3d}' for xv in xs))
# paraboloid z = 4 - x^2 - 2y^2 saddle z = x^2 - y^2
# y= 2 -8 -5 -4 -5 -8 0 -3 -4 -3 0
# y= 1 -2 1 2 1 -2 3 0 -1 0 3
# y= 0 0 3 4 3 0 4 1 0 1 4 row y = 0: the saddle rises in x
# y=-1 -2 1 2 1 -2 3 0 -1 0 3
# y=-2 -8 -5 -4 -5 -8 0 -3 -4 -3 0 middle column x = 0: it falls in y
Key Takeaways
- The graph of $z = f(x, y)$ is a surface in $\mathbb{R}^3$; each point of the domain sits under exactly one point of the surface.
- Sketch surfaces with traces: vertical slices $x = k$ and $y = k$ give profiles, horizontal slices $z = k$ give level curves.
- Recognize the standard shapes: planes, paraboloids, cones, hemispheres, saddles, cylinders, and surfaces of revolution $z = g(\sqrt{x^2 + y^2})$.
🧪 Practice
- Use traces to sketch $z = 1 + x^2 - y$. What kind of surface is it?
- Identify the surface $z = e^{-(x^2 + y^2)}$ by writing it in terms of $r$, and sketch its profile.
- Match each function to its description and explain your reasoning: $f_1 = \lvert x \rvert + \lvert y \rvert$, $f_2 = \sin x$, $f_3 = \sqrt{x^2 + y^2}$; (a) a corrugated sheet, (b) an upside-down square pyramid, (c) a cone.
- Interview-style: You are rendering a terrain from a height map stored as a 2D array. How would you detect cliffs (places where the surface is nearly vertical) from the array alone? Hint: compare neighboring heights. A cliff is where the difference between adjacent cells, divided by the cell spacing, is very large in some direction. This is a discrete version of the partial derivatives of Section 12.3.
Level Curves and Level Surfaces
Theory
A hiking map does not draw the mountain in 3D. It draws contour lines: curves joining points at the same altitude. From these lines alone you can read where the peaks are, where the ground is steep, and where it is flat. This is the second, and often more practical, way to picture a function of two variables.
The level curves (or contour lines) of $f$ are the curves
$$f(x, y) = k$$
for constants $k$ in the range of $f$. Each one is the horizontal trace of the graph at height $k$, pushed down onto the $xy$-plane. A collection of level curves for equally spaced values of $k$ is a contour map.
Reading a contour map.
| Feature of the map | Meaning for the surface |
|---|---|
| Level curves close together | Steep: height changes quickly over a short distance |
| Level curves far apart | Gentle slope or nearly flat |
| Closed curves shrinking to a point, values increasing inward | A peak (local maximum) |
| Closed curves shrinking to a point, values decreasing inward | A pit (local minimum) |
| Two level curves crossing in an X | A saddle point (pass between two hills) |
| Level curves never cross otherwise | One input has only one output |
The spacing is quantitative. If equally spaced levels $k, k + \Delta k, k + 2\Delta k$ have curves a distance $d$ apart, the slope in that region is about $\frac{\Delta k}{d}$. Compare the paraboloid $f = x^2 + y^2$ and the cone $g = \sqrt{x^2 + y^2}$: both have circles as level curves, but the circles of $f$ for $k = 1, 2, 3, 4$ have radii $1, 1.41, 1.73, 2$ and crowd together as you move out (getting steeper), while the circles of $g$ have radii $1, 2, 3, 4$ and are equally spaced (constant slope).
Level surfaces. For a function of three variables $w = f(x, y, z)$ we cannot draw the graph (it lives in four dimensions), but we can draw the level surfaces $f(x, y, z) = k$. For $f = x^2 + y^2 + z^2$ these are spheres. For a temperature function they are isothermal surfaces. The next topic uses them heavily.
Level sets are everywhere in applications:
| Field | Level curves are called | Function |
|---|---|---|
| Cartography | Contour lines | Altitude |
| Meteorology | Isobars, isotherms | Pressure, temperature |
| Economics | Indifference curves, isoquants | Utility, production |
| Physics | Equipotentials | Electric potential |
| Machine learning | Loss contours, decision boundaries | Loss, classifier score |
Examples
- $f(x, y) = x^2 + y^2$: $x^2 + y^2 = k$ is a circle of radius $\sqrt{k}$ for $k > 0$, the origin for $k = 0$, and empty for $k < 0$.
- $f(x, y) = xy$: $xy = k$ is a hyperbola in quadrants I and III for $k > 0$, in quadrants II and IV for $k < 0$, and the two axes for $k = 0$. The axes cross at the origin, which signals a saddle.
- $f(x, y) = 2x + 3y$: the level curves $2x + 3y = k$ are parallel lines, equally spaced, because a plane has constant slope.
The code below draws contour maps in plain text. It evaluates $f$ on a grid, assigns each grid point to the band between consecutive levels, and draws a symbol wherever neighboring points are in different bands (a level curve passes between them). The grid is twice as wide as it is tall because text characters are about twice as tall as they are wide.
import numpy as np
def ascii_contours(f, levels, symbols, lim=2.0, width=49, height=21):
"""Draw the level curves f(x, y) = c for c in `levels`, marking curve c with symbols[c]."""
X, Y = np.meshgrid(np.linspace(-lim, lim, width), np.linspace(lim, -lim, height))
band = np.searchsorted(levels, f(X, Y)) # number of levels below each grid value
for i in range(height):
row = ''
for j in range(width):
nbrs = ([band[i, j + 1]] if j + 1 < width else []) + ([band[i + 1, j]] if i + 1 < height else [])
crossed = [n for n in nbrs if n != band[i, j]]
# a level curve runs between this point and a neighbor in a different band
row += symbols[min(band[i, j], crossed[0])] if crossed else ' '
print(row)
ascii_contours(lambda x, y: x**2 + y**2, [1, 2, 3, 4], '1234') # paraboloid: circles
ascii_contours(lambda x, y: x * y, [-1, 0, 1], '-0+') # saddle z = xy: hyperbolas
Output, with the two maps printed side by side (paraboloid on the left, saddle on the right):
444444444444444444444 - 0 +
4444 333333333333333 44444 - 0 ++
444 33333 22222 333333 4444 - 0 ++
44 33 222222 2222222 333 444 - 0 +++
4 33 2222 1 22222 333 44 -- 0 +++
4 33 2 111111111111111 22 333 44 ---- 0 ++++
4 3 22 11 111 222 33 44 ----- 0 ++++++
4 3 2 1 11 22 33 44 ---- 0 +++++
4 3 2 1 11 2 3 4 0
4 3 2 1 11 2 3 44 0000000000000000000000000
4 3 2 11 1 2 3 4 000000000000000000000000
4 3 2 11 1 2 3 4 0
44 33 22 11 1 2 3 4 +++++ 0 -----
44 33 222 111 11 22 3 4 ++++++ 0 -----
44 333 22 11111111 1111111 2 33 4 +++ 0 ---
44 333 22222 11 2222 33 4 ++++ 0 -
444 333 2222222 222222 33 44 ++ 0 --
4444 333333 222222 33333 444 ++ 0 -
44444 3333333333333333 4444 ++ 0 -
44444444444 4444444444 + 0 -
44 + 0 -
On the left, the gaps between the circles shrink from the center outward: the bowl gets steeper. On the right, the $0$ level is the two axes, crossing at the saddle point, and the $+1$ and $-1$ hyperbolas sit in alternating quadrants.
Key Takeaways
- A level curve $f(x, y) = k$ is the set of inputs giving output $k$; a contour map shows several of them for equally spaced $k$.
- Close spacing means steep, wide spacing means gentle; closed loops mark peaks or pits, and crossing level curves mark saddles.
- For functions of three variables, level surfaces $f(x, y, z) = k$ replace the (four-dimensional) graph.
🧪 Practice
- Sketch the level curves of $f(x, y) = y - x^2$ for $k = -2, -1, 0, 1, 2$. Are they equally spaced? What does that say about the surface?
- Describe the level curves of $f(x, y) = \frac{x^2 + y^2}{1 + x^2 + y^2}$. Which values of $k$ give empty level curves?
- Describe the level surfaces of $f(x, y, z) = x^2 + y^2 - z^2$ for $k > 0$, $k = 0$, and $k < 0$ (Section 11.1 quadric surfaces will help).
- Interview-style: On a contour plot of a model's loss as a function of two weights, the level curves are long, thin ellipses. Why does plain gradient descent zigzag on such a landscape, and what could you do about it? Hint: steepness is very different along the two axes of the ellipse. Section 12.5 shows that the steepest direction is perpendicular to the level curves, which on a thin ellipse points mostly across the valley rather than along it. Rescaling the variables (preconditioning, feature normalization) makes the ellipses rounder.
Functions of Three or More Variables
Theory
Nothing in the definition of a function limits the number of inputs. The temperature in a room is $T(x, y, z)$. The temperature that also changes over time is $T(x, y, z, t)$. The price a model predicts for a house may depend on 50 features. A neural network's loss is a function of all its weights, often billions of them.
A function of $n$ variables assigns a real number $f(x_1, x_2, \ldots, x_n)$ to each point of a domain $D \subseteq \mathbb{R}^n$. Collecting the inputs into a vector $\mathbf{x} = \langle x_1, \ldots, x_n \rangle$ lets us write $f(\mathbf{x})$, and almost every formula in this chapter looks the same for every $n$ once written with vectors. Two families appear constantly:
| Family | Formula | Example |
|---|---|---|
| Linear (affine) function | $f(\mathbf{x}) = \mathbf{c} \cdot \mathbf{x} + b = c_1x_1 + \cdots + c_nx_n + b$ | Total cost of a shopping basket |
| Quadratic form | $f(\mathbf{x}) = \mathbf{x}^{\mathsf{T}} A,\mathbf{x} = \sum_{i,j} a_{ij} x_i x_j$ | Energy of a spring system; the local shape near a minimum (Section 12.6) |
How to visualize. We can see at most three dimensions, so:
- For $f(x, y, z)$, draw level surfaces $f = k$, like layers of an onion.
- Fix all but one or two variables and look at the resulting slice ("hold everything else constant"). This is exactly the idea behind partial derivatives in Section 12.3.
- For large $n$, think in terms of formulas and vectors, not pictures. Most geometric facts from two dimensions (distance, angles, level sets, gradients) carry over unchanged.
Distance in $\mathbb{R}^n$. The distance between $\mathbf{x}$ and $\mathbf{a}$ is $\lvert \mathbf{x} - \mathbf{a} \rvert = \sqrt{(x_1 - a_1)^2 + \cdots + (x_n - a_n)^2}$, the natural extension of Section 11.1. "Close to $\mathbf{a}$" always means "small distance", and that is what limits in Section 12.2 are built on.
Examples
- $f(x, y, z) = \ln(z - x^2 - y^2)$ has domain ${ z > x^2 + y^2 }$, the solid region strictly inside the paraboloid.
- $f(x, y, z) = x^2 + 4y^2 + 9z^2$ has level surfaces $x^2 + 4y^2 + 9z^2 = k$, which are ellipsoids for $k > 0$, nested around the origin.
- A room heater at $H = (1, 2, 0.5)$ gives $T(\mathbf{p}) = 20 + 15e^{-\lvert \mathbf{p} - H \rvert^2}$. $T$ depends on the point only through its distance $d$ to the heater, so every level surface is a sphere around $H$: solving $20 + 15e^{-d^2} = c$ gives radius $d = \sqrt{\ln \frac{15}{c - 20}}$.
import numpy as np
# Temperature (deg C) in a room with a heater at H = (1, 2, 0.5): a function of three variables
H = np.array([1.0, 2.0, 0.5])
def T(p):
d2 = np.sum((p - H)**2, axis=-1) # squared distance to the heater; works on many points
return 20 + 15 * np.exp(-d2)
# The level surface T = c is a sphere around the heater of radius sqrt(ln(15 / (c - 20)))
rng = np.random.default_rng(0)
for c in [22, 26, 30]:
radius = np.sqrt(np.log(15 / (c - 20)))
u = rng.normal(size=(1000, 3))
u /= np.linalg.norm(u, axis=1, keepdims=True) # 1000 random unit directions
temps = T(H + radius * u) # points on the predicted sphere
print(c, round(radius, 4), temps.min().round(10), temps.max().round(10))
# 22 1.4195 22.0 22.0 every point on the sphere has the same temperature
# 26 0.9572 26.0 26.0
# 30 0.6368 30.0 30.0 hotter isotherms are smaller spheres
# Five variables: the mean squared error of a linear model, as a function of its weights w
X = rng.normal(size=(200, 5))
w_true = np.array([1.0, -2.0, 0.5, 0.0, 3.0])
y = X @ w_true
def loss(w):
return np.mean((X @ w - y)**2) # f: R^5 -> R
print(loss(w_true), round(loss(np.zeros(5)), 3), round(loss(w_true + 0.1), 3))
# 0.0 15.956 0.047 minimum value 0 at the true weights
Key Takeaways
- Functions can take any number of inputs; write them as $f(\mathbf{x})$ with $\mathbf{x} \in \mathbb{R}^n$.
- Visualize $f(x, y, z)$ with level surfaces, and higher-dimensional functions with slices and formulas rather than pictures.
- Linear functions $\mathbf{c} \cdot \mathbf{x} + b$ and quadratic forms $\mathbf{x}^{\mathsf{T}} A\mathbf{x}$ are the building blocks of the local approximations used in the rest of the chapter.
- Distance $\lvert \mathbf{x} - \mathbf{a} \rvert$ in $\mathbb{R}^n$ defines what "near" means, for every $n$.
🧪 Practice
- Find the domain of $f(x, y, z) = \sqrt{1 - x^2 - y^2 - z^2} + \ln z$ and describe it geometrically.
- Describe the level surfaces of $f(x, y, z) = x + 2y - z$ and of $g(x, y, z) = x^2 + y^2$.
- Write $f(x, y) = 3x^2 + 4xy - y^2$ as $\mathbf{x}^{\mathsf{T}} A,\mathbf{x}$ with a symmetric matrix $A$.
- Interview-style: You need to find the minimum of a function of 20 variables by evaluating it on a grid with 10 values per variable. Why is this a bad plan, and what does it suggest about how optimization must work in high dimensions? Hint: count the grid points ($10^{20}$). Methods that use local slope information (gradients, Section 12.5) need work that grows roughly linearly with $n$ per step. <a id="122-limits-and-continuity-in-several-variables"></a>
12.2 Limits and Continuity in Several Variables
This subchapter extends limits and continuity to functions of several variables, where the new difficulty is that a point can be approached from infinitely many directions and along curved paths.
Limits Along Paths
Theory
In one variable, $x$ can approach $a$ from only two sides, and the limit exists exactly when the left and right limits agree (Section 2.1). In the plane, a point $(x, y)$ can approach $(a, b)$ along the $x$-axis, along a line of any slope, along a parabola, along a spiral, or by jumping around. The limit has to ignore all of these differences.
Definition. We say $\lim_{(x, y) \to (a, b)} f(x, y) = L$ if $f(x, y)$ gets arbitrarily close to $L$ whenever $(x, y)$ is close enough to $(a, b)$, but not equal to it. Precisely (compare Section 2.5): for every $\varepsilon > 0$ there is a $\delta > 0$ such that
$$0 < \sqrt{(x - a)^2 + (y - b)^2} < \delta \implies \lvert f(x, y) - L \rvert < \varepsilon.$$
The only change from one variable is that $\lvert x - a \rvert$ became the distance between points. The condition talks about a whole punctured disk around $(a, b)$, not about any particular route into it.
One variable: two ways in Two variables: infinitely many ways in
----> a <---- \ | / ~~~~~
\ | / ~~
----(a,b)---- lines, parabolas,
/ | \ ~~ spirals, ...
/ | \ ~~~~~
The path principle. If the limit is $L$, then $f$ approaches $L$ along every path into $(a, b)$, because every path eventually enters every small disk. This gives a powerful one-way test:
| What you find along paths | What you can conclude |
|---|---|
| Two paths give different values | The limit does not exist |
| One path gives no limit at all | The limit does not exist |
| Many paths all give $L$ | Nothing yet. It suggests $L$ but proves nothing |
The last row is the trap. No finite list of paths can prove a limit exists, because there are always more paths you did not check. Proving existence needs an inequality that controls $f$ near $(a, b)$ regardless of direction (Squeeze Theorem from Section 2.2, or polar coordinates later in this subchapter).
Limit laws still hold. Sums, products, quotients (with nonzero denominator limit) and compositions with continuous functions work exactly as in Section 2.2. So limits of "nice" expressions at points where nothing breaks are found by direct substitution: $\lim_{(x, y) \to (1, 2)} (x^2y + e^{xy}) = 2 + e^2$. Paths only matter at points where substitution fails, usually $\frac{0}{0}$.
Examples
$f(x, y) = \frac{x^2 - y^2}{x^2 + y^2}$ near $(0, 0)$:
- Along the $x$-axis ($y = 0$): $f(x, 0) = \frac{x^2}{x^2} = 1$.
- Along the $y$-axis ($x = 0$): $f(0, y) = \frac{-y^2}{y^2} = -1$.
Two paths, two different values, so the limit does not exist.
$g(x, y) = \frac{x^2y}{x^2 + y^2}$ near $(0, 0)$: along the axes, lines, parabolas, and even a spiral, the values go to $0$. That suggests the limit is $0$ but does not prove it. A proof: since $x^2 \le x^2 + y^2$,
$$\lvert g(x, y) \rvert = \frac{x^2}{x^2 + y^2},\lvert y \rvert \le \lvert y \rvert \to 0,$$
so the Squeeze Theorem gives $\lim g = 0$.
import numpy as np
def along_paths(f, paths, ts=(1e-1, 1e-2, 1e-3, 1e-4)):
"""Evaluate f along parametrized paths (x(t), y(t)) -> (0, 0) as t -> 0+."""
for name, (px, py) in paths.items():
vals = [f(px(t), py(t)) for t in ts]
print(f'{name:14s}', ' '.join(f'{v: .6f}' for v in vals))
f = lambda x, y: (x**2 - y**2) / (x**2 + y**2)
along_paths(f, {
'x-axis': (lambda t: t, lambda t: 0 * t),
'y-axis': (lambda t: 0 * t, lambda t: t),
'line y = x': (lambda t: t, lambda t: t),
})
# x-axis 1.000000 1.000000 1.000000 1.000000
# y-axis -1.000000 -1.000000 -1.000000 -1.000000 disagrees with the x-axis: no limit
# line y = x 0.000000 0.000000 0.000000 0.000000
g = lambda x, y: x**2 * y / (x**2 + y**2) # this one does have a limit
along_paths(g, {
'x-axis': (lambda t: t, lambda t: 0 * t),
'line y = 2x':(lambda t: t, lambda t: 2 * t),
'parabola': (lambda t: t, lambda t: t**2),
'spiral': (lambda t: t * np.cos(1 / t), lambda t: t * np.sin(1 / t)),
})
# x-axis 0.000000 0.000000 0.000000 0.000000
# line y = 2x 0.040000 0.004000 0.000400 0.000040
# parabola 0.009901 0.000100 0.000001 0.000000
# spiral -0.038301 -0.003765 0.000262 -0.000028 consistent with 0, but only the squeeze proves it
Key Takeaways
- $\lim_{(x, y) \to (a, b)} f = L$ means $f$ is close to $L$ on a whole small punctured disk around $(a, b)$, whatever the route.
- If the limit exists, every path gives the same value. Two disagreeing paths prove the limit does not exist.
- Agreement along many paths never proves existence; use an inequality (squeeze, polar bounds) for that.
- Away from trouble spots, limit laws and direct substitution work exactly as in one variable.
🧪 Practice
- Find $\lim_{(x, y) \to (2, -1)} \frac{x^2 + xy}{x + y + 1}$ by direct substitution, and explain why substitution is allowed.
- Evaluate $\frac{2xy}{x^2 + 3y^2}$ along the $x$-axis and along the line $y = x$ as $(x, y) \to (0, 0)$. What do you conclude?
- Use the inequality $\lvert xy \rvert \le \frac{1}{2}(x^2 + y^2)$ to prove that $\lim_{(x, y) \to (0, 0)} \frac{xy}{\sqrt{x^2 + y^2}} = 0$.
- Interview-style: A colleague writes a unit test that "verifies" a two-variable limit by checking the function along the two axes and the two diagonals at distance $10^{-8}$ from the point. What can this test miss, and how could you make it more trustworthy? Hint: think about curved paths like $y = x^2$, and about the worst value over a whole small circle rather than a few directions. Also consider floating-point cancellation very close to the point.
Proving Limits Do Not Exist
Theory
The path principle gives the standard method for showing a limit fails, called the two-path test: find two paths into $(a, b)$ along which $f$ approaches different values. The skill is in choosing the paths.
A systematic search. For a limit at the origin:
- Try the axes ($y = 0$, then $x = 0$). They often simplify the expression the most.
- Try all lines at once: substitute $y = mx$ and simplify. If the result depends on $m$, different slopes give different values, and you are done.
- If every line gives the same value, look at the degrees in the expression and try curves $y = mx^k$ (or $x = my^k$) that balance them. Choose $k$ so that the terms in the denominator have the same order of smallness.
- If every path you try agrees, switch to proving the limit exists instead.
Why curves may be needed. In $\frac{x^2y}{x^4 + y^2}$ along $y = mx$, the denominator $x^4 + m^2x^2$ is dominated by $m^2x^2$, and the whole fraction behaves like $\frac{mx^3}{m^2x^2} \to 0$. The $x^4$ term never competes. Along $y = mx^2$, both denominator terms are of order $x^4$, they balance, and the true behavior appears.
| Expression at $(0, 0)$ | Revealing path | Values |
|---|---|---|
| $\frac{xy}{x^2 + y^2}$ | $y = mx$ | $\frac{m}{1 + m^2}$, depends on $m$ |
| $\frac{x^2 - y^2}{x^2 + y^2}$ | The axes | $1$ and $-1$ |
| $\frac{x^2y}{x^4 + y^2}$ | $y = mx^2$ (lines all give $0$) | $\frac{m}{1 + m^2}$ |
| $\frac{xy^2}{x^2 + y^4}$ | $x = my^2$ | $\frac{m}{m^2 + 1}$ |
Examples
- $\lim_{(x, y) \to (0, 0)} \frac{xy}{x^2 + y^2}$. With $y = mx$: $\frac{mx^2}{(1 + m^2)x^2} = \frac{m}{1 + m^2}$. Slope $0$ gives $0$, slope $1$ gives $\frac{1}{2}$. The limit does not exist.
- $\lim_{(x, y) \to (0, 0)} \frac{x^2y}{x^4 + y^2}$. Along every line $y = mx$ the limit is $0$, and along the $y$-axis $f = 0$. But along $y = x^2$: $\frac{x^4}{2x^4} = \frac{1}{2}$. The limit does not exist, even though every straight-line approach agrees. On the surface, there is a ridge of height $\frac{1}{2}$ that follows the parabola into the origin, slipping between all the lines.
Level curves of x^2 y / (x^4 + y^2) near the origin: the parabolas y = m x^2 (value m / (1 + m^2))
y = x^2 (value 1/2)
\ /
'. .' every straight line through the origin
'-. .-' eventually lies BELOW this parabola,
-------'--*--'------- so it sees values near 0
.-' '-.
.' '.
/ \
y = -x^2 (value -1/2)
import sympy as sp
x, m = sp.symbols('x m', real=True)
y = sp.Symbol('y', real=True)
f = x * y / (x**2 + y**2)
along_line = sp.simplify(f.subs(y, m * x)) # restrict to the line y = m x
print(along_line, sp.limit(along_line, x, 0)) # value depends on the slope m
# m/(m**2 + 1) m/(m**2 + 1)
print([along_line.subs(m, v) for v in (0, 1, -1, 2)])
# [0, 1/2, -1/2, 2/5]
g = x**2 * y / (x**4 + y**2)
print(sp.limit(g.subs(y, m * x), x, 0)) # every straight line gives 0 ...
# 0
print(sp.simplify(g.subs(y, x**2))) # ... but the parabola y = x^2 gives 1/2
# 1/2
print(sp.simplify(g.subs(y, m * x**2))) # family of parabolas: value m/(m^2 + 1)
# m/(m**2 + 1)
Key Takeaways
- Two-path test: two paths with different limiting values prove the limit does not exist.
- Substitute $y = mx$ to test all lines at once; a result depending on $m$ settles it.
- If all lines agree, balance the degrees with curves $y = mx^k$; lines can miss behavior that hugs a curve.
- Agreement along every path you try is a signal to switch to proving existence.
🧪 Practice
- Show that $\lim_{(x, y) \to (0, 0)} \frac{3x^2y^2}{x^4 + y^4}$ does not exist.
- Show that $\lim_{(x, y) \to (0, 0)} \frac{xy^2}{x^2 + y^4}$ does not exist, even though it is $0$ along every line.
- For which values of the exponent $p > 0$ does $\lim_{(x, y) \to (0, 0)} \frac{\lvert x \rvert^p,y}{x^2 + y^2}$ exist? Use paths to rule out small $p$ and an inequality to prove existence for large $p$.
- Interview-style: Given a black-box function $f(x, y)$ that you can evaluate but not read, how would you write a heuristic that flags "this limit at $(0, 0)$ probably does not exist"? Hint: evaluate on shrinking circles and track the spread (max minus min) of the values on each circle. A spread that does not shrink signals non-existence. To catch ridges along curves like $y = x^2$, sample many angles, since such ridges occupy a sliver of angles that narrows as the circle shrinks.
Continuity of Multivariable Functions
Theory
The definition is the same as in Section 2.4, with the multivariable limit plugged in. A function $f$ is continuous at $(a, b)$ if
$$\lim_{(x, y) \to (a, b)} f(x, y) = f(a, b).$$
That packs three requirements: $f(a, b)$ is defined, the limit exists, and they are equal. A function is continuous on a set $D$ if it is continuous at every point of $D$. Intuitively, the graph has no tears, holes, or jumps: small moves in the input produce small changes in the output, in every direction at once.
Building continuous functions. The algebra of continuous functions carries over:
| Function type | Continuous where |
|---|---|
| Polynomials, like $x^3y - 4xy^2 + 7$ | Everywhere in $\mathbb{R}^2$ |
| Rational functions $\frac{p(x, y)}{q(x, y)}$ | Wherever $q(x, y) \ne 0$ |
| Sums, products, quotients of continuous functions | Where each part is (and the denominator is nonzero) |
| Compositions $g(f(x, y))$, $f$ continuous, $g$ continuous of one variable | Wherever $g$ is continuous at the value $f(x, y)$ |
So $e^{x - y}\cos(xy)$ is continuous everywhere, $\ln(1 + x^2 + y^2)$ is continuous everywhere, and $\frac{\sin(xy)}{x - y}$ is continuous everywhere except the line $y = x$. The coordinate functions $(x, y) \mapsto x$ and $(x, y) \mapsto y$ are continuous, and everything else is built from them.
Removable discontinuities. If $f$ is undefined at a point but the limit exists there, defining $f$ to equal that limit fills the hole. For example, with $u = x^2 + y^2$, $\frac{\sin(x^2 + y^2)}{x^2 + y^2} = \frac{\sin u}{u} \to 1$ as $(x, y) \to (0, 0)$, because $u \to 0$. Setting the value $1$ at the origin makes it continuous everywhere.
Separate continuity is not enough. A subtle new phenomenon: $f$ can be continuous in $x$ for each fixed $y$, and continuous in $y$ for each fixed $x$, without being continuous. The function
$$g(x, y) = \begin{cases} \frac{xy}{x^2 + y^2} & (x, y) \ne (0, 0) \ 0 & (x, y) = (0, 0) \end{cases}$$
is identically $0$ on both axes, so every horizontal and vertical slice through the origin is continuous. But along $y = x$ it equals $\frac{1}{2}$, so it has no limit at the origin. Checking one variable at a time only checks two directions, and continuity is about all of them.
Why continuity matters. Continuous functions on closed, bounded sets attain a maximum and a minimum (the Extreme Value Theorem, used in Section 12.6), and continuity is the minimum requirement for the derivative ideas in the rest of the chapter.
Examples
- $f(x, y) = \frac{x^2 + y^2}{x^2 - y^2}$ is rational, so it is continuous except where $x^2 = y^2$, the two lines $y = \pm x$.
- $h(x, y) = \begin{cases} \frac{3x^2y}{x^2 + y^2} & (x, y) \ne (0, 0) \ 0 & (x, y) = (0, 0) \end{cases}$ is continuous at the origin, because $\lvert h \rvert \le 3\lvert y \rvert \to 0 = h(0, 0)$. Away from the origin it is rational with a nonzero denominator, so it is continuous everywhere.
- A numerical test for continuity at a point: on a circle of radius $r$ around the point, take the worst difference $\max \lvert f - f(a, b) \rvert$. Continuity means this shrinks to $0$ as $r \to 0$.
import numpy as np
def worst_gap_on_circle(f, value_at_origin, r, n=3600):
"""Largest |f - f(0,0)| over n points on the circle of radius r around the origin."""
th = np.linspace(0, 2 * np.pi, n, endpoint=False)
return np.max(np.abs(f(r * np.cos(th), r * np.sin(th)) - value_at_origin))
f = lambda x, y: np.sin(x**2 + y**2) / (x**2 + y**2) # define f(0, 0) = 1 to fill the hole
g = lambda x, y: x * y / (x**2 + y**2) # define g(0, 0) = 0
for r in [1e-1, 1e-2, 1e-3]:
print(r, f'{worst_gap_on_circle(f, 1.0, r):.2e}', round(worst_gap_on_circle(g, 0.0, r), 6))
# 0.1 1.67e-05 0.5
# 0.01 1.67e-09 0.5 f: the gap shrinks like r^4 / 6, so f is continuous at the origin
# 0.001 1.67e-13 0.5 g: the gap stays at 1/2 (along y = x), so g is not
# g restricted to each axis is identically 0, hence continuous in x alone and in y alone
xs = np.array([-1e-3, -1e-6, 1e-6, 1e-3])
print(g(xs, 0 * xs), g(0 * xs, xs))
# [0. 0. 0. 0.] [0. 0. 0. 0.]
Key Takeaways
- $f$ is continuous at $(a, b)$ when $\lim_{(x, y) \to (a, b)} f(x, y) = f(a, b)$.
- Polynomials are continuous everywhere, rational functions wherever the denominator is nonzero, and sums, products, and compositions of continuous functions are continuous.
- A removable discontinuity is fixed by defining the function to be its limit at that point.
- Continuity in each variable separately does not imply continuity; the limit must work from all directions.
🧪 Practice
- Where is $f(x, y) = \frac{\ln(x + y)}{x^2 + y^2 - 4}$ continuous? Describe the set.
- Can $f(x, y) = \frac{x^4 - y^4}{x^2 + y^2}$ be defined at the origin to make it continuous? If so, with what value?
- Is $f(x, y) = \begin{cases} \frac{x^2y^2}{x^2 + y^2} & (x, y) \ne (0, 0) \ 1 & (x, y) = (0, 0) \end{cases}$ continuous at the origin? What kind of discontinuity, if any, does it have?
- Interview-style: A physics simulation stores a field $f(x, y)$ on a grid and interpolates between grid points bilinearly (separately linear in $x$ and in $y$ within each cell). Is the interpolated function continuous across cell edges, and why does this matter for a particle moving through the field? Hint: on a shared edge, the two neighboring cells use the same two corner values, so they agree. Consider what a jump in the field would do to forces and to an integrator taking small time steps.
Polar Coordinates for Limits
Theory
Paths are good for disproving limits; for proving them we need to control $f$ on whole small circles around the point. Polar coordinates are tailor-made for this. Put
$$x = r\cos\theta, \qquad y = r\sin\theta, \qquad r = \sqrt{x^2 + y^2}.$$
Then "$(x, y) \to (0, 0)$" becomes simply "$r \to 0^+$", for any and every $\theta$. The distance to the origin is now a single variable.
The polar criterion. If you can find a function $g(r)$, not depending on $\theta$, such that
$$\lvert f(r\cos\theta, r\sin\theta) - L \rvert \le g(r) \quad \text{for all } \theta, \qquad g(r) \to 0 \text{ as } r \to 0^+,$$
then $\lim_{(x, y) \to (0, 0)} f(x, y) = L$. This is the Squeeze Theorem with the bound depending only on distance. The usual tools are $\lvert \cos\theta \rvert \le 1$ and $\lvert \sin\theta \rvert \le 1$.
What usually happens. For a ratio of homogeneous polynomials, numerator degree $p$ and denominator degree $q$, the polar form is $r^{p - q}$ times a function of $\theta$ alone. Then:
| Degrees | Polar form | Limit at the origin |
|---|---|---|
| $p > q$, bounded angular part | $r^{p - q} \cdot (\text{bounded in } \theta)$ | $0$ |
| $p = q$ | A function of $\theta$ only | Usually does not exist (different directions, different values) |
| $p < q$ | $r^{p - q} \to \infty$ in most directions | Does not exist (blows up) |
For the "$p > q$" case to work, the angular part must stay bounded. That holds when the denominator is $(x^2 + y^2)^k$, whose polar form is just $r^{2k}$.
The trap: limits for each fixed $\theta$ are not enough. The bound must be uniform in $\theta$. Consider again $f = \frac{x^2y}{x^4 + y^2}$. In polar form,
$$f = \frac{r\cos^2\theta,\sin\theta}{r^2\cos^4\theta + \sin^2\theta}.$$
For each fixed $\theta$ with $\sin\theta \ne 0$, the denominator tends to $\sin^2\theta > 0$ and $f \to 0$. Yet the limit does not exist (previous topic). The catch: when $\theta$ is very small, $\sin^2\theta$ is tiny and the $r^2\cos^4\theta$ term matters. At each radius $r$, the angle $\theta \approx r$ still gives $f \approx \frac{1}{2}$. The bad angles slide toward the axis as $r$ shrinks, exactly following the parabola $y = x^2$.
Other centers. For a limit at $(a, b)$, use $x = a + r\cos\theta$, $y = b + r\sin\theta$.
Examples
- $\lim_{(x, y) \to (0, 0)} \frac{x^3}{x^2 + y^2}$: in polar form this is $\frac{r^3\cos^3\theta}{r^2} = r\cos^3\theta$, and $\lvert r\cos^3\theta \rvert \le r \to 0$. The limit is $0$.
- $\lim_{(x, y) \to (0, 0)} \frac{xy}{\sqrt{x^2 + y^2}}$: polar form $\frac{r^2\cos\theta\sin\theta}{r} = \frac{r}{2}\sin 2\theta$, bounded by $\frac{r}{2}$. The limit is $0$.
- $\lim_{(x, y) \to (0, 0)} (x^2 + y^2)\ln(x^2 + y^2)$: polar form $r^2\ln r^2 = 2r^2\ln r$, which has no $\theta$ at all. One-variable calculus (L'Hopital's Rule, Section 4.5) gives $r^2\ln r \to 0$. The limit is $0$.
- $\lim_{(x, y) \to (0, 0)} \frac{x^2 - y^2}{x^2 + y^2}$: polar form $\cos^2\theta - \sin^2\theta = \cos 2\theta$, which depends on $\theta$ and not on $r$. Different directions give different values, so the limit does not exist.
import numpy as np
import sympy as sp
r, th = sp.symbols('r theta', positive=True)
x, y = r * sp.cos(th), r * sp.sin(th)
for expr in [x**3 / (x**2 + y**2), x * y / sp.sqrt(x**2 + y**2), x**2 * y / (x**4 + y**2)]:
print(sp.simplify(sp.trigsimp(expr))) # the same function written in polar coordinates
# r*cos(theta)**3 bounded by r: limit 0
# r*sin(2*theta)/2 bounded by r/2: limit 0
# r*sin(theta)*cos(theta)**2/(r**2*cos(theta)**4 + sin(theta)**2) no uniform bound: careful
# A theta-dependent limit is not enough: check the WORST theta at each radius
G = sp.lambdify((r, th), x**2 * y / (x**4 + y**2), 'numpy')
thetas = np.linspace(1e-6, np.pi - 1e-6, 200001)
for rv in [1e-1, 1e-2, 1e-3]:
vals = G(rv, thetas)
k = np.argmax(vals)
print(rv, round(vals.max(), 4), f'{thetas[k]:.2e}') # max stays at 1/2, at ever smaller angles
# 0.1 0.5 9.92e-02
# 0.01 0.5 1.00e-02
# 0.001 0.5 1.01e-03 the worst angle is about theta = r: the parabola y = x^2 in disguise
Key Takeaways
- Polar coordinates turn "$(x, y) \to (0, 0)$" into "$r \to 0^+$ for all $\theta$".
- To prove a limit, bound $\lvert f - L \rvert$ by a function of $r$ alone that goes to $0$; $\lvert \cos\theta \rvert, \lvert \sin\theta \rvert \le 1$ are the usual tools.
- If the polar form is a non-constant function of $\theta$ only, the limit does not exist.
- A limit of $0$ for each fixed $\theta$ is not enough: the bound must be uniform in $\theta$.
🧪 Practice
- Use polar coordinates to evaluate $\lim_{(x, y) \to (0, 0)} \frac{x^2y^2}{x^2 + y^2}$ and $\lim_{(x, y) \to (0, 0)} \frac{\sin(x^2 + y^2)}{\sqrt{x^2 + y^2}}$.
- Use polar coordinates to show that $\lim_{(x, y) \to (0, 0)} \frac{x^2 + 3xy}{x^2 + y^2}$ does not exist.
- Evaluate $\lim_{(x, y) \to (1, 2)} \frac{(x - 1)^2(y - 2)}{(x - 1)^2 + (y - 2)^2}$ using shifted polar coordinates.
- Interview-style: You are implementing a function like $\frac{\sin r}{r}$ with $r = \sqrt{x^2 + y^2}$ in a numerical library, and it must be accurate and free of
NaNnear the origin. How would you implement it? Hint: the limit at the origin is $1$, but evaluating the formula at $r = 0$ gives $\frac{0}{0}$. Below a small threshold, switch to the Taylor polynomial $1 - \frac{r^2}{6} + \frac{r^4}{120}$ (Section 9.6), which needs only $r^2 = x^2 + y^2$ and no square root. <a id="123-partial-derivatives"></a>
12.3 Partial Derivatives
This subchapter measures how a function of several variables changes when one input moves and the others are held fixed, and then asks when these one-direction rates combine into a genuine best linear approximation: a tangent plane.
Definition of Partial Derivatives
Theory
How does the pressure of a gas change when you heat it? It depends on whether the container is rigid. To get a clean answer, you hold the volume fixed and change only the temperature. That is the idea of a partial derivative: change one input, freeze all the others, and measure the rate of change of the output. With the other variables frozen, $f$ becomes an ordinary function of one variable, and Chapter 3 applies.
Definition. The partial derivatives of $f(x, y)$ at $(a, b)$ are
$$f_x(a, b) = \lim_{h \to 0} \frac{f(a + h, b) - f(a, b)}{h}, \qquad f_y(a, b) = \lim_{h \to 0} \frac{f(a, b + h) - f(a, b)}{h},$$
when these limits exist. In the first, only $x$ moves; in the second, only $y$. For $f(x_1, \ldots, x_n)$ there is one partial derivative per variable, each nudging one coordinate.
Notation. Several notations are in use, and you should recognize all of them:
| Partial in $x$ | Partial in $y$ | Where you see it |
|---|---|---|
| $f_x$ | $f_y$ | Compact; common in calculus texts |
| $\frac{\partial f}{\partial x}$ | $\frac{\partial f}{\partial y}$ | Leibniz style; physics and engineering |
| $\frac{\partial z}{\partial x}$ | $\frac{\partial z}{\partial y}$ | When $z = f(x, y)$ |
| $D_1 f$, $\partial_1 f$ | $D_2 f$, $\partial_2 f$ | By position; useful when variable names are ambiguous |
The curly $\partial$ (read "partial") is a reminder that other variables exist and are being held constant. Never write $\frac{df}{dx}$ for a function of several variables.
How to compute. Treat every other variable as a constant and differentiate with the usual rules (Sections 3.3 and 3.4). For $f = x^2y^3$, the factor $y^3$ is a constant when differentiating in $x$, just as the $5$ is in $5x^2$:
$$f_x = 2xy^3, \qquad f_y = 3x^2y^2.$$
Units and meaning. $f_x$ has units of "output per unit of $x$". It answers: if $x$ increases by one small unit and nothing else changes, by how much does $f$ change? So $f(a + \Delta x, b) \approx f(a, b) + f_x(a, b),\Delta x$.
Examples
- $f(x, y) = x^2y^3 + \sin(xy)$. Holding $y$ fixed, the chain rule on $\sin(xy)$ gives a factor $y$; holding $x$ fixed, a factor $x$:
$$f_x = 2xy^3 + y\cos(xy), \qquad f_y = 3x^2y^2 + x\cos(xy).$$
At $(1, 0)$: $f_x = 0$ and $f_y = 1$.
-
Three variables: $g(x, y, z) = xe^{yz}$ gives $g_x = e^{yz}$, $g_y = xze^{yz}$, $g_z = xye^{yz}$.
-
The ideal gas law $P = \frac{nRT}{V}$ (1 mole, $R = 8.314$ J/(mol K), volume in liters, pressure in kPa). At $T = 300$ K and $V = 25$ L:
- $\frac{\partial P}{\partial T} = \frac{nR}{V} \approx 0.333$ kPa per kelvin: heating a rigid container by 1 K raises the pressure about $0.33$ kPa.
- $\frac{\partial P}{\partial V} = -\frac{nRT}{V^2} \approx -3.99$ kPa per liter: expanding at fixed temperature lowers the pressure.
The definition is also how a computer estimates a partial derivative when there is no formula: nudge one input and divide.
import numpy as np
import sympy as sp
x, y = sp.symbols('x y', real=True)
f = x**2 * y**3 + sp.sin(x * y)
fx, fy = sp.diff(f, x), sp.diff(f, y) # hold the other variable constant
print(fx, '|', fy)
# 2*x*y**3 + y*cos(x*y) | 3*x**2*y**2 + x*cos(x*y)
print(fx.subs({x: 1, y: 0}), fy.subs({x: 1, y: 0}))
# 0 1
# Difference quotients: nudge ONE input, divide by the nudge
F = sp.lambdify((x, y), f)
a, b, h = 1.0, 2.0, 1e-6
print(round((F(a + h, b) - F(a, b)) / h, 4), round(float(fx.subs({x: a, y: b})), 4))
# 15.1677 15.1677
print(round((F(a, b + h) - F(a, b)) / h, 4), round(float(fy.subs({x: a, y: b})), 4))
# 11.5839 11.5839
# Ideal gas: P = nRT / V. Sensitivities of pressure (kPa) for 1 mol at T = 300 K, V = 25 L
n, R = 1.0, 8.314
P = lambda T, V: n * R * T / V # kPa when V is in liters
T0, V0 = 300.0, 25.0
print(round(P(T0, V0), 3), round(n * R / V0, 4), round(-n * R * T0 / V0**2, 4)) # P, dP/dT, dP/dV
# 99.768 0.3326 -3.9907
Key Takeaways
- $f_x$ is the rate of change of $f$ when only $x$ changes; all other variables are held constant.
- Compute it with one-variable rules, treating the other variables as constants.
- Notations $f_x$, $\frac{\partial f}{\partial x}$, $D_1 f$ all mean the same thing; the $\partial$ signals that other variables exist.
- Numerically, nudge one input: $f_x(a, b) \approx \frac{f(a + h, b) - f(a, b)}{h}$.
🧪 Practice
- Find $f_x$ and $f_y$ for $f(x, y) = \frac{x - y}{x + y}$ and for $f(x, y) = y^x$.
- Find all first partial derivatives of $w = \ln(x + 2y + 3z) + z\arctan(xy)$.
- Use the limit definition to find $f_x(0, 0)$ and $f_y(0, 0)$ for $f(x, y) = \sqrt[3]{x^3 + y^3}$.
- Interview-style: A function in a model has 1000 inputs, and you need all its partial derivatives at one point, but you only have code that evaluates the function. How many function evaluations does a forward-difference estimate take, and how would you pick the step size $h$? Hint: $1000 + 1$ evaluations. Too large an $h$ adds truncation error of order $h$; too small an $h$ amplifies rounding error of order $\frac{\epsilon_{\text{mach}}}{h}$. A step near $\sqrt{\epsilon_{\text{mach}}} \approx 10^{-8}$ (scaled to the size of the input) balances them. Automatic differentiation avoids both problems.
Geometric Interpretation
Theory
The graph of $z = f(x, y)$ is a surface. Slicing it with the vertical plane $y = b$ gives a curve: the profile you would see walking east-west through the point. On that curve, $y$ is frozen at $b$ and only $x$ varies, so the curve is the graph of the one-variable function $x \mapsto f(x, b)$. Its slope at $x = a$ is exactly $f_x(a, b)$.
| Partial | Slice | Curve | Partial derivative is |
|---|---|---|---|
| $f_x(a, b)$ | Plane $y = b$ | $C_1: z = f(x, b)$ | Slope of $C_1$ at $P$: rise in $z$ per unit step east |
| $f_y(a, b)$ | Plane $x = a$ | $C_2: z = f(a, y)$ | Slope of $C_2$ at $P$: rise in $z$ per unit step north |
So at a point $P = (a, b, f(a, b))$ on the surface, the two partial derivatives are the slopes of the surface in the two coordinate directions. A hiker standing at $P$ with $f_x = 2$ and $f_y = -1$ climbs $2$ m per meter walking east and descends $1$ m per meter walking north.
Tangent lines to the traces. The tangent line to $C_1$ at $P$ moves $1$ unit in $x$ and $f_x$ units in $z$, and nothing in $y$, so its direction is $\langle 1, 0, f_x(a, b) \rangle$. Similarly, the tangent line to $C_2$ has direction $\langle 0, 1, f_y(a, b) \rangle$:
$$\ell_1(t) = \langle a + t,\ b,\ f(a, b) + f_x t \rangle, \qquad \ell_2(t) = \langle a,\ b + t,\ f(a, b) + f_y t \rangle.$$
These two lines both pass through $P$ and touch the surface. If the surface has a tangent plane at $P$, it must contain both lines, so its normal vector is their cross product (Section 11.1):
$$\langle 1, 0, f_x \rangle \times \langle 0, 1, f_y \rangle = \langle -f_x,\ -f_y,\ 1 \rangle.$$
This is the bridge to tangent planes later in this subchapter.
What partials do not see. Two partial derivatives describe only two directions. A function can do something completely different along the diagonal (Section 12.2 examples). Section 12.5 handles all directions with the directional derivative.
z Surface z = f(x, y) cut by two vertical planes
| through P:
| C2 (x = a)
| \ slope f_y plane y = b -> curve C1, slope f_x
| \ / plane x = a -> curve C2, slope f_y
| _.--- P ----._
| .-' / \ '-. C1 (y = b)
| / / slope f_x
|/ \
+--------------------------------- y
/ . (a, b)
/
x
Examples
Take $f(x, y) = 4 - x^2 - 2y^2$ at $(1, 1)$, so $P = (1, 1, 1)$.
- $f_x = -2x = -2$: walking in the $+x$ direction, the surface drops $2$ units per unit.
- $f_y = -4y = -4$: walking in the $+y$ direction, it drops twice as fast, matching the steeper $y$-profile seen in Section 12.1.
- Tangent lines: $\ell_1(t) = \langle 1 + t, 1, 1 - 2t \rangle$ and $\ell_2(t) = \langle 1, 1 + t, 1 - 4t \rangle$.
- The cross product of the directions gives the normal $\langle 2, 4, 1 \rangle$.
import sympy as sp
x, y, t = sp.symbols('x y t', real=True)
f = 4 - x**2 - 2 * y**2
a, b = 1, 1
c = f.subs({x: a, y: b})
fx = sp.diff(f, x).subs({x: a, y: b})
fy = sp.diff(f, y).subs({x: a, y: b})
print(c, fx, fy)
# 1 -2 -4
# Tangent lines to the two traces through P = (1, 1, 1)
L1 = sp.Matrix([a + t, b, c + fx * t]) # in the plane y = 1, direction <1, 0, f_x>
L2 = sp.Matrix([a, b + t, c + fy * t]) # in the plane x = 1, direction <0, 1, f_y>
print(L1.T, L2.T)
# Matrix([[t + 1, 1, 1 - 2*t]]) Matrix([[1, t + 1, 1 - 4*t]])
# Each line touches its trace: the gap between curve and line is O(t^2)
trace1 = f.subs({x: a + t, y: b})
print(sp.expand(trace1 - L1[2]), '|', sp.expand(f.subs({x: a, y: b + t}) - L2[2]))
# -t**2 | -2*t**2
# The two direction vectors span the tangent plane; their cross product is its normal
n = sp.Matrix([1, 0, fx]).cross(sp.Matrix([0, 1, fy]))
print(n.T)
# Matrix([[2, 4, 1]])
Key Takeaways
- $f_x(a, b)$ is the slope of the surface's trace in the plane $y = b$; $f_y(a, b)$ is the slope of its trace in the plane $x = a$.
- The tangent lines to the traces have directions $\langle 1, 0, f_x \rangle$ and $\langle 0, 1, f_y \rangle$.
- Their cross product $\langle -f_x, -f_y, 1 \rangle$ is the normal of the tangent plane, when one exists.
- Partial derivatives only see two directions; other directions need the directional derivative.
🧪 Practice
- For $f(x, y) = \sqrt{25 - x^2 - y^2}$ at $(3, 0)$, find $f_x$ and $f_y$ and interpret each as a slope on the hemisphere.
- Find parametric equations of the tangent line to the curve of intersection of $z = x^2 + xy^3$ with the plane $x = 2$, at the point $(2, 1, 6)$.
- On a contour map, the level curves near a point $P$ run exactly north-south. What can you say about $f_x$ and $f_y$ at $P$?
- Interview-style: A robot reads a terrain height map and only computes slopes in the $x$ and $y$ directions. It wants to know whether a straight path at $30^\circ$ from the $x$-axis is too steep. Can it answer with what it has? Hint: if the terrain is differentiable, the slope in direction $\langle \cos\theta, \sin\theta \rangle$ is $f_x\cos\theta + f_y\sin\theta$ (Section 12.5). Without differentiability, two axis slopes cannot determine the others.
Higher-Order Partials and Clairaut's Theorem
Theory
The partial derivatives $f_x$ and $f_y$ are themselves functions of $(x, y)$, so we can differentiate them again. For a function of two variables there are four second-order partial derivatives:
| Notation | Leibniz | Meaning |
|---|---|---|
| $f_{xx} = (f_x)_x$ | $\frac{\partial^2 f}{\partial x^2}$ | How the $x$-slope changes moving in $x$: concavity of the $x$-profile |
| $f_{yy} = (f_y)_y$ | $\frac{\partial^2 f}{\partial y^2}$ | Concavity of the $y$-profile |
| $f_{xy} = (f_x)_y$ | $\frac{\partial^2 f}{\partial y,\partial x}$ | How the $x$-slope changes moving in $y$ |
| $f_{yx} = (f_y)_x$ | $\frac{\partial^2 f}{\partial x,\partial y}$ | How the $y$-slope changes moving in $x$ |
Watch the order: in subscript notation, $f_{xy}$ means "first $x$, then $y$" (read left to right). In Leibniz notation, $\frac{\partial^2 f}{\partial y,\partial x}$ means the same thing (read right to left, innermost first).
The mixed partials $f_{xy}$ and $f_{yx}$ measure twist. On a surface like $z = xy$, the $x$-slope is $y$: it is $0$ along the $x$-axis but grows as you move north. That change of slope with position is $f_{xy} = 1$.
Clairaut's Theorem. In practice, the two mixed partials almost always agree:
If $f_{xy}$ and $f_{yx}$ are both continuous on an open disk around $(a, b)$, then $f_{xy}(a, b) = f_{yx}(a, b)$.
So for all functions built from polynomials, exponentials, trigonometric functions, and so on (away from points where they break), the order of differentiation does not matter. The same holds for higher orders: $f_{xyy} = f_{yxy} = f_{yyx}$ when they are continuous. A useful consequence: a function of $n$ variables has only $\frac{n(n + 1)}{2}$ distinct second partials, not $n^2$, because they form a symmetric matrix (the Hessian, Section 12.6).
When Clairaut fails. The continuity hypothesis is not decoration. For
$$g(x, y) = \begin{cases} \frac{xy(x^2 - y^2)}{x^2 + y^2} & (x, y) \ne (0, 0) \ 0 & (x, y) = (0, 0) \end{cases}$$
the limit definition gives $g_x(0, y) = -y$ and $g_y(x, 0) = x$, so $g_{xy}(0, 0) = -1$ but $g_{yx}(0, 0) = 1$. The mixed partials exist but are not continuous at the origin.
Partial differential equations. Many laws of physics are equations relating partial derivatives:
| Equation | Form | Describes |
|---|---|---|
| Laplace | $u_{xx} + u_{yy} = 0$ | Steady-state temperature, electrostatic potential, ideal fluid flow |
| Heat (diffusion) | $u_t = \alpha,u_{xx}$ | Temperature spreading along a rod over time |
| Wave | $u_{tt} = c^2 u_{xx}$ | Vibrating strings, sound, light |
Checking whether a given function solves one of these is a direct exercise in second partials.
Examples
- $f(x, y) = x^3y^2 - 2xy^4$:
- $f_x = 3x^2y^2 - 2y^4$, $f_y = 2x^3y - 8xy^3$.
- $f_{xx} = 6xy^2$, $f_{yy} = 2x^3 - 24xy^2$.
- $f_{xy} = 6x^2y - 8y^3$ and $f_{yx} = 6x^2y - 8y^3$. Equal, as Clairaut predicts.
- $u = \ln(x^2 + y^2)$ satisfies Laplace's equation: $u_{xx} = \frac{2(y^2 - x^2)}{(x^2 + y^2)^2}$ and $u_{yy} = \frac{2(x^2 - y^2)}{(x^2 + y^2)^2}$, which add to $0$.
- $u = e^{-k^2t}\sin(kx)$ satisfies the heat equation $u_t = u_{xx}$: both sides equal $-k^2e^{-k^2t}\sin(kx)$. Higher frequencies $k$ decay faster, which is why sharp temperature differences smooth out quickly.
import sympy as sp
x, y, h, k, t = sp.symbols('x y h k t', real=True)
f = x**3 * y**2 + sp.exp(x * y) * sp.sin(y)
fxy = sp.diff(f, x, y) # differentiate in x, then in y
fyx = sp.diff(f, y, x)
print(sp.simplify(fxy - fyx)) # Clairaut: the mixed partials agree
# 0
# Counterexample: g = xy(x^2 - y^2)/(x^2 + y^2), g(0, 0) = 0. Work from the limit definition.
def g(X, Y):
return sp.Integer(0) if (X, Y) == (0, 0) else X * Y * (X**2 - Y**2) / (X**2 + Y**2)
gx_on_y_axis = sp.limit((g(h, y) - g(0, y)) / h, h, 0) # g_x(0, y)
gy_on_x_axis = sp.limit((g(x, k) - g(x, 0)) / k, k, 0) # g_y(x, 0)
print(gx_on_y_axis, gy_on_x_axis)
# -y x
print(sp.diff(gx_on_y_axis, y), sp.diff(gy_on_x_axis, x)) # g_xy(0, 0) and g_yx(0, 0) differ
# -1 1
# Classic PDEs checked with second partials
u = sp.log(x**2 + y**2)
print(sp.simplify(sp.diff(u, x, 2) + sp.diff(u, y, 2))) # Laplace: u_xx + u_yy = 0
# 0
w = sp.exp(-k**2 * t) * sp.sin(k * x)
print(sp.simplify(sp.diff(w, t) - sp.diff(w, x, 2))) # heat: u_t = u_xx
# 0
Key Takeaways
- Second partials: $f_{xx}$ and $f_{yy}$ measure concavity along the axes, and the mixed partials $f_{xy}$, $f_{yx}$ measure twist.
- $f_{xy}$ means differentiate in $x$ first, then $y$; in Leibniz notation, $\frac{\partial^2 f}{\partial y,\partial x}$.
- Clairaut's Theorem: if the mixed partials are continuous near a point, they are equal there. This covers almost every function you will meet.
- Partial differential equations (Laplace, heat, wave) are relations among partial derivatives; verifying a solution is direct differentiation.
🧪 Practice
- Find all second partial derivatives of $f(x, y) = e^{x}\cos y + x^2y$ and verify that $f_{xy} = f_{yx}$.
- For $f(x, y, z) = x^2y^3z^4$, find $f_{xyz}$ and $f_{zyx}$. How many distinct third-order partials does a smooth function of three variables have?
- Show that $u(x, t) = \sin(x - ct) + (x + ct)^2$ satisfies the wave equation $u_{tt} = c^2u_{xx}$.
- Interview-style: A second-order optimizer needs the Hessian matrix of a function of $n = 10{,}000$ variables. Using Clairaut's Theorem, how many distinct entries must be stored, and is that still practical? Hint: the matrix is symmetric, so $\frac{n(n + 1)}{2} \approx 5 \times 10^7$ entries. That fits in memory but is expensive to compute; methods that only need Hessian-vector products (or approximate the Hessian, like L-BFGS) avoid forming it.
Differentiability and Tangent Planes
Theory
In one variable, having a derivative at $a$ means the graph has a tangent line there, and differentiability implies continuity (Section 3.2). In several variables, the partial derivatives alone do not give this. The function $\frac{xy}{x^2 + y^2}$ (with value $0$ at the origin) has $f_x(0, 0) = f_y(0, 0) = 0$ because it is $0$ on both axes, yet it is not even continuous there. The partials only look along two lines; a surface can misbehave in between.
So we need a stronger notion. The right idea, which generalizes the tangent line, is: near the point, the function is well approximated by a plane.
Definition. $f$ is differentiable at $(a, b)$ if there is a linear function
$$L(x, y) = f(a, b) + f_x(a, b)(x - a) + f_y(a, b)(y - b)$$
such that the error is small compared to the distance from $(a, b)$:
$$\lim_{(x, y) \to (a, b)} \frac{f(x, y) - L(x, y)}{\sqrt{(x - a)^2 + (y - b)^2}} = 0.$$
Every continuous function has error going to $0$; differentiability demands that the error shrink faster than the distance, from every direction. Zoom in on the graph near a point of differentiability and it looks more and more like a flat plane.
The tangent plane. When $f$ is differentiable at $(a, b)$, the graph of $L$ is the tangent plane to the surface at $P = (a, b, f(a, b))$:
$$\boxed{z = f(a, b) + f_x(a, b)(x - a) + f_y(a, b)(y - b)}$$
Its normal vector is $\langle f_x, f_y, -1 \rangle$, or equivalently $\langle -f_x, -f_y, 1 \rangle$, matching the cross product from the geometric interpretation topic.
How to check differentiability in practice.
| Fact | Consequence |
|---|---|
| If $f_x$ and $f_y$ exist near $(a, b)$ and are continuous at $(a, b)$, then $f$ is differentiable at $(a, b)$ | The everyday test: polynomials, exponentials, trig, and their combinations are differentiable wherever their partials are continuous |
| Differentiable implies continuous | A discontinuous function is never differentiable |
| Differentiable implies all directional derivatives exist (Section 12.5) | |
| Partials exist does not imply differentiable, or even continuous | Always check the stronger condition at suspicious points |
A useful hierarchy:
$$\text{continuous partials} \implies \text{differentiable} \implies \text{continuous}, \qquad \text{differentiable} \implies \text{partials exist}.$$
None of these arrows can be reversed.
A continuous, non-differentiable example. $g(x, y) = \sqrt{\lvert xy \rvert}$ is continuous and is $0$ on both axes, so $g_x(0, 0) = g_y(0, 0) = 0$ and the candidate tangent plane is $z = 0$. But along the diagonal, $g(t, t) = \lvert t \rvert$, and the error divided by the distance $\sqrt{2}\lvert t \rvert$ is $\frac{1}{\sqrt{2}}$, not $0$. The surface has a crease along the diagonals, like a folded sheet, and no plane fits it.
Examples
Tangent plane to $z = x^2 + y^2$ at $(1, 2, 5)$: $f_x = 2x = 2$ and $f_y = 2y = 4$, so
$$z = 5 + 2(x - 1) + 4(y - 2), \qquad \text{that is,} \qquad 2x + 4y - z = 5.$$
At $(1 + d, 2 + d)$ the surface and plane differ by exactly $2d^2$: the gap shrinks much faster than the distance $\sqrt{2},d$.
The code below measures the ratio $\frac{\lvert f - L \rvert}{\text{distance}}$ along several directions. For a differentiable function it goes to $0$ in every direction. For $\sqrt{\lvert xy \rvert}$ it goes to $0$ along the axis but stays at $0.707$ along the diagonal.
import numpy as np
def linearization_ratio(f, fx, fy, a, b, direction, hs=(1e-1, 1e-2, 1e-3, 1e-4)):
"""|f(P + h u) - L(P + h u)| / |h u| for shrinking steps along one direction u."""
u = np.asarray(direction, float) / np.linalg.norm(direction)
out = []
for s in hs:
dx, dy = s * u
L = f(a, b) + fx * dx + fy * dy # the candidate tangent plane
out.append(abs(f(a + dx, b + dy) - L) / s)
return ' '.join(f'{v:.2e}' for v in out)
# Smooth: f = x^2 y + e^y at (1, 0) has f_x = 0, f_y = 2. Error/distance -> 0 in every direction
f = lambda x, y: x**2 * y + np.exp(y)
for d in [(1, 0), (1, 1), (-2, 1)]:
print('smooth ', d, linearization_ratio(f, 0.0, 2.0, 1.0, 0.0, d))
# smooth (1, 0) 0.00e+00 0.00e+00 0.00e+00 0.00e+00
# smooth (1, 1) 1.29e-01 1.25e-02 1.25e-03 1.25e-04 shrinks like h: differentiable
# smooth (-2, 1) 6.63e-02 6.96e-03 7.00e-04 7.00e-05
# g = sqrt(|xy|): continuous, g_x(0,0) = g_y(0,0) = 0, yet NOT differentiable at the origin
g = lambda x, y: np.sqrt(np.abs(x * y))
for d in [(1, 0), (1, 1)]:
print('sqrt|xy| ', d, linearization_ratio(g, 0.0, 0.0, 0.0, 0.0, d))
# sqrt|xy| (1, 0) 0.00e+00 0.00e+00 0.00e+00 0.00e+00
# sqrt|xy| (1, 1) 7.07e-01 7.07e-01 7.07e-01 7.07e-01 never shrinks: no tangent plane
# Tangent plane to z = x^2 + y^2 at (1, 2): z = 5 + 2(x - 1) + 4(y - 2)
z = lambda x, y: x**2 + y**2
T = lambda x, y: 5 + 2 * (x - 1) + 4 * (y - 2)
for d in [0.5, 0.1, 0.01]:
zz, tt = z(1 + d, 2 + d), T(1 + d, 2 + d)
print(d, round(zz, 6), round(tt, 6), round(zz - tt, 6)) # gap = 2 d^2, shrinks faster than d
# 0.5 8.5 8.0 0.5
# 0.1 5.62 5.6 0.02
# 0.01 5.0602 5.06 0.0002
Key Takeaways
- Differentiable at $(a, b)$ means $f(x, y) - L(x, y)$ is small compared to the distance to $(a, b)$, from every direction.
- The tangent plane is $z = f(a, b) + f_x(a, b)(x - a) + f_y(a, b)(y - b)$, with normal $\langle f_x, f_y, -1 \rangle$.
- Continuous partial derivatives guarantee differentiability; existence of partials alone does not even guarantee continuity.
- Differentiable implies continuous, so a discontinuous function is never differentiable.
🧪 Practice
- Find the tangent plane to $z = e^{x}\ln y$ at the point $(0, 1, 0)$.
- Find the tangent plane to $z = \sqrt{x^2 + y^2}$ at $(3, 4, 5)$. Is the cone differentiable at the origin? Explain using the definition.
- Explain why $f(x, y) = x^2 + \sin(xy) + e^{y}$ is differentiable everywhere without computing any limits.
- Interview-style: Neural networks often use the activation $\max(0, u)$, which is not differentiable at $0$, yet training uses gradients anyway. Why does this rarely cause trouble, and what do frameworks do at the kink? Hint: the set of inputs landing exactly on a kink has probability essentially zero for continuous data, and frameworks pick a fixed value (a "subgradient", usually $0$) there. Compare with $\sqrt{\lvert xy \rvert}$, whose crease lies along whole lines rather than isolated points.
Linear Approximation and Total Differential
Theory
The tangent plane is the multivariable version of the tangent line approximation from Section 4.5: near a point where $f$ is differentiable,
$$f(x, y) \approx L(x, y) = f(a, b) + f_x(a, b)(x - a) + f_y(a, b)(y - b).$$
The function $L$ is the linearization of $f$ at $(a, b)$. In words: the change in output is approximately the sum of the changes caused by each input separately. Each input contributes "its sensitivity times its change", and the contributions add.
The total differential. Writing the changes as differentials $dx = \Delta x$, $dy = \Delta y$, the approximate change in $z = f(x, y)$ is
$$dz = f_x,dx + f_y,dy = \frac{\partial z}{\partial x}dx + \frac{\partial z}{\partial y}dy.$$
Compare with the actual change $\Delta z = f(a + dx, b + dy) - f(a, b)$. Differentiability means $\Delta z - dz$ is small compared to $\sqrt{dx^2 + dy^2}$. For $n$ variables: $dw = f_{x_1}dx_1 + \cdots + f_{x_n}dx_n$.
Error propagation. The main practical use. If you measure $x$ and $y$ with possible errors $\lvert dx \rvert \le \delta_x$ and $\lvert dy \rvert \le \delta_y$, the computed $z$ has error about
$$\lvert dz \rvert \le \lvert f_x \rvert,\delta_x + \lvert f_y \rvert,\delta_y \quad \text{(worst case)}, \qquad \sigma_z \approx \sqrt{(f_x\sigma_x)^2 + (f_y\sigma_y)^2} \quad \text{(independent random errors)}.$$
The worst case assumes all errors push the same way. The root-sum-square formula, standard in laboratory work, assumes independent random errors that partly cancel, and gives a smaller, more typical figure.
Relative error and products. For a product of powers $z = Cx^py^q$, the differential becomes especially clean after dividing by $z$:
$$\frac{dz}{z} = p\frac{dx}{x} + q\frac{dy}{y}.$$
Percentage errors are multiplied by the exponents and added. Squaring a measurement doubles its percentage error.
Examples
- Estimate $\sqrt{3.02^2 + 3.97^2}$. Let $f = \sqrt{x^2 + y^2}$ near $(3, 4)$, where $f = 5$, $f_x = \frac{x}{f} = 0.6$, and $f_y = \frac{y}{f} = 0.8$:
$$f(3.02, 3.97) \approx 5 + 0.6(0.02) + 0.8(-0.03) = 4.988.$$
The true value is $4.98812$.
- A cylinder has radius $r = 10 \pm 0.1$ cm and height $h = 25 \pm 0.1$ cm. With $V = \pi r^2h$:
$$dV = 2\pi rh,dr + \pi r^2,dh = 2\pi(250)(0.1) + \pi(100)(0.1) = 60\pi \approx 188.5 \text{ cm}^3,$$
about $2.4%$ of $V = 2500\pi \approx 7854$ cm³. The relative form gives the same thing at a glance: $\frac{dV}{V} = 2\frac{dr}{r} + \frac{dh}{h} = 2(1%) + 0.4% = 2.4%$. The radius error dominates, because it is squared and because $1%$ is a bigger fraction than $0.4%$. To improve the measurement, measure the radius better.
import numpy as np
# Linear approximation of f(x, y) = sqrt(x^2 + y^2) near (3, 4)
f = lambda x, y: np.hypot(x, y)
L = lambda x, y: 5 + 0.6 * (x - 3) + 0.8 * (y - 4) # f_x = x/f = 3/5, f_y = y/f = 4/5
print(round(f(3.02, 3.97), 6), round(L(3.02, 3.97), 6))
# 4.988116 4.988
# Error propagation: cylinder V = pi r^2 h with r = 10 +/- 0.1 cm, h = 25 +/- 0.1 cm
r, h, dr, dh = 10.0, 25.0, 0.1, 0.1
V = np.pi * r**2 * h
dV = 2 * np.pi * r * h * dr + np.pi * r**2 * dh # total differential = worst-case linear bound
true_worst = np.pi * (r + dr)**2 * (h + dh) - V
print(round(V, 1), round(dV, 2), round(true_worst, 2), f'{dV / V:.2%}')
# 7854.0 188.5 189.91 2.40%
# Monte Carlo: independent errors rarely all push the same way
rng = np.random.default_rng(1)
R = r + rng.uniform(-dr, dr, 1_000_000)
H = h + rng.uniform(-dh, dh, 1_000_000)
dev = np.pi * R**2 * H - V
rss = np.sqrt((2 * np.pi * r * h * dr)**2 + (np.pi * r**2 * dh)**2) # root-sum-square combination
print(round(np.abs(dev).max(), 2), round(np.percentile(np.abs(dev), 95), 2), round(rss, 2))
# 189.82 157.19 160.19 worst case is rare; 95% of errors are within about the RSS figure
Key Takeaways
- Linearization: $f(x, y) \approx f(a, b) + f_x,(x - a) + f_y,(y - b)$ near a point of differentiability.
- Total differential: $dz = f_x,dx + f_y,dy$. Each input contributes its sensitivity times its change, and the contributions add.
- Worst-case error bound: $\lvert f_x \rvert\delta_x + \lvert f_y \rvert\delta_y$; for independent random errors, combine in quadrature.
- For products of powers, relative errors scale by the exponents and add: $\frac{dz}{z} = p\frac{dx}{x} + q\frac{dy}{y}$.
🧪 Practice
- Use a linear approximation to estimate $(1.98)^3\sqrt{4.03}$, and compare with the exact value.
- The period of a pendulum is $T = 2\pi\sqrt{\frac{L}{g}}$. If $L$ is measured with a $1%$ error and $g$ is known with a $0.2%$ error, estimate the maximum percentage error in $T$.
- The sides of a rectangular box are measured as $75$, $60$, and $40$ cm, each with a possible error of $0.2$ cm. Use differentials to estimate the maximum error in the computed volume and in the computed surface area.
- Interview-style: You are building a sensor-fusion system that computes a robot's distance $d = \sqrt{x^2 + y^2}$ from noisy position estimates with standard deviations $\sigma_x$ and $\sigma_y$. How would you estimate $\sigma_d$, and when does this linear estimate break down? Hint: $\sigma_d \approx \sqrt{(\frac{x}{d}\sigma_x)^2 + (\frac{y}{d}\sigma_y)^2}$. Near $d = 0$ the function is not differentiable (a cone tip), and the linearization is meaningless there; a Monte Carlo or unscented approach is safer. <a id="124-chain-rules-and-implicit-functions"></a>
12.4 Chain Rules and Implicit Functions
This subchapter shows how rates of change pass through chains of dependent variables, organizes the bookkeeping with tree diagrams, and uses the result to differentiate equations that define functions implicitly.
Chain Rule for Several Variables
Theory
A drone flies through a room where the temperature is $T(x, y)$. Its position changes with time, $x = x(t)$, $y = y(t)$, so the temperature it feels, $T(x(t), y(t))$, also changes with time. How fast? Two things are happening at once: the drone moves in $x$, which changes $T$ at rate $T_x$ per unit of $x$, and it moves in $y$, which changes $T$ at rate $T_y$ per unit of $y$. For a short time $dt$, the linear approximation (Section 12.3) says
$$dT \approx T_x,dx + T_y,dy = T_x,x'(t),dt + T_y,y'(t),dt.$$
Dividing by $dt$ gives the chain rule. The rates through each input simply add.
Case 1: one parameter. If $z = f(x, y)$ is differentiable, and $x = x(t)$, $y = y(t)$ are differentiable, then
$$\frac{dz}{dt} = \frac{\partial z}{\partial x}\frac{dx}{dt} + \frac{\partial z}{\partial y}\frac{dy}{dt}.$$
Case 2: several parameters. If $x = x(s, t)$ and $y = y(s, t)$, then $z$ depends on $s$ and $t$, and each partial is computed the same way, with the other parameter held fixed:
$$\frac{\partial z}{\partial s} = \frac{\partial z}{\partial x}\frac{\partial x}{\partial s} + \frac{\partial z}{\partial y}\frac{\partial y}{\partial s}, \qquad \frac{\partial z}{\partial t} = \frac{\partial z}{\partial x}\frac{\partial x}{\partial t} + \frac{\partial z}{\partial y}\frac{\partial y}{\partial t}.$$
General form. With $n$ intermediate variables $x_1, \ldots, x_n$, each depending on parameters $t_1, \ldots, t_m$:
$$\frac{\partial z}{\partial t_j} = \sum_{i=1}^{n} \frac{\partial z}{\partial x_i}\frac{\partial x_i}{\partial t_j}.$$
Matrix form. This sum is a matrix product. Put the partials of the outer function in a row and the partials of the inner functions in a matrix, the Jacobian matrix of the inner map:
$$\begin{bmatrix} z_s & z_t \end{bmatrix} = \begin{bmatrix} z_x & z_y \end{bmatrix} \begin{bmatrix} x_s & x_t \ y_s & y_t \end{bmatrix}.$$
This is the true generalization of the one-variable rule $(f \circ g)' = f'(g),g'$: the derivative of a composition is the product of the derivatives, where derivatives are now matrices. It is also exactly what backpropagation computes in neural networks, layer by layer.
Hypothesis. The outer function must be differentiable, not merely have partial derivatives. For the non-differentiable $f = \frac{xy}{x^2 + y^2}$ (Section 12.3) with $x = t$, $y = t$, the formula predicts $f_x\cdot 1 + f_y\cdot 1 = 0$ at the origin, but $f(t, t) = \frac{1}{2}$ for $t \ne 0$ and $0$ at $t = 0$, which is not even continuous.
Examples
- $z = x^2y + 3xy^4$ with $x = \sin 2t$, $y = \cos t$. Then
$$\frac{dz}{dt} = (2xy + 3y^4)(2\cos 2t) + (x^2 + 12xy^3)(-\sin t).$$
At $t = 0$: $x = 0$, $y = 1$, so $\frac{dz}{dt} = (3)(2) + (0)(0) = 6$.
- Polar coordinates. If $u = f(x, y)$ and $x = r\cos\theta$, $y = r\sin\theta$, then
$$u_r = u_x\cos\theta + u_y\sin\theta, \qquad u_\theta = -u_x,r\sin\theta + u_y,r\cos\theta.$$
These formulas convert derivatives between coordinate systems and are used to rewrite PDEs like Laplace's equation in polar form.
- $z = e^x\sin y$ with $x = st^2$, $y = s^2t$, at $(s, t) = (1, 2)$: the matrix form gives $[z_s, z_t] \approx [107.70, 175.86]$, and finite differences agree.
import numpy as np
import sympy as sp
t, s = sp.symbols('t s', real=True)
X, Y = sp.symbols('x y', real=True)
f = X**2 * Y + 3 * X * Y**4
# Case 1: one parameter. x = sin 2t, y = cos t
xt, yt = sp.sin(2 * t), sp.cos(t)
chain = (sp.diff(f, X) * sp.diff(xt, t) + sp.diff(f, Y) * sp.diff(yt, t)).subs({X: xt, Y: yt})
direct = sp.diff(f.subs({X: xt, Y: yt}), t) # substitute first, then differentiate
print(sp.simplify(chain - direct), chain.subs(t, 0))
# 0 6
# Case 2: two parameters. z = e^x sin y with x = s t^2, y = s^2 t
g = sp.exp(X) * sp.sin(Y)
xs, ys = s * t**2, s**2 * t
z_s = (sp.diff(g, X) * sp.diff(xs, s) + sp.diff(g, Y) * sp.diff(ys, s)).subs({X: xs, Y: ys})
print(sp.simplify(z_s - sp.diff(g.subs({X: xs, Y: ys}), s)))
# 0
# Matrix form: derivative of the composition = (row of f's partials) @ (Jacobian of the inner map)
s0, t0 = 1.0, 2.0
x0, y0 = s0 * t0**2, s0**2 * t0
grad_g = np.array([np.exp(x0) * np.sin(y0), np.exp(x0) * np.cos(y0)]) # [g_x, g_y]
J = np.array([[t0**2, 2 * s0 * t0], # [[x_s, x_t],
[2 * s0 * t0, s0**2]]) # [y_s, y_t]]
print(np.round(grad_g @ J, 4)) # [z_s, z_t]
# [107.7004 175.863 ]
h = 1e-6
zf = lambda s_, t_: np.exp(s_ * t_**2) * np.sin(s_**2 * t_)
print(np.round([(zf(s0 + h, t0) - zf(s0 - h, t0)) / (2 * h), (zf(s0, t0 + h) - zf(s0, t0 - h)) / (2 * h)], 4))
# [107.7004 175.863 ]
Key Takeaways
- When $z$ depends on $t$ through several intermediate variables, $\frac{dz}{dt}$ is the sum of the rates through each one: $z_x x' + z_y y'$.
- With several parameters, compute each partial the same way, holding the other parameters fixed.
- In matrix form, the derivative of a composition is the product of the outer row of partials and the inner Jacobian matrix.
- The outer function must be differentiable; existence of partials is not enough.
🧪 Practice
- Use the chain rule to find $\frac{dw}{dt}$ for $w = xe^{y/z}$ with $x = t^2$, $y = 1 - t$, $z = 1 + 2t$, at $t = 0$.
- Let $z = f(x, y)$ with $x = r\cos\theta$, $y = r\sin\theta$. Show that $(z_x)^2 + (z_y)^2 = (z_r)^2 + \frac{1}{r^2}(z_\theta)^2$.
- The radius of a cylinder grows at $0.2$ cm/s while its height shrinks at $0.5$ cm/s. How fast is the volume changing when $r = 4$ cm and $h = 10$ cm?
- Interview-style: Backpropagation through a network with layers $\mathbf{x} \to \mathbf{h}_1 \to \mathbf{h}_2 \to L$ computes $\frac{\partial L}{\partial \mathbf{x}}$ as a product of Jacobian matrices. Why do frameworks multiply these matrices starting from the loss end rather than from the input end? Hint: the loss is a single number, so its row of partials is $1 \times n$. Multiplying from the left keeps every intermediate result a single row (vector-matrix products), while multiplying from the right builds full matrix-matrix products.
Tree Diagrams for Dependencies
Theory
When variables depend on other variables through several layers, it is easy to forget a term in the chain rule. A tree diagram makes the bookkeeping mechanical. Put the final variable at the top, draw a branch to each variable it depends on directly, and continue down to the independent variables. Label each branch with the partial derivative of the upper variable with respect to the lower one.
The rule. To find $\frac{\partial (\text{top})}{\partial (\text{bottom variable})}$:
- Find every path from the top down to that bottom variable.
- Multiply the branch labels along each path.
- Add the products over all paths.
For $w = f(x, y, z)$ with $x$, $y$, $z$ each depending on $s$ and $t$:
w
w_x / | w_y \ w_z
x y z
x_s / \ x_t y_s / \ y_t z_s / \ z_t
s t s t s t
dw/ds: three paths ending in s -> w_x x_s + w_y y_s + w_z z_s
dw/dt: three paths ending in t -> w_x x_t + w_y y_t + w_z z_t
The tree diagram is exactly a computational graph. Automatic differentiation (as in PyTorch or JAX) stores this graph while computing a value, then walks it from the top down, multiplying along edges and adding across paths. That walk is backpropagation.
Pitfall: a variable appearing at two levels. Suppose $w = f(x, y, t)$ where $x$ and $y$ also depend on $t$. Then $t$ influences $w$ directly and through $x$ and $y$:
w
w_x / |w_y \ w_t dw/dt = w_x x'(t) + w_y y'(t) + w_t
x y t \___ indirect ___/ \_ direct
x' | | y'
t t
Now "$\frac{\partial w}{\partial t}$" is ambiguous: does it mean the direct partial (holding $x$ and $y$ fixed) or the total rate? Write $\frac{dw}{dt}$ for the total rate and $\frac{\partial f}{\partial t}$ for the direct partial. Physics uses this constantly: for a temperature field $T(x, y, z, t)$ measured by a moving probe, the probe sees
$$\frac{dT}{dt} = \frac{\partial T}{\partial t} + T_x x' + T_y y' + T_z z',$$
the local change of the field plus the change due to motion through it (the material derivative in fluid mechanics). Thermodynamics avoids the ambiguity with subscripts, writing $\left(\frac{\partial U}{\partial T}\right)_V$ for "with $V$ held fixed".
Examples
-
For $w = f(x, t) = xt$ with $x = t^2$: the direct partial is $\frac{\partial f}{\partial t} = x = t^2$, but the total derivative is $\frac{dw}{dt} = f_x,x' + f_t = t\cdot 2t + t^2 = 3t^2$, which matches $w = t^3$ after substitution.
-
$w = xy + yz + zx$ with $x = st$, $y = e^{st}$, $z = t^2$, at $(s, t) = (1, 2)$. The tree has three paths to each of $s$ and $t$ (the path through $z$ to $s$ contributes $0$, since $z$ does not depend on $s$). The code below builds the same tree as a computation graph, a tiny reverse-mode automatic differentiator, and checks it against the hand formula.
import math
class Node:
"""A value in a computation graph that remembers which nodes it was built from."""
def __init__(self, value, parents=()):
self.value = value
self.parents = parents # list of (parent node, local partial derivative)
self.grad = 0.0
def __add__(self, other):
return Node(self.value + other.value, [(self, 1.0), (other, 1.0)])
def __mul__(self, other):
return Node(self.value * other.value, [(self, other.value), (other, self.value)])
def exp(n):
v = math.exp(n.value)
return Node(v, [(n, v)])
def backward(output):
"""Sum over all paths from output to each input of the product of local partials."""
order, seen = [], set()
def visit(n): # topological order: parents before children
if id(n) not in seen:
seen.add(id(n))
for p, _ in n.parents:
visit(p)
order.append(n)
visit(output)
output.grad = 1.0
for n in reversed(order): # walk the tree diagram from the top down
for p, local in n.parents:
p.grad += n.grad * local # multiply along the edge, add across paths
# w = xy + yz + zx with x = st, y = e^(st), z = t^2, at s = 1, t = 2
s, t = Node(1.0), Node(2.0)
x = s * t
y = exp(s * t)
z = t * t
w = x * y + y * z + z * x
backward(w)
print(round(s.grad, 4), round(t.grad, 4))
# 111.4468 93.2796
# Check with the tree-diagram formula dw/ds = w_x x_s + w_y y_s + w_z z_s
S, T = 1.0, 2.0
X, Y, Z = S * T, math.exp(S * T), T**2
w_x, w_y, w_z = Y + Z, X + Z, Y + X
print(round(w_x * T + w_y * T * Y + w_z * 0, 4), round(w_x * S + w_y * S * Y + w_z * 2 * T, 4))
# 111.4468 93.2796
Key Takeaways
- Draw the dependency tree; for each path from the top to the target variable, multiply the branch derivatives, then add over all paths.
- A variable that is not reached by any path contributes nothing; a variable reached by several paths contributes one term per path.
- When a variable affects the output both directly and indirectly, distinguish the total derivative $\frac{dw}{dt}$ from the direct partial $\frac{\partial f}{\partial t}$.
- The tree diagram is a computation graph, and the path rule is what reverse-mode automatic differentiation (backpropagation) implements.
🧪 Practice
- Draw the tree diagram and write the chain rule for $\frac{\partial u}{\partial t}$ when $u = f(x, y, z)$, $x = x(r, s, t)$, $y = y(r, s, t)$, and $z = z(r, s, t)$.
- Let $w = f(u, v)$ with $u = x^2 - y^2$ and $v = 2xy$. Find $w_x$ and $w_y$ in terms of $f_u$ and $f_v$, then compute $w_{xx}$ (draw a new tree for $f_u$ and $f_v$, which also depend on $u$ and $v$).
- A boat moves through a lake whose temperature is $T(x, y, t) = 15 + 0.1x - 0.05y + 0.2t$ (degrees, with $x$, $y$ in km and $t$ in hours). The boat's position is $(3t, 4t)$. Find the rate of change of temperature the boat measures, and separate it into the direct and the motion parts.
- Interview-style: In an automatic-differentiation engine, a node's value is used by several later nodes (for example, $y$ above feeds both $xy$ and $yz$). What must the backward pass do at such a node, and what bug appears if gradients are overwritten instead of accumulated? Hint: each use is a separate path in the tree diagram, so the contributions must be added (
+=). Overwriting keeps only the last path, giving silently wrong gradients that a finite-difference check would catch.
Implicit Differentiation in Several Variables
Theory
Section 3.5 differentiated equations like $x^2 + y^2 = 25$ without solving for $y$, by differentiating both sides and solving for $y'$. The multivariable chain rule turns this procedure into a formula.
Curves. Suppose the equation $F(x, y) = 0$ defines $y$ as a differentiable function of $x$. Then $F(x, y(x)) = 0$ for all $x$ in an interval, so its derivative with respect to $x$ is zero. By the chain rule (tree: $F$ depends on $x$ directly and on $x$ through $y$),
$$F_x \cdot 1 + F_y \cdot \frac{dy}{dx} = 0 \qquad \implies \qquad \boxed{\frac{dy}{dx} = -\frac{F_x}{F_y}} \quad (F_y \ne 0).$$
Surfaces. If $F(x, y, z) = 0$ defines $z$ as a function of $x$ and $y$, differentiate in $x$ with $y$ held fixed, then in $y$ with $x$ held fixed:
$$\frac{\partial z}{\partial x} = -\frac{F_x}{F_z}, \qquad \frac{\partial z}{\partial y} = -\frac{F_y}{F_z} \qquad (F_z \ne 0).$$
The minus sign is easy to forget, and it is not a typo. The formula is not "$\frac{\partial F/\partial x}{\partial F/\partial y}$ cancels the $\partial F$"; the $\partial$ symbols do not cancel like fractions.
Why it works geometrically. Moving along the curve $F = 0$ keeps $F$ constant, so the change from moving in $x$ ($F_x,dx$) must be exactly canceled by the change from the accompanying move in $y$ ($F_y,dy$). The ratio of the moves is forced: $\frac{dy}{dx} = -\frac{F_x}{F_y}$.
The denominator. The formula fails where $F_y = 0$ (or $F_z = 0$). These are exactly the points where the curve has a vertical tangent, or where $y$ stops being a function of $x$ nearby. The next topic makes this precise.
Examples
- Folium of Descartes, $x^3 + y^3 = 6xy$. With $F = x^3 + y^3 - 6xy$:
$$\frac{dy}{dx} = -\frac{3x^2 - 6y}{3y^2 - 6x} = \frac{x^2 - 2y}{2x - y^2}.$$
At $(3, 3)$: $\frac{9 - 6}{6 - 9} = -1$.
- The surface $x^3 + y^3 + z^3 + 6xyz = 1$. With $F = x^3 + y^3 + z^3 + 6xyz - 1$:
$$\frac{\partial z}{\partial x} = -\frac{3x^2 + 6yz}{3z^2 + 6xy} = -\frac{x^2 + 2yz}{z^2 + 2xy}, \qquad \frac{\partial z}{\partial y} = -\frac{y^2 + 2xz}{z^2 + 2xy}.$$
- The ideal gas law $PV = nRT$ links three variables. Holding $n$ fixed and using the formula on $F = PV - nRT$ gives the classic identity $\left(\frac{\partial P}{\partial V}\right)\left(\frac{\partial V}{\partial T}\right)\left(\frac{\partial T}{\partial P}\right) = -1$. The product of the three partials is $-1$, not $+1$, a famous warning against cancelling $\partial$'s.
import sympy as sp
from scipy.optimize import brentq
x, y, z = sp.symbols('x y z', real=True)
# Curve: folium of Descartes x^3 + y^3 = 6xy
F = x**3 + y**3 - 6 * x * y
dydx = -sp.diff(F, x) / sp.diff(F, y) # dy/dx = -F_x / F_y
print(sp.simplify(dydx), dydx.subs({x: 3, y: 3}))
# (x**2 - 2*y)/(2*x - y**2) -1
print(sp.idiff(F, y, x).subs({x: 3, y: 3})) # sympy's built-in implicit differentiation agrees
# -1
# Surface: x^3 + y^3 + z^3 + 6xyz = 1, find dz/dx and dz/dy
G = x**3 + y**3 + z**3 + 6 * x * y * z - 1
zx = -sp.diff(G, x) / sp.diff(G, z)
zy = -sp.diff(G, y) / sp.diff(G, z)
print(sp.simplify(zx), '|', sp.simplify(zy))
# (-x**2 - 2*y*z)/(2*x*y + z**2) | (-2*x*z - y**2)/(2*x*y + z**2)
# Numerical check: solve G = 0 for z near a point, nudge x, and difference
Gf = sp.lambdify((x, y, z), G)
x0, y0 = 0.5, 0.2
z_of = lambda xv, yv: brentq(lambda zv: Gf(xv, yv, zv), -3, 3)
z0, h = z_of(x0, y0), 1e-6
print(round(z0, 6), round((z_of(x0 + h, y0) - z_of(x0 - h, y0)) / (2 * h), 6), round(float(zx.subs({x: x0, y: y0, z: z0})), 6))
# 0.747864 -0.723225 -0.723225 z, finite-difference slope, formula slope
Key Takeaways
- For a curve $F(x, y) = 0$: $\frac{dy}{dx} = -\frac{F_x}{F_y}$ wherever $F_y \ne 0$.
- For a surface $F(x, y, z) = 0$: $z_x = -\frac{F_x}{F_z}$ and $z_y = -\frac{F_y}{F_z}$ wherever $F_z \ne 0$.
- These come from applying the chain rule to $F(x, y(x)) = 0$; the minus sign is essential, and $\partial$'s do not cancel.
- The formulas break down exactly where the denominator vanishes: vertical tangents and points where the implicit function fails to exist.
🧪 Practice
- Find $\frac{dy}{dx}$ for $\cos(xy) = 1 + \sin y$ using $-\frac{F_x}{F_y}$.
- Find $\frac{\partial z}{\partial x}$ and $\frac{\partial z}{\partial y}$ for $e^z = xyz$ and for $yz + x\ln y = z^2$.
- Verify the identity $\left(\frac{\partial P}{\partial V}\right)_T\left(\frac{\partial V}{\partial T}\right)_P\left(\frac{\partial T}{\partial P}\right)_V = -1$ for the ideal gas law, by computing each factor.
- Interview-style: A circuit simulator needs $\frac{dV}{dI}$ for a diode, whose voltage and current are related by an equation $F(V, I) = 0$ that has no closed-form solution for $V$. How would you get the derivative at an operating point, and what condition makes it well defined? Hint: solve $F = 0$ numerically for the operating point, then use $\frac{dV}{dI} = -\frac{F_I}{F_V}$, which needs only partial derivatives of $F$ at that point and requires $F_V \ne 0$.
Implicit Function Theorem
Theory
Implicit differentiation assumes that the equation $F(x, y) = 0$ really defines $y$ as a differentiable function of $x$. Does it? The circle $x^2 + y^2 = 1$ is not the graph of one function, but near most of its points it is: near $(0.6, 0.8)$ it is $y = \sqrt{1 - x^2}$, and near $(0.6, -0.8)$ it is $y = -\sqrt{1 - x^2}$. Near $(1, 0)$, however, every vertical line just left of $x = 1$ meets the circle twice and every vertical line just right meets it not at all. No single function $y(x)$ describes the circle near $(1, 0)$. There, $F_y = 2y = 0$.
The theorem says the derivative of $F$ is what decides this.
Implicit Function Theorem (two variables). Suppose $F$ has continuous partial derivatives near $(a, b)$, with $F(a, b) = 0$ and $F_y(a, b) \ne 0$. Then there is an open interval $I$ around $a$ and a unique differentiable function $y = g(x)$ on $I$ with $g(a) = b$ and $F(x, g(x)) = 0$, and $g'(x) = -\frac{F_x}{F_y}$.
Why $F_y \ne 0$ is the right condition. Near $(a, b)$, $F$ is approximately its linearization: $F_x(a, b)(x - a) + F_y(a, b)(y - b) \approx 0$. If $F_y \ne 0$, this linear equation can be solved for $y$ uniquely, and the theorem promises the curved equation can be too, at least locally. If $F_y = 0$, the linear equation says nothing about $y$, and the curve may turn back on itself, cross itself, or have a vertical tangent.
| At a point on $F = 0$ | Conclusion |
|---|---|
| $F_y \ne 0$ | Locally $y = g(x)$, a smooth graph over $x$ |
| $F_y = 0$, $F_x \ne 0$ | Locally $x = h(y)$ instead: the curve has a vertical tangent |
| $F_x = F_y = 0$ | Singular point; anything can happen (a crossing, a cusp, an isolated point) |
The general theorem. For a system of $m$ equations $\mathbf{F}(\mathbf{x}, \mathbf{y}) = \mathbf{0}$ in $n + m$ unknowns, $\mathbf{y}$ can be solved locally as a smooth function of $\mathbf{x}$ when the $m \times m$ Jacobian matrix $\frac{\partial \mathbf{F}}{\partial \mathbf{y}}$ is invertible. The derivative is then the matrix version of $-\frac{F_x}{F_y}$:
$$\frac{\partial \mathbf{y}}{\partial \mathbf{x}} = -\left(\frac{\partial \mathbf{F}}{\partial \mathbf{y}}\right)^{-1}\frac{\partial \mathbf{F}}{\partial \mathbf{x}}.$$
This justifies treating constraint surfaces as smooth (Lagrange multipliers, Section 12.6), describes how equilibria of a system move when parameters change (comparative statics in economics), and underlies the inverse function theorem.
Computational meaning. The condition $F_y \ne 0$ is also exactly what Newton's method in $y$ (Section 4.5) needs: the iteration $y \leftarrow y - \frac{F(x, y)}{F_y(x, y)}$ divides by $F_y$. Tracking a solution $y(x)$ as $x$ changes, using the last solution as the next starting guess, is called numerical continuation. It works smoothly until $F_y$ approaches $0$. There the slope $-\frac{F_x}{F_y}$ blows up, Newton slows down, and past that point the branch simply ends.
Examples
- Circle $F = x^2 + y^2 - 1$: $F_y = 2y \ne 0$ except at $(\pm 1, 0)$. Everywhere else, $y$ is locally a function of $x$, with $y' = -\frac{x}{y}$. At $(\pm 1, 0)$, $F_x = \pm 2 \ne 0$, so $x$ is locally a function of $y$ there instead.
- Folium $x^3 + y^3 = 6xy$: $F_y = 3y^2 - 6x$. Solving $F = 0$ and $F_y = 0$ together gives $(0, 0)$ (where also $F_x = 0$: the curve crosses itself) and $(2^{5/3}, 2^{4/3}) \approx (3.175, 2.520)$, where the loop has a vertical tangent.
Folium near its rightmost point (3.175, 2.520)
y | ,--.
| / \ F_y != 0: the branch through (3, 3)
| / ) <- vertical tangent, F_y = 0
|/ _.-' beyond x = 3.175 this branch does not exist
+------------- x
The code continues the solution from $(3, 3)$ toward the vertical tangent. Watch $F_y$ shrink and the slope explode; just past the turning point, Newton's method converges to a completely different branch of the curve.
import numpy as np
# Follow the branch of the folium x^3 + y^3 = 6xy through (3, 3) as x increases
F = lambda x, y: x**3 + y**3 - 6 * x * y
Fx = lambda x, y: 3 * x**2 - 6 * y
Fy = lambda x, y: 3 * y**2 - 6 * x
def solve_y(x, y_guess, tol=1e-12):
"""Newton's method in y alone: possible exactly when F_y != 0 (the implicit function theorem)."""
y = y_guess
for it in range(1, 100):
step = F(x, y) / Fy(x, y)
y -= step
if abs(step) < tol:
return y, it
return y, None
y = 3.0
for x in [3.0, 3.1, 3.15, 3.17, 3.174, 3.1748]:
y, iters = solve_y(x, y) # previous solution is the next starting guess
slope = -Fx(x, y) / Fy(x, y) # dy/dx from implicit differentiation
print(f'{x:7.4f} {y:9.6f} {Fy(x, y):10.6f} {slope:12.3f} {iters}')
# 3.0000 3.000000 9.000000 -1.000 1
# 3.1000 2.861905 5.971499 -1.952 5
# 3.1500 2.728464 3.433553 -3.902 6
# 3.1700 2.615220 1.498121 -9.649 6
# 3.1740 2.559465 0.608577 -24.428 6
# 3.1748 2.521892 0.031022 -486.964 9 F_y -> 0, slope -> infinity
x_star = 2 ** (5 / 3) # where F = 0 and F_y = 0 together
print(round(x_star, 6), round(2 ** (4 / 3), 6), F(x_star, 2 ** (4 / 3)), Fy(x_star, 2 ** (4 / 3)))
# 3.174802 2.519842 0.0 0.0
print(solve_y(3.18, 2.52)) # past the turning point: no solution on this branch
# (-5.045184463757283, 38) Newton wanders to another branch of the curve
Key Takeaways
- If $F(a, b) = 0$ and $F_y(a, b) \ne 0$ (with continuous partials), then near $(a, b)$ the equation $F = 0$ defines a unique smooth function $y = g(x)$ with $g' = -\frac{F_x}{F_y}$.
- Where $F_y = 0$ the curve may have a vertical tangent, turn back, or cross itself; try solving for $x$ instead, or treat it as a singular point.
- In general, $\mathbf{F}(\mathbf{x}, \mathbf{y}) = \mathbf{0}$ defines $\mathbf{y}(\mathbf{x})$ locally when the Jacobian $\frac{\partial \mathbf{F}}{\partial \mathbf{y}}$ is invertible.
- The same condition is what Newton's method in $\mathbf{y}$ needs; numerical continuation breaks down where it fails.
🧪 Practice
- At which points of the ellipse $x^2 + 4y^2 = 4$ does the Implicit Function Theorem fail to give $y$ as a function of $x$? Near those points, can $x$ be written as a function of $y$?
- Show that the equation $x^2 + y^3 - xy = 1$ defines $y$ as a differentiable function of $x$ near $(1, 1)$, and find $y'(1)$.
- Near which points of the surface $x^2 + y^2 - z^2 = 1$ can $z$ be written as a smooth function of $(x, y)$? Describe the set where it cannot.
- Interview-style: You track the equilibrium price $p$ of a market as a tax rate $\tau$ slowly changes, by re-solving $F(p, \tau) = 0$ with Newton's method from the previous solution. At some value of $\tau$ the solver starts taking many iterations and then jumps to a wildly different price. What is likely happening, and how would you detect it in advance? Hint: monitor $F_p$ along the path. As it approaches $0$, the branch is reaching a fold (a turning point), and beyond it the old equilibrium no longer exists. Arc-length continuation, which parametrizes by distance along the curve instead of by $\tau$, can follow the curve around the fold. <a id="125-directional-derivatives-and-gradients"></a>
12.5 Directional Derivatives and Gradients
This subchapter measures the rate of change of a function in any direction, packages all of these rates into a single vector (the gradient), and shows that the gradient points uphill and is perpendicular to level curves and level surfaces.
Directional Derivative Definition
Theory
Partial derivatives give the slope of a surface walking due east ($f_x$) or due north ($f_y$). A hiker can walk in any direction, though, and a weather balloon drifts wherever the wind takes it. We want the rate of change of $f$ per unit distance in an arbitrary direction.
A direction in the plane is a unit vector $\mathbf{u} = \langle a, b \rangle$ with $a^2 + b^2 = 1$, for example $\mathbf{u} = \langle \cos\theta, \sin\theta \rangle$. Starting at $(x_0, y_0)$ and moving a distance $h$ in direction $\mathbf{u}$ brings you to $(x_0 + ha, y_0 + hb)$.
Definition. The directional derivative of $f$ at $(x_0, y_0)$ in the direction of the unit vector $\mathbf{u}$ is
$$D_{\mathbf{u}} f(x_0, y_0) = \lim_{h \to 0} \frac{f(x_0 + ha,\ y_0 + hb) - f(x_0, y_0)}{h},$$
if the limit exists. Geometrically, slice the surface $z = f(x, y)$ with the vertical plane through $(x_0, y_0)$ in direction $\mathbf{u}$; $D_{\mathbf{u}} f$ is the slope of the resulting curve.
The partial derivatives are special cases: $f_x = D_{\mathbf{i}} f$ with $\mathbf{u} = \langle 1, 0 \rangle$, and $f_y = D_{\mathbf{j}} f$ with $\mathbf{u} = \langle 0, 1 \rangle$. Reversing the direction flips the sign: $D_{-\mathbf{u}} f = -D_{\mathbf{u}} f$.
Why a unit vector. The definition divides by $h$, which is the distance traveled only when $\lvert \mathbf{u} \rvert = 1$. With a vector of length $2$, a step of $h$ covers distance $2h$, and the quotient doubles. If you are given a direction as any nonzero vector $\mathbf{v}$, first normalize: $\mathbf{u} = \frac{\mathbf{v}}{\lvert \mathbf{v} \rvert}$.
The computing formula. Applying the one-parameter chain rule (Section 12.4) to $g(h) = f(x_0 + ha, y_0 + hb)$ gives $g'(0) = f_x,a + f_y,b$. So, if $f$ is differentiable,
$$D_{\mathbf{u}} f(x_0, y_0) = f_x(x_0, y_0),a + f_y(x_0, y_0),b.$$
Two partial derivatives determine the slope in every direction. For $n$ variables, $D_{\mathbf{u}} f = \sum_i f_{x_i} u_i$.
The differentiability requirement matters. For $g = \frac{x^2y}{x^2 + y^2}$ with $g(0, 0) = 0$, the definition gives $D_{\mathbf{u}} g(0, 0) = a^2b$ in every direction (since $g(ha, hb) = h,a^2b$). Both partials are $0$, so the formula would predict $0$ in every direction; along $\langle \frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}} \rangle$ the true value is $\frac{1}{2\sqrt{2}} \approx 0.354$. All directional derivatives exist, but the function is not differentiable at the origin, so they are not linear in $\mathbf{u}$.
Examples
$f(x, y) = x^3 - 3xy + 4y^2$ at $(1, 2)$, in the direction making angle $\frac{\pi}{6}$ with the $x$-axis, $\mathbf{u} = \langle \frac{\sqrt{3}}{2}, \frac{1}{2} \rangle$:
- $f_x = 3x^2 - 3y = -3$ and $f_y = -3x + 8y = 13$ at $(1, 2)$.
- $D_{\mathbf{u}} f = -3\cdot\frac{\sqrt{3}}{2} + 13\cdot\frac{1}{2} = \frac{13 - 3\sqrt{3}}{2} \approx 3.902$.
Walking in that direction, the surface rises about $3.9$ units per unit distance.
import numpy as np
f = lambda x, y: x**3 - 3 * x * y + 4 * y**2
a, b = 1.0, 2.0
u = np.array([np.cos(np.pi / 6), np.sin(np.pi / 6)]) # unit vector at 30 degrees
# Definition: rate of change per unit distance along u
for h in [1e-1, 1e-3, 1e-5]:
print(h, round((f(a + h * u[0], b + h * u[1]) - f(a, b)) / h, 6))
# 0.1 4.103515
# 0.001 3.903875
# 1e-05 3.901943
fx, fy = 3 * a**2 - 3 * b, -3 * a + 8 * b # partials at (1, 2): -3 and 13
print(round(fx * u[0] + fy * u[1], 6), round((13 - 3 * np.sqrt(3)) / 2, 6))
# 3.901924 3.901924
# Pitfall: a non-unit vector scales the answer by its length
v = 2 * u
h = 1e-6
print(round((f(a + h * v[0], b + h * v[1]) - f(a, b)) / h, 4)) # twice the true rate
# 7.8039
# All directional derivatives can exist and still not obey the formula: f = x^2 y / (x^2 + y^2)
g = lambda x, y: x**2 * y / (x**2 + y**2) if (x, y) != (0, 0) else 0.0
for deg in [0, 45, 90, 135]:
w = np.array([np.cos(np.radians(deg)), np.sin(np.radians(deg))])
print(deg, round((g(*(1e-8 * w)) - 0) / 1e-8, 6), round(w[0]**2 * w[1], 6)) # g_x = g_y = 0 would predict 0
# 0 0.0 0.0
# 45 0.353553 0.353553
# 90 0.0 0.0
# 135 0.353553 0.353553
Key Takeaways
- $D_{\mathbf{u}} f$ is the rate of change of $f$ per unit distance in the direction of the unit vector $\mathbf{u}$; always normalize a given direction first.
- Partial derivatives are directional derivatives along $\mathbf{i}$ and $\mathbf{j}$.
- If $f$ is differentiable, $D_{\mathbf{u}} f = f_x a + f_y b$: two partials determine every direction.
- Without differentiability, directional derivatives may all exist and still not follow this formula.
🧪 Practice
- Find the directional derivative of $f(x, y) = x^2y^3 - 4y$ at $(2, -1)$ in the direction of $\mathbf{v} = \langle 2, 5 \rangle$.
- Find $D_{\mathbf{u}} f(0, 0)$ for $f(x, y) = e^{x}\cos(y)$ in every direction $\mathbf{u} = \langle \cos\theta, \sin\theta \rangle$. For which $\theta$ is it zero?
- Use the definition directly to find $D_{\mathbf{u}} f(0, 0)$ for $f(x, y) = \sqrt{x^2 + y^2}$. Show it exists in no direction, and explain what that says about the cone's tip.
- Interview-style: A function of $n = 10^6$ variables can be evaluated cheaply, and you need its directional derivative along one given direction $\mathbf{v}$, but not the full gradient. How many function evaluations do you need, and why? Hint: the directional derivative is a single one-variable derivative of $g(h) = f(\mathbf{x} + h\mathbf{u})$, so two evaluations suffice for a central difference. Forward-mode automatic differentiation computes the exact value at about the cost of one evaluation.
The Gradient Vector
Theory
The formula $D_{\mathbf{u}} f = f_x a + f_y b$ is a dot product: of $\langle f_x, f_y \rangle$ with $\langle a, b \rangle = \mathbf{u}$. The first vector depends only on $f$ and the point, not on the direction. It deserves a name.
Definition. The gradient of $f$ is the vector of its partial derivatives:
$$\nabla f(x, y) = \langle f_x(x, y),\ f_y(x, y) \rangle, \qquad \nabla f(x, y, z) = \langle f_x,\ f_y,\ f_z \rangle,$$
and in $n$ variables $\nabla f = \langle f_{x_1}, \ldots, f_{x_n} \rangle$. The symbol $\nabla$ is read "del" or "nabla"; some texts write $\operatorname{grad} f$.
With it, the directional derivative is a dot product:
$$\boxed{D_{\mathbf{u}} f = \nabla f \cdot \mathbf{u}}$$
So the gradient is a compact summary of how $f$ changes in all directions at once. Knowing this one vector, you know the slope in every direction. This is the multivariable counterpart of "the derivative is a number": for a differentiable function of several variables, the derivative at a point is the gradient vector (or, as a linear map, $\mathbf{h} \mapsto \nabla f \cdot \mathbf{h}$).
The linearization of Section 12.3 also becomes cleaner:
$$f(\mathbf{x}) \approx f(\mathbf{a}) + \nabla f(\mathbf{a}) \cdot (\mathbf{x} - \mathbf{a}).$$
This has exactly the form of the one-variable tangent line $f(a) + f'(a)(x - a)$, with a dot product in place of a product.
Properties. The gradient obeys the familiar rules, applied component by component:
| Rule | Formula |
|---|---|
| Linearity | $\nabla(af + bg) = a\nabla f + b\nabla g$ |
| Product rule | $\nabla(fg) = f\nabla g + g\nabla f$ |
| Quotient rule | $\nabla!\left(\frac{f}{g}\right) = \frac{g\nabla f - f\nabla g}{g^2}$ |
| Chain rule | $\nabla(h(f)) = h'(f),\nabla f$ for a one-variable function $h$ |
| Useful special cases | $\nabla(\mathbf{c} \cdot \mathbf{x}) = \mathbf{c}$, $\nabla(\lvert \mathbf{x} \rvert^2) = 2\mathbf{x}$, $\nabla \lvert \mathbf{x} \rvert = \frac{\mathbf{x}}{\lvert \mathbf{x} \rvert}$ |
A vector field. The gradient assigns a vector to every point, so $\nabla f$ is a vector field. Gravitational and electric forces are gradients of potential functions (up to sign). Chapter 14 studies such fields.
Examples
- $f(x, y) = \sin x + e^{xy}$: $\nabla f = \langle \cos x + ye^{xy},\ xe^{xy} \rangle$, and $\nabla f(0, 1) = \langle 2, 0 \rangle$.
- $f(x, y, z) = x\sin(yz)$: $\nabla f = \langle \sin(yz),\ xz\cos(yz),\ xy\cos(yz) \rangle$. At $(1, 3, 0)$ this is $\langle 0, 0, 3 \rangle$. In the direction of $\mathbf{v} = \langle 1, 2, -1 \rangle$, $\mathbf{u} = \frac{1}{\sqrt{6}}\langle 1, 2, -1 \rangle$, so
$$D_{\mathbf{u}} f(1, 3, 0) = \langle 0, 0, 3 \rangle \cdot \tfrac{1}{\sqrt{6}}\langle 1, 2, -1 \rangle = -\tfrac{3}{\sqrt{6}} = -\tfrac{\sqrt{6}}{2} \approx -1.225.$$
- The squared distance $f(\mathbf{x}) = \lvert \mathbf{x} - \mathbf{a} \rvert^2$ has $\nabla f = 2(\mathbf{x} - \mathbf{a})$: it points straight away from $\mathbf{a}$, the direction in which distance grows fastest.
import numpy as np
import sympy as sp
x, y, z = sp.symbols('x y z', real=True)
f = x * sp.sin(y * z)
grad = sp.Matrix([sp.diff(f, var) for var in (x, y, z)]) # vector of partials
print(grad.T)
# Matrix([[sin(y*z), x*z*cos(y*z), x*y*cos(y*z)]])
P = {x: 1, y: 3, z: 0}
gP = grad.subs(P)
v = sp.Matrix([1, 2, -1])
print(gP.T, gP.dot(v / v.norm()), sp.N(gP.dot(v / v.norm()), 6)) # D_u f = grad . u
# Matrix([[0, 0, 3]]) -sqrt(6)/2 -1.22474
def numerical_gradient(F, p, h=1e-6):
"""Central differences, one coordinate at a time: 2n evaluations for n variables."""
p = np.asarray(p, float)
g = np.zeros_like(p)
for i in range(p.size):
e = np.zeros_like(p)
e[i] = h
g[i] = (F(p + e) - F(p - e)) / (2 * h)
return g
F = sp.lambdify([(x, y, z)], f)
print(np.round(numerical_gradient(F, [1, 3, 0]), 6))
# [0. 0. 3.]
p = [2.0, 0.5, 1.5]
print(np.round(numerical_gradient(F, p), 6), np.round(np.array(grad.subs(dict(zip((x, y, z), p))), float).ravel(), 6))
# [0.681639 2.195067 0.731689] [0.681639 2.195067 0.731689]
Key Takeaways
- $\nabla f = \langle f_x, f_y, \ldots \rangle$, the vector of all partial derivatives.
- $D_{\mathbf{u}} f = \nabla f \cdot \mathbf{u}$: the gradient encodes the rate of change in every direction.
- Linearization in vector form: $f(\mathbf{x}) \approx f(\mathbf{a}) + \nabla f(\mathbf{a}) \cdot (\mathbf{x} - \mathbf{a})$.
- The gradient follows linearity, product, quotient, and chain rules; $\nabla \lvert \mathbf{x} \rvert^2 = 2\mathbf{x}$ and $\nabla(\mathbf{c} \cdot \mathbf{x}) = \mathbf{c}$ are worth memorizing.
🧪 Practice
- Find $\nabla f$ for $f(x, y, z) = \frac{x}{y} + \ln(yz)$, and evaluate it at $(1, 1, 1)$.
- Show that $\nabla \lvert \mathbf{x} \rvert = \frac{\mathbf{x}}{\lvert \mathbf{x} \rvert}$ for $\mathbf{x} \ne \mathbf{0}$ in $\mathbb{R}^3$, and interpret it geometrically.
- For $f(\mathbf{x}) = \mathbf{x}^{\mathsf{T}} A,\mathbf{x}$ with a $2 \times 2$ matrix $A$, show by expanding that $\nabla f = (A + A^{\mathsf{T}})\mathbf{x}$, which equals $2A\mathbf{x}$ when $A$ is symmetric.
- Interview-style: Before trusting a hand-written gradient function in a training loop, how would you test it, and what tolerance would you use? Hint: a gradient check: compare against central differences $\frac{f(\mathbf{x} + h\mathbf{e}i) - f(\mathbf{x} - h\mathbf{e}i)}{2h}$ at a few random points, using a relative error $\frac{\lvert g{\text{code}} - g{\text{num}} \rvert}{\max(\lvert g_{\text{code}} \rvert, \lvert g_{\text{num}} \rvert, \text{tiny})}$. With double precision and $h \approx 10^{-5}$, relative errors around $10^{-7}$ or smaller are healthy; $10^{-2}$ signals a bug. Avoid points near kinks such as $\max(0, u)$.
Maximizing the Rate of Change
Theory
Standing on a hillside in fog, which way should you walk to climb fastest? The gradient answers this directly. Write the dot product with the angle $\theta$ between $\nabla f$ and $\mathbf{u}$ (Section 11.1), and use $\lvert \mathbf{u} \rvert = 1$:
$$D_{\mathbf{u}} f = \nabla f \cdot \mathbf{u} = \lvert \nabla f \rvert,\lvert \mathbf{u} \rvert\cos\theta = \lvert \nabla f \rvert\cos\theta.$$
As $\mathbf{u}$ turns, only $\cos\theta$ changes, and everything follows from its range $[-1, 1]$:
| Direction $\mathbf{u}$ | Angle $\theta$ | Rate $D_{\mathbf{u}} f$ | Meaning |
|---|---|---|---|
| Same as $\nabla f$ | $0$ | $\lvert \nabla f \rvert$ (maximum) | Steepest ascent |
| Opposite to $\nabla f$ | $\pi$ | $-\lvert \nabla f \rvert$ (minimum) | Steepest descent |
| Perpendicular to $\nabla f$ | $\frac{\pi}{2}$ | $0$ | No change: along the level curve |
| In between | $\theta$ | $\lvert \nabla f \rvert\cos\theta$ | Partial climb |
Theorem. If $f$ is differentiable and $\nabla f \ne \mathbf{0}$ at a point, then the maximum rate of increase of $f$ is $\lvert \nabla f \rvert$, attained in the direction of $\nabla f$. The maximum rate of decrease is $\lvert \nabla f \rvert$, in the direction of $-\nabla f$.
So the gradient has a meaning independent of coordinates: its direction is the uphill direction, and its length is the steepness in that direction. When $\nabla f = \mathbf{0}$, no direction gives a first-order change; such points are the critical points of Section 12.6.
Gradient descent. To minimize a function, repeatedly take a small step downhill:
$$\mathbf{x}_{k+1} = \mathbf{x}_k - \alpha,\nabla f(\mathbf{x}_k),$$
where $\alpha > 0$ is the step size (learning rate). This is how most machine-learning models are trained, with $f$ the loss and $\mathbf{x}$ the parameters. The step size is critical:
- Too small, and progress is slow.
- Too large, and steps overshoot the valley floor and can diverge. For $f = x^2 + 10y^2$, the $y$-coordinate is multiplied by $1 - 20\alpha$ at each step, which grows in size when $\alpha > 0.1$.
- When the function is much steeper in some directions than others (long thin valleys, Section 12.1 level curves), the safe step is limited by the steepest direction, and progress along the shallow direction is slow.
Steepest descent is a local rule. It points downhill here and now, not toward the global minimum, which is why it zigzags in narrow valleys.
Examples
A hill has height $h(x, y) = 1000 - 0.01x^2 - 0.02y^2$ meters. A hiker stands at $(100, 50)$.
- $\nabla h = \langle -0.02x, -0.04y \rangle = \langle -2, -2 \rangle$.
- Steepest ascent direction: $\frac{1}{\sqrt{2}}\langle -1, -1 \rangle$, roughly toward the summit, at a slope of $\lvert \nabla h \rvert = 2\sqrt{2} \approx 2.83$ (meters up per meter walked).
- Walking in direction $\frac{1}{\sqrt{2}}\langle 1, -1 \rangle$, perpendicular to the gradient, the hiker stays at the same height.
Note that the steepest direction is not straight toward the summit at $(0, 0)$ (which would be direction $\langle -2, -1 \rangle$), because the hill is steeper in $y$.
import numpy as np
# Height of a hill: h(x, y) = 1000 - 0.01 x^2 - 0.02 y^2 (meters; x, y in meters)
grad_h = lambda x, y: np.array([-0.02 * x, -0.04 * y])
g = grad_h(100.0, 50.0) # hiker at (100, 50)
print(g, round(np.linalg.norm(g), 4), np.round(g / np.linalg.norm(g), 4))
# [-2. -2.] 2.8284 [-0.7071 -0.7071]
# Brute force over 3600 directions: the best one is the gradient direction, the best rate is |grad|
th = np.linspace(0, 2 * np.pi, 3600, endpoint=False)
U = np.stack([np.cos(th), np.sin(th)], axis=1)
rates = U @ g # D_u h = grad . u for every direction at once
k = np.argmax(rates)
print(round(rates.max(), 4), np.round(U[k], 4), round(rates.min(), 4))
# 2.8284 [-0.7071 -0.7071] -2.8284
# Gradient descent on f(x, y) = x^2 + 10 y^2 (minimum at the origin)
grad_f = lambda p: np.array([2 * p[0], 20 * p[1]])
for lr in [0.01, 0.09, 0.11]:
p = np.array([5.0, 1.0])
for _ in range(50):
p = p - lr * grad_f(p) # step downhill, against the gradient
print(lr, f'{p[0]:.4f} {p[1]:.4f}')
# 0.01 1.8208 0.0000 too timid: x has barely moved after 50 steps
# 0.09 0.0002 0.0000 well chosen: both coordinates converge
# 0.11 0.0000 9100.4382 too bold: |1 - 20 * 0.11| > 1, so y explodes
Key Takeaways
- $D_{\mathbf{u}} f = \lvert \nabla f \rvert\cos\theta$: the gradient points in the direction of steepest ascent, and its length is that steepest rate.
- $-\nabla f$ is the direction of steepest descent; directions perpendicular to $\nabla f$ give zero rate of change.
- Gradient descent $\mathbf{x} \leftarrow \mathbf{x} - \alpha\nabla f$ minimizes functions by repeated steepest-descent steps; the step size must respect the steepest direction.
- Steepest descent is local and can zigzag in elongated valleys.
🧪 Practice
- For $f(x, y) = 4 - x^2 - y^2$ at $(1, 2)$, find the direction of fastest increase, the maximum rate, and a direction of zero change.
- The temperature in a metal ball is $T(x, y, z) = \frac{200}{1 + x^2 + 2y^2 + 3z^2}$. In which direction does the temperature increase fastest at $(1, 1, -2)$? What is the maximum rate?
- Show that for gradient descent on $f(x) = \frac{c}{2}x^2$, each step multiplies $x$ by $1 - \alpha c$. For which step sizes does it converge? Apply this to each coordinate of $f(x, y) = x^2 + 10y^2$.
- Interview-style: A model's training loss decreases quickly for a few hundred steps, then oscillates without improving. A colleague suggests the gradient computation is buggy. What else could explain it, and what would you try first? Hint: a step size that is fine at first can be too large near the bottom of a narrow valley, causing bouncing across it. Try decaying the learning rate, rescaling the features, or using a method with per-coordinate step sizes (momentum, Adam); verify the gradient with a finite-difference check before blaming it.
Gradients Normal to Level Sets
Theory
The table above showed that directions perpendicular to $\nabla f$ give zero rate of change. But the directions of zero change are exactly the directions along a level curve, where $f$ stays constant. This gives one of the most useful facts in the chapter:
At each point where $\nabla f \ne \mathbf{0}$, the gradient is perpendicular (normal) to the level curve (or level surface) through that point.
Proof with the chain rule. Let $\mathbf{r}(t)$ be any differentiable curve lying in the level set $f = k$, passing through the point at $t = t_0$. Then $f(\mathbf{r}(t)) = k$ for all $t$. Differentiating with the chain rule (Section 12.4):
$$\nabla f(\mathbf{r}(t)) \cdot \mathbf{r}'(t) = 0.$$
So $\nabla f$ is perpendicular to the tangent vector $\mathbf{r}'(t_0)$ of every curve in the level set through the point. In two dimensions there is only one tangent direction; in three dimensions this means $\nabla f$ is perpendicular to the whole tangent plane.
On a contour map. The gradient at a point is an arrow crossing the contour line at a right angle, pointing toward higher values, and longer where the contours are crowded (steep). The paths of steepest ascent are curves that cross every contour at right angles. Water running downhill follows $-\nabla f$, which is why streams on a topographic map cross contour lines perpendicularly.
Level curves of f = x^2 + y^2 (circles) and gradient arrows grad f = <2x, 2y>
^
|
.------+------.
.-' '-.
/ .---------. \
| | * | +---> at (r, 0): grad f = <2r, 0>
\ '---------' /
'-. .-' every arrow crosses its circle at a right angle,
'------+------' points toward larger f, and has length 2r
|
v
Tangent line to a level curve. At a point $P = (a, b)$ on $f(x, y) = k$, the gradient is a normal vector, so the tangent line is
$$f_x(a, b)(x - a) + f_y(a, b)(y - b) = 0.$$
This is often easier than implicit differentiation, and it handles vertical tangents ($f_y = 0$) without special cases.
Examples
- $f = x^2 + 4y^2$, level curve $x^2 + 4y^2 = 4$ (an ellipse, parametrized as $\langle 2\cos t, \sin t \rangle$). The gradient $\langle 2x, 8y \rangle$ has zero dot product with the tangent $\langle -2\sin t, \cos t \rangle$ at every $t$: $-8\cos t\sin t + 8\sin t\cos t = 0$.
- Tangent line to that ellipse at $P = (\sqrt{2}, \frac{\sqrt{2}}{2})$: $\nabla f(P) = \langle 2\sqrt{2}, 4\sqrt{2} \rangle$, so the line is $2\sqrt{2}(x - \sqrt{2}) + 4\sqrt{2}(y - \frac{\sqrt{2}}{2}) = 0$, which simplifies to $x + 2y = 2\sqrt{2}$.
- At $(2, 0)$ on the same ellipse, $\nabla f = \langle 4, 0 \rangle$, and the tangent line is $4(x - 2) = 0$, the vertical line $x = 2$. Implicit differentiation, $y' = -\frac{x}{4y}$, would divide by zero there.
import numpy as np
# Level curve x^2 + 4y^2 = 4 (an ellipse), parametrized as r(t) = (2 cos t, sin t)
grad = lambda x, y: np.array([2 * x, 8 * y])
for t in np.linspace(0, 2 * np.pi, 7)[:-1]:
x, y = 2 * np.cos(t), np.sin(t)
tangent = np.array([-2 * np.sin(t), np.cos(t)]) # r'(t), tangent to the level curve
gx, gy = grad(x, y)
print(f'{t:5.3f} grad=({gx:+.3f}, {gy:+.3f}) grad . tangent = {grad(x, y) @ tangent:+.1e}')
# 0.000 grad=(+4.000, +0.000) grad . tangent = +0.0e+00
# 1.047 grad=(+2.000, +6.928) grad . tangent = +0.0e+00
# 2.094 grad=(-2.000, +6.928) grad . tangent = +0.0e+00
# 3.142 grad=(-4.000, +0.000) grad . tangent = +0.0e+00
# 4.189 grad=(-2.000, -6.928) grad . tangent = +0.0e+00
# 5.236 grad=(+2.000, -6.928) grad . tangent = +0.0e+00
# Tangent line to the level curve at (sqrt 2, sqrt 2 / 2): grad . (X - P) = 0
P = np.array([np.sqrt(2), np.sqrt(2) / 2])
n = grad(*P)
print(np.round(n, 4), np.round(n / n[0], 4), round(n @ P / n[0], 6)) # line: x + 2y = 2 sqrt 2
# [2.8284 5.6569] [1. 2.] 2.828427
Key Takeaways
- Where $\nabla f \ne \mathbf{0}$, the gradient is perpendicular to the level curve or level surface through the point.
- Proof: for any curve $\mathbf{r}(t)$ in the level set, $\frac{d}{dt}f(\mathbf{r}(t)) = \nabla f \cdot \mathbf{r}' = 0$.
- Steepest-ascent paths and streams of water cross contour lines at right angles.
- Tangent line to $f(x, y) = k$ at $(a, b)$: $f_x(a, b)(x - a) + f_y(a, b)(y - b) = 0$, with no trouble at vertical tangents.
🧪 Practice
- Find the tangent line to the hyperbola $x^2 - y^2 = 3$ at $(2, 1)$ using the gradient, and check it with implicit differentiation.
- Sketch the level curves of $f(x, y) = xy$ for $k = \pm 1, \pm 4$ and draw the gradient vector at $(1, 1)$, $(2, 2)$, and $(-1, 1)$. Verify that each is perpendicular to its level curve.
- A ball is released at $(1, 2)$ on the surface $z = x^2 + 3y^2$. In which direction in the $xy$-plane does it start to roll? Find the path it follows if it always moves in the direction of $-\nabla z$. (Solve $\mathbf{r}' = -\nabla z(\mathbf{r})$, two separate linear differential equations from Chapter 8.)
- Interview-style: In a 2D game, you have a signed distance function $d(x, y)$ for a wall (negative inside, positive outside, zero on the wall). A character has penetrated the wall slightly. How do you push it out along the shortest route, and how far? Hint: the wall is the level set $d = 0$, so $\nabla d$ is normal to it, and for a true distance function $\lvert \nabla d \rvert = 1$. Move the character by $-d(x, y),\nabla d(x, y)$.
Tangent Planes to Level Surfaces
Theory
Many surfaces are not graphs of a single function $z = f(x, y)$: spheres, ellipsoids, and closed shapes in general. They are naturally described as level surfaces $F(x, y, z) = k$. The previous topic showed that $\nabla F$ is perpendicular to every curve in the surface through a point. That makes it the normal vector of the tangent plane.
Tangent plane and normal line. At a point $P = (x_0, y_0, z_0)$ on $F(x, y, z) = k$ with $\nabla F(P) \ne \mathbf{0}$:
$$\text{Tangent plane:} \quad F_x(P)(x - x_0) + F_y(P)(y - y_0) + F_z(P)(z - z_0) = 0,$$
$$\text{Normal line:} \quad \mathbf{r}(t) = \langle x_0, y_0, z_0 \rangle + t,\nabla F(P).$$
Graphs are a special case. The graph $z = f(x, y)$ is the level surface $G = f(x, y) - z = 0$. Then $\nabla G = \langle f_x, f_y, -1 \rangle$, and the tangent plane formula reduces exactly to the one from Section 12.3: $z - z_0 = f_x(x - x_0) + f_y(y - y_0)$. The level-surface version is more general and treats all three variables symmetrically.
Applications.
| Problem | Tool |
|---|---|
| Tangent plane to a level surface | Normal vector $\nabla F(P)$ |
| Angle between two surfaces at a common point | Angle between their normals $\nabla F$ and $\nabla G$ |
| Surfaces tangent to each other at $P$ | $\nabla F(P)$ parallel to $\nabla G(P)$ |
| Surfaces orthogonal at $P$ | $\nabla F(P) \cdot \nabla G(P) = 0$ |
| Tangent line to the curve where $F = k_1$ and $G = k_2$ meet | Direction $\nabla F \times \nabla G$, perpendicular to both normals |
| Surface normals for lighting in computer graphics | $\frac{\nabla F}{\lvert \nabla F \rvert}$ for implicit surfaces |
The fourth and fifth rows return in Section 12.6: Lagrange multipliers are exactly the condition that a level curve of the objective is tangent to the constraint.
Examples
- Ellipsoid $\frac{x^2}{4} + y^2 + \frac{z^2}{9} = 3$ at $P = (-2, 1, -3)$. $\nabla F = \langle \frac{x}{2}, 2y, \frac{2z}{9} \rangle$, so $\nabla F(P) = \langle -1, 2, -\frac{2}{3} \rangle$. The tangent plane is
$$-(x + 2) + 2(y - 1) - \tfrac{2}{3}(z + 3) = 0 \quad \iff \quad 3x - 6y + 2z + 18 = 0,$$
and the normal line is $\langle -2 - t,\ 1 + 2t,\ -3 - \frac{2}{3}t \rangle$.
- The sphere $x^2 + y^2 + z^2 = 6$ and the paraboloid $z = x^2 + y^2$ meet at $(1, 1, 2)$. Their normals there are $\langle 2, 2, 4 \rangle$ and $\langle 2, 2, -1 \rangle$. The intersection curve has tangent direction $\langle 2, 2, 4 \rangle \times \langle 2, 2, -1 \rangle = \langle -10, 10, 0 \rangle$: horizontal, as expected, since the curve is the circle $x^2 + y^2 = 2$ at height $z = 2$. The surfaces meet at an angle of about $74.2^\circ$.
import numpy as np
import sympy as sp
x, y, z, t = sp.symbols('x y z t', real=True)
# Tangent plane and normal line to the ellipsoid x^2/4 + y^2 + z^2/9 = 3 at (-2, 1, -3)
F = x**2 / 4 + y**2 + z**2 / 9
P = {x: -2, y: 1, z: -3}
n = sp.Matrix([sp.diff(F, v) for v in (x, y, z)]).subs(P)
plane = sp.expand(3 * n.dot(sp.Matrix([x + 2, y - 1, z + 3]))) # scaled by 3 to clear fractions
print(n.T, plane, '= 0')
# Matrix([[-1, 2, -2/3]]) -3*x + 6*y - 2*z - 18 = 0
print((sp.Matrix([-2, 1, -3]) + t * n).T) # normal line
# Matrix([[-t - 2, 2*t + 1, -2*t/3 - 3]])
# Graph z = f(x, y) as a level surface G = f(x, y) - z = 0: same plane as Section 12.3
G = x**2 + y**2 - z
print(sp.Matrix([sp.diff(G, v) for v in (x, y, z)]).subs({x: 1, y: 2, z: 5}).T)
# Matrix([[2, 4, -1]]) normal <f_x, f_y, -1>
# Two surfaces: sphere x^2 + y^2 + z^2 = 6 and paraboloid z = x^2 + y^2 meet at (1, 1, 2)
g1 = np.array([2.0, 2.0, 4.0]) # grad of x^2 + y^2 + z^2
g2 = np.array([2.0, 2.0, -1.0]) # grad of x^2 + y^2 - z
print(np.cross(g1, g2)) # tangent direction of the intersection curve
# [-10. 10. 0.]
cos_angle = g1 @ g2 / (np.linalg.norm(g1) * np.linalg.norm(g2))
print(round(np.degrees(np.arccos(cos_angle)), 2)) # angle between the surfaces (between normals)
# 74.21
Key Takeaways
- For a level surface $F(x, y, z) = k$, $\nabla F(P)$ is normal to the tangent plane at $P$: $\nabla F(P) \cdot (\mathbf{x} - P) = 0$.
- The normal line is $P + t,\nabla F(P)$.
- The graph $z = f(x, y)$ is the level surface $f(x, y) - z = 0$, with normal $\langle f_x, f_y, -1 \rangle$.
- Two surfaces meet at the angle between their gradients; the curve where they intersect has tangent direction $\nabla F \times \nabla G$.
🧪 Practice
- Find the tangent plane and normal line to $x^2 + 2y^2 - 3z^2 = 3$ at $(2, -1, 1)$.
- Find the points on the ellipsoid $x^2 + 2y^2 + 3z^2 = 1$ where the tangent plane is parallel to the plane $3x - y + 3z = 1$.
- Show that the sphere $x^2 + y^2 + z^2 = 4$ and the cone $z^2 = x^2 + y^2$ meet at right angles at every point where they intersect.
- Interview-style: A ray tracer renders an implicit surface $F(x, y, z) = 0$ (for example, a blend of spheres) and needs the surface normal at each hit point for shading, but only has code that evaluates $F$. How would you compute the normals, and what can go wrong? Hint: the normal is $\frac{\nabla F}{\lvert \nabla F \rvert}$, which you can estimate by central differences in $x$, $y$, $z$ (six evaluations, or four with a tetrahedral stencil). Where $\nabla F \approx \mathbf{0}$ (creases, singular points of the surface), the normal is undefined and shading artifacts appear. <a id="126-multivariable-optimization"></a>
12.6 Multivariable Optimization
This subchapter finds maxima and minima of functions of several variables: locating critical points, classifying them with second derivatives, finding absolute extremes on closed regions, and optimizing under constraints with Lagrange multipliers.
Critical Points of Functions of Two Variables
Theory
At the top of a smooth hill, the ground is level in every direction. Walk east or north (or any other way) and, to first order, you neither climb nor descend. The same is true at the bottom of a valley. This is the multivariable version of Fermat's theorem (Section 4.2), and it is how we locate candidates for extremes.
Definitions. $f$ has a local maximum at $(a, b)$ if $f(x, y) \le f(a, b)$ for all $(x, y)$ in some disk around $(a, b)$, and a local minimum if $f(x, y) \ge f(a, b)$ there. If the inequality holds on the whole domain, the extremum is absolute (global).
Fermat's theorem in several variables. If $f$ has a local extremum at $(a, b)$ and its partial derivatives exist there, then
$$f_x(a, b) = 0 \quad \text{and} \quad f_y(a, b) = 0, \qquad \text{that is,} \qquad \nabla f(a, b) = \mathbf{0}.$$
Reason: the one-variable function $g(x) = f(x, b)$ has a local extremum at $x = a$, so $g'(a) = f_x(a, b) = 0$; likewise for $y$. If $f$ is differentiable, the tangent plane at an extremum is horizontal: $z = f(a, b)$.
Critical points. A point $(a, b)$ in the domain is a critical point if $\nabla f(a, b) = \mathbf{0}$, or if a partial derivative does not exist there. Local extremes can occur only at critical points, but not every critical point is an extreme: at a saddle point the gradient vanishes while $f$ rises in some directions and falls in others.
Finding them. Solve the system $f_x = 0$, $f_y = 0$ simultaneously. Unlike in one variable, this is a system of (usually nonlinear) equations. Useful strategies:
- Solve one equation for one variable and substitute into the other.
- Factor: an equation like $x(y - 1) = 0$ splits into the cases $x = 0$ or $y = 1$, and each case must be pursued.
- Look for symmetry: if swapping $x$ and $y$ leaves $f$ unchanged, try $y = x$ first, but do not assume all solutions have this form.
- Complete the square when $f$ is quadratic.
Examples
- $f(x, y) = x^2 + y^2 - 2x - 6y + 14$. Then $f_x = 2x - 2 = 0$ and $f_y = 2y - 6 = 0$ give the single critical point $(1, 3)$. Completing the square, $f = (x - 1)^2 + (y - 3)^2 + 4 \ge 4$, so $(1, 3)$ is an absolute minimum with value $4$.
- $f(x, y) = x^4 + y^4 - 4xy + 1$. The equations are $4x^3 - 4y = 0$ and $4y^3 - 4x = 0$. The first gives $y = x^3$; substituting, $x^9 = x$, so $x(x^8 - 1) = 0$, giving $x = 0, 1, -1$. The critical points are $(0, 0)$, $(1, 1)$, $(-1, -1)$. Sampling $f$ on a small circle around each shows that $(\pm 1, \pm 1)$ look like minima (values rise in every direction) but $(0, 0)$ is not an extreme (values both rise and fall).
- $f(x, y) = \sqrt{x^2 + y^2}$ has no point where $\nabla f = \mathbf{0}$, but its partials do not exist at the origin. The origin is a critical point, and it is the absolute minimum (the tip of a cone).
import numpy as np
import sympy as sp
x, y = sp.symbols('x y', real=True)
f = x**4 + y**4 - 4 * x * y + 1
grad = [sp.diff(f, v) for v in (x, y)]
crit = sp.solve(grad, [x, y], dict=True) # solve f_x = 0 and f_y = 0 together
real_crit = [c for c in crit if all(v.is_real for v in c.values())]
print([(c[x], c[y], f.subs(c)) for c in real_crit])
# [(-1, -1, -1), (0, 0, 1), (1, 1, -1)]
# What does f do near each critical point? Sample a small circle around it
F = sp.lambdify((x, y), f)
th = np.linspace(0, 2 * np.pi, 360, endpoint=False)
for c in real_crit:
a, b = float(c[x]), float(c[y])
ring = F(a + 0.05 * np.cos(th), b + 0.05 * np.sin(th)) - F(a, b)
print((a, b), f'{ring.min():+.5f} {ring.max():+.5f}') # both signs => neither max nor min
# (-1.0, -1.0) +0.00965 +0.02003 always higher nearby: looks like a local min
# (0.0, 0.0) -0.00500 +0.00500 higher and lower nearby: not an extremum
# (1.0, 1.0) +0.00965 +0.02003
Key Takeaways
- At a local extremum where the partials exist, $\nabla f = \mathbf{0}$: the tangent plane is horizontal.
- Critical points are where $\nabla f = \mathbf{0}$ or a partial derivative fails to exist; they are the only candidates for local extremes.
- Finding them means solving a system; factor, substitute, and follow every case.
- A critical point need not be an extremum: saddle points have $\nabla f = \mathbf{0}$ too.
🧪 Practice
- Find all critical points of $f(x, y) = x^2 + xy + y^2 + 3x - 3y + 4$.
- Find all critical points of $f(x, y) = xy(1 - x - y)$. (There are four.)
- Find the critical points of $f(x, y) = e^{-(x^2 + y^2)}(x^2 + 2y^2)$. Use symmetry and factoring.
- Interview-style: You need the critical points of a function of 50 variables, and solving $\nabla f = \mathbf{0}$ symbolically is hopeless. What numerical approaches are available, and what do they actually find? Hint: Newton's method applied to the system $\nabla f = \mathbf{0}$ (using the Hessian as its Jacobian) converges quickly to some nearby critical point, including saddles. Gradient descent tends to find minima instead. Neither guarantees finding all critical points or the global minimum; multiple random starts help.
Second Derivative Test and the Hessian
Theory
In one variable, the second derivative tells a minimum ($f'' > 0$, concave up) from a maximum ($f'' < 0$) at a critical point (Section 4.3). In two variables, the surface can curve up in one direction and down in another, so we need all the second partials together.
Why the second partials decide. Near a critical point $(a, b)$, where $\nabla f = \mathbf{0}$, the second-order Taylor expansion (Section 9.6, extended to two variables) has no linear term:
$$f(a + h, b + k) \approx f(a, b) + \tfrac{1}{2}\left(f_{xx}h^2 + 2f_{xy}hk + f_{yy}k^2\right).$$
The behavior is governed by the quadratic form $Q(h, k) = Ah^2 + 2Bhk + Ck^2$, with $A = f_{xx}$, $B = f_{xy}$, $C = f_{yy}$ at the point. Completing the square (when $A \ne 0$):
$$Q = A\left(h + \tfrac{B}{A}k\right)^2 + \frac{AC - B^2}{A}k^2.$$
If $AC - B^2 > 0$, both terms have the sign of $A$ and $Q$ never changes sign. If $AC - B^2 < 0$, the two terms have opposite signs and $Q$ takes both signs.
The Second Derivative Test. Suppose the second partials are continuous near $(a, b)$ and $\nabla f(a, b) = \mathbf{0}$. Let
$$D = f_{xx}(a, b),f_{yy}(a, b) - [f_{xy}(a, b)]^2.$$
| Condition | Conclusion | Picture |
|---|---|---|
| $D > 0$ and $f_{xx} > 0$ | Local minimum | Bowl: curves up in every direction |
| $D > 0$ and $f_{xx} < 0$ | Local maximum | Dome: curves down in every direction |
| $D < 0$ | Saddle point | Up in some directions, down in others |
| $D = 0$ | Inconclusive | Second-order terms do not decide; look at higher order |
When $D > 0$, $f_{xx}$ and $f_{yy}$ have the same sign, so either can be used.
The Hessian matrix. $D$ is the determinant of the Hessian, the matrix of second partials:
$$H_f = \begin{bmatrix} f_{xx} & f_{xy} \ f_{yx} & f_{yy} \end{bmatrix}, \qquad D = \det H_f.$$
By Clairaut's Theorem (Section 12.3) it is symmetric, so it has real eigenvalues $\lambda_1, \lambda_2$, with $D = \lambda_1\lambda_2$. The eigenvalues are the curvatures along the two principal directions (the eigenvectors), and the eigenvalue version of the test works in any number of variables:
| Eigenvalues of $H_f$ at a critical point | Name of $H_f$ | Conclusion |
|---|---|---|
| All positive | Positive definite | Local minimum |
| All negative | Negative definite | Local maximum |
| Some positive, some negative | Indefinite | Saddle point |
| Some zero, the rest of one sign | Semidefinite | Inconclusive |
The quadratic approximation $f(\mathbf{a} + \mathbf{h}) \approx f(\mathbf{a}) + \nabla f \cdot \mathbf{h} + \frac{1}{2}\mathbf{h}^{\mathsf{T}}H_f,\mathbf{h}$ is also the model behind Newton's method for optimization, which jumps to the critical point of this quadratic: $\mathbf{h} = -H_f^{-1}\nabla f$.
When $D = 0$. The test gives no information, and the point can be anything. For $x^4 + y^4$, $-(x^4 + y^4)$, and $x^4 - y^4$, every second partial vanishes at the origin, yet these are a minimum, a maximum, and a saddle respectively. Examine $f$ directly along lines and curves, or use higher-order terms.
Examples
- $f = x^4 + y^4 - 4xy + 1$ (previous topic): $f_{xx} = 12x^2$, $f_{yy} = 12y^2$, $f_{xy} = -4$.
- At $(0, 0)$: $D = 0 - 16 = -16 < 0$, a saddle point.
- At $(\pm 1, \pm 1)$: $D = 144 - 16 = 128 > 0$ and $f_{xx} = 12 > 0$, local minima with value $-1$.
- $f = x^3 + y^3 - 3xy$: critical points from $3x^2 = 3y$ and $3y^2 = 3x$ are $(0, 0)$ and $(1, 1)$. With $f_{xx} = 6x$, $f_{yy} = 6y$, $f_{xy} = -3$: $D(0, 0) = -9$ (saddle) and $D(1, 1) = 27$ with $f_{xx} = 6 > 0$ (local minimum, $f = -1$). It is not an absolute minimum, since $f(x, 0) = x^3 \to -\infty$.
import numpy as np
import sympy as sp
x, y = sp.symbols('x y', real=True)
def classify(f):
"""Second derivative test at every real critical point, plus the Hessian's eigenvalues."""
H = sp.hessian(f, (x, y)) # [[f_xx, f_xy], [f_yx, f_yy]]
for c in sp.solve([sp.diff(f, x), sp.diff(f, y)], [x, y], dict=True):
if not all(v.is_real for v in c.values()):
continue
Hc = H.subs(c)
D = Hc.det() # f_xx f_yy - f_xy^2
if D > 0:
kind = 'local min' if Hc[0, 0] > 0 else 'local max'
elif D < 0:
kind = 'saddle'
else:
kind = 'inconclusive'
eig = np.linalg.eigvalsh(np.array(Hc, dtype=float)).round(3).tolist()
print((c[x], c[y]), 'D =', D, kind, eig)
classify(x**4 + y**4 - 4 * x * y + 1)
# (-1, -1) D = 128 local min [8.0, 16.0]
# (0, 0) D = -16 saddle [-4.0, 4.0] eigenvalues of opposite signs
# (1, 1) D = 128 local min [8.0, 16.0] both eigenvalues positive
classify(x**3 + y**3 - 3 * x * y)
# (0, 0) D = -9 saddle [-3.0, 3.0]
# (1, 1) D = 27 local min [3.0, 9.0]
for g in [x**4 + y**4, -(x**4 + y**4), x**4 - y**4]: # D = 0 at the origin for all three
classify(g)
# (0, 0) D = 0 inconclusive [0.0, 0.0] actually a minimum
# (0, 0) D = 0 inconclusive [0.0, 0.0] actually a maximum
# (0, 0) D = 0 inconclusive [0.0, 0.0] actually a saddle
Key Takeaways
- At a critical point, $D = f_{xx}f_{yy} - f_{xy}^2$ decides: $D > 0$ with $f_{xx} > 0$ is a minimum, $D > 0$ with $f_{xx} < 0$ is a maximum, $D < 0$ is a saddle, $D = 0$ is inconclusive.
- $D$ is the determinant of the Hessian; the signs of the Hessian's eigenvalues give the same test in any dimension.
- The test comes from the second-order Taylor expansion, whose quadratic part dominates near a critical point.
- When $D = 0$, examine the function directly; anything is possible.
🧪 Practice
- Find and classify the critical points of $f(x, y) = x^3 - 12xy + 8y^3$.
- Find and classify the critical points of $f(x, y) = (x^2 + y^2)e^{y^2 - x^2}$.
- Show that $f(x, y) = x^2 + kxy + y^2$ has a minimum at the origin when $\lvert k \rvert < 2$ and a saddle when $\lvert k \rvert > 2$. What happens at $\lvert k \rvert = 2$?
- Interview-style: After training, you want to check whether a model's parameters sit at a true local minimum rather than a saddle point, but the Hessian has $10^8$ rows. How could you estimate its smallest eigenvalue without forming it? Hint: Hessian-vector products $H\mathbf{v}$ can be computed at about the cost of two gradient evaluations (by differentiating $\nabla f \cdot \mathbf{v}$, or by finite differences of gradients). Iterative methods such as the Lanczos algorithm or power iteration on a shifted matrix need only these products to estimate extreme eigenvalues.
Saddle Points
Theory
A mountain pass is the classic picture: walking along the trail that crosses the ridge, the pass is the highest point; walking along the ridge from one peak to the other, it is the lowest point. At the pass itself the ground is level, so $\nabla f = \mathbf{0}$, but it is neither a maximum nor a minimum. The surface is shaped like a horse saddle, which gives the name.
Definition. A critical point $(a, b)$ with $\nabla f(a, b) = \mathbf{0}$ is a saddle point if every disk around it contains points where $f > f(a, b)$ and points where $f < f(a, b)$.
Recognizing saddles.
| Evidence | Why it shows a saddle |
|---|---|
| $D < 0$ in the second derivative test | The quadratic part takes both signs |
| Hessian has a positive and a negative eigenvalue | $f$ curves up along one eigenvector and down along the other |
| Two lines (or curves) through the point, $f$ has a min along one and a max along the other | Directly from the definition |
| Level curves cross at the point (an X in the contour map, Section 12.1) | The level curve at the critical value splits the neighborhood into regions above and below |
The standard saddles.
- $z = x^2 - y^2$: up along the $x$-axis, down along the $y$-axis. Level curves are hyperbolas, and the level curve $z = 0$ is the pair of lines $y = \pm x$, crossing at the saddle.
- $z = xy$: the same surface rotated by $45^\circ$. Up along $y = x$, down along $y = -x$.
- Monkey saddle $z = x^3 - 3xy^2$: three ups and three downs around the origin (room for two legs and a tail). Its Hessian at the origin is zero, so $D = 0$ and the test is inconclusive; the saddle is "degenerate" and is identified from $f$ in polar form, $r^3\cos 3\theta$.
Saddle z = x^2 - y^2 near the origin
y-axis: z = -y^2 (a maximum along this line)
|
- - | - - + : f > 0
- | - - : f < 0
+ - | - +
+ \|/ + lines y = +x and y = -x:
------+--------*--------+------ x the level curve f = 0, crossing in an X
+ /|\ +
+ - | - + x-axis: z = x^2 (a minimum along this line)
- | -
- - | - -
Saddles in high dimensions. For a function of $n$ variables, a critical point is a minimum only if all $n$ Hessian eigenvalues are positive. For a "random" critical point, the chance that all $n$ signs come out positive shrinks rapidly with $n$, so most critical points of high-dimensional functions are saddles. This is the situation in training neural networks: the danger is less "getting stuck in a bad local minimum" and more "slowing down near saddles", where the gradient is tiny.
Escaping saddles. A saddle is an unstable equilibrium for gradient descent. Started exactly on the "uphill" axis, descent converges to the saddle, but any tiny perturbation along a negative-curvature direction grows geometrically, and the iterate eventually slides off. Random noise (as in stochastic gradient descent) helps.
Examples
- $f(x, y) = x^2 - y^2$: $\nabla f = \langle 2x, -2y \rangle = \mathbf{0}$ only at the origin; $f_{xx} = 2$, $f_{yy} = -2$, $f_{xy} = 0$, so $D = -4 < 0$. Saddle. The eigenvectors of the Hessian are the axes: the directions of steepest upward and downward curvature.
- $f(x, y) = xy$: $D = 0\cdot 0 - 1^2 = -1 < 0$. Saddle. Hessian eigenvalues $\pm 1$ with eigenvectors along $y = \pm x$.
- Gradient descent on $x^2 - y^2$ started at $(1, 10^{-8})$: the $x$-coordinate shrinks by a factor $0.8$ per step (toward the saddle), while the $y$-coordinate grows by a factor $1.2$ per step (away from it). After about 100 steps, the tiny perturbation has grown to order $1$ and the iterate has left the saddle.
import numpy as np
# Saddle z = x^2 - y^2: gradient descent started almost exactly at the saddle point
grad = lambda p: np.array([2 * p[0], -2 * p[1]])
p = np.array([1.0, 1e-8]) # tiny offset in the "downhill" y direction
for k in range(1, 101):
p = p - 0.1 * grad(p)
if k in (10, 50, 100):
print(k, f'{p[0]:.2e} {p[1]:.2e}') # x shrinks to 0, then y grows: it escapes
# 10 1.07e-01 6.19e-08
# 50 1.43e-05 9.10e-05 lingering near the saddle: gradient is tiny, progress is slow
# 100 2.04e-10 8.28e-01 escaped along the negative-curvature direction
# In many dimensions, a "random" critical point is almost never a minimum:
# fraction of random symmetric Hessians whose eigenvalues are all positive
rng = np.random.default_rng(0)
for n in [1, 2, 3, 5, 8]:
trials, hits = 20000, 0
for _ in range(trials):
A = rng.normal(size=(n, n))
H = (A + A.T) / 2 # random symmetric matrix
hits += np.all(np.linalg.eigvalsh(H) > 0)
print(n, hits / trials)
# 1 0.498
# 2 0.145
# 3 0.02715
# 5 0.0001
# 8 0.0 already at n = 8, essentially every random critical point is a saddle
Key Takeaways
- A saddle point is a critical point where $f$ takes both larger and smaller values arbitrarily close by.
- Detect it with $D < 0$, with Hessian eigenvalues of both signs, or with two directions showing a min and a max; level curves cross there.
- Degenerate saddles like the monkey saddle $x^3 - 3xy^2$ have $D = 0$ and must be analyzed directly.
- In high dimensions most critical points are saddles; gradient descent slows near them but escapes along negative-curvature directions.
🧪 Practice
- Show that $f(x, y) = x^2 + 4xy + y^2$ has a saddle point at the origin, and find the two lines through the origin along which $f$ is identically $0$.
- Write the monkey saddle $x^3 - 3xy^2$ in polar coordinates and use the result to describe how many times $f$ changes sign on a small circle around the origin.
- For $f(x, y) = x^2 - y^2$, find all starting points from which gradient descent with step size $0.1$ converges to the saddle point.
- Interview-style: A colleague claims that their optimizer has converged because the gradient norm dropped below $10^{-6}$. Why is that not enough to claim they found a minimum, and what cheap checks would you add? Hint: a tiny gradient also occurs near saddles and on flat plateaus. Check the loss along a few random directions (or along a negative-curvature direction found with Hessian-vector products) to see whether it decreases, and restart from perturbed points to see whether the optimizer returns.
Absolute Extrema on Closed Regions
Theory
In one variable, a continuous function on a closed interval $[a, b]$ attains an absolute maximum and minimum, and you find them by comparing the values at critical points and endpoints (the Closed Interval Method, Section 4.2). In two variables, the interval becomes a region, and the two endpoints become an entire boundary curve.
When are extremes guaranteed? A set $D \subseteq \mathbb{R}^2$ is
- closed if it contains all its boundary points (a disk with its circle, a rectangle with its edges);
- bounded if it fits inside some disk (it does not run off to infinity).
Extreme Value Theorem. If $f$ is continuous on a closed, bounded set $D$, then $f$ attains an absolute maximum value and an absolute minimum value at some points of $D$.
Drop either hypothesis and the guarantee fails: $f = x$ on the open disk has no maximum (it approaches $1$ but never reaches it), and $f = x$ on the whole plane is unbounded.
The closed region method. For continuous $f$ on a closed, bounded region $D$:
- Find the values of $f$ at the critical points in the interior of $D$. (No need to classify them.)
- Find the extreme values of $f$ on the boundary of $D$. This is usually a one-variable problem: parametrize each piece of the boundary and use the Closed Interval Method on it, including the endpoints of each piece (the corners).
- The largest value from steps 1 and 2 is the absolute maximum; the smallest is the absolute minimum.
| Boundary type | How to reduce it to one variable |
|---|---|
| Straight edge, like $y = 0$, $0 \le x \le 3$ | Substitute: $g(x) = f(x, 0)$ on $[0, 3]$ |
| Circle $x^2 + y^2 = r^2$ | Parametrize $x = r\cos t$, $y = r\sin t$, $t \in [0, 2\pi]$; or substitute $y^2 = r^2 - x^2$ |
| Curve $y = h(x)$ | Substitute $g(x) = f(x, h(x))$ |
| Smooth curve $g(x, y) = c$ | Lagrange multipliers (next topic) |
The most common mistake is forgetting the boundary, where the extremes very often occur (for example, for any non-constant linear function, both extremes lie on the boundary).
Examples
Find the absolute extremes of $f(x, y) = x^2 - 2xy + 2y$ on the rectangle $D = {0 \le x \le 3,\ 0 \le y \le 2}$.
- Interior: $f_x = 2x - 2y = 0$ and $f_y = -2x + 2 = 0$ give $(1, 1)$, with $f(1, 1) = 1$.
- Boundary, edge by edge:
| Edge | Restriction | Range on the edge |
|---|---|---|
| $y = 0$, $0 \le x \le 3$ | $f = x^2$ | min $0$ at $(0, 0)$, max $9$ at $(3, 0)$ |
| $x = 3$, $0 \le y \le 2$ | $f = 9 - 4y$ | min $1$ at $(3, 2)$, max $9$ at $(3, 0)$ |
| $y = 2$, $0 \le x \le 3$ | $f = x^2 - 4x + 4 = (x - 2)^2$ | min $0$ at $(2, 2)$, max $4$ at $(0, 2)$ |
| $x = 0$, $0 \le y \le 2$ | $f = 2y$ | min $0$ at $(0, 0)$, max $4$ at $(0, 2)$ |
- Compare: the absolute maximum is $9$ at $(3, 0)$; the absolute minimum is $0$, attained at both $(0, 0)$ and $(2, 2)$. The interior critical point (value $1$) is neither.
y
2 +-----------*-----------+ * (2, 2): min 0, inside an edge
| (0,2): 4 | (3,2): 1
| |
1 + o (1, 1): 1 | o interior critical point (a saddle)
| |
| |
0 *-----------------------*-- x * (0, 0): min 0 * (3, 0): MAX 9
0 2 3
import numpy as np
import sympy as sp
x, y = sp.symbols('x y', real=True)
f = x**2 - 2 * x * y + 2 * y # on the rectangle 0 <= x <= 3, 0 <= y <= 2
cands = []
# 1. Interior critical points
for c in sp.solve([sp.diff(f, x), sp.diff(f, y)], [x, y], dict=True):
if 0 < c[x] < 3 and 0 < c[y] < 2:
cands.append(('interior', c[x], c[y]))
# 2. Each edge becomes a one-variable problem: critical points of the restriction
edges = [('y = 0', {y: 0}, x, (0, 3)), ('y = 2', {y: 2}, x, (0, 3)),
('x = 0', {x: 0}, y, (0, 2)), ('x = 3', {x: 3}, y, (0, 2))]
for name, fix, var, (lo, hi) in edges:
g = f.subs(fix)
for r in sp.solve(sp.diff(g, var), var):
if lo < r < hi:
pt = {**fix, var: r}
cands.append((name, pt[x], pt[y]))
# 3. Corners
cands += [('corner', a, b) for a in (0, 3) for b in (0, 2)]
for kind, a, b in cands:
print(f'{kind:9s} ({a}, {b}) f = {f.subs({x: a, y: b})}')
# interior (1, 1) f = 1
# y = 2 (2, 2) f = 0
# corner (0, 0) f = 0
# corner (0, 2) f = 4
# corner (3, 0) f = 9
# corner (3, 2) f = 1
# Brute-force check on a fine grid
F = sp.lambdify((x, y), f)
X, Y = np.meshgrid(np.linspace(0, 3, 601), np.linspace(0, 2, 401))
Z = F(X, Y)
i = Z.argmax()
print(Z.max(), X.flat[i], Y.flat[i], Z.min()) # max value, where it occurs, min value
# 9.0 3.0 0.0 0.0
Key Takeaways
- A continuous function on a closed and bounded set always attains an absolute maximum and minimum.
- Candidates: interior critical points, critical points of $f$ restricted to each boundary piece, and corners.
- Reduce each boundary piece to a one-variable problem by substitution or parametrization.
- Evaluate $f$ at every candidate and compare; no classification is needed.
🧪 Practice
- Find the absolute maximum and minimum of $f(x, y) = x^2 + y^2 - 2x$ on the closed triangle with vertices $(2, 0)$, $(0, 2)$, and $(0, -2)$.
- Find the absolute extremes of $f(x, y) = x^2 + 2y^2 - x$ on the closed disk $x^2 + y^2 \le 1$, parametrizing the boundary circle.
- Explain why $f(x, y) = \frac{1}{x^2 + y^2}$ has no absolute maximum on the set $0 < x^2 + y^2 \le 1$. Which hypothesis of the Extreme Value Theorem fails?
- Interview-style: A warehouse robot's battery drain is a smooth function $E(x, y)$ of its position on a rectangular floor. Your teammate finds the minimum drain by running gradient descent from the center of the floor and clipping positions to stay inside the rectangle. When can this miss the true minimum, and how would you make the search reliable? Hint: gradient descent finds one local minimum; the true minimum may be in a different basin or on an edge or corner. Check the interior critical points, the minimum along each edge (a one-variable problem), and the corners, or run the descent from many starting points.
Lagrange Multipliers
Theory
Many optimization problems have a constraint: maximize profit within a budget, minimize material for a fixed volume, find the point on a curve closest to the origin. Sometimes you can solve the constraint for one variable and substitute, but often that is messy or impossible. Lagrange multipliers handle the constraint without solving it.
The geometric idea. Maximize $f(x, y)$ subject to $g(x, y) = k$. Picture the contour map of $f$ with the constraint curve $g = k$ drawn on top. Walk along the constraint curve and watch the values of $f$. Where the curve crosses a level curve of $f$, you can keep going in one direction to increase $f$, so that point is not a maximum. At a maximum (or minimum), the constraint curve must touch a level curve of $f$ without crossing it: the two curves are tangent.
Maximize f = x on a closed constraint curve g = k: level curves of f are vertical lines
f = 3 f = 4 f = 5
: : :
: : :
---------:---. : : the constraint curve crosses f = 3,
: '--. : : so moving right along it still raises f ...
: '-P : ... until P, where it only touches f = 4:
: .--': : the two curves share a tangent line there,
: .-' : : so grad f and grad g are parallel at P
---------:--' : :
: : :
Two curves are tangent where their normal vectors are parallel. The normal to a level curve of $f$ is $\nabla f$, and the normal to $g = k$ is $\nabla g$ (Section 12.5). So:
Method of Lagrange multipliers. If $f$ has an extreme value on the curve $g(x, y) = k$ at a point where $\nabla g \ne \mathbf{0}$, then at that point
$$\nabla f = \lambda,\nabla g \quad \text{for some number } \lambda.$$
To find the candidates, solve the system of three equations in three unknowns $x, y, \lambda$:
$$f_x = \lambda g_x, \qquad f_y = \lambda g_y, \qquad g(x, y) = k.$$
Then evaluate $f$ at every solution. If the constraint set is closed and bounded (a circle, an ellipse), the Extreme Value Theorem guarantees that the largest value is the maximum and the smallest is the minimum. The number $\lambda$ is the Lagrange multiplier. The same method works in three or more variables: $\nabla f = \lambda\nabla g$ with $g(x, y, z) = k$ gives four equations.
An equivalent packaging is the Lagrangian $\mathcal{L}(x, y, \lambda) = f(x, y) - \lambda,(g(x, y) - k)$. Setting all its partial derivatives to zero reproduces the system above, which turns a constrained problem into an unconstrained critical point problem.
The meaning of $\lambda$. If $M(k)$ is the optimal value of $f$ for constraint level $k$, then
$$\frac{dM}{dk} = \lambda.$$
So $\lambda$ measures how much the optimum improves per unit of loosening the constraint. Economists call it the shadow price or marginal value of the constrained resource: if the budget is the constraint, $\lambda$ is the extra output from one more dollar.
Pitfalls.
- Points where $\nabla g = \mathbf{0}$ are not covered by the method; check them separately.
- Solving the system is the hard part. Typical moves: divide equations to eliminate $\lambda$ (watch for division by zero), or consider cases like $x = 0$ versus $x \ne 0$.
- On an unbounded constraint curve, a solution of the system need not be a maximum or a minimum.
Examples
-
Extremes of $f = x^2 + 2y^2$ on the unit circle $x^2 + y^2 = 1$. The equations are $2x = 2\lambda x$, $4y = 2\lambda y$, $x^2 + y^2 = 1$. From the first, $x = 0$ or $\lambda = 1$.
- $x = 0$: then $y = \pm 1$, $\lambda = 2$, and $f = 2$.
- $\lambda = 1$: then $4y = 2y$, so $y = 0$, $x = \pm 1$, and $f = 1$.
The maximum is $2$ at $(0, \pm 1)$ and the minimum is $1$ at $(\pm 1, 0)$. Geometrically, the level curves of $f$ are ellipses, and the unit circle touches the ellipse $f = 1$ at its ends and the ellipse $f = 2$ at its top and bottom.
-
A factory's output is $P = 100x^{3/4}y^{1/4}$ (labor $x$ costs $150$ per unit, capital $y$ costs $250$ per unit), with budget $150x + 250y = B$. Lagrange's conditions give $x = \frac{B}{200}$ and $y = \frac{B}{1000}$: spend $75%$ of the budget on labor and $25%$ on capital, matching the exponents. For $B = 50{,}000$: $x = 250$, $y = 50$, $P \approx 16{,}718.5$, and $\lambda \approx 0.334$. One more dollar of budget yields about $0.334$ more units of output, as the code confirms.
import sympy as sp
x, y, lam, B = sp.symbols('x y lambda B', real=True)
# Extremes of f = x^2 + 2y^2 on the unit circle g = x^2 + y^2 = 1
f = x**2 + 2 * y**2
g = x**2 + y**2
eqs = [sp.diff(f, x) - lam * sp.diff(g, x), # grad f = lambda grad g ...
sp.diff(f, y) - lam * sp.diff(g, y),
g - 1] # ... and stay on the constraint
for s in sp.solve(eqs, [x, y, lam], dict=True):
print((s[x], s[y]), 'lambda =', s[lam], 'f =', f.subs(s))
# (-1, 0) lambda = 1 f = 1
# (1, 0) lambda = 1 f = 1
# (0, -1) lambda = 2 f = 2
# (0, 1) lambda = 2 f = 2
# Cobb-Douglas production P = 100 x^(3/4) y^(1/4) with budget 150 x + 250 y = B
P = 100 * x**sp.Rational(3, 4) * y**sp.Rational(1, 4)
budget = 150 * x + 250 * y
sol = sp.solve([sp.diff(P, x) - lam * 150, sp.diff(P, y) - lam * 250, budget - B], [x, y, lam], dict=True)[0]
print(sp.simplify(sol[x]), sp.simplify(sol[y]))
# B/200 B/1000
B0 = 50000
best = lambda b: P.subs({x: sol[x], y: sol[y]}).subs(B, b)
print(sp.N(best(B0), 7), sp.N(sol[lam].subs(B, B0), 6))
# 16718.51 0.334370
print(sp.N(best(B0 + 1) - best(B0), 6)) # one more dollar of budget buys about lambda more units
# 0.334370
Key Takeaways
- At a constrained extremum of $f$ on $g = k$ (with $\nabla g \ne \mathbf{0}$), the level curve of $f$ is tangent to the constraint: $\nabla f = \lambda\nabla g$.
- Solve $\nabla f = \lambda\nabla g$ together with $g = k$, then compare the values of $f$ at all solutions.
- On a closed, bounded constraint set, the largest and smallest values are the absolute maximum and minimum.
- $\lambda = \frac{dM}{dk}$: the rate at which the optimal value changes as the constraint is relaxed (the shadow price).
🧪 Practice
- Find the maximum and minimum of $f(x, y) = xy$ on the ellipse $x^2 + 4y^2 = 8$.
- Find the point on the plane $x + 2y + 3z = 14$ closest to the origin, by minimizing the squared distance with a Lagrange multiplier. Check your answer with the normal vector of the plane.
- A closed rectangular box must have volume $32$ m³. Use Lagrange multipliers to find the dimensions that minimize its surface area, then interpret $\lambda$.
- Interview-style: In a regularized model, you can either minimize the loss subject to $\lvert \mathbf{w} \rvert^2 \le c$ or minimize the loss plus $\alpha\lvert \mathbf{w} \rvert^2$ with no constraint. Why are these two formulations related, and how does $\alpha$ correspond to $c$? Hint: at a solution where the constraint is active, $\nabla(\text{loss}) = -2\lambda\mathbf{w}$, which is exactly the critical point condition of loss $+ \lambda\lvert \mathbf{w} \rvert^2$. So $\alpha$ plays the role of the multiplier, and a tighter budget $c$ corresponds to a larger $\alpha$.
Multiple Constraints
Theory
With two constraints in three variables, $g(x, y, z) = k_1$ and $h(x, y, z) = k_2$, the feasible set is (usually) a curve: the intersection of two surfaces. Think of a plane slicing a cylinder to give an ellipse. We want the extremes of $f$ along this curve.
The geometric idea. At an extreme point $P$ on the curve, $f$ cannot change to first order as you move along the curve, so $\nabla f$ is perpendicular to the curve's tangent direction $\mathbf{T}$. The tangent direction is perpendicular to both $\nabla g$ and $\nabla h$ (it lies in both surfaces; Section 12.5 gave $\mathbf{T} \parallel \nabla g \times \nabla h$). The vectors perpendicular to $\mathbf{T}$ form the plane spanned by $\nabla g$ and $\nabla h$. Therefore:
$$\nabla f = \lambda,\nabla g + \mu,\nabla h,$$
for some numbers $\lambda$ and $\mu$, provided $\nabla g$ and $\nabla h$ are not parallel at $P$. Together with the two constraints, this gives five equations in the five unknowns $x, y, z, \lambda, \mu$:
$$f_x = \lambda g_x + \mu h_x, \quad f_y = \lambda g_y + \mu h_y, \quad f_z = \lambda g_z + \mu h_z, \quad g = k_1, \quad h = k_2.$$
General pattern. For $m$ equality constraints in $n$ variables ($m < n$), at an extremum $\nabla f = \lambda_1\nabla g_1 + \cdots + \lambda_m\nabla g_m$, provided the constraint gradients are linearly independent. The Lagrangian becomes $\mathcal{L} = f - \sum_i \lambda_i(g_i - k_i)$, and each multiplier is the sensitivity of the optimum to its own constraint, $\lambda_i = \frac{\partial M}{\partial k_i}$.
| Setting | Feasible set | Condition | Unknowns |
|---|---|---|---|
| 2 variables, 1 constraint | Curve | $\nabla f = \lambda\nabla g$ | $x, y, \lambda$ |
| 3 variables, 1 constraint | Surface | $\nabla f = \lambda\nabla g$ | $x, y, z, \lambda$ |
| 3 variables, 2 constraints | Curve | $\nabla f = \lambda\nabla g + \mu\nabla h$ | $x, y, z, \lambda, \mu$ |
| $n$ variables, $m$ constraints | $(n - m)$-dimensional set | $\nabla f = \sum \lambda_i\nabla g_i$ | $n + m$ |
Inequality constraints (a preview). Real problems often have constraints like $g(\mathbf{x}) \le k$ (a budget you do not have to spend fully). At the optimum, either the constraint is inactive ($g < k$, and the point is an ordinary critical point, multiplier $0$) or active ($g = k$, and the Lagrange condition holds with a multiplier of the correct sign). These are the Karush-Kuhn-Tucker (KKT) conditions, the foundation of constrained numerical optimization and of support vector machines. The closed region method from earlier in this subchapter is the hands-on version: interior critical points (inactive) plus boundary (active).
Examples
Find the maximum of $f(x, y, z) = x + 2y + 3z$ on the curve where the plane $x - y + z = 1$ meets the cylinder $x^2 + y^2 = 1$.
With $g = x - y + z$ and $h = x^2 + y^2$, the condition $\nabla f = \lambda\nabla g + \mu\nabla h$ reads
$$1 = \lambda + 2\mu x, \qquad 2 = -\lambda + 2\mu y, \qquad 3 = \lambda.$$
So $\lambda = 3$, then $x = -\frac{1}{\mu}$ and $y = \frac{5}{2\mu}$. The cylinder gives $\frac{1}{\mu^2} + \frac{25}{4\mu^2} = 1$, so $\mu = \pm\frac{\sqrt{29}}{2}$, and $x = \mp\frac{2}{\sqrt{29}}$, $y = \pm\frac{5}{\sqrt{29}}$, $z = 1 - x + y = 1 \pm \frac{7}{\sqrt{29}}$. Then $f = 3 \pm \sqrt{29}$. The curve is a closed ellipse, so the maximum is $3 + \sqrt{29} \approx 8.385$ and the minimum is $3 - \sqrt{29} \approx -2.385$.
import numpy as np
import sympy as sp
from scipy.optimize import minimize
x, y, z, lam, mu = sp.symbols('x y z lambda mu', real=True)
# Extremes of f = x + 2y + 3z on the curve where the plane x - y + z = 1 meets the cylinder x^2 + y^2 = 1
f = x + 2 * y + 3 * z
g = x - y + z - 1
h = x**2 + y**2 - 1
L = f - lam * g - mu * h # Lagrangian: grad f = lambda grad g + mu grad h
eqs = [sp.diff(L, v) for v in (x, y, z)] + [g, h]
for s in sp.solve(eqs, [x, y, z, lam, mu], dict=True):
print([sp.nsimplify(s[v]) for v in (x, y, z)], 'f =', sp.simplify(f.subs(s)), '~', round(float(f.subs(s)), 4))
# [2*sqrt(29)/29, -5*sqrt(29)/29, 1 - 7*sqrt(29)/29] f = 3 - sqrt(29) ~ -2.3852
# [-2*sqrt(29)/29, 5*sqrt(29)/29, 1 + 7*sqrt(29)/29] f = 3 + sqrt(29) ~ 8.3852
# The same maximum from a general-purpose constrained optimizer
res = minimize(lambda p: -(p[0] + 2 * p[1] + 3 * p[2]), x0=[1.0, 0.0, 0.0], method='SLSQP',
constraints=[{'type': 'eq', 'fun': lambda p: p[0] - p[1] + p[2] - 1},
{'type': 'eq', 'fun': lambda p: p[0]**2 + p[1]**2 - 1}])
print(np.round(res.x, 4), round(-res.fun, 4), round(3 + np.sqrt(29), 4))
# [-0.3714 0.9285 2.2999] 8.3852 8.3852
Key Takeaways
- With two constraints, the extremes of $f$ on the intersection curve satisfy $\nabla f = \lambda\nabla g + \mu\nabla h$, plus both constraint equations.
- Geometrically, $\nabla f$ must be perpendicular to the curve, so it lies in the plane spanned by the two constraint normals.
- With $m$ constraints, add one multiplier per constraint; each multiplier is the sensitivity of the optimum to its constraint.
- Inequality constraints lead to the KKT conditions: each constraint is either inactive (multiplier $0$) or active (a Lagrange condition).
🧪 Practice
- Find the extreme values of $f(x, y, z) = z$ on the curve of intersection of the cone $x^2 + y^2 = z^2$ and the plane $x + y + z = 24$.
- Find the highest and lowest points on the ellipse where the plane $x + y + 2z = 2$ meets the paraboloid $z = x^2 + y^2$.
- Maximize $f(x, y, z) = x + y + z$ subject to $x^2 + y^2 + z^2 \le 1$ and $z \le \frac{1}{2}$. Decide which constraints are active at the optimum.
- Interview-style: A portfolio optimizer minimizes risk $\mathbf{w}^{\mathsf{T}}\Sigma\mathbf{w}$ subject to two constraints: the weights sum to $1$, and the expected return $\mathbf{r} \cdot \mathbf{w}$ equals a target $R$. How would you solve it, and what do the two multipliers tell a portfolio manager? Hint: the Lagrange conditions $2\Sigma\mathbf{w} = \lambda\mathbf{1} + \mu\mathbf{r}$ plus the two constraints form a linear system, solvable directly. $\mu$ is the marginal risk of demanding one more unit of return (the slope of the efficient frontier), and $\lambda$ is the marginal risk of investing one more unit of capital.
<a id="13-multiple-integrals"></a>