3
$\begingroup$

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?

$\endgroup$
1
  • 1
    $\begingroup$ The best way is probably to exploit the fact that the summands decay quite rapidly. If I set lim = 20 the summation is 100 times as fast and has still a precision of 15 digits... $\endgroup$ Dec 3 '20 at 17:33
3
$\begingroup$

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
*)
$\endgroup$
1
  • $\begingroup$ Thanks a lot! Indeed, Henrik's comment is useful for this toy sum, but Michael's method is useful for a general hypergeometric series. $\endgroup$
    – Epsilon
    Dec 3 '20 at 19:12
2
$\begingroup$

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}

$\endgroup$
2
  • $\begingroup$ Thanks a lot for this very nice approach. $\endgroup$
    – Epsilon
    Dec 3 '20 at 21:18
  • $\begingroup$ But for the more general case where we have at least 4-5 Euler gamma functions in the numerator, this approach might not be so efficient as we will have to perform a multi-fold integration. $\endgroup$
    – Epsilon
    Dec 3 '20 at 21:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.