Efficient Dyson series implementation

Question

I'm trying to implement a Dyson series \begin{array}{lcl} U(x,x_0) & = & 1 + \int_{x_0}^{x}{dy_1V(y_1)}+\int_{x_0}^x{dy_1\int_{x_0}^{y_1}{dy_2V(y_1)V(y_2)}}+\cdots \\ & &{} + \int_{x_0}^x{dy_1\int_{x_0}^{y_1}{dy_2 \cdots \int_{x_0}^{y_{n-1}}{dy_nV(y_1)V(y_2) \cdots V(y_n)}}} +\cdots \end{array} where $V(x)$ is a large matrix. (I'm aware of the solution for symmetric functions, it doesn't work in my case.)

Probably the main bottleneck is the integration but, since it involves multiple integrals, I don't see any easy way of dealing with it. Compile seems to be out of reach as well. I'll greatly appreciate any kind of advice.

Here's the code I've written so far:

myInt[exp_, x1_, rat_] := 
 Module[{mulInt, len, int, sum, s1, s2, rar, x, x2 = rat x1, r2, 
   r1},
  SetAttributes[x, NHoldFirst];
  (*to prevent switching into numerical values*)

  mulInt[y1_, y2_, n_Integer /; n > 0] := Module[{tab},
    tab = Table[{x[i], y1, x[i - 1]}, {i, 1, n - 1, 1}];
    Return@Prepend[tab, {x[0], y1, y2}]
    ];
  (*mulInt creates lists of variables to be integrated*)

  int[n_] := 
   Integrate[Dot @@ (#1 /@ #2[[;; , 1]]), Sequence @@ #2] &[exp[#] &, 
    mulInt[x1, x2, n]];
  (*maps first elements of mulInt's lists as variables of exp, 
  creates a matrix product out of all exps*)

  sum[n_] := Module[{set},
    set = Sum[int[m], {m, 1, n}];
    len = Length[set[[1]]];
    Return[IdentityMatrix[len] + set]
    ];
  (*Dyson series*)

  s1 = sum[3];
  rar = Total[Abs@Flatten[#]] &;
  r1 = rar[s1];
  (*a function for convergence check*)
  Do[s2 = s1 + int[n];
   r2 = rar[s2];
    Print[{n, r2 - r1}];
    If[Abs[r2 - r1] <= 10^-3. len^2, Break[]]; s1 = s2; 
   r1 = r2;, {n, 4, 15}];
  (*convergence test*)

  Return[s2]
  ]

Here's a simple example with the 3rd Pauli matrix:

MatrixExp[Integrate[PauliMatrix[3], {x, 1., 2.}]]
myInt[PauliMatrix[3] &, 1., 2.]
(*results:
{{2.718282, 0.}, {0., 0.3678794}}
{{2.716667, 0}, {0, 0.3666667}}
*)

Of course, MatrixExp[Integrate[function, {x, 1., 2.}]] does not work in more general cases.

The operators I work with are sparse arrays (code below) and it take ages to get the result. Typical values of M: 10 to 20 (M=2 can be a good working example). To run the code, type for instance myInt[toMT[ #, 2]& , 10^-7, 2.].

toMT[x_, M_] := 
 Module[{nMax = 8 M + 4, sw, ret, v = 0.1, t = 5 10^8., ntol, mrul},
  ntol[nn_, Mn_] := Mn - (Ceiling[nn/4] - 1.);
mrul[n_, m_, y_] := {n, m} -> y;

  sw[n_] := 
   Switch[{Mod[n, 4], n + 12 <= nMax, n + 15 <= nMax, n - 12 >= 1},
    {1, True, True, _},
    {mrul[n, n, -ntol[n, M] /x],
     mrul[n, n + 12, I v  t],
     mrul[n, n + 15, 2 v   (ntol[n, M] - 2) /x]},

    {1, True, False, _},
    {mrul[n, n, -ntol[n, M] /x],
     mrul[n, n + 12, I v  t]},

    {1, False, _, _},
    {mrul[n, n, -ntol[n, M] /x]},

    {2, _, _, _},
    {mrul[n, n, (ntol[n, M] - 1) /x],
     mrul[n, n + 1, -I t ]},

    {3, _, _, True},
    {mrul[n, n, ntol[n, M] /x],
     mrul[n, n - 12, I t v  ],
     mrul[n, n - 13, -2 v  (ntol[n, M] + 2) /x]},

    {3, _, _, False},
    {mrul[n, n, ntol[n, M] /x]},

    {0, _, _, _},
    {mrul[n, n, -(ntol[n, M] + 1) /x],
     mrul[n, n - 3, -I t ]}
    ];
  ret = Map[sw[#] &, Range[1, 8 M + 4]];
  Return[SparseArray@Flatten[ret]]]

Answer 1

One can decrease the difficulty of the problem by reducing the Dyson series to a matrix ODE. Let's start from the definition

$$ U(x,x_0) = 1 + \int_{x_0}^{x}{dy_1V(y_1)}+\int_{x_0}^x{dy_1\int_{x_0}^{y_1}{dy_2V(y_1)V(y_2)}}+\ldots $$

and take the derivative with respect to $x$

$$ \frac{\partial}{\partial x}U(x,x_0) = V(x)\left(1+\int_{x_0}^{x}{dy_2V(y_2)}+\ldots\right). $$

It has the same tail as $U(x,x_0)$, therefore $$ \frac{\partial}{\partial x}U(x,x_0) = V(x)U(x,x_0). $$

Now it is very easy to solve it with NDSolve

σ3 = PauliMatrix[3];
u0 = IdentityMatrix@Length@σ3;
{x1, x2} = {1, 2};
NDSolve[{u'[x] == σ3.u[x], u[x1] == u0}, u, {x, x1, x2}][[1, 1, 2]][x2]

{{2.71828, 0.}, {0., 0.367879}}

m = 20;
u0 = IdentityMatrix@Length@toMT[x, m];
{x1, x2} = 10^-7 {1, 2};
res = NDSolve[{u'[x] == toMT[x, m].u[x], u[x1] == u0}, 
      u, {x, x1, x2}][[1, 1, 2]][x2]; // AbsoluteTiming

{5.583331, Null}

It takes about 5 sec to obtain res for M=20!