# 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? May 8, 2017 at 13:22
• Sorry, my fault. Edited my post. The values are of no importance. It is more a programming style question. May 8, 2017 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. May 8, 2017 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. May 8, 2017 at 13:48
• Good to know, thank you! May 8, 2017 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). May 8, 2017 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}    *)


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 // AbsoluteTiming


I get about a 3x speed up this way.

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.