[proofplan]
We reduce the multivariate Taylor expansion to the one-dimensional case by parameterising the line segment from $x$ to $x+h$ via the auxiliary function $g(t) = f(x+th)$. We apply the one-dimensional Taylor theorem with integral remainder to $g$, then compute the higher-order derivatives $g^{(j)}(t)$ using the chain rule and the multinomial theorem. Substituting these expressions and simplifying the factorials yields the polynomial part and the integral remainder $R_k(h)$ as stated.
[/proofplan]
[step:Parameterise the line segment and apply the one-dimensional Taylor theorem]
Define the auxiliary function
\begin{align*}
g: [0,1] &\to \mathbb{R}^m, \\
t &\mapsto f(x + th).
\end{align*}
Since $f$ is $C^{k+1}$ and $t \mapsto x + th$ is smooth, $g$ is of class $C^{k+1}$ on $[0,1]$. By the one-dimensional Taylor theorem with integral remainder applied component-wise:
\begin{align*}
g(1) = \sum_{j=0}^k \frac{g^{(j)}(0)}{j!} + \int_0^1 \frac{(1-t)^k}{k!}\, g^{(k+1)}(t) \, d\mathcal{L}^1(t).
\end{align*}
Since $g(1) = f(x+h)$, it remains to express $g^{(j)}(t)$ in multi-index notation.
[/step]
[step:Compute $g^{(j)}(t)$ via the chain rule and the multinomial theorem]
By the chain rule, the first derivative is
\begin{align*}
g'(t) = \sum_{i=1}^n h_i\, \partial_{x_i} f(x + th).
\end{align*}
By induction, the $j$-th derivative is obtained by applying the directional derivative operator $(h \cdot \nabla) = \sum_{i=1}^n h_i\, \partial_{x_i}$ a total of $j$ times:
\begin{align*}
g^{(j)}(t) = \left(\sum_{i=1}^n h_i\, \partial_{x_i}\right)^j f(x+th).
\end{align*}
By the multinomial theorem, the operator power expands as
\begin{align*}
\left(\sum_{i=1}^n h_i\, \partial_{x_i}\right)^j = \sum_{|\alpha|=j} \frac{j!}{\alpha!}\, h^\alpha\, D^\alpha.
\end{align*}
Therefore
\begin{align*}
g^{(j)}(t) = \sum_{|\alpha|=j} \frac{j!}{\alpha!}\, h^\alpha\, D^\alpha f(x+th).
\end{align*}
[guided]
The key computation is to express $g^{(j)}(t)$ in terms of the partial derivatives of $f$. Each differentiation of $g(t) = f(x+th)$ with respect to $t$ produces a factor of $h \cdot \nabla$ by the chain rule. Repeating $j$ times gives the operator $(h \cdot \nabla)^j$ applied to $f$.
To expand $(h \cdot \nabla)^j = \left(\sum_{i=1}^n h_i \partial_{x_i}\right)^j$, we use the multinomial theorem: for commuting operators $\partial_{x_1}, \ldots, \partial_{x_n}$ (Schwarz's theorem guarantees this since $f \in C^{k+1}$),
\begin{align*}
\left(\sum_{i=1}^n h_i \partial_{x_i}\right)^j = \sum_{\substack{\alpha \in \mathbb{N}_0^n \\ |\alpha| = j}} \frac{j!}{\alpha!} \prod_{i=1}^n h_i^{\alpha_i} \prod_{i=1}^n \partial_{x_i}^{\alpha_i} = \sum_{|\alpha|=j} \frac{j!}{\alpha!}\, h^\alpha\, D^\alpha.
\end{align*}
Here $\alpha! = \alpha_1! \cdots \alpha_n!$ and $h^\alpha = h_1^{\alpha_1} \cdots h_n^{\alpha_n}$.
Evaluating at $x + th$:
\begin{align*}
g^{(j)}(t) = \sum_{|\alpha|=j} \frac{j!}{\alpha!}\, h^\alpha\, D^\alpha f(x+th).
\end{align*}
[/guided]
[/step]
[step:Evaluate the polynomial terms at $t = 0$]
At $t = 0$:
\begin{align*}
g^{(j)}(0) = \sum_{|\alpha|=j} \frac{j!}{\alpha!}\, h^\alpha\, D^\alpha f(x).
\end{align*}
Substituting into the polynomial part of the one-dimensional expansion:
\begin{align*}
\sum_{j=0}^k \frac{g^{(j)}(0)}{j!} &= \sum_{j=0}^k \frac{1}{j!} \sum_{|\alpha|=j} \frac{j!}{\alpha!}\, h^\alpha\, D^\alpha f(x) \\
&= \sum_{j=0}^k \sum_{|\alpha|=j} \frac{h^\alpha}{\alpha!}\, D^\alpha f(x) \\
&= \sum_{|\alpha| \le k} \frac{D^\alpha f(x)}{\alpha!}\, h^\alpha.
\end{align*}
This matches the polynomial part in the theorem statement.
[/step]
[step:Substitute into the integral remainder and simplify]
The integral remainder from the one-dimensional expansion is
\begin{align*}
\int_0^1 \frac{(1-t)^k}{k!}\, g^{(k+1)}(t)\, d\mathcal{L}^1(t).
\end{align*}
Substituting the formula for $g^{(k+1)}(t)$:
\begin{align*}
&= \int_0^1 \frac{(1-t)^k}{k!} \sum_{|\alpha|=k+1} \frac{(k+1)!}{\alpha!}\, h^\alpha\, D^\alpha f(x+th)\, d\mathcal{L}^1(t).
\end{align*}
Exchanging the finite sum and integral by linearity:
\begin{align*}
&= \sum_{|\alpha|=k+1} \frac{(k+1)!}{k!\, \alpha!}\, h^\alpha \int_0^1 (1-t)^k\, D^\alpha f(x+th)\, d\mathcal{L}^1(t) \\
&= (k+1) \sum_{|\alpha|=k+1} \frac{h^\alpha}{\alpha!} \int_0^1 (1-t)^k\, D^\alpha f(x+th)\, d\mathcal{L}^1(t).
\end{align*}
This equals $R_k(h)$ as defined in the theorem statement.
[/step]
