Vectorization of multifold summation to speedup

Question

I searched this website but didn't find any suitable answer describing how one can speed up summation in Mathematica using vectorization techniques and other techniques.

I often have to numerically sum over a multi-fold series of the hypergeometric type in my research work. One toy example is

lim = 150;
Sum[
  Gamma[1 + n1 + n2 + n3]/(n1! n2! n3!) (0.1)^n1 (0.1)^n2 (0.1)^
   n3, {n1, 0, lim}, {n2, 0, lim}, {n3, 0, lim}] // AbsoluteTiming

which takes about 42 sec on my laptop. The only way I know to speed-up is by using ParallelSum instead of Sum, which takes 9 sec, thanks to my 8 core processor. I want to know if there are any tricks or techniques to speed-up?

Answer 1

This may be useful as an idea, but @Henrik's observation suggests a superior approach on the case at hand. The following takes a single-core run of the OP's code from around 40s down to a little over 1s on 4 cores. Using Map instead of ParallelMap takes about three times as long. The code uses the highly efficient Outer[List,...] to construct the array of inputs and LogGamma instead of Gamma to keep the numerics in the machine-number realm.

lim = 150;
log10 = Log[0.1];
ParallelMap[
   Apply@Function[{n1, n2, n3}, 
     Total[
       LogGamma[1 + n1 + n2 + n3] - LogGamma[1 + n1] - 
         LogGamma[1 + n2] - LogGamma[1 + n3] + (n1 + n2 + n3) log10 //
         Exp, 2]],
   Transpose[
    Outer[List, Range[0., lim], Range[0., lim], Range[0., lim]],
    {1, 3, 4, 2}]
   ] // Total // AbsoluteTiming

(*  {1.14043, 1.42857}  *)

Demonstration of Henrik's comment:

lim = 20;
log10 = Log[0.1];
res2 = ParallelMap[
    Apply@
     Function[{n1, n2, n3}, 
      Total[
        LogGamma[1 + n1 + n2 + n3] - LogGamma[1 + n1] - 
          LogGamma[1 + n2] - LogGamma[1 + n3] + (n1 + n2 + n3) log10 //
          Exp, 2]],
    Transpose[
     Outer[List, Range[0., lim], Range[0., lim], Range[0., lim]],
     {1, 3, 4, 2}]
    ] // Total // AbsoluteTiming
res - res2 // Last

(*
  {0.01348, 1.42857}
  -6.66134*10^-16
*)

Answer 2

Your toy example admits another approach that may be helpful. First, note that:

$$\Gamma (n+1)=\int_0^{\infty } t^n e^{-t} \, dt$$

So, your original expression can be rewritten as (assuming the order of integration and summation can be interchanged):

$$\int_0^{\infty } e^{-t} \left(\sum _{n=0}^L \frac{a^n t^n}{n!}\right){}^3 \, \ dt$$

The sum can be computed:

Sum[(a^n t^n)/n!, {n, 0, L}]

(E^(a t) Gamma[1 + L, a t])/Gamma[1 + L]

So, your original expression is equivalent to the following integral:

$$\int_0^{\infty } e^{-t} \left(\frac{e^{a t} \Gamma (L+1,a t)}{\Gamma (L+1)}\right)^3 \, dt$$

A function that does this computation:

s[a_, L_, opts:OptionsPattern[NIntegrate]] := NIntegrate[
    Exp[-t] (Exp[a t] Gamma[L+1, a t]/Gamma[L+1])^3,
    {t, 0, Infinity},
    opts
]

Your example:

s[1/10, 150, WorkingPrecision -> 50] //AbsoluteTiming

{0.144731, 1.4285714285714285714285714285714285714285714285714}