Taylor's theorem and the symmetric group Anytime I see an $n!$ in some formula, my instinct is to look for the symmetric group on $n$ letters coming in somewhere. I have never done this seriously with the $n!$ in Taylor's theorem.
Question: Is there some way to see the $n!$ in Taylor's theorem coming naturally from a symmetry group?
Possible lead:
Here is a definition of $f^{(n)}(a)$ which does not depend on finding earlier derivatives: Let $g: \mathbb{R}^n \rightarrow \mathbb{R}$ be defined by $g(x_1,x_2,x_3, ..., x_n)$ is the lead coefficient of the unique $n^{th}$ degree polynomial passing through $(a, f(a)), (x_1, f(x_1)),(x_2, f(x_2)),...,(x_n, f(x_n))$. Then $f^{(n)}(a)$ is $1/n!$ times the limit of $g(x_1,x_2,x_3, ..., x_n)$ as $(x_1,x_2,x_3, ..., x_n)$ approaches $(a,a,a,...,a)$. Could the $1/n!$ be related to the symmetry of $g$ under exchange of coordinates?
 A: There must be many ways to think of this. Here's one: 
The symmetric group is involved with homogeneous polynomials of degree $n$ because they correspond to symmetric multilinear functions of $n$ variables, and division by $n$ factorial appears when recovering the former from the latter.
For example, a homogeneous polynomial map $f:V\to W$ of degree $3$ determines a map $F:V\times V\times V\to W$ by
$$
F(x,y,z)=f(x+y+z)-f(x+y)-f(x+z)-f(y+z)
+f(x)+f(y)+f(z)-f(0),
$$
and $f(x)$ can be recovered as 
$$\frac{F(x,x,x)}{6}.
$$
The same formula for $F(x,y,z)$ can be applied to other functions $f(x)$. Applied to a polynomial of degree $<3$ it gives zero. Applied to a polynomial of degree $\le3$ it gives the multilinear function that corresponds to the purely cubic part of $f$. Applied to any function at all it gives a symmetric function of $(x,y,z)$ that vanishes when $x$ or $y$ or $z$ is zero.
The foundation of all differential calculus is the process that takes $f$, subtracts $f(0)$, and then makes the best linear approximation of the result. If you apply this process in all three variables to $f(x+y+z)$, you are first making $F(x,y,z)$ and then (assuming that $f$ was sufficiently smooth) making best linear approximation in all variables. If you then set $x=y=z$ and divide by $3$ factorial, you get the third term in the Taylor series of $f$.
Of course, if you perform the linearizations in the three variables one after another then you see this trilinear map as a derivative of a  derivative of a derivative. 
The reason I look at it this way is that a homotopical categorical analogue of this plays a big role in my "functor calculus". The analogue there of division by $n$ factorial is a homotopy orbit spectrum for an action of the symmetric group. 
In a little more detail: 
Let $V$ and $W$ be model categories in which filtered homotopy colimits commute with finite homotopy limits. Assume that $W$ is stable, meaning that the final object is equivalent to initial object and that homotopy pushout squares are the same as homotopy pullback squares. For example, $W$ might be the category of spectra and $V$ might be spaces or spectra. Consider functors $f:V\to W$ that are homotopy-invariant (preserve weak equivalences). Say that $f$ has degree $\le 1$ if it preserves homotopy pushout squares. Say that it has degree $\le d$ if it takes those $(d+1)$-dimensional square diagrams all of whose $2$-dimensional faces are homotopy pushouts to homotopy pushout cubes, those in which the last object is the homotopy colimit of the others. Call $f$ reduced if it takes the trivial object $\star$ to itself (up to equivalence). Call it linear if it has degree $\le 1$ and is reduced. There is a linearization process that takes a reduced functor and makes the universal linear functor under it: basically $hocolim \Omega^nf\Sigma^n$. This can be easily generalized to make something called $P_1$ (first polynomial approximation) which takes any functor $f$, reduced or not, and makes the universal degree $\le 1$ functor $P_1f$ under it. (If we first reduce $f$ by taking the homotopy fiber of $f\to f(\star)$ and then linearize, it's the same as first doing $P_1$ and then reducing.) It can be further generalized to make something called $P_d$ ($d$th polynomial approximation) which takes any functor $f$ and makes the universal degree $\le d$ functor $P_df$ under it. The homotopy fiber of the canonical map $P_df\to P_{d-1}f$ is always a degree $\le d$ functor such that $P_{d-1}$ of it is trivial. Call such functors homogeneous of degree $d$. When $d=1$ this is the same as linear. 
There is the functor of $d$ variables $(X_1,\dots,X_d)\mapsto f(X_1+\dots +X_d)$, where $+$ means (derived) coproduct. This is symmetric in the sense that there are isomorphisms when you permute the inputs, satisfying the obvious identities. Reduce it in all variables simultaneously, by first making a cubical diagram consisting of the objects $f(X_S)$, coproduct of all the $X_i$ for $i\in S$, and then looking at the homotopy fiber of the canonical  map from $f(X_1+\dots +X_d)$ to the holim of all the rest. I call this reduced symmetric functor of $d$ variables the $d$th crosseffect of $f$. If we simultaneously linearize it in all variables, we get a symmetric multilinear functor.
If $f$ has degree $\le (d-1)$ then its $d$th crosseffect is contractible. If $f$ is homogeneous of degree $d$ then its $d$th crosseffect is already multilinear, and this way of making a symmetric multilinear functor from a homogeneous degree $d$ functor can be inverted: $f(X)$ is the homotopy orbit object for the action of the $d$th symmetric group on $F(X,\dots ,X)$, where the action on $F(X,\dots ,X)$ is the one that you get from the symmetry on $F(X_1,\dots ,X_d)$.
For general $f$ the multilinearization of the crosseffect is the same symmetric multilinear functor that corresponds in this way to the $d$th homogeneous part of $f$, i.e. to the homotopy fiber of $P_df\to P_{d-1}f$.
Again, by linearizing in one variable at a time you can work out a sense in which the multilinear functor corresponding to the $d$th homogeneous part is in fact an iterated derivative. But in this context it is important to have the option of performing linearization in all variables simultaneously, so to speak, in order that the result should be symmetric in the sense that is required for reconstructing the corresponding homogeneous functor. (Here symmetry is a structure, not a property.)
A: Seems to me the question is just a mighty complicated way to ask "how does $n!$ creep into the $n$-th derivative (analytic or formal) of $x^n$", to which the answer is ob course as Will Sawin noted "through repeated application of the product rule".
I think the question should say Taylor's formula, not Taylor's theorem (and this is true for many uses of Taylor's theorem elsewhere). Taylor did not have any theorem, and what is nowadays called Taylor's theorem (any of the many forms) can be stated without mentioning $n!$ at all, after subtracting off the uninteresting polynomial part (once we've learned about linearity and how to differentiate monomials): "if $f$ is $n$ times differentiable at $a$ with its value and all $n$ derivatives zero at $a$ (and possible additional regularity hypothesis in a neighbourhood of $a$) then (the error term) $f(x)$ satisfies (you favorite estimate or formula) for $x$ in some (specific) neighbourhood of $a$".
And if the question is about relating the coefficients of the Taylor polynomial to the derivatives, then it really boils down to understanding the differentiation of monomials.
A: *

*Don't you mean to say that $f^{(n)}(a)$ is $n!$ times the lead coefficient? For instance, take $f=x^n$.

*The functions $f_k(n)=\frac{n!}{(n-k)!k!}$ satisfy a discrete version of the derivative. $f_k(n)-f_k(n-1)=f_{k-1}(n)$. If we view Taylor's theorem as being about approximating with a very fine version of these functions, we just have to see where and why $k!$ arises in their combinatorics.

*Fundamentally, we can break Taylor's theorem up into two parts - that functions are approximable by polynomials, which with your definition of $g$ could create an $n!$-free version, and that the $n$th derivative of $x^n$ is $n!$. The first part has no $n!$ in it. The second one doesn't seem very deep, but is very clearly about $S_n$. Take $x \cdot x \cdot x \cdot ... \cdot x$, and differentiate it $n$ times. By the product rule, we have to choose an order to differentiate the $x$s in. Because the second derivative of $x$ is zero, we must avoid repetition. Clearly, this is where your permutations come from.
A: One way is to use a combinatorial definition of the derivative. Let $A(z) = \sum a_n z^n$ be a power series. In combinatorics, where $A$ is likely to be an ordinary generating function, $a_n$ is likely to count the number of possible combinatorial structures of some kind on an ordered set with $n$ elements. For example when $a_n = 1$ the structure is "being an ordered set," when $a_n = a^n$ the structure is "being an ordered set together with a choice, for each element of the order, of a letter from an alphabet of size $a$," and so forth.
Then a combinatorial definition of the derivative $A'(z) = \sum na_n z^{n-1}$ is as follows: the derivative is an operation that given an $A$-structure produces a new type of structure, an $A'$-structure. An $A'$-structure is "being an ordered set together with an extra element, and an $A$-structure on the new order." This is because for an ordered set with $n$ elements there are $n+1$ possibilities for the extra element and $a_{n+1}$ possibilities for the $A$-structure on the new order. 
Then after applying the derivative $k$ times we have added $k$ new elements. If we start from the empty set then there are $k!$ ways to do this, and the set of such ways to do this is naturally a torsor for the symmetric group $S_k$. (A simple example of one benefit of categorification: when you replace numbers by sets, they can support more complicated structures such as group actions.) 
One nice property of this definition is that it offers a conceptual interpretation of the Leibniz rule. First, recall that if $A, B$ are two generating functions for $A$-structures and $B$-structures, then an $AB$-structure is a partition of an ordered set into a first segment and a second segment together with an $A$-structure on the first segment and a $B$-structure on the second segment, and the corresponding generating function is the product $AB$. Now, the derivative $(AB)'$ counts the number of ways to add one element to an ordered set and then put an $AB$-structure on the new order. The new element may be either in the first segment or the second segment, and this gives the two terms $A B'$ and $A' B$ in the Leibniz rule.

There is a second, possibly more satisfying, explanation using the theory of combinatorial species and groupoid cardinality. I seem to recall that this explanation was written up by John Baez somewhere in This Week's Finds, but I can't currently find it.
A: I think the most obvious way to see the symmetric group appearing is the following.  To calculate the difference between $f(x)$ and the constant $f(0)$
$f(x) - f(0) = \int_0^x f'(s_1) ds_1$
Often we do this move because $x$ is small, so we figure that a small change in the input to $f$ will result in the output changing depending on the derivative of $f$.  But the same logic goes for $f'$, so applying the same formula to $f'$, we get:
$f(x) = f(0) + f'(0) x + \int_{0 < s_2 < s_1 < x} f''(s_2) ds_2 ds_1$
Note $f'$ may not change significantly from its initial value if $f''$ is under control, but there can be oscillations at the level of the next derivative responsible for $f'$ experiencing a small change -- in this case Taylor expansion usually is not helping us understand the problem.  In any case, here we have an integral over the region $ \{ 0 < s_2 < s_1 < x \} $.  This is a fundamental domain for the action of $S_2$ on the square $\{ 0 < s_1, s_2 < x\}$, so if $f''$ is a constant $f''(0)$, you get the volume $x^2 / 2!$ times $f''(0)$.
Even if $f''$ is not constant, you can certainly write $f''(s_2) = f''(0) + \int_0^{s_2} f''(s_3) ds_3 $ and get $f''(0) x^2 / 2!$ plus an integral over $\{ 0 < s_3 < s_2 < s_1 < x \}$, which is again a fundamental domain for the action of $S_3$ therefore has volume $x^3/3!$.  Similar considerations apply to the higher order terms.
Normally we replace the use of Fubini's theorem here with an equivalent integration by parts so that we never see iterated integrals.  (Fubini's theorem is one way to prove integration by parts, so it really is the same move.)
A: Here is an observation related to Taylor's Formula. Unfortunately, it doesn't involve the symmetric group. 
If $\mathcal A$ is the category of commutative $\mathbb Q$-algebras and $\mathcal D$ the category of commutative differential $\mathbb Q$-algebras, then the forgetful functor $\mathcal D\to\mathcal A$ has a right adjoint, and the adjunction is given by Taylor's Formula. 
More precisely, the right adjoint attaches $(A[[X]],d/dX)$ to $A$, and the adjunction maps $\varphi:D\to A$ to $\Phi:D\to A[[X]]$ with
$$
\Phi(f)=\sum_{n\ge0}\ \frac{\varphi(f^{(n)})}{n!}\ X^n.
$$ 
(The map $\varphi\mapsto\Phi$ is of course the evaluation at $0$.)
