How best to write an exponential of differential operators?

Question

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 *)

Answer 1

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] /. 
    Thread[vars -> 0]]

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".

Answer 2

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]] /. 
    Thread[vars -> 0]
   ];

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] /. 
    Thread[vars -> 0]
   ];

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] /. 
    Thread[vars -> 0]
   ],
  Module[{k = Array[K, Length[vars]]}, 
   FourierTransform[
     Exp[-(1/2) k.ainverse.k] InverseFourierTransform[f, vars, k], k, 
     vars] /. Thread[vars -> 0]]
  ]