Speed up double summation

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.

• Does tau2 have a value? What about dt? – Szabolcs May 8 '17 at 13:22
• Sorry, my fault. Edited my post. The values are of no importance. It is more a programming style question. – Display Name May 8 '17 at 13:25
• You can reduce it by a factor of about two by noting that factor[x,y] == factor[y,x], so you need only iterate for y > x. – rogerl May 8 '17 at 13:28

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.

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

• There is some additional information. Exp is Listable, otherwise, it won't work. If Exp is not a Listable Function and Exp can't accept a list as parameter, then you can use SetAttributes[f,{Listable}] to make it Listable. – UnchartedWorks May 8 '17 at 13:48
• Good to know, thank you! – Display Name May 8 '17 at 13:58
• @DisplayName That's a bit misleading. The Listable attribute on its own brings no improvements in performance. Many basic functions like Exp, Sin, etc. can operate on arrays, and use an efficient parallelized implementation under the hood. If you define your own function, then make it Listable explicitly, you will lose this performance. i.stack.imgur.com/V1GRC.png That's because you're back to evaluating the function for each element of the array separately (instead of using true vectorization). – Szabolcs May 8 '17 at 14:03

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