1
$\begingroup$

I have to make a sum over 4 variables. My code is very very slow. I want to know how to speed up this code. This problem is related to but different from one previous problem. Any help or suggestion will be highly appreciated! The code is shown below:

data = Table[   Exp[-((i + j - 20.)/5)^2] Exp[-((i - j)/5)^2], {i, 20}, {j, 20}];
data = Chop[data, 0.00001];
data = data/Sqrt[Sum[(data[[i, j]])^2, {i, 1, 20}, {j, 1, 20}]];
ListDensityPlot[data, InterpolationOrder -> 0, Mesh -> All,  PlotRange -> All, ColorFunction -> (Blend[{Hue[2/3], Hue[0]}, #] &)]
c = 3*10^8; 
Δ = 0.5; 
λ0 = 1500;
CC1[i_, j_, k_, l_, t_] := (data[[i, l]] data[[j, 
  k]] Cos[π*(c/(λ0 - 10 + i*Δ - 
      0.5 Δ) + 
     c/(λ0 - 10 + l*Δ - 
      0.5 Δ)) t] Cos[π*(c/(λ0 - 10 + 
      j*Δ - 0.5 Δ) + 
     c/(λ0 - 10 + k*Δ - 
      0.5 Δ)) t] + 
data[[i, k]] data[[j, 
  l]] Cos[π*(c/(λ0 - 10 + i*Δ - 
      0.5 Δ) + 
     c/(λ0 - 10 + k*Δ - 
      0.5 Δ)) t] Cos[π*(c/(λ0 - 10 + 
      j*Δ - 0.5 Δ) + 
     c/(λ0 - 10 + l*Δ - 
      0.5 Δ)) t] - 
data[[i, j]] data[[k, 
  l]] Sin[π*(c/(λ0 - 10 + i*Δ - 
      0.5 Δ) + 
     c/(λ0 - 10 + j*Δ - 
      0.5 Δ)) t] Sin[π*(c/(λ0 - 10 + 
      k*Δ - 0.5 Δ) + 
     c/(λ0 - 10 + l*Δ - 
      0.5 Δ)) t])^2;
CC2[t_] := \!\(\*UnderoverscriptBox[\(∑\), \(i = 1\), \(20\)]\(\*UnderoverscriptBox[\(∑\), \(j = 1\), \(20\)]\(\*UnderoverscriptBox[\(∑\), \(k = 1\), \(20\)]\(\*UnderoverscriptBox[\(∑\), \(l = 1\), \(20\)]CC1[i, j, k, l,t]\)\)\)\);
ListPlot[Table[{i, CC2[i*0.001]}, {i, -10, 10, 1}], Joined -> True, Axes -> None, PlotRange -> All, Frame -> True, ImageSize -> {400, 250}]
ListPlot[Table[{i, CC2[i*0.001]}, {i, -10, 10, 0.001}],  Joined -> True, Axes -> None, PlotRange -> All, Frame -> True,  ImageSize -> {400, 250}]
$\endgroup$
9
$\begingroup$

The problem is that it reevaluates the sum every single time you call it, recomputing every 20^4 term again and again. You just need to compile the function CC2 so that it performs the summation only once. Using the code you have, it takes my machine about 6 seconds to compute a single data point:

CC2[0.003] // AbsoluteTiming

(*  {6.069311, 1.49893}  *)

But if I compile it first,

CC2comp = Compile[{{t, _Real}},
  Evaluate[
     Sum[CC1[i, j, k, l, t], {i, 20}, {j, 20}, {k, 20}, {l, 20}]
     ]];

I now get the answer in milliseconds

CC2comp[0.003] // AbsoluteTiming

{0.002888, 1.49893}

Now you can make the plot you were going for in about a minute.

ListPlot[Table[{i, CC2comp[i*0.001]}, {i, -10, 10, 0.001}],  
  Joined -> True, Axes -> None, PlotRange -> All, Frame -> True,  
  ImageSize -> {400, 250}] // AbsoluteTiming

enter image description here

$\endgroup$
  • $\begingroup$ Thank you so much for this good suggestion! My final goal to use this code for a 100*100 data. So I hope the speed can be faster. The compile time for 20*20 data is 13 seconds; but the compile time for 100*100 is more than one hour. I hope to simplify the equation use the command of "Total", as it was used in last example [link] (mathematica.stackexchange.com/questions/48640/…). @Simon Woods $\endgroup$ – user14634 Aug 24 '15 at 11:10
  • $\begingroup$ The compile time for 100*100 data data = Table[ Exp[-((i + j - 100.)/10)^2] Exp[-((i - j)/10)^2], {i, 100}, {j, 100}]; data = Chop[data, 0.00001]; data = data/Sqrt[Sum[(data[[i, j]])^2, {i, 1, 100}, {j, 1, 100}]]; ListDensityPlot[data, InterpolationOrder -> 0, Mesh -> All, PlotRange -> All, ColorFunction -> (Blend[{Hue[2/3], Hue[0]}, #] &)] $\endgroup$ – user14634 Aug 24 '15 at 11:15
  • $\begingroup$ @user14634 Well it seems to me that Jason has given you an excellent example of how to solve your speed problem as presented in your own question. Using Total as you suggest may also be a good idea, which I encourage you to pursue. However, I'd like to caution you: changing the nature of the question after it has been answered is frowned upon in this forum. $\endgroup$ – MarcoB Aug 24 '15 at 13:04
  • $\begingroup$ @user14634, I do agree with MarcoB that if you have a different problem, it might be better as a new question. Compiling a function that is a sum of 20^4 = 160,000 terms is not a big deal on a modern computer, but a function with 100^4, or a hundred million terms, is a different story. I can't say whether you would gain something by writing your function in terms of matrix multiplication, or using Total. $\endgroup$ – Jason B. Aug 24 '15 at 13:21
  • 1
    $\begingroup$ But you could gain some time by taking advantage of the sparseness of your data matrix. Before applying Chop there are 10,000 nonzero elements in data, but after applying Chop there are only 1,789. By my reckoning, of the 100 million potential terms in the summation, only a little less than 4 million aren't equal to zero for all values of t. Try doing the sum only over this smaller list of elements. $\endgroup$ – Jason B. Aug 24 '15 at 13:24
0
$\begingroup$

Okay, so I tried to write the function like you were saying, taking advantage of the Total built-in function, and there is some improvement.

cc2Comp = Compile[{{t, _Real}},
Evaluate@
Module[{ndim, data, c, Δ, λ0, costable, 
  sintable, table1, table2, table3, table4},
  ndim = 20;
  data = 
  Array[Exp[-((#1 + #2 - 20.)/5)^2] Exp[-((#1 - #2)/5)^2] &
    , {ndim, ndim}];
  data = Chop[data, 0.00001];
  data = data/Norm[data];
  data = SparseArray[data];
  c = 3*10^8;
  Δ = 0.5;
  λ0 = 1500;

  costable = 
  Table[Cos[π t (c/(-10 - 0.5` Δ + 
    j Δ + λ0) + c/(-10 - 0.5` Δ + k Δ + λ0))], {j, ndim}, {k,ndim}];
  sintable = 
  Table[Sin[π t (c/(-10 - 0.5` Δ + 
    i Δ + λ0) + c/(-10 - 0.5` Δ + j Δ + λ0))] , {i, ndim}, {j, ndim}];
  table1 = 
  Transpose[
    Outer[Times, costable data, costable data], {1, 4, 2, 3}];
  table2 = 
  Transpose[
    Outer[Times, costable data, costable data], {1, 3, 2, 4}];
  table3 = Outer[Times, sintable data, sintable data];


  (* This next part can be modified, changing the amount of time needed to 
    compile the code versus the time needed to execute the code  *)
  Total[
  table1^2 + 2 table1 table2 + table2^2 - 2 table1 table3 - 2 table2 table3 + table3^2, 4]
  ]
]; // AbsoluteTiming


cc2Comp[0.003] // AbsoluteTiming
ListPlot[Table[{i, cc2Comp[i*0.001]}, {i, -10, 10, 0.001}], 
  Joined -> True, Axes -> None, PlotRange -> All, Frame -> True, 
  ImageSize -> {400, 250}] // AbsoluteTiming

enter image description here

As you can see, on my machine, this took about 10 seconds to compile (similar to the method in my first answer), about .026 seconds to compute a single data point (ten times longer than above), and 16 seconds to make the plot (one fourth of the time above).

If I change the last line of the Module to read

table4 = SparseArray[table1 + table2 - table3];

Total[(table4)^2, 4]

Then the result is it takes less time to compile but more time to make the plot, and a little more time altogether.

Just straightforwardly applying this method to larger arrays, as in your other question, is problematic. I don't have the RAM to work do this method when the dimension is raised to 100 (since we are talking about 100 million matrix elements in the summation), and it crashes the kernel.

enter image description here

$\endgroup$

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.