The Chain Rule in Several Variables

Suppose a factory's output depends on capital and labor. Both capital and labor drift with time, following some economic policy. How fast does output change with time? In one-variable calculus, the answer is a single multiplication of derivatives. In several variables, the answer is a matrix product. Understanding why requires taking seriously what a derivative in several variables actually is. It is not a number, but a linear map. This map best approximates a function near a point. The chain rule in several variables is a key theorem. It states that composition of differentiable maps is differentiable. It also states that the derivative of the composition is the composition of the derivatives. Equivalently, this composition is the matrix product of the derivatives. It is the single most-used computational tool in multivariable calculus. It justifies implicit differentiation. It underlies change of variables in partial differential equations. It powers backpropagation in neural networks. It is also the technical engine behind the Implicit Function Theorem and Taylor's Theorem in Several Variables. The subtlety is this. In one variable, "derivative" is a scalar. The chain rule is (g∘f)′(x)=g′(f(x))f′(x)(g\circ f)'(x) = g'(f(x))f'(x), ordinary multiplication. In Rn\mathbb{R}^n, the derivative Df(a)Df(a) of ff at a point aa is different. It is a linear transformation. Once bases are fixed, it is represented by the Jacobian Matrix. Composing two differentiable functions composes their best linear approximations. Composition of linear maps corresponds to matrix multiplication. This is not an analogy; it is a theorem. Its proof is a careful piece of ε\varepsilon–δ\delta bookkeeping. It is built directly on the definition of the Total Derivative. This article states the precise hypotheses needed. These require differentiability at a point, not merely existence of partial derivatives. It proves the theorem in full using the total-derivative definition. It works two genuinely multivariable examples. It shows exactly what fails if we relax "differentiable" to "partial derivatives exist."

Key takeaways
  • The multivariable chain rule says D(g∘f)(a)=Dg(f(a))∘Df(a)D(g\circ f)(a) = Dg(f(a)) \circ Df(a), a composition of linear maps.
  • In coordinates this composition is the matrix product of the Jacobian matrices: Jg∘f(a)=Jg(f(a)) Jf(a)J_{g\circ f}(a) = J_g(f(a))\, J_f(a).
  • Hypothesis is differentiability (existence of a good linear approximation), not merely existence of partial derivatives.
  • The scalar chain rule and all "multivariable chain rule" formulas from calculus courses are special cases of one matrix identity.
  • The proof is a direct ε\varepsilon–δ\delta estimate: it bounds the composition's remainder using differentiability of both ff and gg, plus continuity of ff.
  • Key idea: differentiation linearizes; composing differentiable maps composes their linearizations, and composing linear maps is matrix multiplication.

Where this lives today

The chain rule in matrix form is the computational backbone of backpropagation in machine learning: a neural network is literally a long composition of differentiable maps (affine layers composed with differentiable activation functions), and training requires the gradient of the loss with respect to every parameter, computed by repeatedly multiplying Jacobian matrices (or their transposes acting on vectors) layer by layer — reverse-mode automatic differentiation is precisely a computationally efficient application of this theorem. In physics and engineering, the chain rule underlies change of variables in partial differential equations: converting the heat or wave equation from Cartesian to polar, cylindrical, or curvilinear coordinates requires expressing derivatives with respect to new coordinates in terms of the old ones, exactly via the Jacobian product formula. In continuum mechanics, the deformation gradient of a material (itself a Jacobian matrix) composes according to this rule when successive deformations are applied. In pure mathematics, the theorem is the direct technical ingredient in the proof of the Inverse Function Theorem and the Implicit Function Theorem — both of which analyze compositions of a function with its (candidate) local inverse or with parametrizations of a level set, and both of which are indispensable in differential geometry, dynamical systems, and constrained optimization (Lagrange multipliers). In probability and statistics, the multivariate change-of-variables formula for densities under a smooth transformation uses the Jacobian determinant, and correctly deriving how a transformed random vector's density behaves under composition of transformations again invokes this chain rule.

Ask the tutor: Can you show me how the chain rule proof changes if ff maps into an infinite-dimensional Banach space instead of Rm\mathbb{R}^m?

DefinitionComposition, differentiability, and the derivative as a linear map

Let U⊆RnU \subseteq \mathbb{R}^n be open and let f:U→Rmf: U \to \mathbb{R}^m. We say ff is differentiable at a∈Ua \in U if there exists a linear map Df(a):Rn→RmDf(a): \mathbb{R}^n \to \mathbb{R}^m such that

lim⁡h→0∥f(a+h)−f(a)−Df(a)h∥∥h∥=0.\lim_{h \to 0} \frac{\lVert f(a+h) - f(a) - Df(a)h \rVert}{\lVert h \rVert} = 0.

The linear map Df(a)Df(a), when it exists, is unique, and is called the (total) derivative of ff at aa. With respect to the standard bases, Df(a)Df(a) is represented by the Jacobian matrix

Jf(a)=[∂fi∂xj(a)]1≤i≤m, 1≤j≤n.J_f(a) = \left[ \frac{\partial f_i}{\partial x_j}(a) \right]_{1 \le i \le m,\ 1 \le j \le n}.

Now let V⊆RmV \subseteq \mathbb{R}^m be open, let g:V→Rpg: V \to \mathbb{R}^p, and suppose f(U)⊆Vf(U) \subseteq V, so the composition

g∘f:U→Rp,(g∘f)(x)=g(f(x))g \circ f : U \to \mathbb{R}^p, \qquad (g\circ f)(x) = g(f(x))

is well-defined. If ff is differentiable at a∈Ua \in U and gg is differentiable at b:=f(a)∈Vb := f(a) \in V, the chain rule asserts that g∘fg \circ f is differentiable at aa with derivative

D(g∘f)(a)=Dg(b)∘Df(a),D(g \circ f)(a) = Dg(b) \circ Df(a),

the composite linear map Rn→Df(a)Rm→Dg(b)Rp\mathbb{R}^n \xrightarrow{Df(a)} \mathbb{R}^m \xrightarrow{Dg(b)} \mathbb{R}^p. In matrix form,

Jg∘f(a)=Jg(f(a)) Jf(a),J_{g\circ f}(a) = J_g(f(a)) \, J_f(a),

an ordinary p×mp \times m by m×nm \times n matrix product, giving a p×np \times n matrix.

Intuition

The one-variable chain rule says: to first order, zooming in on g(f(x))g(f(x)) near aa, you see two magnifications applied in sequence — stretch by f′(a)f'(a), then stretch by g′(f(a))g'(f(a)) — and sequential stretches multiply. In several variables the "stretch" at a point is not a single number but a linear map that can rotate, shear, and scale different directions differently. Composing two such maps means applying one linear transformation and then the other — and the matrix that represents "first apply AA, then apply BB" is exactly the product BABA (applied to a column vector vv: B(Av)B(Av)). So differentiability packages the intuitive idea "near aa, ff looks like an affine function f(a)+Df(a)hf(a) + Df(a)h" into a rigorous statement, and the chain rule packages "the best linear approximation to a composition is the composition of the best linear approximations." The content of the theorem is that these local linear pictures glue correctly under composition — that the error terms from approximating ff and approximating gg don't conspire to produce a non-negligible error in the composition. That is exactly what the proof below controls. A concrete visual: think of f:R2→R2f: \mathbb{R}^2 \to \mathbb{R}^2 as a coordinate change (say, from polar-like parameters to Cartesian) and g:R2→Rg: \mathbb{R}^2 \to \mathbb{R} as a temperature field. Near a point, ff's Jacobian tells you how a little rectangle of parameter-space gets stretched and sheared into a little parallelogram in physical space; gg's Jacobian (a gradient, as a 1×21\times 2 matrix) tells you the rate of temperature change per unit displacement in physical space. Multiplying the matrices tells you the rate of temperature change per unit displacement in parameter space — exactly the composite sensitivity.

Figure
Composition of two linear maps: derivative matrices multiply
Parameter space ℝ² →(Df(a))→ ℝ² →(Dg(b))→ ℝ (output) u = a + h unit square h → Df(a) h x = f(u) stretched parallelogram → Dg(b)(·) g(f(a)) temperature g(x) Composite derivative: D(g∘f)(a) = Dg(b)∘Df(a) In coordinates: J_{g∘f}(a) = J_g(b) · J_f(a) (matrix product)
TheoremThe Chain Rule for Differentiable Maps

Let U⊆RnU \subseteq \mathbb{R}^n and V⊆RmV \subseteq \mathbb{R}^m be open sets. Let

f:U→Rm,g:V→Rp,f: U \to \mathbb{R}^m, \qquad g: V \to \mathbb{R}^p,

with f(U)⊆Vf(U) \subseteq V, so g∘f:U→Rpg \circ f: U \to \mathbb{R}^p is defined. Let a∈Ua \in U and set b=f(a)∈Vb = f(a) \in V. Suppose: 1. ff is differentiable at aa, with derivative Df(a):Rn→RmDf(a): \mathbb{R}^n \to \mathbb{R}^m; 2. gg is differentiable at bb, with derivative Dg(b):Rm→RpDg(b): \mathbb{R}^m \to \mathbb{R}^p. Then g∘fg \circ f is differentiable at aa, and

D(g∘f)(a)=Dg(b)∘Df(a).D(g\circ f)(a) = Dg(b) \circ Df(a).

Equivalently, writing Jacobian matrices with respect to the standard bases,

Jg∘f(a)=Jg(b) Jf(a).J_{g \circ f}(a) = J_g(b)\, J_f(a).
Proof

Technique. Direct ε\varepsilon–δ\delta estimation from the definition of differentiability (the "little-o" definition), using continuity of ff at aa (itself a consequence of differentiability of ff at aa) to control the size of the increment fed into gg's remainder term.

Step 1: Set up remainder functions. Since ff is differentiable at aa, define for h∈Rnh \in \mathbb{R}^n with a+h∈Ua + h \in U:

r(h):=f(a+h)−f(a)−Df(a)h.r(h) := f(a+h) - f(a) - Df(a)h.

By hypothesis (1),

∥r(h)∥∥h∥→0as h→0,r(0)=0.(1)\frac{\lVert r(h) \rVert}{\lVert h \rVert} \to 0 \quad \text{as } h \to 0, \qquad r(0) = 0. \tag{1}

Since gg is differentiable at bb, define for k∈Rmk \in \mathbb{R}^m with b+k∈Vb+k \in V:

s(k):=g(b+k)−g(b)−Dg(b)k.s(k) := g(b+k) - g(b) - Dg(b)k.

By hypothesis (2),

∥s(k)∥∥k∥→0as k→0,s(0)=0.(2)\frac{\lVert s(k) \rVert}{\lVert k \rVert} \to 0 \quad \text{as } k \to 0, \qquad s(0)=0. \tag{2}

It is convenient to extend ss continuously at 00: define η(k)=∥s(k)∥/∥k∥\eta(k) = \lVert s(k)\rVert / \lVert k \rVert for k≠0k \neq 0 and η(0)=0\eta(0) = 0; by (2), η\eta is continuous at 00, and

∥s(k)∥=η(k) ∥k∥for all k (including k=0).(2’)\lVert s(k) \rVert = \eta(k)\,\lVert k \rVert \quad \text{for all } k \text{ (including } k=0\text{)}. \tag{2'}

Step 2: Continuity of ff at aa. We first record that differentiability of ff at aa implies ff is continuous at aa. Indeed, by the triangle inequality and the definition of Df(a)hDf(a)h as a linear (hence bounded — every linear map on Rn\mathbb{R}^n is bounded) map, there is a constant M≥0M \geq 0 with ∥Df(a)h∥≤M∥h∥\lVert Df(a)h \rVert \leq M \lVert h \rVert for all hh. Then

∥f(a+h)−f(a)∥≤∥Df(a)h∥+∥r(h)∥≤M∥h∥+∥r(h)∥→0\lVert f(a+h) - f(a) \rVert \leq \lVert Df(a) h \rVert + \lVert r(h) \rVert \leq M\lVert h \rVert + \lVert r(h)\rVert \to 0

as h→0h \to 0, using (1). Hence f(a+h)→f(a)f(a+h) \to f(a) as h→0h \to 0, i.e. ff is continuous at aa. This step is essential: it lets us pass h→0h \to 0 and conclude the corresponding increment k=k(h)k = k(h) in gg's domain also tends to 00, which is what allows us to invoke (2).

Step 3: Express the increment of g∘fg \circ f. For hh small enough that a+h∈Ua + h \in U (so that f(a+h)f(a+h) is defined and, by continuity established in Step 2 together with openness of VV and b∈Vb \in V, lies in VV for all sufficiently small hh — since VV open means some ball Bδ(b)⊆VB_\delta(b) \subseteq V, and continuity of ff at aa gives hh small enough that f(a+h)∈Bδ(b)f(a+h) \in B_\delta(b)), define

k(h):=f(a+h)−f(a)=Df(a)h+r(h).k(h) := f(a+h) - f(a) = Df(a) h + r(h).

Then f(a+h)=b+k(h)f(a+h) = b + k(h), and

(g∘f)(a+h)−(g∘f)(a)=g(b+k(h))−g(b)=Dg(b) k(h)+s(k(h)),(g\circ f)(a+h) - (g\circ f)(a) = g(b + k(h)) - g(b) = Dg(b)\,k(h) + s(k(h)),

using the definition of ss from Step 1. Substituting k(h)=Df(a)h+r(h)k(h) = Df(a)h + r(h) and using linearity of Dg(b)Dg(b):

(g∘f)(a+h)−(g∘f)(a)=Dg(b)(Df(a)h)+Dg(b)(r(h))+s(k(h)).(g\circ f)(a+h) - (g\circ f)(a) = Dg(b)\big(Df(a) h\big) + Dg(b)\big(r(h)\big) + s(k(h)).

Since Dg(b)∘Df(a)Dg(b) \circ Df(a) is linear (composition of linear maps is linear), we may write this as

(g∘f)(a+h)−(g∘f)(a)−(Dg(b)∘Df(a))h  =  Dg(b)(r(h))⏟(A)+s(k(h))⏟(B).(3)(g\circ f)(a+h) - (g\circ f)(a) - \big(Dg(b)\circ Df(a)\big) h \;=\; \underbrace{Dg(b)\big(r(h)\big)}_{(A)} + \underbrace{s(k(h))}_{(B)}. \tag{3}

To prove the theorem it suffices to show the right side of (3), divided by ∥h∥\lVert h \rVert, tends to 00 as h→0h \to 0.

Step 4: Estimate term (A). Since Dg(b)Dg(b) is linear on Rm\mathbb{R}^m, it is bounded: there is N≥0N \geq 0 with ∥Dg(b)v∥≤N∥v∥\lVert Dg(b)v \rVert \leq N \lVert v \rVert for all v∈Rmv \in \mathbb{R}^m. Hence

∥Dg(b)(r(h))∥∥h∥≤N∥r(h)∥∥h∥→h→00\frac{\lVert Dg(b)(r(h)) \rVert}{\lVert h \rVert} \leq N \frac{\lVert r(h) \rVert}{\lVert h \rVert} \xrightarrow{h \to 0} 0

by (1). So term (A)/∥h∥→0\lVert h\rVert \to 0.

Step 5: Estimate term (B). By (2'), ∥s(k(h))∥=η(k(h)) ∥k(h)∥\lVert s(k(h)) \rVert = \eta(k(h)) \, \lVert k(h) \rVert. Now

∥k(h)∥=∥Df(a)h+r(h)∥≤M∥h∥+∥r(h)∥.\lVert k(h) \rVert = \lVert Df(a) h + r(h) \rVert \leq M \lVert h \rVert + \lVert r(h) \rVert.

Therefore

∥s(k(h))∥∥h∥≤η(k(h))(M+∥r(h)∥∥h∥).\frac{\lVert s(k(h)) \rVert}{\lVert h \rVert} \leq \eta(k(h)) \left( M + \frac{\lVert r(h) \rVert}{\lVert h \rVert} \right).

As h→0h \to 0: by Step 2 (continuity of ff at aa, equivalently k(h)→0k(h) \to 0) and continuity of η\eta at 00 (Step 1), η(k(h))→η(0)=0\eta(k(h)) \to \eta(0) = 0. Also ∥r(h)∥/∥h∥→0\lVert r(h)\rVert/\lVert h \rVert \to 0 by (1), so the bracketed factor stays bounded (in fact tends to MM). A product of a factor tending to 00 and a bounded factor tends to 00. Hence term (B)/∥h∥→0\lVert h \rVert \to 0.

Step 6: Conclude. Combining Steps 4 and 5 with (3),

∥(g∘f)(a+h)−(g∘f)(a)−(Dg(b)∘Df(a))h∥∥h∥≤∥(A)∥∥h∥+∥(B)∥∥h∥→h→00.\frac{\big\lVert (g\circ f)(a+h) - (g\circ f)(a) - (Dg(b)\circ Df(a))h \big\rVert}{\lVert h \rVert} \leq \frac{\lVert (A)\rVert}{\lVert h\rVert} + \frac{\lVert (B) \rVert}{\lVert h \rVert} \xrightarrow{h\to 0} 0.

This is precisely the statement that g∘fg \circ f is differentiable at aa with derivative Dg(b)∘Df(a)Dg(b) \circ Df(a). By uniqueness of the derivative (a standard consequence of the definition — if two linear maps both satisfy the defining limit condition, their difference LL satisfies ∥Lh∥/∥h∥→0\lVert Lh\rVert/\lVert h\rVert \to 0, and testing along h=tvh = tv for fixed vv and t→0t \to 0 forces Lv=0Lv = 0 for every vv, so L=0L = 0), this is the derivative D(g∘f)(a)D(g\circ f)(a).

Step 7: Matrix form. The matrix of a composition of linear maps, with respect to standard bases, is the matrix product of the matrices of the two maps (this is a standard fact from linear algebra: if S:Rn→RmS: \mathbb{R}^n \to \mathbb{R}^m has matrix AA and T:Rm→RpT:\mathbb{R}^m \to \mathbb{R}^p has matrix BB, then T∘ST\circ S has matrix BABA, since (T∘S)(ej)=T(Aej)=T(∑iAijei)=∑iAijT(ei)=∑iAij∑kBkiek=∑k(∑iBkiAij)ek(T\circ S)(e_j) = T(Ae_j) = T\left(\sum_i A_{ij} e_i\right) = \sum_i A_{ij} T(e_i) = \sum_i A_{ij} \sum_k B_{ki} e_k = \sum_k \left(\sum_i B_{ki}A_{ij}\right)e_k, and the coefficient ∑iBkiAij\sum_i B_{ki}A_{ij} is exactly the (k,j)(k,j) entry of BABA). Applying this with S=Df(a)S = Df(a) (matrix Jf(a)J_f(a)) and T=Dg(b)T = Dg(b) (matrix Jg(b)J_g(b)) gives

Jg∘f(a)=Jg(b) Jf(a)=Jg(f(a)) Jf(a).J_{g\circ f}(a) = J_g(b)\, J_f(a) = J_g(f(a))\,J_f(a).

■\blacksquare

∎
CorollaryEntrywise / partial-derivative form

Under the hypotheses of the theorem, write y=f(x)∈Rmy = f(x) \in \mathbb{R}^m with components y1,…,ymy_1,\dots,y_m, and z=g(y)∈Rpz = g(y) \in \mathbb{R}^p with components z1,…,zpz_1, \dots, z_p. Then for each 1≤k≤p1 \le k \le p and 1≤j≤n1 \le j \le n,

∂zk∂xj(a)=∑i=1m∂zk∂yi(b) ∂yi∂xj(a).\frac{\partial z_k}{\partial x_j}(a) = \sum_{i=1}^m \frac{\partial z_k}{\partial y_i}(b) \, \frac{\partial y_i}{\partial x_j}(a).

This is simply the (k,j)(k,j) entry of the matrix identity Jg∘f=JgJfJ_{g\circ f} = J_g J_f, expanded via the row-times-column rule for matrix multiplication.

Worked examples

Example 1 (asymmetric coordinate change feeding a scalar field). Let f:R2→R3f: \mathbb{R}^2 \to \mathbb{R}^3 be given by

f(u,v)=(u2v,  u+3v2,  euv),f(u,v) = \big(u^2 v,\ \ u + 3v^2,\ \ e^{uv}\big),

and let g:R3→Rg: \mathbb{R}^3 \to \mathbb{R} be given by

g(x,y,z)=xy2−ln⁡(z+2),z>−2.g(x,y,z) = x y^2 - \ln(z + 2), \qquad z > -2.

We compute D(g∘f)D(g\circ f) at the point a=(u,v)=(1,2)a = (u,v) = (1,2) using the chain rule, and check it against direct substitution. Step 1: evaluate ff at aa. With u=1,v=2u=1, v=2:

x=u2v=1⋅2=2,y=u+3v2=1+12=13,z=euv=e2.x = u^2v = 1\cdot 2 = 2, \qquad y = u+3v^2 = 1 + 12 = 13, \qquad z = e^{uv} = e^{2}.

So b=f(1,2)=(2,13,e2)b = f(1,2) = (2, 13, e^2). Step 2: Jacobian of ff at (1,2)(1,2).

Jf(u,v)=(2uvu216vveuvueuv)⟹Jf(1,2)=(411122e2e2).J_f(u,v) = \begin{pmatrix} 2uv & u^2 \\ 1 & 6v \\ v e^{uv} & u e^{uv} \end{pmatrix} \quad\Longrightarrow\quad J_f(1,2) = \begin{pmatrix} 4 & 1 \\ 1 & 12 \\ 2e^2 & e^2 \end{pmatrix}.

(Row 1: ∂(u2v)/∂u=2uv=4\partial(u^2v)/\partial u = 2uv = 4, ∂(u2v)/∂v=u2=1\partial(u^2v)/\partial v = u^2 = 1. Row 2: ∂(u+3v2)/∂u=1\partial(u+3v^2)/\partial u = 1, /∂v=6v=12/\partial v = 6v = 12. Row 3: ∂(euv)/∂u=veuv=2e2\partial(e^{uv})/\partial u = v e^{uv} = 2e^2, /∂v=ueuv=e2/\partial v = u e^{uv} = e^2.) Step 3: Jacobian of gg at b=(2,13,e2)b = (2,13,e^2).

Jg(x,y,z)=(y22xy−1z+2)⟹Jg(2,13,e2)=(16952−1e2+2).J_g(x,y,z) = \begin{pmatrix} y^2 & 2xy & -\dfrac{1}{z+2} \end{pmatrix} \quad\Longrightarrow\quad J_g(2,13,e^2) = \begin{pmatrix} 169 & 52 & -\dfrac{1}{e^2+2} \end{pmatrix}.

(Since ∂g/∂x=y2=169\partial g/\partial x = y^2 = 169; ∂g/∂y=2xy=2⋅2⋅13=52\partial g/\partial y = 2xy = 2\cdot2\cdot13 = 52; ∂g/∂z=−1/(z+2)=−1/(e2+2)\partial g/\partial z = -1/(z+2) = -1/(e^2+2).) Step 4: multiply. By the theorem,

Jg∘f(1,2)=Jg(b) Jf(1,2)=(16952−1e2+2)(411122e2e2).J_{g\circ f}(1,2) = J_g(b)\,J_f(1,2) = \begin{pmatrix} 169 & 52 & -\dfrac{1}{e^2+2} \end{pmatrix} \begin{pmatrix} 4 & 1 \\ 1 & 12 \\ 2e^2 & e^2 \end{pmatrix}.

Entry 1 (the ∂/∂u\partial/\partial u component):

169(4)+52(1)+(−1e2+2)(2e2)=676+52−2e2e2+2=728−2e2e2+2.169(4) + 52(1) + \left(-\frac{1}{e^2+2}\right)(2e^2) = 676 + 52 - \frac{2e^2}{e^2+2} = 728 - \frac{2e^2}{e^2+2}.

Entry 2 (the ∂/∂v\partial/\partial v component):

169(1)+52(12)+(−1e2+2)(e2)=169+624−e2e2+2=793−e2e2+2.169(1) + 52(12) + \left(-\frac{1}{e^2+2}\right)(e^2) = 169 + 624 - \frac{e^2}{e^2+2} = 793 - \frac{e^2}{e^2+2}.

So

Jg∘f(1,2)=(728−2e2e2+2,  793−e2e2+2).J_{g\circ f}(1,2) = \left( 728 - \frac{2e^2}{e^2+2},\ \ 793 - \frac{e^2}{e^2+2} \right).

Check by direct substitution. (g∘f)(u,v)=(u2v)(u+3v2)2−ln⁡(euv+2)(g\circ f)(u,v) = (u^2v)(u+3v^2)^2 - \ln(e^{uv}+2). Differentiating directly with respect to uu:

∂u[(u2v)(u+3v2)2]=2uv(u+3v2)2+u2v⋅2(u+3v2),∂u[−ln⁡(euv+2)]=−veuveuv+2.\partial_u\big[(u^2v)(u+3v^2)^2\big] = 2uv(u+3v^2)^2 + u^2v\cdot 2(u+3v^2), \qquad \partial_u\big[-\ln(e^{uv}+2)\big] = -\frac{v e^{uv}}{e^{uv}+2}.

At (1,2)(1,2): (u+3v2)=13(u+3v^2) = 13, so first term is 2(1)(2)(169)+(1)(2)(2)(13)=676+52=7282(1)(2)(169) + (1)(2)(2)(13) = 676 + 52 = 728, matching. Second term: −2e2/(e2+2)-2e^2/(e^2+2), matching exactly. The ∂/∂v\partial/\partial v computation checks similarly. The chain-rule matrix product reproduces the direct computation, as it must.

Example 2 (a map to itself, second-order check via composition of two nontrivial maps). Let f:R2→R2f: \mathbb{R}^2 \to \mathbb{R}^2, f(x,y)=(x2−y, xy)f(x,y) = (x^2 - y,\ xy), and let g:R2→R2g: \mathbb{R}^2 \to \mathbb{R}^2, g(s,t)=(s+t2, st)g(s,t) = (s+t^2,\ s t). We compute the Jacobian of g∘fg \circ f at a=(2,1)a = (2,1). Step 1: f(2,1)f(2,1). x2−y=4−1=3x^2 - y = 4-1=3, xy=2xy = 2. So b=f(2,1)=(3,2)b = f(2,1) = (3,2). Step 2: JfJ_f.

Jf(x,y)=(2x−1yx),Jf(2,1)=(4−112).J_f(x,y) = \begin{pmatrix} 2x & -1 \\ y & x \end{pmatrix}, \qquad J_f(2,1) = \begin{pmatrix} 4 & -1 \\ 1 & 2 \end{pmatrix}.

Step 3: JgJ_g.

Jg(s,t)=(12tts),Jg(3,2)=(1423).J_g(s,t) = \begin{pmatrix} 1 & 2t \\ t & s \end{pmatrix}, \qquad J_g(3,2) = \begin{pmatrix} 1 & 4 \\ 2 & 3 \end{pmatrix}.

Step 4: product.

Jg∘f(2,1)=Jg(3,2) Jf(2,1)=(1423)(4−112)=(1⋅4+4⋅11⋅(−1)+4⋅22⋅4+3⋅12⋅(−1)+3⋅2)=(87114).J_{g\circ f}(2,1) = J_g(3,2)\,J_f(2,1) = \begin{pmatrix} 1 & 4 \\ 2 & 3 \end{pmatrix}\begin{pmatrix} 4 & -1 \\ 1 & 2 \end{pmatrix} = \begin{pmatrix} 1\cdot4+4\cdot1 & 1\cdot(-1)+4\cdot2 \\ 2\cdot4+3\cdot1 & 2\cdot(-1)+3\cdot2 \end{pmatrix} = \begin{pmatrix} 8 & 7 \\ 11 & 4 \end{pmatrix}.

Check by direct substitution. (g∘f)(x,y)=((x2−y)+(xy)2, (x2−y)(xy))(g\circ f)(x,y) = \big((x^2-y) + (xy)^2,\ (x^2-y)(xy)\big). Let P=x2−yP = x^2-y, Q=xyQ = xy; the output is (P+Q2,PQ)(P+Q^2, PQ). First component u1=P+Q2u_1 = P + Q^2: ∂u1/∂x=2x+2QQx=2x+2xy⋅y=2x+2xy2\partial u_1/\partial x = 2x + 2QQ_x = 2x + 2xy\cdot y = 2x+2xy^2; at (2,1)(2,1): 4+2⋅2⋅1=4+4=84 + 2\cdot2\cdot1 = 4+4=8. Matches entry (1,1)=8(1,1)=8. ∂u1/∂y=−1+2QQy=−1+2xy⋅x=−1+2x2y\partial u_1/\partial y = -1 + 2Q Q_y = -1+2xy\cdot x = -1+2x^2y; at (2,1)(2,1): −1+2(4)(1)=−1+8=7-1+2(4)(1) = -1+8=7. Matches entry (1,2)=7(1,2)=7. Second component u2=PQu_2 = PQ: ∂u2/∂x=PxQ+PQx=2x⋅xy+(x2−y)y=2x2y+x2y−y2=3x2y−y2\partial u_2/\partial x = P_x Q + P Q_x = 2x\cdot xy + (x^2-y)y = 2x^2y + x^2y - y^2 = 3x^2y - y^2; at (2,1)(2,1): 3(4)(1)−1=12−1=113(4)(1) - 1 = 12-1=11. Matches entry (2,1)=11(2,1)=11. ∂u2/∂y=PyQ+PQy=(−1)(xy)+(x2−y)(x)=−xy+x3−xy=x3−2xy\partial u_2/\partial y = P_y Q + P Q_y = (-1)(xy) + (x^2-y)(x) = -xy + x^3 - xy = x^3 - 2xy; at (2,1)(2,1): 8−4=48 - 4 = 4. Matches entry (2,2)=4(2,2)=4. Both entries agree exactly with the matrix-product computation, confirming Jg∘f(2,1)=(87114)J_{g\circ f}(2,1) = \begin{pmatrix} 8 & 7\\ 11 & 4\end{pmatrix}.

CounterexampleWhy "differentiable," not merely "partial derivatives exist," is needed

The theorem requires ff differentiable at aa and gg differentiable at b=f(a)b = f(a) — the full total-derivative condition, not just existence of the partial derivatives (which is a strictly weaker condition in several variables). Dropping this to "partial derivatives exist" breaks the conclusion. Let g:R2→Rg: \mathbb{R}^2 \to \mathbb{R} be the classical example

g(x,y)={xyx2+y2(x,y)≠(0,0)0(x,y)=(0,0).g(x,y) = \begin{cases} \dfrac{xy}{x^2+y^2} & (x,y) \ne (0,0) \\ 0 & (x,y) = (0,0). \end{cases}

Both partial derivatives of gg exist at the origin: along y=0y=0, g≡0g\equiv 0, so gx(0,0)=0g_x(0,0)=0; similarly gy(0,0)=0g_y(0,0) = 0. But gg is not differentiable at (0,0)(0,0) — indeed gg is not even continuous there, since along the line y=xy=x, g(x,x)=x2/(2x2)=1/2g(x,x) = x^2/(2x^2) = 1/2 for all x≠0x\ne 0, while g(0,0)=0g(0,0)=0. Now let f:R→R2f: \mathbb{R} \to \mathbb{R}^2, f(t)=(t,t)f(t) = (t, t), which is differentiable everywhere (it's linear), with f(0)=(0,0)=:a′f(0) = (0,0) =: a' image point, i.e. take the composite variable a=0a=0, b=f(0)=(0,0)b = f(0) = (0,0). Consider h(t):=g(f(t))=g(t,t)=t⋅tt2+t2=12h(t) := g(f(t)) = g(t,t) = \dfrac{t\cdot t}{t^2+t^2} = \dfrac12 for t≠0t \ne 0, and h(0)=g(0,0)=0h(0) = g(0,0) = 0. So hh is not even continuous at t=0t=0 (it jumps from the constant value 1/21/2 to 00), hence certainly not differentiable there. If one attempted to apply a "chain rule" using only gx(0,0)=0g_x(0,0)=0 and gy(0,0)=0g_y(0,0)=0 formally, the naive formula would predict

h′(0)gx(0,0)f1′(0)+gy(0,0)f2′(0)=0⋅1+0⋅1=0,h'(0) \overset{?}{=} g_x(0,0) f_1'(0) + g_y(0,0) f_2'(0) = 0\cdot 1 + 0 \cdot 1 = 0,

suggesting hh is even differentiable with derivative 00 — but hh is not even continuous at 00, a contradiction. The failure is entirely attributable to gg failing to be differentiable (in the total-derivative sense) at (0,0)(0,0), even though both its partial derivatives exist there. This shows the hypothesis "differentiable at the point," not merely "partials exist," is exactly what the theorem needs and cannot be weakened.

Scope, edge cases & common traps

SituationWhat to checkCommon pitfall
f,gf,g both C1C^1 (continuously differentiable) on open setsC1⇒C^1 \Rightarrow differentiable at every point, chain rule applies freelyNone — this is the safe, standard case used in most calculus courses
Partial derivatives of ff or gg exist but are not known to be continuousMust verify differentiability directly (via the limit definition), not just existence of partialsAssuming existence of partials suffices; see the counterexample above
Domains don't compose correctly (f(U)⊈Vf(U) \not\subseteq V)g∘fg\circ f may not even be defined near aaSilently plugging in formulas without checking the image lies in gg's domain
Dimensions mismatched in matrix productJgJ_g is p×mp\times m, JfJ_f is m×nm\times n; inner dimensions (mm) must matchWriting JfJgJ_f J_g instead of JgJfJ_g J_f — matrix multiplication is not commutative, and only one order is even conformable in general
Chain rule along a curve (implicit differentiation setup)f:R→Rmf:\mathbb{R}\to\mathbb{R}^m, g:Rm→Rg:\mathbb{R}^m\to\mathbb{R}: result is ddtg(f(t))=∇g(f(t))⋅f′(t)\frac{d}{dt}g(f(t)) = \nabla g(f(t))\cdot f'(t), a dot product, which is the 1×11\times1 case of the matrix productForgetting this is literally the same theorem, not a separate rule to memorize

The mathematician behind it

Origins · 1 minCarl Gustav Jacob Jacobi and how this was discovered
CG
Carl Gustav Jacob Jacobi
1804–1851 · Prussia, University of Königsberg

Jacobi was a prodigiously productive 19th-century mathematician working in analysis, algebra, and mechanics. He systematized the study of the determinant of partial derivatives that now bears his name — the Jacobian — in his work on transformations of multiple integrals and systems of differential equations, and his functional-determinant calculus underlies the modern matrix formulation of the multivariable chain rule, though the fully rigorous ε\varepsilon–δ\delta formulation of differentiability used in this article is a later, Weierstrass-era refinement.

History
- 1841 — Jacobi publishes his major memoir on functional determinants ("De determinantibus functionalibus"), formalizing the object now called the Jacobian and its role in transforming multiple integrals.
- 1870s–1880s — Karl Weierstrass and collaborators rigorize the ε\varepsilon–δ\delta foundations of analysis, giving the precise "linear approximation" definition of differentiability used to state and prove the chain rule cleanly.
- 1900s–1930s — Maurice Fréchet and others generalize differentiability to normed vector spaces and Banach spaces, showing the chain rule as proved here is the finite-dimensional case of a much more general theorem (the Fréchet derivative chain rule).
- 1950s–1960s — The Implicit and Inverse Function Theorems, whose proofs rest crucially on the chain rule and matrix Jacobians, become standard core material in analysis courses following texts such as Rudin's Principles of Mathematical Analysis.
- 1980s–present — The chain rule for compositions of differentiable maps becomes the mathematical basis of automatic differentiation and the backpropagation algorithm central to training neural networks.
Why is the multivariable chain rule a matrix product and not just a sum of products?

It is a sum of products entrywise — the corollary's formula ∂zk/∂xj=∑i(∂zk/∂yi)(∂yi/∂xj)\partial z_k/\partial x_j = \sum_i (\partial z_k/\partial y_i)(\partial y_i/\partial x_j) is exactly what you get from expanding matrix multiplication entry by entry. The point of stating it as a matrix product JgJfJ_g J_f rather than as a list of summation formulas is conceptual and practical: it makes clear that the derivative of a composition is the composition of linear maps (this is what "chain" really means — chaining two linear approximations together), it generalizes uniformly to any input/output dimensions without re-deriving index formulas each time, and it connects directly to the definition of the total derivative as a linear map rather than a table of numbers.

Does the chain rule require ff and gg to be continuously differentiable (C1C^1), or just differentiable at a point?

Only differentiability at the relevant single points is required: ff differentiable at aa, and gg differentiable at b=f(a)b=f(a). The proof given uses only the defining limit of differentiability at those two points (plus the automatic continuity of ff at aa that differentiability there implies). C1C^1 is a stronger, global condition (differentiable at every point of an open set with continuously varying derivative) that is often assumed in practice because it's easy to verify from formulas, but the chain rule itself is a pointwise statement and holds under the weaker pointwise hypothesis.

How does this relate to the one-variable chain rule I already know?

It specializes to it exactly. If n=m=p=1n=m=p=1, then Df(a)Df(a) and Dg(b)Dg(b) are 1×11\times1 matrices — literally the numbers f′(a)f'(a) and g′(b)g'(b) — and matrix multiplication of 1×11\times1 matrices is ordinary multiplication of real numbers. So Jg∘f(a)=Jg(b)Jf(a)J_{g\circ f}(a) = J_g(b)J_f(a) reads (g∘f)′(a)=g′(f(a))f′(a)(g\circ f)'(a) = g'(f(a))f'(a), the familiar formula. Every version of the chain rule taught in a calculus sequence (single variable, multivariable "tree diagram" versions, directional derivative versions) is a special case of this one matrix statement.

What exactly goes wrong if I only assume the partial derivatives of ff and gg exist, rather than full differentiability?

Existence of partial derivatives only guarantees good linear approximation along the coordinate axis directions; it says nothing about behavior along other directions, and in particular does not even guarantee continuity of the function at the point, as the counterexample g(x,y)=xy/(x2+y2)g(x,y)=xy/(x^2+y^2) shows. Since the chain rule proof crucially uses continuity of the inner function (to know the increment fed into the outer function's remainder term also shrinks) and the full defining limit of differentiability for both maps, weakening either hypothesis to mere existence of partials removes exactly the ingredients the proof needs, and the conclusion can fail outright, including failure of the composite to be continuous, let alone differentiable.

Is the chain rule the same thing as implicit differentiation?

Implicit differentiation is a direct application of the chain rule, not a separate rule. If F(x,y)=0F(x,y)=0 implicitly defines yy as a differentiable function of xx, differentiating the composition x↦F(x,y(x))x \mapsto F(x, y(x)) using the chain rule gives Fx+Fy y′(x)=0F_x + F_y \, y'(x) = 0, hence y′(x)=−Fx/Fyy'(x) = -F_x/F_y where Fy≠0F_y \ne 0. This is exactly the scalar chain rule for the composite map x↦(x,y(x))↦F(x,y(x))x \mapsto (x, y(x)) \mapsto F(x,y(x)), and its rigorous justification (that y(x)y(x) actually is differentiable) is precisely the content of the Implicit Function Theorem, whose proof in turn uses this chain rule.

Does the chain rule generalize beyond Rn\mathbb{R}^n, e.g. to infinite dimensions or manifolds?

Yes. In normed vector spaces (Banach spaces), the same statement and essentially the same proof (with "linear map" meaning "bounded linear operator" and norms replacing Euclidean norms) yield the chain rule for the Fréchet derivative, which is the natural generalization used in the calculus of variations and nonlinear functional analysis. On smooth manifolds, the chain rule for the differential of a smooth map between manifolds is stated and proved using local coordinate charts, reducing again to compositions of Euclidean-space maps and this exact theorem — it's the local model in every more general setting.

How does this connect to the Jacobian determinant and change of variables in multiple integrals?

When n=m=pn=m=p (same input and output dimension at every stage), the Jacobian matrices are square, and one can take determinants. Since det⁡(AB)=det⁡(A)det⁡(B)\det(AB) = \det(A)\det(B) for square matrices, the chain rule immediately gives det⁡Jg∘f(a)=det⁡Jg(b)⋅det⁡Jf(a)\det J_{g\circ f}(a) = \det J_g(b) \cdot \det J_f(a) — the determinants of composed coordinate changes multiply. This is exactly the fact used to verify that the change-of-variables formula for multiple integrals composes correctly when you chain two coordinate transformations (e.g., Cartesian → cylindrical → a further curvilinear system).

Before this

Where this leads