0
$\begingroup$

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]}, #] &)]

enter image description here

$\endgroup$
  • $\begingroup$ Look into Compile. $\endgroup$ – Ajasja Jul 15 '14 at 7:34
4
$\begingroup$

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

Mathematica graphics

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

Mathematica graphics

$\endgroup$
  • $\begingroup$ @xzczd Thank you so much! Your answer is really helpful. $\endgroup$ – user14634 Jul 15 '14 at 9:37
  • $\begingroup$ @blochwave, thank you for your kind encouragement. $\endgroup$ – user14634 Jul 15 '14 at 9:43
  • $\begingroup$ @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. $\endgroup$ – user14634 Jul 15 '14 at 11:10
  • $\begingroup$ @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. $\endgroup$ – xzczd Jul 15 '14 at 11:48
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Many thanks to Teake Nutma and Ajasja. The Compile function works well. $\endgroup$ – user14634 Jul 15 '14 at 9:21
  • 1
    $\begingroup$ You're welcome, but @xzczd's answer is a whole lot better :). $\endgroup$ – Teake Nutma Jul 15 '14 at 9:23
  • $\begingroup$ @ 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? $\endgroup$ – user14634 Jul 15 '14 at 9:34
  • $\begingroup$ @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. $\endgroup$ – Teake Nutma Jul 15 '14 at 9:51
  • $\begingroup$ I am using Mathematica 8.0; The Compile function can be used, but the warnings appear. $\endgroup$ – user14634 Jul 15 '14 at 9:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.