# How can I speed up this code with 4 variables sum

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}]


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


• 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 Aug 24, 2015 at 11:10
• 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]}, #] &)] Aug 24, 2015 at 11:15
• @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. Aug 24, 2015 at 13:04
• @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. Aug 24, 2015 at 13:21
• 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. Aug 24, 2015 at 13:24

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


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.