Suppose you know a function's value at a single point. You also know its rate of change in every direction. You know how those rates of change are themselves changing. How well can you predict the function's value nearby? In one variable, Taylor's theorem answers this. It builds a polynomial from derivatives at a point, plus a controlled remainder. The multivariable version answers the same question in Rn. Here, "derivative" now means a gradient vector and a Hessian matrix. The bookkeeping of higher partial derivatives threatens to become unmanageable without the right notation. The idea matters because almost every local analysis in higher-dimensional calculus reduces to one thing. This includes classifying critical points, justifying Newton's method, deriving error bounds in numerical analysis, and expanding energies in physics near equilibrium. All of these reduce to approximating a smooth function by a quadratic model near a point. The quadratic model is exactly the second-order Taylor polynomial. The Hessian matrix that appears in it governs whether a critical point is a max, a min, or a saddle. None of this can be stated cleanly without first taming the notation for higher-order partial derivatives. The big picture: multi-index notation converts sums like ∑i,j∂xi∂xj∂2fhihj into single symbols Dαf. This exactly mirrors the one-variable formula ∑kk!f(k)(a)(x−a)k. Once this notation is in place, the multivariable Taylor theorem is proved by reducing it to the single-variable case. Restrict f to the line segment joining the two points. Apply ordinary one-variable Taylor theorem to the resulting function of one real variable. Then translate back using the chain rule. This article builds multi-index notation from scratch. It states and proves the second-order Taylor theorem with Lagrange remainder in Rn. It works through two genuinely asymmetric numerical examples. It also shows exactly why each hypothesis is needed — open domain, continuity of second partials, convexity of the segment — by exhibiting what fails without it.
Key takeaways
A multi-indexα=(α1,…,αn) packages partial derivative orders into one symbol, with ∣α∣=∑αi and α!=∏αi!.
The second-order Taylor polynomial of f at a is f(a)+∇f(a)⋅h+21hTHf(a)h, where Hf(a) is the Hessian matrix.
The exact Lagrange remainder involves third partial derivatives evaluated at some point on the segment from a to a+h — its existence, not its location, is what we prove.
The proof works by restricting f to a line through a, reducing an n-variable problem to one-variable Taylor's theorem with remainder.
Clairaut's Theorem on Mixed Partials is essential: it lets the Hessian be symmetric and lets us collapse cross terms like hihj∂i∂jf without ambiguity of order.
This expansion is the algebraic engine behind the second-derivative test for local extrema, foreshadowing The Hessian Matrix and Critical Points and Extrema.
Key idea: every multivariable Taylor expansion is secretly a single-variable Taylor expansion of g(t)=f(a+th), unpacked by the chain rule.
Where this lives today
The second-order Taylor expansion is the algebraic core of the second-derivative test for classifying critical points in Critical Points and Extrema: at a critical point a where ∇f(a)=0, the expansion collapses to f(a+h)=f(a)+21hTHf(a)h+o(∥h∥2), so the local behavior of f near a is governed entirely by the definiteness of The Hessian Matrix — positive definite gives a local minimum, negative definite a local maximum, indefinite a saddle. In numerical optimization, Newton's method for minimizing a function repeatedly builds the quadratic model P2 at the current iterate and jumps to its minimizer, exploiting exactly this expansion; the remainder bound O(∥h∥3) explains the method's famous quadratic local convergence rate. In numerical analysis more broadly, finite-difference schemes for partial differential equations are derived and error-bounded using multivariable Taylor expansions with explicit remainder terms. In physics, expanding a potential energy function to second order near an equilibrium point produces the harmonic oscillator approximation, with the Hessian playing the role of a stiffness (or spring-constant) matrix whose eigenvalues give normal-mode frequencies. In statistics and machine learning, second-order (quadratic) approximations of loss functions underlie the Laplace approximation for Bayesian posteriors and quasi-Newton optimization algorithms such as BFGS.
:::faq[Can the theorem fail if U is open but not convex, even if the specific segment [a,a+h] stays inside U?] No — convexity of U is never actually used in the proof. The only requirement is that the particular segment from a to a+h lies in U; the rest of U's shape is irrelevant. So on a non-convex open set, Taylor's theorem still holds for pairs (a,h) whose connecting segment happens to avoid the "missing" region, but may fail (or simply not apply) for pairs whose segment exits the domain, exactly as illustrated in the counterexample with the slit plane.
Ask the tutor: How does the Hessian's definiteness at a critical point actually determine whether it's a max, min, or saddle?
DefinitionMulti-indices, and Taylor polynomials in Rn
A multi-index is an n-tuple α=(α1,…,αn) of nonnegative integers. For a multi-index α and a vector h=(h1,…,hn)∈Rn, define
If f is a function on an open set U⊆Rn possessing all partial derivatives of order ∣α∣ at a point, write
Dαf=∂x1α1∂x2α2⋯∂xnαn∂∣α∣f.
Let U⊆Rn be open, let a∈U, and let f:U→R be a function whose partial derivatives up to order k+1 all exist and are continuous on U (i.e. f∈Ck+1(U)). For h∈Rn such that the segment [a,a+h]⊂U, the k-th order Taylor polynomial of f at a is
where ∇f(a)=(∂1f(a),…,∂nf(a)) is the gradient and Hf(a)=[∂i∂jf(a)]i,j=1n is the Hessian matrix, a symmetric n×n matrix.
Intuition
In one variable, Taylor's theorem says a smooth function looks, near a, like a polynomial in (x−a): a constant term (the value), a linear term (the slope), a quadratic term (the curvature), and so on. In several variables the same idea holds, but "slope" becomes a direction-dependent number — the directional derivative ∇f(a)⋅h — and "curvature" becomes a quadratic formhTHf(a)h that can be positive in some directions and negative in others. The Hessian matrix is precisely the object that records this direction-dependent bending. The cleanest way to see why the formula must involve ∇f(a)⋅h and hTHf(a)h is to freeze the direction: fix h, and look at the single real function g(t)=f(a+th) for t near 0. This is an ordinary one-variable function, so its own Taylor expansion at t=0 is
g(t)=g(0)+g′(0)t+21g′′(0)t2+⋯,
and the chain rule converts g′(0) into ∇f(a)⋅h and g′′(0) into hTHf(a)h. Setting t=1 recovers the multivariable expansion. This restriction-to-a-line trick is not just an analogy; it is the actual proof technique used below.
Figure
The quadratic approximation of f(x,y)=excosy near a=(0,0): level curves of f (solid) versus level curves of its Taylor quadratic P2 (dashed) agree to second order
Central theorem
TheoremSecond-Order Taylor's Theorem with Lagrange Remainder, in Rn
Let U⊆Rn be open and let f:U→R be of class C3 on U (all partial derivatives up to third order exist and are continuous). Let a∈U and let h∈Rn be such that the closed segment
[a,a+h]={a+th:t∈[0,1]}
is entirely contained in U. Then there exists θ∈(0,1) such that
Proof of the Lemma, by repeated application of The Chain Rule in Several Variables. Since [a,a+h]⊂U and U is open, and t↦a+th is a smooth (affine, hence C∞) map from [0,1] into U, the composite g(t)=f(a+th) is well defined on [0,1]; because f∈C3(U), the chain rule guarantees g∈C3([0,1]) (composition of C3 with an affine — in particular C∞ — map is C3). By the chain rule, writing γ(t)=a+th so γ′(t)=h for all t,
This exhibits g′ as a new function of t, namely t↦∑ihi∂if(a+th), a finite linear combination (with constant coefficients hi) of functions of the form t↦∂if(a+th). Differentiating again, applying the chain rule to each term ∂if(a+th) exactly as before (since ∂if∈C2(U), being a component of f∈C3),
Note that this double sum is exactly the quadratic form associated to the matrix Hf(a+th)=[∂i∂jf(a+th)]; Clairaut's Theorem on Mixed Partials guarantees ∂i∂jf=∂j∂if throughout U (since f∈C2(U) at least), so Hf is symmetric — though symmetry is not even needed to evaluate the sum, only to justify calling it "the" Hessian without ambiguity of index order. Differentiating once more, term by term, using ∂i∂jf∈C1(U) (a consequence of f∈C3):
where continuity of third partials (again from f∈C3) lets us reorder the indices in each triple product freely by repeated use of Clairaut's theorem. This proves the Lemma. ■
Proof of the Theorem, using one-variable Taylor's theorem with Lagrange remainder (assumed known from single-variable calculus) applied to g. By the Lemma, g∈C3([0,1]), so ordinary one-variable Taylor's theorem with Lagrange remainder applies to g at the base point t=0, evaluated at t=1: there exists θ∈(0,1) such that
g(1)=g(0)+g′(0)⋅1+2g′′(0)⋅12+6g′′′(θ)⋅13.
(This is the single-variable Taylor theorem with n=2, i.e., quadratic Taylor polynomial plus Lagrange remainder for the third derivative, applied to the interval [0,1]; its proof — via the Cauchy Mean Value Theorem or Rolle's theorem applied to an auxiliary function — is standard one-variable material and is not reproduced here, only invoked.) Now translate each term back to f using the Lemma:
g(1)=f(a+h) and g(0)=f(a), directly from the definition g(t)=f(a+th).
g′(0)=∇f(a)⋅h, by the Lemma's first formula at t=0.
g′′(0)=hTHf(a)h, by the Lemma's second formula at t=0.
g′′′(θ)=∑i,j,k∂i∂j∂kf(a+θh)hihjhk, by the Lemma's third formula at t=θ.
Substituting these four identities into the displayed one-variable Taylor expansion of g gives exactly
which is the claimed formula for R2(h) (the coefficient α!3! in the multi-index form arises from collecting repeated indices i,j,k into a single multi-index α with ∣α∣=3, a standard multinomial bookkeeping identity: the number of ordered triples (i,j,k) collapsing to a given α is α!3!). For the bound: since each of the n3 terms ∂i∂j∂kf(a+θh)hihjhk satisfies
∣∂i∂j∂kf(a+θh)hihjhk∣≤M3∥h∥∞3
(using ∣hi∣,∣hj∣,∣hk∣≤∥h∥∞ and the definition of M3 as a supremum over the segment, which contains a+θh since θ∈(0,1)), the triangle inequality applied to the triple sum over i,j,k∈{1,…,n} gives
Since f∈C3(U) and [a,a+h] is a fixed compact segment, M3<∞ automatically (continuous functions on compact sets are bounded), so this bound is always available for h small enough that [a,a+h]⊂U. Finally, since ∥h∥∞≤∥h∥ (Euclidean norm) up to the fixed constant n, we get ∣R2(h)∣=O(∥h∥3)=o(∥h∥2) as h→0, because ∥h∥3/∥h∥2=∥h∥→0. ■
∎
CorollarySecond-order approximation
Under the hypotheses of the theorem, as h→0,
f(a+h)=f(a)+∇f(a)⋅h+21hTHf(a)h+o(∥h∥2).
Worked examples
Example 1: A non-polynomial function at an asymmetric point. Let f(x,y)=excosy (globally C∞ on all of R2, in particular C3), and expand about the deliberately asymmetric base point a=(1,π/3), i.e. compute P2(h) for h=(h1,h2) near 0. First compute the needed partials:
Numerical check. Take h=(0.1,−0.2), i.e. approximate f(1.1,π/3−0.2). Exact value: x=1.1, y=π/3−0.2≈1.04720−0.2=0.84720. Compute e1.1≈3.00417 and cos(0.84720)≈0.66235. So f(a+h)≈3.00417×0.66235≈1.99039. Taylor estimate: plug h1=0.1, h2=−0.2 into P2(h) with e/2≈1.35914:
So P2(h)≈1.35914×1.46605≈1.99277. The two agree to about 0.0024 absolute error, and the theorem guarantees this discrepancy is O(∥h∥3) — with ∥h∥∞=0.2, indeed 0.23=0.008 is the right order of magnitude, comfortably bounding the observed error once the constant 6n3M3 is included.
Example 2: Verifying the remainder bound explicitly for a cubic. Let f(x,y)=x3+2x2y−xy2+y3, expanded at a=(2,−1) (again asymmetric, with mixed signs), evaluated at h=(1,1) (so target point (3,0)). Since f is a cubic polynomial, all fourth-order partials vanish, so the Lagrange remainder R2(h) should reproduce the exact cubic terms with no error left over (a good sanity check on the formula, since the theorem's remainder must exactly capture whatever the theorem's second-order part misses). Compute all needed derivatives:
fx=3x2+4xy−y2,fy=2x2−2xy+3y2,
fxx=6x+4y,fyy=−2x+6y,fxy=4x−2y,
fxxx=6,fxxy=4,fxyy=−2,fyyy=6.
All third partials are constant, so evaluating them at a+θh for any θ gives the same values — meaning the remainder is exact for any choice of θ∈(0,1), consistent with the theorem (which only asserts existence of someθ, but here every θ works because the third derivatives don't vary). At a=(2,−1):
Remainder, using constant third partials fxxx=6,fxxy=4,fxyy=−2,fyyy=6 and multinomial coefficients (3,0,03)=1,(2,1,03)=3,(1,2,03)=3,(0,3,03)=1 for the pure/mixed terms in h=(1,1):
Exact value: f(3,0)=27+2(9)(0)−3(0)+0=27. And indeed P2(h)+R2(h)=24+3=27. The identity is exact, confirming the formula's bookkeeping (multinomial coefficients, index collapsing) is correct.
:::counterexample[Why the segment [a,a+h] must lie inside U] Drop the hypothesis that [a,a+h]⊂U, keeping f∈C3 on its domain U. Let
U=R2∖{(0,y):y≥0}
(the plane with the non-negative y-axis removed), and let f(x,y)=arg(x+iy)∈(−π,π) be a continuous branch of the argument function on U (well-defined and C∞ there, since removing a ray from the plane makes the multivalued argument single-valued and smooth). Take a=(1,−1), so f(a)=arg(1−i)=−π/4, and take h=(−2,2), so a+h=(−1,1) and f(a+h)=arg(−1+i)=3π/4. The straight segment from (1,−1) to (−1,1) is the set {(t,−t):t∈[−1,1]} — wait, more carefully, parametrize a+th=(1−2t,−1+2t) for t∈[0,1]; at t=1/2 this is (0,0), which lies on the boundary x=0 but is not in U (since y=0≥0 is exactly the removed ray's endpoint), so the segment actually exits U near t=1/2 (for t slightly less than 1/2, x=1−2t>0 small and y=−1+2t<0 close to the boundary, but for t slightly more than 1/2, x<0 — the path crosses arbitrarily close to the removed ray at points with x=0,y>0 only in the limit; to make the violation unambiguous, perturb h slightly). Sharper version: take a=(1,−ϵ) and a+h=(−1,ϵ) for small ϵ>0, so the segment passes through (0,0)∈/U exactly at t=1/2 regardless of ϵ, hence [a,a+h]⊂U. As ϵ→0+, f(a)=arg(1−iϵ)→0− while f(a+h)=arg(−1+iϵ)→π−. Every partial derivative of f (gradient, Hessian, third derivatives) stays uniformly bounded near a as ϵ→0 (they depend smoothly on a, which stays near (1,0), away from the branch cut), so the right-hand side P2(h)+R2(h), if the theorem's formula were blindly applied using derivatives only at a, would predict f(a+h)−f(a) to be a small, bounded quantity — yet the true difference f(a+h)−f(a)→π as ϵ→0, a jump of order π, not small. The naive formula fails outright because the straight path from a to a+h leaves the domain of smoothness of f: there is no valid θ∈(0,1) with a+θh∈U at the crossing point, so the Lagrange-form remainder cannot even be evaluated there, and no true equality holds. This shows the segment hypothesis is not a technicality: it is what allows the one-variable reduction g(t)=f(a+th) to be defined and differentiable on all of [0,1] in the first place.
Scope, edge cases & common traps
Situation
What happens
f∈C2 only (not C3)
Second-order polynomial P2 still exists and equals the stated quadratic; but Lagrange remainder (with a specificθ) is not guaranteed — only the weaker Peano form R2(h)=o(∥h∥2) holds, via a more delicate argument.
U not convex, but segment [a,a+h]⊂U anyway
Theorem still applies verbatim — convexity of U is never used, only that the one particular segment needed lies in U.
n=1
Recovers exactly the classical single-variable Taylor theorem with Lagrange remainder; the Hessian becomes the scalar f′′(a).
Higher order k>2
Identical proof strategy (reduce to g(t)=f(a+th), apply 1-variable Taylor of order k) gives the general multi-index formula f(a+h)=∑∣α∣≤kα!Dαf(a)hα+Rk(h).
Mistaking the Hessian for asymmetric
If f∈/C2 (mixed partials fail to be continuous), ∂i∂jf can differ from ∂j∂if at isolated points, and "the" Hessian is ill-defined without specifying an order — see Clairaut's Theorem on Mixed Partials.
The mathematician behind it
Origins · 1 minBrook Taylor and how this was discovered
BT
Brook Taylor
1685–1731 · Early 18th-century England, Newtonian calculus circle
Brook Taylor was an English mathematician associated with Isaac Newton's circle at the Royal Society, active during the calculus priority disputes between Newton and Leibniz. His 1715 treatise Methodus Incrementorum Directa et Inversa introduced what is now called Taylor's theorem for functions of a single variable, expressed at the time in the language of finite differences and fluxions rather than modern limits. Multivariable generalizations, multi-index notation, and rigorous remainder estimates came much later, developed over the eighteenth and nineteenth centuries as mathematical analysis was put on a fully rigorous footing by Cauchy, Lagrange (who gave the remainder form used here), and their successors.
History
- 1715 — Brook Taylor publishes Methodus Incrementorum, stating the single-variable expansion in the language of finite differences.
- 1772 — Joseph-Louis Lagrange gives the mean-value form of the remainder, Rk=(k+1)!f(k+1)(ξ)(x−a)k+1, sharpening Taylor's original statement.
- 1797 — Lagrange's Théorie des fonctions analytiques further systematizes power series expansions and remainder estimates within a rigorous framework for its era.
- 1860s — Rigorous ϵ–δ analysis (Weierstrass and collaborators) supplies the modern proof technique of reducing multivariable statements to one-variable ones via restriction to line segments.
- 20th century — Multi-index notation, formalized for the systematic study of partial differential equations and distributions, becomes the standard language for stating Taylor's theorem, Sobolev embeddings, and related expansions in Rn.
Why do we need multi-index notation at all — couldn't we just write out the sums?
For low orders like k=2 in low dimension, writing sums explicitly (as in ∑i,j∂i∂jfhihj) is perfectly manageable. But for order k in dimension n, the number of distinct terms if you don't collapse repeated mixed partials grows like nk, while multi-index notation groups all index tuples giving the same partial derivative into one coefficient Dαf/α!, so the formula ∑∣α∣≤kα!Dαf(a)hα has a fixed, dimension-independent shape mirroring the one-variable case exactly. It also makes the multinomial coefficients k!/α! transparent, which is essential once you go beyond second order or work with several variables simultaneously, as in PDE theory.
Does the Hessian have to be symmetric for the formula to make sense?
The quadratic form hTHf(a)h=∑i,j∂i∂jf(a)hihj makes sense as a number regardless of whether ∂i∂jf=∂j∂if. But calling the matrix "the Hessian" and treating it as a well-defined symmetric object relies on Clairaut's Theorem on Mixed Partials, which requires f∈C2 near a. Without that hypothesis, ∂i∂jf(a) and ∂j∂if(a) could genuinely differ, and the "Hessian matrix" would depend on an arbitrary choice of differentiation order — a genuinely different, non-symmetric matrix for each ordering, undermining later uses like eigenvalue-based classification of critical points.
Is the Taylor polynomial P2 the unique quadratic polynomial approximating f to second order at a?
Yes. If Q(h) is any polynomial of degree ≤2 in h satisfying f(a+h)−Q(h)=o(∥h∥2) as h→0, then Q=P2. This follows by matching Taylor coefficients: subtracting two such polynomials gives a degree-≤2 polynomial that is o(∥h∥2), and a nonzero polynomial of degree ≤2 cannot be o(∥h∥2) (evaluate along suitable rays and use that a nonzero polynomial's lowest-degree homogeneous part dominates as h→0 along a generic direction), forcing the two polynomials to be identical. This uniqueness is why the Hessian is the correct object encoding second-order behavior, not merely one convenient encoding among many.
How does this connect to the classification of critical points?
At a critical point a (where ∇f(a)=0), the linear term of P2 vanishes, leaving f(a+h)−f(a)≈21hTHf(a)h up to o(∥h∥2). The sign behavior of this quadratic form — determined by the eigenvalues of the Hessian — dictates whether a is a local min (all eigenvalues positive), local max (all negative), saddle (mixed signs), or an inconclusive degenerate case (a zero eigenvalue). This is the entire content of the second-derivative test, developed fully in Critical Points and Extrema and built on the object introduced here, The Hessian Matrix.
What changes if we want a Taylor expansion of order higher than 2?
Nothing conceptually: the same proof — restrict to g(t)=f(a+th), apply the one-variable Taylor theorem of the corresponding order, translate derivatives of g back into partial derivatives of f via repeated chain rule — produces, for f∈Ck+1(U),
for some θ∈(0,1). The bookkeeping of multi-indices becomes more valuable as k grows, since the number of distinct mixed partials of order k in n variables grows combinatorially.