How can I speed up this code with multiple sum？

There are 4 variables in this multiple sum, therefore it may take a long time. I have run this program for 12 hours, but no result untill now. I want to know how to speed up this code. Any help or suggestion will be highly appreciated!

the code is bellow:

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


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[t_] := 1/4*\!$$\*UnderoverscriptBox[\(\[Sum]$$, $$i = 1$$, $$100$$]$$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$100$$]$$\*UnderoverscriptBox[\(\[Sum]$$, $$k = 1$$, $$100$$]$$\*UnderoverscriptBox[\(\[Sum]$$, $$m = 1$$, $$100$$]S1[i, j, k, m, t]\)\)\)\);
ListPlot[Table[S2[i], {i, -0.01, 0.01, 0.0001}], Joined -> True,
PlotRange -> All, Frame -> True]

-
I guess you could speed up your code by using SparseArray and ParallelSum. –  Öskå May 27 at 12:38
I'd suggest precomputing S1 and S2, instead of putting them into functions, and I'd do it by breaking it into pieces. For example, data[[i, k]] data[[j, m]] is identical to data[[j, k]] data[[i, m]] up to a permutation of indices, so you only have to calculate it once. I'll try and get back to this later today, but that should give you a direction, at least. –  rcollyer May 27 at 14:25
Some general tips based on a quick glance at your code: 1. Try to use vectorization where you can. Instead of indexing over arrays and processing elements at a time, operate on the whole array. 2. You can use data /= Norm@Flatten[data] or data /= Sqrt@Total[data^2, 2] (this is not an answer, just some tips) –  Szabolcs May 27 at 16:13

Your sum can be computed more efficiently as follows:

data = DeveloperToPackedArray[data];
n = Length@data;

f[i_, j_] = Cos[(2 π*(3*10^8)/(1584 - 5 + j*0.1 - 0.05) -
2 π*(3*10^8)/(1584 - 5 + i*0.1 - 0.05)) t];

S3[t_] = (Total[data^2, -1]^2 - Total[(data.Transpose[data])^2 Array[f, {n, n}], -1])/2;

ListPlot[Table[S3[i], {i, -0.01, 0.01, 0.0001}], Joined -> True,
PlotRange -> All, Frame -> True]


The total time is about 1 second. Note that the term Total[data^2, -1]^2 is actually equal to 1 for your data, because of the normalisation.

How it works

The first term of the sum looks like:

$$\sum _{i,j,k,m} {\left(d_{ik} d_{jm}\right)}^2$$

Since the square of a product is the product of the squares, this can be written as

$$\sum _{i,j,k,m} a_{ik} a_{jm}$$

where $a=d^2$. This takes the product of every possible pair of elements in $a$ and sums them, which is equivalent to summing all the elements of $a$ and then squaring the result. So this gives us

$${\left(\sum _{i,k} a_{ik}\right)}^2$$

The second term is exactly the same (it's just an interchange of indices). So the first two terms of the sum are given by 2 Total[data^2, -1]^2.

The third term looks like (ignoring the factor of -2):

$$\sum _{i,j,k,m} f_{ij} \left(d_{ik} d_{jk}\right) \left(d_{im} d_{jm}\right)$$

where the cosine term has been expressed as $f_{ij}$. I have also re-ordered the terms to collect the $k$ and $m$ indices together. Again, we have a sum of products which can be expressed as the square of the sum:

$$\sum _{i,j} f_{ij} {\left(\sum _k d_{ik} d_{jk}\right)}^2$$

If we express the term $d_{jk}$ as $(d^T)_{kj}$ we get

$$\sum _{i,j} f_{ij} {\left(\sum _k d_{ik} (d^T)_{kj}\right)}^2$$

where the sum over $k$ is now just an inner product, so we have

$$\sum _{i,j} f_{ij} {\left((d.d^T)^2\right)}_{ij}$$

So the third term of your sum is given by -2 Total[(data.Transpose[data])^2 Array[f,{n,n}], -1]

Expressing the sum as I have done makes us of the Listable attribute of Power and Times, and the very fast implementation of Dot for numerical arrays. As Szabolcs commented, you will generally get much better performance by acting on a whole array at once instead of on each element individually.

-
Thank you for this answer! It works well. But I don't understand why the 4-variable sum can be replaced by Outer[Dot, data, data, 1]^2 Array[f, {n, n}]. Why can't the 4-variable sum cannot be replaced by data^4 Array[f, {n, n}]`? –  user14634 May 28 at 7:33
@user14634, I added some explanation. –  Simon Woods May 28 at 11:48
Now I can understand this detailed explanation. Thank you so much for your kind help. –  user14634 May 29 at 2:59