How can I speed up this code with multiple sum？

Question

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]

Answer 1

Your sum can be computed more efficiently as follows:

data = Developer`ToPackedArray[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.