Speed up double summation

Question

Is there any way how to speed up this double summation?

(*Define Function*)
a = -0.07;
b = -0.45;
c = 2.;
d = 3.;
tau = 2.;
dt=10.^-2;
e = c*tau;
ft[t1_] = a*Exp[-t1/d] - b*Exp[-t1/e];
g[t1_, t2_, t3_] = ft[t1]*ft[t2]*ft[t3];

factor[x_, y_] := factor[x, y] = g[0., x*dt, y*dt]
sumfactor[h_, g_] := Sum[factor[x, y], {x, 0, h}, {y, 0, g}]

Module[{k = 3000}, sumfactor[k, k]]

For arguments like 3000 in the function sumfactor this is pretty slow unfortunately.

Szabolcs · Accepted Answer

Use vectorization.

factor[x_, y_] = g[0., x*dt, y*dt]

Assuming h=g=3000:

xx = ConstantArray[N@Range[0, 3000], 3001];
yy = Transpose[xx];

Total[factor[xx, yy], 2] // AbsoluteTiming
(* {0.375605, 253819.} *)

The key is to pass in a vector (or matrix) containing all possible x, y values, and operate on the entire array at once. This works well if you only use scalar operations, like here.  If you use any vector operations, like Dot, you must be careful to ensure that the operation does what you expect, even on a matrix.

Your code for comparison:

factor[x_, y_] := factor[x, y] = g[0., x*dt, y*dt]
sumfactor[h_, g_] := Sum[factor[x, y], {x, 0, h}, {y, 0, g}]

Module[{k = 3000}, sumfactor[k, k]] // AbsoluteTiming
(* {88.9004, 253819.} *)

Akku14 · Answer

There is a much faster way. Since the sumand is quite simple, you can get an analytic expression of your sum.

There are two ways, either with Rationalize or without.

su[h_, g_] =Sum[Rationalize[factor[x, y], 0], {x, 0, h}, {y, 0, g}]

su2[h_, g_] = Sum[factor[x, y], {x, 0, h}, {y, 0, g}]

Results not schown here, because to long.

Then you don't need to sum up, but just evaluate the calculated functions su or su2

N[su[2000, 2000], 30] // Timing

(*     {2.47025*10^-15, 9486.97325178896101045672676230}    *)

N[su2[2000, 2000], 30] // Timing

(*    {8.32667*10^-17, 9486.97}     *)

Exactly the same as your result

N[sumfactor[2000, 2000], 30] // Timing

(*    {17.688, 9486.97}    *)

Arnoud Buzing · Answer

In version 12.1, you can speed up the accepted answer from Szabolcs with FunctionCompile:
cf = FunctionCompile[
  Function[{Typed[k, "MachineInteger"]},
   Module[{a, b, c, d, tau, dt, e, ft, g, factor, xx, yy},
    a = -0.07;
    b = -0.45;
    c = 2.0;
    d = 3.0;
    tau = 2.0;
    dt = 0.01;
    e = c*tau;
    xx = ConstantArray[N@Range[0, k], k + 1]; yy = Transpose[xx];
    Typed[
      KernelFunction[
       Total], {TypeSpecifier["PackedArray"]["Real64", 2], 
        "MachineInteger"} -> 
       "Real64"][(a*Exp[0.0] - b*Exp[0.0])*(a*Exp[-xx*dt/d] - 
        b*Exp[-xx*dt/e])*(a*Exp[-yy*dt/d] - b*Exp[-yy*dt/e]), 2]
    ]]]

cf[3000] // AbsoluteTiming

I get about a 3x speed up this way.

Roman · Answer

I'd recommend doing the whole calculation with symbols and substituting only at the end:
ft[t1_] = a*Exp[-t1/d] - b*Exp[-t1/e];
g[t1_, t2_, t3_] = ft[t1]*ft[t2]*ft[t3];
factor[x_, y_] = factor[x, y] = g[0, x*dt, y*dt];
sumfactor[h_, g_] = Sum[factor[x, y], {x, 0, h}, {y, 0, g}] // FullSimplify
(*    (1/((-1 + E^(dt/d))^2 (-1 + E^(dt/e))^2))(a - b) E^(-((dt (d + e) (g + h))/(d e))) (a E^((dt g)/e) (-1 + E^(dt/e)) (-1 + E^((dt (1 + g))/d)) - b E^((dt g)/d) (-1 + E^(dt/d)) (-1 + E^((dt (1 + g))/e))) (a E^((dt h)/e) (-1 + E^(dt/e)) (-1 + E^((dt (1 + h))/d)) - b E^((dt h)/d) (-1 + E^(dt/d)) (-1 + E^((dt (1 + h))/e)))    *)

sumfactor[3000, 3000] //. {a -> -0.07, b -> -0.45, c -> 2, d -> 3, tau -> 2, dt -> 10^-2, e -> c*tau}
(*    9617.85    *)

In this way the calculation is accurate and immediate.

Speed up double summation

4 Answers

Add your own answers!

Ask a Question