Taylor's Theorem in Several Variables

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\mathbb{R}^n. 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∂2f∂xi∂xjhihj\sum_{i,j} \frac{\partial^2 f}{\partial x_i \partial x_j} h_i h_j into single symbols DαfD^\alpha f. This exactly mirrors the one-variable formula ∑kf(k)(a)k!(x−a)k\sum_k \frac{f^{(k)}(a)}{k!}(x-a)^k. Once this notation is in place, the multivariable Taylor theorem is proved by reducing it to the single-variable case. Restrict ff 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\mathbb{R}^n. 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)\alpha = (\alpha_1,\dots,\alpha_n) packages partial derivative orders into one symbol, with ∣α∣=∑αi|\alpha| = \sum \alpha_i and α!=∏αi!\alpha! = \prod \alpha_i!.
  • The second-order Taylor polynomial of ff at aa is f(a)+∇f(a)⋅h+12hTHf(a)hf(a) + \nabla f(a)\cdot h + \tfrac12 h^T H_f(a) h, where Hf(a)H_f(a) is the Hessian matrix.
  • The exact Lagrange remainder involves third partial derivatives evaluated at some point on the segment from aa to a+ha+h — its existence, not its location, is what we prove.
  • The proof works by restricting ff to a line through aa, reducing an nn-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∂jfh_ih_j\partial_i\partial_j f 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)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 aa where ∇f(a)=0\nabla f(a) = 0, the expansion collapses to f(a+h)=f(a)+12hTHf(a)h+o(∥h∥2)f(a+h) = f(a) + \tfrac12 h^T H_f(a) h + o(\|h\|^2), so the local behavior of ff near aa 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 P2P_2 at the current iterate and jumps to its minimizer, exploiting exactly this expansion; the remainder bound O(∥h∥3)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 UU is open but not convex, even if the specific segment [a,a+h][a,a+h] stays inside UU?] No — convexity of UU is never actually used in the proof. The only requirement is that the particular segment from aa to a+ha+h lies in UU; the rest of UU's shape is irrelevant. So on a non-convex open set, Taylor's theorem still holds for pairs (a,h)(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\mathbb{R}^n

A multi-index is an nn-tuple α=(α1,…,αn)\alpha = (\alpha_1, \dots, \alpha_n) of nonnegative integers. For a multi-index α\alpha and a vector h=(h1,…,hn)∈Rnh = (h_1,\dots,h_n) \in \mathbb{R}^n, define

∣α∣=α1+⋯+αn,α!=α1! α2!⋯αn!,hα=h1α1h2α2⋯hnαn.|\alpha| = \alpha_1 + \cdots + \alpha_n, \qquad \alpha! = \alpha_1!\,\alpha_2!\cdots\alpha_n!, \qquad h^\alpha = h_1^{\alpha_1} h_2^{\alpha_2}\cdots h_n^{\alpha_n}.

If ff is a function on an open set U⊆RnU \subseteq \mathbb{R}^n possessing all partial derivatives of order ∣α∣|\alpha| at a point, write

Dαf=∂∣α∣f∂x1α1∂x2α2⋯∂xnαn.D^\alpha f = \frac{\partial^{|\alpha|} f}{\partial x_1^{\alpha_1} \partial x_2^{\alpha_2} \cdots \partial x_n^{\alpha_n}}.

Let U⊆RnU \subseteq \mathbb{R}^n be open, let a∈Ua \in U, and let f:U→Rf : U \to \mathbb{R} be a function whose partial derivatives up to order k+1k+1 all exist and are continuous on UU (i.e. f∈Ck+1(U)f \in C^{k+1}(U)). For h∈Rnh \in \mathbb{R}^n such that the segment [a,a+h]⊂U[a, a+h] \subset U, the kk-th order Taylor polynomial of ff at aa is

Pk(h)=∑∣α∣≤kDαf(a)α! hα.P_k(h) = \sum_{|\alpha| \le k} \frac{D^\alpha f(a)}{\alpha!}\, h^\alpha.

For k=2k=2 this expands, using symmetry of second partials from Clairaut's Theorem on Mixed Partials, to

P2(h)=f(a)+∇f(a)⋅h+12hTHf(a) h,P_2(h) = f(a) + \nabla f(a) \cdot h + \frac{1}{2} h^T H_f(a)\, h,

where ∇f(a)=(∂1f(a),…,∂nf(a))\nabla f(a) = \left(\partial_1 f(a), \dots, \partial_n f(a)\right) is the gradient and Hf(a)=[∂i∂jf(a)]i,j=1nH_f(a) = \left[\partial_i \partial_j f(a)\right]_{i,j=1}^n is the Hessian matrix, a symmetric n×nn\times n matrix.

Intuition

In one variable, Taylor's theorem says a smooth function looks, near aa, like a polynomial in (x−a)(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\nabla f(a)\cdot h — and "curvature" becomes a quadratic form hTHf(a)hh^T H_f(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\nabla f(a)\cdot h and hTHf(a)hh^T H_f(a) h is to freeze the direction: fix hh, and look at the single real function g(t)=f(a+th)g(t) = f(a+th) for tt near 00. This is an ordinary one-variable function, so its own Taylor expansion at t=0t=0 is

g(t)=g(0)+g′(0)t+12g′′(0)t2+⋯ ,g(t) = g(0) + g'(0) t + \frac{1}{2} g''(0) t^2 + \cdots,

and the chain rule converts g′(0)g'(0) into ∇f(a)⋅h\nabla f(a)\cdot h and g′′(0)g''(0) into hTHf(a)hh^T H_f(a) h. Setting t=1t=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)=excos⁡yf(x,y)=e^{x}\cos y near a=(0,0)a=(0,0): level curves of ff (solid) versus level curves of its Taylor quadratic P2P_2 (dashed) agree to second order
x y a = (0,0) level curves of f level curves of P₂ (Taylor quadratic) Near a, dashed and solid curves nearly coincide; they separate visibly farther from a — the remainder growing.

Central theorem

TheoremSecond-Order Taylor's Theorem with Lagrange Remainder, in Rn\mathbb{R}^n

Let U⊆RnU \subseteq \mathbb{R}^n be open and let f:U→Rf: U \to \mathbb{R} be of class C3C^3 on UU (all partial derivatives up to third order exist and are continuous). Let a∈Ua \in U and let h∈Rnh \in \mathbb{R}^n be such that the closed segment

[a,a+h]={a+th:t∈[0,1]}[a, a+h] = \{ a + th : t \in [0,1] \}

is entirely contained in UU. Then there exists θ∈(0,1)\theta \in (0,1) such that

f(a+h)=f(a)+∇f(a)⋅h+12hTHf(a)h+R2(h),f(a+h) = f(a) + \nabla f(a) \cdot h + \frac{1}{2} h^T H_f(a) h + R_2(h),

where the remainder is

R2(h)=16∑∣α∣=33!α!Dαf(a+θh) hα=∑i,j,k=1n16 ∂i∂j∂kf(a+θh) hihjhk.R_2(h) = \frac{1}{6} \sum_{|\alpha| = 3} \frac{3!}{\alpha!} D^\alpha f(a+\theta h)\, h^\alpha = \sum_{i,j,k=1}^n \frac{1}{6}\, \partial_i \partial_j \partial_k f(a + \theta h)\, h_i h_j h_k.

Consequently, if M3=sup⁡x∈[a,a+h]max⁡i,j,k∣∂i∂j∂kf(x)∣M_3 = \sup_{x \in [a,a+h]} \max_{i,j,k} |\partial_i\partial_j\partial_k f(x)| is finite, then

∣R2(h)∣≤n3M36∥h∥∞3,|R_2(h)| \le \frac{n^3 M_3}{6} \|h\|_\infty^3,

and in particular R2(h)=o(∥h∥2)R_2(h) = o(\|h\|^2) as h→0h \to 0.

LemmaChain rule reduction to one variable

Under the hypotheses of the theorem, define g:[0,1]→Rg: [0,1] \to \mathbb{R} by g(t)=f(a+th)g(t) = f(a+th). Then g∈C3([0,1])g \in C^3([0,1]) and for every t∈[0,1]t \in [0,1],

g′(t)=∇f(a+th)⋅h,g′′(t)=hTHf(a+th) h=∑i,j=1n∂i∂jf(a+th) hihj,g'(t) = \nabla f(a+th)\cdot h, \qquad g''(t) = h^T H_f(a+th)\, h = \sum_{i,j=1}^n \partial_i\partial_j f(a+th)\, h_i h_j,
g′′′(t)=∑i,j,k=1n∂i∂j∂kf(a+th) hihjhk.g'''(t) = \sum_{i,j,k=1}^n \partial_i\partial_j\partial_k f(a+th)\, h_i h_j h_k.
Proof

Proof of the Lemma, by repeated application of The Chain Rule in Several Variables. Since [a,a+h]⊂U[a,a+h] \subset U and UU is open, and t↦a+tht \mapsto a+th is a smooth (affine, hence C∞C^\infty) map from [0,1][0,1] into UU, the composite g(t)=f(a+th)g(t) = f(a+th) is well defined on [0,1][0,1]; because f∈C3(U)f \in C^3(U), the chain rule guarantees g∈C3([0,1])g \in C^3([0,1]) (composition of C3C^3 with an affine — in particular C∞C^\infty — map is C3C^3). By the chain rule, writing γ(t)=a+th\gamma(t) = a + th so γ′(t)=h\gamma'(t) = h for all tt,

g′(t)=∇f(γ(t))⋅γ′(t)=∑i=1n∂if(a+th) hi=∇f(a+th)⋅h.g'(t) = \nabla f(\gamma(t)) \cdot \gamma'(t) = \sum_{i=1}^n \partial_i f(a+th)\, h_i = \nabla f(a+th)\cdot h.

This exhibits g′g' as a new function of tt, namely t↦∑ihi ∂if(a+th)t \mapsto \sum_i h_i\, \partial_i f(a+th), a finite linear combination (with constant coefficients hih_i) of functions of the form t↦∂if(a+th)t \mapsto \partial_i f(a+th). Differentiating again, applying the chain rule to each term ∂if(a+th)\partial_i f(a+th) exactly as before (since ∂if∈C2(U)\partial_i f \in C^2(U), being a component of f∈C3f \in C^3),

ddt[∂if(a+th)]=∑j=1n∂j∂if(a+th) hj.\frac{d}{dt}\big[\partial_i f(a+th)\big] = \sum_{j=1}^n \partial_j \partial_i f(a+th)\, h_j.

Hence

g′′(t)=∑i=1nhi∑j=1n∂j∂if(a+th) hj=∑i,j=1n∂i∂jf(a+th) hihj=hTHf(a+th) h.g''(t) = \sum_{i=1}^n h_i \sum_{j=1}^n \partial_j\partial_i f(a+th)\, h_j = \sum_{i,j=1}^n \partial_i \partial_j f(a+th)\, h_i h_j = h^T H_f(a+th)\, h.

Note that this double sum is exactly the quadratic form associated to the matrix Hf(a+th)=[∂i∂jf(a+th)]H_f(a+th) = [\partial_i \partial_j f(a+th)]; Clairaut's Theorem on Mixed Partials guarantees ∂i∂jf=∂j∂if\partial_i \partial_j f = \partial_j \partial_i f throughout UU (since f∈C2(U)f\in C^2(U) at least), so HfH_f 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)\partial_i\partial_j f \in C^1(U) (a consequence of f∈C3f \in C^3):

ddt[∂i∂jf(a+th)]=∑k=1n∂k∂i∂jf(a+th) hk,\frac{d}{dt}\big[\partial_i\partial_j f(a+th)\big] = \sum_{k=1}^n \partial_k\partial_i\partial_j f(a+th)\, h_k,

so

g′′′(t)=∑i,j=1nhihj∑k=1n∂k∂i∂jf(a+th) hk=∑i,j,k=1n∂i∂j∂kf(a+th) hihjhk,g'''(t) = \sum_{i,j=1}^n h_i h_j \sum_{k=1}^n \partial_k \partial_i \partial_j f(a+th)\, h_k = \sum_{i,j,k=1}^n \partial_i\partial_j\partial_k f(a+th)\, h_i h_j h_k,

where continuity of third partials (again from f∈C3f \in C^3) lets us reorder the indices in each triple product freely by repeated use of Clairaut's theorem. This proves the Lemma. ■\blacksquare

Proof of the Theorem, using one-variable Taylor's theorem with Lagrange remainder (assumed known from single-variable calculus) applied to gg. By the Lemma, g∈C3([0,1])g \in C^3([0,1]), so ordinary one-variable Taylor's theorem with Lagrange remainder applies to gg at the base point t=0t=0, evaluated at t=1t=1: there exists θ∈(0,1)\theta \in (0,1) such that

g(1)=g(0)+g′(0)⋅1+g′′(0)2⋅12+g′′′(θ)6⋅13.g(1) = g(0) + g'(0)\cdot 1 + \frac{g''(0)}{2}\cdot 1^2 + \frac{g'''(\theta)}{6}\cdot 1^3.

(This is the single-variable Taylor theorem with n=2n=2, i.e., quadratic Taylor polynomial plus Lagrange remainder for the third derivative, applied to the interval [0,1][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 ff using the Lemma:

  • g(1)=f(a+h)g(1) = f(a+h) and g(0)=f(a)g(0) = f(a), directly from the definition g(t)=f(a+th)g(t) = f(a+th).
  • g′(0)=∇f(a)⋅hg'(0) = \nabla f(a) \cdot h, by the Lemma's first formula at t=0t=0.
  • g′′(0)=hTHf(a)hg''(0) = h^T H_f(a) h, by the Lemma's second formula at t=0t=0.
  • g′′′(θ)=∑i,j,k∂i∂j∂kf(a+θh) hihjhkg'''(\theta) = \sum_{i,j,k} \partial_i\partial_j\partial_k f(a+\theta h)\, h_i h_j h_k, by the Lemma's third formula at t=θt=\theta.

Substituting these four identities into the displayed one-variable Taylor expansion of gg gives exactly

f(a+h)=f(a)+∇f(a)⋅h+12hTHf(a)h+16∑i,j,k=1n∂i∂j∂kf(a+θh) hihjhk⏟R2(h),f(a+h) = f(a) + \nabla f(a)\cdot h + \frac{1}{2} h^T H_f(a) h + \underbrace{\frac{1}{6}\sum_{i,j,k=1}^n \partial_i\partial_j\partial_k f(a+\theta h)\, h_i h_j h_k}_{R_2(h)},

which is the claimed formula for R2(h)R_2(h) (the coefficient 3!α!\tfrac{3!}{\alpha!} in the multi-index form arises from collecting repeated indices i,j,ki,j,k into a single multi-index α\alpha with ∣α∣=3|\alpha|=3, a standard multinomial bookkeeping identity: the number of ordered triples (i,j,k)(i,j,k) collapsing to a given α\alpha is 3!α!\tfrac{3!}{\alpha!}). For the bound: since each of the n3n^3 terms ∂i∂j∂kf(a+θh) hihjhk\partial_i\partial_j\partial_k f(a+\theta h)\, h_i h_j h_k satisfies

∣∂i∂j∂kf(a+θh) hihjhk∣≤M3∥h∥∞3\left|\partial_i\partial_j\partial_k f(a+\theta h)\, h_i h_j h_k\right| \le M_3 \|h\|_\infty^3

(using ∣hi∣,∣hj∣,∣hk∣≤∥h∥∞|h_i|,|h_j|,|h_k| \le \|h\|_\infty and the definition of M3M_3 as a supremum over the segment, which contains a+θha+\theta h since θ∈(0,1)\theta\in(0,1)), the triangle inequality applied to the triple sum over i,j,k∈{1,…,n}i,j,k \in \{1,\dots,n\} gives

∣R2(h)∣≤16∑i,j,k=1n∣∂i∂j∂kf(a+θh)∣∣hi∣∣hj∣∣hk∣≤16⋅n3⋅M3⋅∥h∥∞3=n3M36∥h∥∞3.|R_2(h)| \le \frac{1}{6}\sum_{i,j,k=1}^n \left|\partial_i\partial_j\partial_k f(a+\theta h)\right| |h_i||h_j||h_k| \le \frac{1}{6} \cdot n^3 \cdot M_3 \cdot \|h\|_\infty^3 = \frac{n^3 M_3}{6}\|h\|_\infty^3.

Since f∈C3(U)f \in C^3(U) and [a,a+h][a,a+h] is a fixed compact segment, M3<∞M_3 < \infty automatically (continuous functions on compact sets are bounded), so this bound is always available for hh small enough that [a,a+h]⊂U[a,a+h]\subset U. Finally, since ∥h∥∞≤∥h∥\|h\|_\infty \le \|h\| (Euclidean norm) up to the fixed constant n\sqrt n, we get ∣R2(h)∣=O(∥h∥3)=o(∥h∥2)|R_2(h)| = O(\|h\|^3) = o(\|h\|^2) as h→0h \to 0, because ∥h∥3/∥h∥2=∥h∥→0\|h\|^3/\|h\|^2 = \|h\| \to 0. ■\blacksquare

∎
CorollarySecond-order approximation

Under the hypotheses of the theorem, as h→0h \to 0,

f(a+h)=f(a)+∇f(a)⋅h+12hTHf(a)h+o(∥h∥2).f(a+h) = f(a) + \nabla f(a)\cdot h + \frac{1}{2}h^T H_f(a) h + o(\|h\|^2).

Worked examples

Example 1: A non-polynomial function at an asymmetric point. Let f(x,y)=excos⁡yf(x,y) = e^x \cos y (globally C∞C^\infty on all of R2\mathbb{R}^2, in particular C3C^3), and expand about the deliberately asymmetric base point a=(1,π/3)a = (1, \pi/3), i.e. compute P2(h)P_2(h) for h=(h1,h2)h=(h_1,h_2) near 00. First compute the needed partials:

fx=excos⁡y,fy=−exsin⁡y,fxx=excos⁡y,fyy=−excos⁡y,fxy=−exsin⁡y.f_x = e^x\cos y, \quad f_y = -e^x \sin y, \quad f_{xx} = e^x\cos y, \quad f_{yy} = -e^x\cos y, \quad f_{xy} = -e^x \sin y.

At a=(1,π/3)a=(1,\pi/3): cos⁡(π/3)=12\cos(\pi/3) = \tfrac12, sin⁡(π/3)=32\sin(\pi/3) = \tfrac{\sqrt3}{2}, so

f(a)=e⋅12=e2,fx(a)=e2,fy(a)=−e⋅32=−e32,f(a) = e\cdot\tfrac12 = \tfrac{e}{2}, \qquad f_x(a) = \tfrac{e}{2}, \qquad f_y(a) = -e\cdot\tfrac{\sqrt3}{2} = -\tfrac{e\sqrt3}{2},
fxx(a)=e2,fyy(a)=−e2,fxy(a)=−e32.f_{xx}(a) = \tfrac{e}{2}, \qquad f_{yy}(a) = -\tfrac{e}{2}, \qquad f_{xy}(a) = -\tfrac{e\sqrt3}{2}.

So the gradient is ∇f(a)=(e2,−e32)\nabla f(a) = \left(\tfrac{e}{2}, -\tfrac{e\sqrt3}{2}\right) and the Hessian is

Hf(a)=(e2−e32−e32−e2).H_f(a) = \begin{pmatrix} \tfrac{e}{2} & -\tfrac{e\sqrt3}{2} \\[2pt] -\tfrac{e\sqrt3}{2} & -\tfrac{e}{2} \end{pmatrix}.

By the theorem, writing h=(h1,h2)h=(h_1,h_2),

P2(h)=e2+e2h1−e32h2+12(e2h12−e3 h1h2−e2h22).P_2(h) = \frac{e}{2} + \frac{e}{2}h_1 - \frac{e\sqrt3}{2} h_2 + \frac12\left(\frac{e}{2}h_1^2 - e\sqrt3\, h_1 h_2 - \frac{e}{2}h_2^2\right).

Simplify:

P2(h)=e2[1+h1−3 h2+h12−h222−3 h1h2].P_2(h) = \frac{e}{2}\left[1 + h_1 - \sqrt3\, h_2 + \frac{h_1^2 - h_2^2}{2} - \sqrt3\, h_1 h_2\right].

Numerical check. Take h=(0.1,−0.2)h = (0.1, -0.2), i.e. approximate f(1.1,π/3−0.2)f(1.1, \pi/3 - 0.2). Exact value: x=1.1x=1.1, y=π/3−0.2≈1.04720−0.2=0.84720y = \pi/3 - 0.2 \approx 1.04720 - 0.2 = 0.84720. Compute e1.1≈3.00417e^{1.1} \approx 3.00417 and cos⁡(0.84720)≈0.66235\cos(0.84720) \approx 0.66235. So f(a+h)≈3.00417×0.66235≈1.99039f(a+h) \approx 3.00417 \times 0.66235 \approx 1.99039. Taylor estimate: plug h1=0.1h_1 = 0.1, h2=−0.2h_2=-0.2 into P2(h)P_2(h) with e/2≈1.35914e/2 \approx 1.35914:

1+0.1−3(−0.2)+0.01−0.042−3(0.1)(−0.2)=1+0.1+0.34641−0.015+0.034641=1.46605.1 + 0.1 - \sqrt3(-0.2) + \frac{0.01 - 0.04}{2} - \sqrt3(0.1)(-0.2) = 1 + 0.1 + 0.34641 - 0.015 + 0.034641 = 1.46605.

So P2(h)≈1.35914×1.46605≈1.99277P_2(h) \approx 1.35914 \times 1.46605 \approx 1.99277. The two agree to about 0.00240.0024 absolute error, and the theorem guarantees this discrepancy is O(∥h∥3)O(\|h\|^3) — with ∥h∥∞=0.2\|h\|_\infty = 0.2, indeed 0.23=0.0080.2^3 = 0.008 is the right order of magnitude, comfortably bounding the observed error once the constant n3M36\tfrac{n^3M_3}{6} is included.

Example 2: Verifying the remainder bound explicitly for a cubic. Let f(x,y)=x3+2x2y−xy2+y3f(x,y) = x^3 + 2x^2y - xy^2 + y^3, expanded at a=(2,−1)a = (2,-1) (again asymmetric, with mixed signs), evaluated at h=(1,1)h=(1,1) (so target point (3,0)(3,0)). Since ff is a cubic polynomial, all fourth-order partials vanish, so the Lagrange remainder R2(h)R_2(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,f_x = 3x^2 + 4xy - y^2, \quad f_y = 2x^2 - 2xy + 3y^2,
fxx=6x+4y,fyy=−2x+6y,fxy=4x−2y,f_{xx} = 6x + 4y, \quad f_{yy} = -2x+6y, \quad f_{xy} = 4x - 2y,
fxxx=6,fxxy=4,fxyy=−2,fyyy=6.f_{xxx} = 6, \quad f_{xxy} = 4, \quad f_{xyy} = -2, \quad f_{yyy} = 6.

All third partials are constant, so evaluating them at a+θha+\theta h for any θ\theta gives the same values — meaning the remainder is exact for any choice of θ∈(0,1)\theta \in (0,1), consistent with the theorem (which only asserts existence of some θ\theta, but here every θ\theta works because the third derivatives don't vary). At a=(2,−1)a=(2,-1):

f(a)=8+2(4)(−1)−2(1)+(−1)=8−8−2−1=−3.f(a) = 8 + 2(4)(-1) - 2(1) + (-1) = 8 - 8 - 2 - 1 = -3.
fx(a)=3(4)+4(2)(−1)−1=12−8−1=3,fy(a)=2(4)−2(2)(−1)+3(1)=8+4+3=15.f_x(a) = 3(4) + 4(2)(-1) - 1 = 12 - 8 - 1 = 3, \qquad f_y(a) = 2(4) - 2(2)(-1) + 3(1) = 8+4+3=15.
fxx(a)=12−4=8,fyy(a)=−4−6=−10,fxy(a)=8+2=10.f_{xx}(a) = 12 - 4 = 8, \qquad f_{yy}(a) = -4-6=-10, \qquad f_{xy}(a) = 8+2 = 10.

So

P2(h)=−3+(3)(1)+(15)(1)+12(8(1)2+2(10)(1)(1)+(−10)(1)2)=−3+3+15+12(8+20−10)=15+9=24.P_2(h) = -3 + (3)(1)+(15)(1) + \frac12\left(8(1)^2 + 2(10)(1)(1) + (-10)(1)^2\right) = -3+3+15+\frac12(8+20-10) = 15 + 9 = 24.

Remainder, using constant third partials fxxx=6, fxxy=4, fxyy=−2, fyyy=6f_{xxx}=6,\,f_{xxy}=4,\,f_{xyy}=-2,\,f_{yyy}=6 and multinomial coefficients (33,0,0)=1,(32,1,0)=3,(31,2,0)=3,(30,3,0)=1\binom{3}{3,0,0}=1,\binom{3}{2,1,0}=3,\binom{3}{1,2,0}=3,\binom{3}{0,3,0}=1 for the pure/mixed terms in h=(1,1)h=(1,1):

R2(h)=16[6⋅13+3(4)(1)2(1)+3(−2)(1)(1)2+6⋅13]=16[6+12−6+6]=186=3.R_2(h) = \frac{1}{6}\left[6\cdot1^3 + 3(4)(1)^2(1) + 3(-2)(1)(1)^2 + 6\cdot1^3\right] = \frac16\left[6+12-6+6\right] = \frac{18}{6}=3.

Exact value: f(3,0)=27+2(9)(0)−3(0)+0=27f(3,0) = 27 + 2(9)(0) - 3(0) + 0 = 27. And indeed P2(h)+R2(h)=24+3=27P_2(h) + R_2(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][a,a+h] must lie inside UU] Drop the hypothesis that [a,a+h]⊂U[a,a+h] \subset U, keeping f∈C3f \in C^3 on its domain UU. Let

U=R2∖{(0,y):y≥0}U = \mathbb{R}^2 \setminus \{(0,y) : y \ge 0\}

(the plane with the non-negative yy-axis removed), and let f(x,y)=arg⁡(x+iy)∈(−π,π)f(x,y) = \arg(x+iy) \in (-\pi, \pi) be a continuous branch of the argument function on UU (well-defined and C∞C^\infty there, since removing a ray from the plane makes the multivalued argument single-valued and smooth). Take a=(1,−1)a = (1,-1), so f(a)=arg⁡(1−i)=−π/4f(a) = \arg(1-i) = -\pi/4, and take h=(−2,2)h = (-2, 2), so a+h=(−1,1)a+h = (-1,1) and f(a+h)=arg⁡(−1+i)=3π/4f(a+h) = \arg(-1+i) = 3\pi/4. The straight segment from (1,−1)(1,-1) to (−1,1)(-1,1) is the set {(t,−t):t∈[−1,1]}\{(t,-t): t\in[-1,1]\} — wait, more carefully, parametrize a+th=(1−2t, −1+2t)a+th = (1-2t,\,-1+2t) for t∈[0,1]t\in[0,1]; at t=1/2t=1/2 this is (0,0)(0,0), which lies on the boundary x=0x=0 but is not in UU (since y=0≥0y=0\ge 0 is exactly the removed ray's endpoint), so the segment actually exits UU near t=1/2t=1/2 (for tt slightly less than 1/21/2, x=1−2t>0x = 1-2t>0 small and y=−1+2t<0y=-1+2t<0 close to the boundary, but for tt slightly more than 1/21/2, x<0x<0 — the path crosses arbitrarily close to the removed ray at points with x=0,y>0x=0, y>0 only in the limit; to make the violation unambiguous, perturb hh slightly). Sharper version: take a=(1,−ϵ)a = (1,-\epsilon) and a+h=(−1,ϵ)a+h = (-1,\epsilon) for small ϵ>0\epsilon>0, so the segment passes through (0,0)∉U(0,0)\notin U exactly at t=1/2t=1/2 regardless of ϵ\epsilon, hence [a,a+h]⊄U[a,a+h]\not\subset U. As ϵ→0+\epsilon \to 0^+, f(a)=arg⁡(1−iϵ)→0−f(a) = \arg(1-i\epsilon) \to 0^- while f(a+h)=arg⁡(−1+iϵ)→π−f(a+h) = \arg(-1+i\epsilon) \to \pi^-. Every partial derivative of ff (gradient, Hessian, third derivatives) stays uniformly bounded near aa as ϵ→0\epsilon\to0 (they depend smoothly on aa, which stays near (1,0)(1,0), away from the branch cut), so the right-hand side P2(h)+R2(h)P_2(h) + R_2(h), if the theorem's formula were blindly applied using derivatives only at aa, would predict f(a+h)−f(a)f(a+h) - f(a) to be a small, bounded quantity — yet the true difference f(a+h)−f(a)→πf(a+h)-f(a) \to \pi as ϵ→0\epsilon\to 0, a jump of order π\pi, not small. The naive formula fails outright because the straight path from aa to a+ha+h leaves the domain of smoothness of ff: there is no valid θ∈(0,1)\theta \in (0,1) with a+θh∈Ua+\theta h \in 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)g(t)=f(a+th) to be defined and differentiable on all of [0,1][0,1] in the first place.

Scope, edge cases & common traps

SituationWhat happens
f∈C2f \in C^2 only (not C3C^3)Second-order polynomial P2P_2 still exists and equals the stated quadratic; but Lagrange remainder (with a specific θ\theta) is not guaranteed — only the weaker Peano form R2(h)=o(∥h∥2)R_2(h)=o(\|h\|^2) holds, via a more delicate argument.
UU not convex, but segment [a,a+h]⊂U[a,a+h]\subset U anywayTheorem still applies verbatim — convexity of UU is never used, only that the one particular segment needed lies in UU.
n=1n=1Recovers exactly the classical single-variable Taylor theorem with Lagrange remainder; the Hessian becomes the scalar f′′(a)f''(a).
Higher order k>2k>2Identical proof strategy (reduce to g(t)=f(a+th)g(t)=f(a+th), apply 1-variable Taylor of order kk) gives the general multi-index formula f(a+h)=∑∣α∣≤kDαf(a)α!hα+Rk(h)f(a+h)=\sum_{|\alpha|\le k} \frac{D^\alpha f(a)}{\alpha!}h^\alpha + R_k(h).
Mistaking the Hessian for asymmetricIf f∉C2f\notin C^2 (mixed partials fail to be continuous), ∂i∂jf\partial_i\partial_j f can differ from ∂j∂if\partial_j\partial_i f 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=f(k+1)(ξ)(k+1)!(x−a)k+1R_k = \frac{f^{(k+1)}(\xi)}{(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 ϵ\epsilon–δ\delta 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\mathbb{R}^n.
Why do we need multi-index notation at all — couldn't we just write out the sums?

For low orders like k=2k=2 in low dimension, writing sums explicitly (as in ∑i,j∂i∂jf hihj\sum_{i,j}\partial_i\partial_jf\,h_ih_j) is perfectly manageable. But for order kk in dimension nn, the number of distinct terms if you don't collapse repeated mixed partials grows like nkn^k, while multi-index notation groups all index tuples giving the same partial derivative into one coefficient Dαf/α!D^\alpha f/\alpha!, so the formula ∑∣α∣≤kDαf(a)α!hα\sum_{|\alpha|\le k} \frac{D^\alpha f(a)}{\alpha!}h^\alpha has a fixed, dimension-independent shape mirroring the one-variable case exactly. It also makes the multinomial coefficients k!/α!k!/\alpha! 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)hihjh^T H_f(a) h = \sum_{i,j} \partial_i\partial_j f(a) h_i h_j makes sense as a number regardless of whether ∂i∂jf=∂j∂if\partial_i\partial_jf = \partial_j\partial_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∈C2f\in C^2 near aa. Without that hypothesis, ∂i∂jf(a)\partial_i\partial_jf(a) and ∂j∂if(a)\partial_j\partial_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 P2P_2 the unique quadratic polynomial approximating ff to second order at aa?

Yes. If Q(h)Q(h) is any polynomial of degree ≤2\le 2 in hh satisfying f(a+h)−Q(h)=o(∥h∥2)f(a+h) - Q(h) = o(\|h\|^2) as h→0h\to 0, then Q=P2Q = P_2. This follows by matching Taylor coefficients: subtracting two such polynomials gives a degree-≤2\le2 polynomial that is o(∥h∥2)o(\|h\|^2), and a nonzero polynomial of degree ≤2\le 2 cannot be o(∥h∥2)o(\|h\|^2) (evaluate along suitable rays and use that a nonzero polynomial's lowest-degree homogeneous part dominates as h→0h\to0 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 aa (where ∇f(a)=0\nabla f(a)=0), the linear term of P2P_2 vanishes, leaving f(a+h)−f(a)≈12hTHf(a)hf(a+h) - f(a) \approx \tfrac12 h^T H_f(a) h up to o(∥h∥2)o(\|h\|^2). The sign behavior of this quadratic form — determined by the eigenvalues of the Hessian — dictates whether aa 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)g(t)=f(a+th), apply the one-variable Taylor theorem of the corresponding order, translate derivatives of gg back into partial derivatives of ff via repeated chain rule — produces, for f∈Ck+1(U)f\in C^{k+1}(U),

f(a+h)=∑∣α∣≤kDαf(a)α!hα+∑∣α∣=k+1k+1α!Dαf(a+θh) hαf(a+h) = \sum_{|\alpha|\le k}\frac{D^\alpha f(a)}{\alpha!}h^\alpha + \sum_{|\alpha|=k+1}\frac{k+1}{\alpha!}D^\alpha f(a+\theta h)\,h^\alpha

for some θ∈(0,1)\theta\in(0,1). The bookkeeping of multi-indices becomes more valuable as kk grows, since the number of distinct mixed partials of order kk in nn variables grows combinatorially.

Before this

Where this leads