# How best to write an exponential of differential operators?

I want to evaluate a term like this,

$$\left.\exp\left({1\over 2}\sum_{i,j=1}^{n}(A^{-1})_{ij}{\partial \over \partial x_i}{\partial \over \partial x_j}\right) f(\vec{x})\right|_{\vec{x}=0}$$

and I want to make a function that can evaluate this term given an arbitrary polynomial function $f$, with an arbitrary number of variables $x_i$.

The trouble I'm having, is compactly writing out all the derivatives that will survive the series expansion. Since the function $f$ is polynomial, we only need to take the series expansion out to a finite number of terms, equal to $\lfloor d/2\rfloor$ where $d$ is the degree of the polynomial.

$$f + \frac{1}{2} \sum_{i,j=1}^{n}(A^{-1})_{ij}{\partial \over \partial x_i}{\partial \over \partial x_j}f + \frac{1}{8}\left( \sum_{i,j=1}^{n}(A^{-1})_{ij}{\partial \over \partial x_i}{\partial \over \partial x_j}\right)^2f + ....$$

Now the mixed partial derivatives of $f$ will be symmetric and since the matrix $A$ is also symmetric, the $(A)^{-1})_{ij}=(A)^{-1})_{ji}$, but I haven't taken advantage of that in my code below.

expD[f_, vars_, ainverse_] :=
Module[{combinations, nvars, fdegree},

fdegree =
Exponent[# /. ((# -> \[FormalX] RandomReal[]) & /@
Variables[#]), \[FormalX]] &@f;
nvars = Length@vars;

(f + Sum[combinations = Tuples[Range@Length@vars, 2 m];
1/(2^m m!)
Sum[Times @@ (Extract[ainverse, #] & /@
Partition[n, 2]) (Fold[D, f, vars[[#]] & /@ n]), {n,
combinations}],
{m, 1, Floor[fdegree/2]}]) /.
Inner[Rule, vars, ConstantArray[0, nvars], List]
]


Thanks to klgr for the code to find the degree of the polynomial. You can verify that it works by using a test function and a test matrix,

testf = 4 + 8 x^4 + 12 x^3 y - y^2 x;
testainv = {{7/2, -1}, {-1, 1/2}};


and comparing the output of this function to what you get if you manually write it out,

expD[testf, {x, y}, testainv]
(testf +
1/2 (testainv[[1, 1]] D[testf, {x, 2}, {y, 0}] +
2 testainv[[1, 2]] D[testf, {x, 1}, {y, 1}] +
testainv[[2, 2]] D[testf, {x, 0}, {y, 2}]) +
1/8 (testainv[[1, 1]] testainv[[1, 1]] D[testf, {x, 4}, {y, 0}] +
4 testainv[[1, 1]] testainv[[1, 2]] D[
testf, {x, 3}, {y, 1}])) /. {x -> 0, y -> 0}
(* 172 *)
(* 172 *)

• What about iteratively applying the operator to the previous partial result? Can be done with NestList for example. Jan 28, 2016 at 16:06
• @DanielLichtblau - you mean do the terms in the sum one by one, and use the result from one step as the start of the next? I will think on that, thank you for the suggestion! Jan 28, 2016 at 16:31
• Yes, that's what I had in mind. Jan 28, 2016 at 16:33
• Closely related, possible duplicate: Exponential of a Differential Operator. If that Q&A isn't what you need, could you please point out in the question what is different from the linked one?
– Jens
Jan 28, 2016 at 20:30

Here is the shortest form of the operator I could come up with. It works for any function that can be Fourier-transformed, not just polynomials:

expD[f_, vars_, ai_] :=
Module[{k = Array[K, Length[vars]]},
FourierTransform[
Exp[-(1/2) k.ai.k] InverseFourierTransform[f, vars, k], k, vars] /.


Here is the test data from the question:

testf = 4 + 8 x^4 + 12 x^3 y - y^2 x;
testainv = {{7/2, -1}, {-1, 1/2}};
expD[testf, {x, y}, testainv]

(* ==> 172 *)


The result is the same as in the original approach. The calculation is slower (for this polynomial) because of the transforms which contain delta functions. But the advantage is that the function need not be a polynomial, it can essentially be arbitrary.

Here, I followed the general idea described in approach (1) of my answer to "Exponential of a Differential Operator".

• Jens, thanks for this great answer! I didn't get your earlier comment until this morning, I agree this could be marked as a duplicate, but if it is you should post this answer over there - assuming it can be written in terms of a general differential operator instead of just my specific operator. Jan 29, 2016 at 12:38
• I also just wanted to say that when I was at UO for grad school, I found your web pages on Mathematica to be extremely helpful. Jan 29, 2016 at 12:39

Jens's answer is great, especially how it can be applied to a general function, I just want to post here a couple solutions that are strictly only valid for polynomials.

The first takes Jens's approach (2) from this answer, which uses Fold to apply the differential operator,

expD2[f_, vars_, ainverse_] := Module[{nvars, fdegree, dop},
fdegree =
Exponent[# /. ((# -> \[FormalX] RandomReal[]) & /@
Variables[#]), \[FormalX]] &@f;
nvars = Length@vars;
dop = Function[{f2},
1/2 Sum[ainverse[[i, j]] D[
f2, {vars[[i]], 1}, {vars[[j]], 1}], {i, nvars}, {j, nvars}]];
Fold[f + dop[#1]/#2 &, f, Reverse@Range[fdegree]] /.
];


The second approach uses the method suggested by @DanielLichtblau, using NestList to repeatedly apply the differential operator to the previous result.

expD3[f_, vars_, ainverse_] := Module[{nvars, fdegree, dop},
fdegree =
Exponent[# /. ((# -> \[FormalX] RandomReal[]) & /@
Variables[#]), \[FormalX]] &@f;
nvars = Length@vars;
dop = Function[{f2},
1/2 Sum[ainverse[[i, j]] D[
f2, {vars[[i]], 1}, {vars[[j]], 1}], {i, nvars}, {j, nvars}]];
(1/Range[0, fdegree]!).NestList[dop, f, fdegree] /.
];


When I took a randomly generated large 2-variable polynomial with total degree 18 and applied Jens's Fourier method, the Fold method, and the NestList method, the timings were 113s, 5.8s, and 0.19s, respectively. So clearly there is an advantage to determining whether the function is a polynomial in the given variables prior to deciding how to proceed.

expD[f_, vars_, ainverse_] := If[PolynomialQ[f, vars],
Module[{nvars, fdegree, dop},
fdegree =
Exponent[# /. ((# -> \[FormalX] RandomReal[]) & /@
Variables[#]), \[FormalX]] &@f;
nvars = Length@vars;
dop = Function[{f2},
1/2 Sum[ainverse[[i, j]] D[
f2, {vars[[i]], 1}, {vars[[j]], 1}], {i, nvars}, {j, nvars}]];
(1/Range[0, fdegree]!).NestList[dop, f, fdegree] /.