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.
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.} *)
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).
$\endgroup$
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[3000] // 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.
tau2have a value? What aboutdt? $\endgroup$factor[x,y] == factor[y,x], so you need only iterate fory > x. $\endgroup$