# How to simplify the equation and speed up the codes for four 100*100 matrices

I have to make a sum over 4 variables. Each variable corresponds to a 100*100 matrix. I want to know how to speed up this code. This problem is related to the previous problem 1 and 2. But now I need higher speed than the case in the previous two problems. I hope to simplify the equation using the command of “Total”， which works well in 1. But simplifying the equation correctly is not easy for me. Any help or suggestion will be highly appreciated! The code is shown below：

data = Table[ 10 Exp[-((i + j - 100.)/10)^2] Exp[-((i - j)/10)^2] -
2 Exp[-((i + j - 100.)/10)^2] Exp[-((i - j - 30)/8)^2] -
2 Exp[-((i + j - 100.)/10)^2] Exp[-((i - j + 30)/8)^2] +
Exp[-((i + j - 100.)/10)^2] Exp[-((i - j - 60)/6)^2] +
Exp[-((i + j - 100.)/10)^2] Exp[-((i - j + 60)/6)^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]}, #] &)]

c = 3*10^8; Δ = 0.1; λ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[\(\[Sum]$$, $$i = 1$$, $$100$$]$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$100$$]$$\*UnderoverscriptBox[\(\[Sum]$$, $$k = 1$$, $$100$$]$$\*UnderoverscriptBox[\(\[Sum]$$, $$l = 1$$, $$100$$]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}]

• This is the same exact code you posted earlier, but in that case the sum was over a smaller 20x20 grid. Nevertheless, JasonB showed you how to speed up this code 3000x by compiling the function involved in the sum. have you tried applying that approach to this form of the problem on your own? If you did, where did you get stuck exactly? Commented Aug 25, 2015 at 11:51
• @MarcoB - It's because in this problem he's doing a four-fold summation where each index goes up to 100, for a total of 10^8 terms. Trying to create a simple compiled function with that many terms in the sum would require a lot of RAM, so it won't work in this case. Commented Aug 25, 2015 at 14:34
• Many thanks to Jason B! Yes, if I try 100*100 data, the speed is very slow. Commented Aug 26, 2015 at 4:18

Simply taking the answer from the 20x20 case and applying it directly to this case will isn't really applicable (it crashes my machine, requiring a hard restart). But I think we can improve upon your code and get the answer (about 700 times faster on my machine).

First, I can say that to get the value of your function CC2 for a single time point will take about an hour on my machine. How did I figure this out? I sure didn't waste an hour trying it, I did a little detective work. I tried it on a smaller sample first, and got the timing and extrapolated. Like this:

Monitor[
temp = Sum[
CC1[i, j, k, l, t3], {i, 3}, {j, 100}, {k, 100}, {l, 100}]; //AbsoluteTiming
, {i, j, k, l}]


Note the function Monitor, it is your friend when trying to watch these kinds of code. So doing the above sum over just 3 values of i took my machine 112.83 seconds, so to do all 100 would take 63 minutes.

Therefore, just trying the code you have for the second plot would take 8.7 days. That's why you need to monitor your code execution, and figure out how long it will take before you try generating a large dataset, rather than just writing the code and twiddling your thumbs wondering when your plot is going to pop up on the screen.

The following function should be able to get a single data point in a much shorter time:

ndim = 100;
data = Table[
10 Exp[-((i + j - 100.)/10)^2] Exp[-((i - j)/10)^2] -
2 Exp[-((i + j - 100.)/10)^2] Exp[-((i - j - 30)/8)^2] -
2 Exp[-((i + j - 100.)/10)^2] Exp[-((i - j + 30)/8)^2] +
Exp[-((i + j - 100.)/10)^2] Exp[-((i - j - 60)/6)^2] +
Exp[-((i + j - 100.)/10)^2] Exp[-((i - j + 60)/6)^2], {i, ndim}, {j, ndim}];
data = Chop[data, 0.00001];
data = data/Sqrt[Sum[(data[[i, j]])^2, {i, ndim}, {j, ndim}]];
data = SparseArray[data];
c = 3*10^8;
Δ = 0.1;
λ0 = 1500;

cc2func[t_] :=
Module[{costable, sintable, table1, table2, table3, table4},
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];
table4 =
table1^2 + 2 table1 table2 + table2^2 - 2 table1 table3 -
2 table2 table3 + table3^2;
Total[table4, 4]
];


When I execute cc2func[0.0008] it takes about 5.5 seconds to return an answer, so it is about 700 times faster.

Now you have to ask what you want to have in the end? Do you want a single plot? If so, we can run this code for a few hours and get it. I don't see how to make it much faster than that.

I think you are running your Table at the end with too fine a step size. First, you should save the data as a list, then plot that list. Otherwise you won't be able to plot it again without recalculating the data.

Monitor[
table2plot = Table[{t, cc2func[t*.001]}, {t, -.5, .5, .01}];
table2plotB = Table[{t, cc2func[t*.001]}, {t, -.05, .05, .005}];
, t];
Show[ListLinePlot[table2plot],ListLinePlot[table2plotB,PlotStyle->Red]]


gives this result

So there really isn't any need to have the stepsize smaller than 0.01. Also, I see your function is symmetric with respect to t=0 so why calculate both positive and negative timepoints? If each call to cc2func takes 5.5 seconds, you should be able to make your plot in about an hour and a half.

• Thank you so much for the new answer. I have tried to execute cc2func[0.0008] using your program, but the speed is very very slow. The result doesnot appear in 10 minutes. What is reason for this? Commented Aug 26, 2015 at 4:23
• @user14634, perhaps you don't have enough RAM on your machine? When I copy the code from the first block and execute it, then try cc2func[0.0008], 5.74 seconds later it returns the value 1.8319. Within the function, there are several tables created before the final call to Total; try putting PrintTemporary["costable created"], etc, in between each line so that you can see where you are getting hung up. There are many tools like Monitor, Print, PrintTemporary that I find useful when evaluating my code. Commented Aug 26, 2015 at 6:53
• Thanks a lot for the good comments! Yes, I donot have enough RAM on my computer. After I close all the other programs and only run Mathematica, I get In[19]:= cc2func[0.0008] // AbsoluteTiming Out[19]:= {1063.9924969, 1.8319}. It takes 18 min on my computer for one point. Too long time for me. Maybe the best strategy is to simplify the equation first. One clue is that the four variables are separable, e.g., data[[i, j]] data[[k, l]]Sin[i,j]Sin[k,l]=(data[[i, j]]Sin[i,j])(data[[k, l]]Sin[k,l]). Commented Aug 26, 2015 at 11:43
• I have tried to simply the equation in the following way. The speed is very fast, only 0.93 seconds for one point. But the result is not correct. Commented Aug 26, 2015 at 11:57
• n = 100; ff1[i_, j_] = Cos[[Pi]*(c/([Lambda]0 - 10 + i*[CapitalDelta] - 0.5 [CapitalDelta]) + c/([Lambda]0 - 10 + j*[CapitalDelta] - 0.5 [CapitalDelta])) t]; ff2[i_, j_] = Sin[[Pi]*(c/([Lambda]0 - 10 + i*[CapitalDelta] - 0.5 [CapitalDelta]) + c/([Lambda]0 - 10 + j*[CapitalDelta] - 0.5 [CapitalDelta])) t]; Commented Aug 26, 2015 at 11:57