# Fast way to the Taylor series expansion coefficients of multivariable function?

Is there a fast method to get the coefficients of Taylor series expansion of function $$f(x_1,x_2,...,x_d)$$ with maximal summed partial derivative up to $$n$$, where $$d,n$$ can be relatively large (for example, $$d=10,n=10$$)?

For example, in mathematica, the series expansion

n = 1; Series[f[x, y], {x, 0, n}, {y, 0, n}]


$$\left(f(0,0)+y f^{(0,1)}(0,0)+O\left(y^2\right)\right)+x \left(f^{(1,0)}(0,0)+f^{(1,1)}(0,0) y+O\left(y^2\right)\right)+O\left(x^2\right)$$

I need the coefficient list $$(1, x, y)$$ for function list: ($$f(0,0), f^{(1,0)}(0,0), f^{(0,1)}(0,0)$$)

The conventional Series function is very slow, for example,

In[655]:=n = 1; Series[f[x1, x2, x3, x4, x5, x6, x7], {x1, 0, n}, {x2, 0, n}, {x3, 0, n}, {x4, 0, n}, {x5, 0, n}, {x6, 0, n}, {x7, 0, n}]; // AbsoluteTiming

Out[655]= {0.000259252, Null}

In[656]:= n = 2; Series[f[x1, x2, x3, x4, x5, x6, x7], {x1, 0, n}, {x2, 0, n}, {x3, 0, n}, {x4, 0, n}, {x5, 0, n}, {x6, 0, n}, {x7, 0, n}]; // AbsoluteTiming

Out[656]= {0.442849, Null}

In[657]:= n = 3; Series[f[x1, x2, x3, x4, x5, x6, x7], {x1, 0, n}, {x2, 0, n}, {x3, 0, n}, {x4, 0, n}, {x5, 0, n}, {x6, 0, n}, {x7, 0, n}]; // AbsoluteTiming

Out[657]= {2.95541, Null}

In[665]:= n = 5; Series[f[x1, x2, x3, x4, x5, x6, x7], {x1, 0, n}, {x2, 0, n}, {x3, 0, n}, {x4, 0, n}, {x5, 0, n}, {x6, 0, n}, {x7, 0, n}]; // AbsoluteTiming

Out[665]= {42.6434, Null}


The series function seems to be inadequate for obtaining the coefficient list and derivative list of function with large $$(n,d)$$.

I have found a solution in https://arxiv.org/abs/1905.02809

Mathematica code for Taylor series in several variables

Based on the multi-index, the Taylor series expansion of a multi-variable scalar function $$u(x_{1},...,x_{d})$$ at $$\mathbf {a} =(0,...,0)$$ can be written as

$$u(x_{1},...,x_{d})=\sum _{(n_{1},...,n_{d})\in \alpha _{d}^{n}}{\frac {x_{1}^{n_{1}}...x_{d}^{n_{d}}}{n_{1}!...n_{d}!}}u_{,n_{1}...n_{d}}+O(r^{n+1})$$ ,

with $$r={\sqrt {x_{1}^{2}+...+x_{d}^{2}}}$$ , $$u_{,n_{1}...n_{d}}={\frac {\partial ^{n_{1}+...+n_{d}}u}{\partial x_{1}^{n_{1}}...\partial x_{d}^{n_{d}}}}|_{x_{1}=0,...,x_{d}=0}$$ and $$\alpha _{d}^{n}$$ is the set of multi-indexes in $$d$$ dimensions with maximal single index up to $$n$$ , i.e. $$\mathbf {\alpha } _{d}^{n}=\{(n_{1},...,n_{d})|0\leq \sum _{i=1}^{d}n_{i}\leq n,\,n_{i}\in \mathbb {N} ^{0},1\leq i\leq d\},$$

which can be obtained by an efficient Mathematica code as

MultiIndexList[d_,n_]:=Block[{a,b,c},a=Subsets[Range[d+n],{d}];
Do[ c=a[[i]];b=c-1;b[[2;;]]-=c[[1;;-2]];a[[i]]=b,{i,Length[a]}];a];
(*note: d=number of spatial dimensions, n=maximal order of derivative*)

(*
(*Some Performance test*)
MultiIndexList[5, 20] // Length // AbsoluteTiming
MultiIndexList[5, 40] // Length // AbsoluteTiming
MultiIndexList[10, 10] // Length // AbsoluteTiming
MultiIndexList[200, 3] // Length // AbsoluteTiming
MultiIndexList[200, 1] // Length // AbsoluteTiming
(*
{0.5810 seconds,53130 terms}
{9.777 seconds,1221759 terms}
{1.583 seconds,184756 terms}
{89.013 seconds,1373701 terms}
{0.01451 seconds,201 terms}
*)
*)


For each multi-index $$(n_{1},...,n_{d})\in \alpha _{d}^{n}$$ , the polynomial and partial derivative are

$${\frac {x_{1}^{n_{1}}...x_{d}^{n_{d}}}{n_{1}!...n_{d}!}},\,\,u_{,n_{1}...n_{d}},\quad \forall (n_{1},...,n_{d})\in \alpha _{d}^{n}$$

When the multi-index is written explicitly, the Taylor series expansion of the multi-variable function is straightforward.

• Other than the function term, isn't the coefficient of a derivative term given by term /. Derivative[nn__][f][xx__] :> Multinomial[nn] Apply[Times, {xx}^{nn}]? Jan 23, 2020 at 16:39
• Yes, the coefficient is right. But this needs the list of nn at the first step Jan 23, 2020 at 21:34
• But you know the list nn: i.stack.imgur.com/puvUO.png Jan 26, 2020 at 21:36
• Can you explain why Jens's answer in the other thread did not suit your needs? With[{n = 5}, Normal[Series[Apply[f, {x1, x2, x3, x4, x5, x6, x7} t], {t, 0, n}]] /. t -> 1] Jan 28, 2020 at 13:22
• Additionally, you can adapt the code from this answer to generate your required multi-indices: MultiIndexList[d_, n_] := Flatten[Table[Permutations[ip] - 1, {k, n}, {ip, IntegerPartitions[k + d, {d}]}], 2] Jan 28, 2020 at 13:56

## 2 Answers

With some discussion with @J.M.willbebacksoon and @MichaelE2, I found the key to the Taylor series expansion of multi-variable function is the generation of the multi-index list. Two solutions for the multi-indice are

 MultiIndexList0[d_,n_]:=Block[{a,b,c},a=Subsets[Range[d+n],{d}];
Do[ c=a[[i]];b=c-1;b[[2;;]]-=c[[1;;-2]];a[[i]]=b,{i,Length[a]}];a];


and

 MultiIndexList1[d_, n_] := Flatten[Table[Permutations[ip] - 1, {k, n}, {ip,
IntegerPartitions[k + d, {d}]}], 2]


MultiIndexList1 is proposed by @J.M.willbebacksoon. Both MultiIndexList0 and MultiIndexList1 are very concise and fast. MultiIndexList1 is more than 10 times faster than MultiIndexList0.

The mathematica functions based on two multi-indices functions for the Taylor series expansion of a one variable or multi-variable function can be written as

TaylorSeries0[f_, {X_, X0_, n_}] :=
Block[{vars = Flatten[{X}], vars0 = Flatten[{X0}], d, alist, xlist,
dflist, xx}, d = Length[vars]; xx = vars - vars0;
alist = MultiIndexList0[d, n];
xlist = Table[1/Apply[Times, Factorial[nn]] Apply[Times, xx^nn], {nn, alist}];
dflist = Table[D[f, Sequence @@ Transpose[{vars, nn}]]/. (Rule @@@ ({vars, vars0}\[Transpose])), {nn,alist}];{xlist, dflist}];
TaylorSeries1[f_, {X_, X0_, n_}] :=
Block[{vars = Flatten[{X}], vars0 = Flatten[{X0}], d, alist, xlist,
dflist, xx}, d = Length[vars]; xx = vars - vars0;
alist = MultiIndexList1[d, n];
xlist = Table[ 1/Apply[Times, Factorial[nn]] Apply[Times, xx^nn], {nn, alist}];
dflist =Table[D[f, Sequence @@ Transpose[{vars, nn}]]/. (Rule @@@ ({vars, vars0}\[Transpose])), {nn,alist}]; {xlist, dflist}];


The coefficients and derivatives are listed seperately.

Verfication:

Single variable function:

 In[84]:= TaylorSeries0[f[x, y], {x, 1, 5}] // TeXForm

Out[84]//TeXForm=


$$\left( \begin{array}{cccccc} 1 & x-1 & \frac{1}{2} (x-1)^2 & \frac{1}{6} (x-1)^3 & \frac{1}{24} (x-1)^4 & \frac{1}{120} (x-1)^5 \\ f(1,y) & f^{(1,0)}(1,y) & f^{(2,0)}(1,y) & f^{(3,0)}(1,y) & f^{(4,0)}(1,y) & f^{(5,0)}(1,y) \\ \end{array} \right)$$

Multi-variable function:

 In[82]:= TaylorSeries0[f[x, y], {{x, y}, {1, 1}, 2}] // TeXForm

Out[82]//TeXForm=


$$\left( \begin{array}{cccccc} 1 & y-1 & \frac{1}{2} (y-1)^2 & x-1 & (x-1) (y-1) & \frac{1}{2} (x-1)^2 \\ f(1,1) & f^{(0,1)}(1,1) & f^{(0,2)}(1,1) & f^{(1,0)}(1,1) & f^{(1,1)}(1,1) & f^{(2,0)}(1,1) \\ \end{array} \right)$$

Some performance tests:

In[73]:= tf = TaylorSeries0[ f[x1, x2, x3, x4, x5, x6, x7],
{{x1, x2, x3, x4, x5, x6, x7}, Table[0, {7}], 4}]; // AbsoluteTiming

Out[73]= {0.027772, Null}

In[76]:= tf =TaylorSeries0[f[x1, x2, x3, x4, x5, x6, x7],
{{x1, x2, x3, x4, x5, x6, x7}, Table[0, {7}], 10}]; // AbsoluteTiming

Out[76]= {1.4329, Null}

In[77]:= tf =TaylorSeries1[f[x1, x2, x3, x4, x5, x6, x7],
{{x1, x2, x3, x4, x5, x6, x7},Table[0, {7}], 10}]; // AbsoluteTiming

Out[77]= {1.59547, Null}

In[75]:= tf =TaylorSeries1[f[x1, x2, x3, x4, x5, x6, x7],
{{x1, x2, x3, x4, x5, x6, x7},Table[0, {7}], 4}]; // AbsoluteTiming

Out[75]= {0.0153569, Null}

In[55]:=tf = TaylorSeries0[f[x1, x2, x3, x4, x5, x6, x7],
{{x1, x2, x3, x4, x5, x6, x7},Table[0, {7}], 20}]; // AbsoluteTiming

Out[55]= {84.3272, Null}

In[78]:= tf = TaylorSeries1[f[x1, x2, x3, x4, x5, x6, x7],
{{x1, x2, x3, x4, x5, x6, x7},Table[0, {7}], 20}]; // AbsoluteTiming

Out[78]= {116.976, Null}


The TaylorSeries0 based on MultiIndexList0 seems to be slightly faster than TaylorSeries1.

A simple answer:

(*ref/Grad  Maclaurin*)
orderedForm[poly_, var_List] :=
HoldForm[+##] & @@
MonomialList[poly,
var][[Ordering[-Total[#] & @@@ CoefficientRules[poly, var], All,
GreaterEqual]]];
(*f[x_,y_,z_] :=g[x,y,z]*)
grads = NestList[Grad[#, {x, y, z}] &, f[x, y, z], 6] /. {x -> 0,
y -> 0, z -> 0};
\[DoubleStruckCapitalX] = {x, y, z};

orderedForm[
Sum[1/(k - 1)! Dot[grads[[k]], ##] & @@
Table[\[DoubleStruckCapitalX], {k - 1}] , {k,
3}], \[DoubleStruckCapitalX]]

• Thanks. You are the expert in mathematica. But I am new to mathematica. It would be great if some examples are added. Jul 8, 2020 at 14:27