How can I speed up this code？

There are 4 variables in this multiple sum, and only 2 are used. This question has the similar scenario as the previous one How can I speed up this code with multiple sum？, but the calculation process is different. This code takes about 1 min 38 sec to obtain the results, but I hope the time can be shortened to several seconds. Any help or suggestion will be highly appreciated!

the code is bellow :

data = Table[Exp[-((i + j - 50.)/5)^2] Exp[-((i - j)/5)^2], {i, 50}, {j, 50}];
data = Chop[data, 0.00001];
data = data/Sqrt[Sum[(data[[i, j]])^2, {i, 1, 50}, {j, 1, 50}]];
ListDensityPlot[data, InterpolationOrder -> 0, Mesh -> All,  PlotRange -> All, ColorFunction -> (Blend[{Hue[2/3], Hue[0]}, #] &)]

S1[i_, j_, k_, m_, t_] := (data[[i, k]] data[[j, m]])^2 + (data[[j, k]] data[[i,
m]])^2 - 2 data[[i, k]] data[[j, m]] data[[j, k]] data[[i,
m]] Cos[(2 \[Pi]*(3*10^8)/(1584 - 5 + j*0.1 - 0.05) -
2 \[Pi]*(3*10^8)/(1584 - 5 + i*0.1 - 0.05)) t];
S2[i_, j_, t_] := 1/4*\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$k = 1$$, $$50$$]$$\*UnderoverscriptBox[\(\[Sum]$$, $$m = 1$$, $$50$$]S1[i, j, k, m,
t]\)\);
data2 = Table[S2[i, j, 4*0.001], {i, 1, 50}, {j, 1, 50}];
ListDensityPlot[data2, InterpolationOrder -> 0, Mesh -> All, PlotRange -> All, ColorFunction -> (Blend[{Hue[2/3], Hue[0]}, #] &)]


• Look into Compile. Commented Jul 15, 2014 at 7:34

You should have read Simon Woods's answer more carefully. The difference between this question and the previous one is just you don't need to sum over $j$ and $k$ now, with almost same analysis, you can figure out that data2 is equal to

$$2 {\left((d.d^T)^2\right)}_{ij}(1 - f_{ij})/4$$

So here is the code:

a = AbsoluteTiming;
data0 = Table[Exp[-((i + j - 50.)/5)^2] Exp[-((i - j)/5)^2], {i, 50}, {j, 50}];
data = Partition[Normalize@Flatten@Chop[data0, 0.00001], 50]; // a
ArrayPlot[data, Mesh -> All, MeshStyle -> Black,
ColorFunction -> (Blend[{Hue[2/3], Hue[0]}, #] &)] // a


Notice my modification for data. (As mentioned by Simon in the previous post, you will generally get much better performance by acting on a whole array at once instead of on each element individually.) Also, I chose ArrayPlot instead of ListDensityPlot because it's faster and has the same effect as your ListDensityPlot with the option InterpolatingOrder -> 0.

The rest part is nearly the same as Simon's:

f[i_, j_, t_] = Cos[(2 Pi (3 10^8)/(1584 - 5 + 0.1 j - 0.05) -
2 Pi (3 10^8)/(1584 - 5 + 0.1 i - 0.05)) t];
data2 = 2 (data.Transpose[data])^2 (1 - Array[f[##, 4 0.001] &, {50, 50}]) / 4; // a
ArrayPlot[data2, Mesh -> All, MeshStyle -> Black,
ColorFunction -> (Blend[{Hue[2/3], Hue[0]}, #1] &)] // a


• @xzczd Thank you so much! Your answer is really helpful. Commented Jul 15, 2014 at 9:37
• @blochwave, thank you for your kind encouragement. Commented Jul 15, 2014 at 9:43
• @xzczd How to prove the first term? I mean sum\limits_ {k, m} (d_ {ik} d_ {jm})^2 = (d_. d^T)^2 _ {ij} ? Can we prove: sum\limits_ {k, m} (d_ {ik} d_ {jm})^2 =sum\limits_ {k, k} (d_ {ik} d_ {jk})^2 ? It seems difficult. Commented Jul 15, 2014 at 11:10
• @user14634 This is exactly what Simon proves in the last 4 steps of his explanation. By the way, the justification can be even easier with tensor notation. Commented Jul 15, 2014 at 11:48

Like Ajasja said in a comment, Compile is your best friend here (short of completely rewriting S1 and S2 like in your other question). Simply compiling everything together:

With[{d = data},
S2 = Compile[
{{i, _Integer}, {j, _Integer}, {t, _Real}},
1/4*Sum[
(d[[i, k]] d[[j, m]])^2
+ (d[[j, k]] d[[i, m]])^2
- 2 d[[i, k]] d[[j, m]] d[[j, k]] d[[i, m]] Cos[
(2 \[Pi]*(3*10^8)/(1584 - 5 + j*0.1 - 0.05) -
2 \[Pi]*(3*10^8)/(1584 - 5 + i*0.1 - 0.05)) t
],
{k, 1, 50}, {m, 1, 50}
],
CompilationTarget -> "C", RuntimeOptions -> "Speed"
]
];


gives a considerably speed increase:

(data2 = Table[S2[i, j, 4*0.001], {i, 1, 50}, {j, 1, 50}];) // Timing


{0.915097, Null}

Note that inserting data into the compiled function with With is key to getting a big speed increase. My guess is that the compiled function then doesn't have to go to the main evaluation loop to acces data each time, which would severely slow performance.

• Many thanks to Teake Nutma and Ajasja. The Compile function works well. Commented Jul 15, 2014 at 9:21
• You're welcome, but @xzczd's answer is a whole lot better :). Commented Jul 15, 2014 at 9:23
• @ Teake Nutma, I met the error of "A C compiler cannot be found on your system. Please consult the documentation to learn how to set up suitable compilers." Should I install some compilers first? Commented Jul 15, 2014 at 9:34
• @user14634 That shouldn't be necessary. What Mathematica version are you on? The options to Compile I gave (CompilationTarget, RuntimeOptions) only work for Mathematica 8 or newer; if you have 7 or older you can simply omit them. Note that doing so on Mathematica 7 will slow down performance to around 10 seconds. Commented Jul 15, 2014 at 9:51
• I am using Mathematica 8.0; The Compile function can be used, but the warnings appear. Commented Jul 15, 2014 at 9:57