Sum of combination of elements of a matrix

Question

Give a $N \times N$ matrix $M$ in MMA, I am interested in finding a very particular sum using its matrix elements:

$ \sum_{i \neq j \neq k \neq l}^{N} M_{ij}M_{jk}M_{kl}M_{li}$,

which is of some use in studies in Physics. The only problem is that the calculation quickly becomes very tedious, even with small values of $N$, as the number of terms $\sim N^{4}$.

I am approaching the problem in a very straightforward manner as follows. Another thing to note is that my original matrix is Hermitian (or Symmetric for the real case), so maybe one can try to use the structure, but I am not sure how to do that. We use a random GOE matrix as a minimal example.

(*Defining a matrix using seed for reproducibility*)
n = 5; 
SeedRandom[1234];
mat = RandomVariate[GaussianOrthogonalMatrixDistribution[2^n]];

Then I am using ParallelSum and explicitly defining the constraint on the sum, that none of the indices should be the same, using Boole (there are 6 in all )

ParallelSum[
 mat[[i, j]] mat[[j, k]] mat[[k, l]]  mat[[l, i]]  
Boole[i =!= j ] Boole[j =!= k ] Boole[k =!= l ] Boole[l =!= i ] Boole[j =!= l ] Boole[k =!= i], 
{i, 1, 2^n} , {j, 1, 2^n}, {k, 1, 2^n}, {l, 1, 2^n} ]
(*-767.583*)

I would like to know if there is any way to scale the calculation for large $N$. It is already very time-consuming for me to run anything more than $2^6$. For eg for $2^7$ takes around 10 minutes for me, and the time quickly escalates after that. My actual matrices are of the size $2^{16}$ for which I will use a powerful workstation, but still the calculation scales badly with $N$, and it would be useless before doing any further optimisation.

Another solution I tried was using Map to map all possible indices (pre-generated), but I quickly ran out of memory when I pre-generated all the possible indices.

I tried using ParallelTable to generate the list and do the sum, but it is slow as well.

Any help or insight on this is appreciated.

Answer 1

n = 5;
M = RandomVariate[GaussianOrthogonalMatrixDistribution[2^n]];

sum = Sum[
   Times[
    M[[i, j]] M[[j, k]] M[[k, l]] M[[l, i]] ,
    Boole[i =!= j] Boole[j =!= k] Boole[k =!= l] Boole[l =!= i], 
    Boole[j =!= l] Boole[k =!= i]
    ]
   , {i, 1, 2^n}, {j, 1, 2^n}, {k, 1, 2^n}, {l, 1, 2^n}];

My first observation was that the 4-fold sum looks strikingly similar to Tr[M.M.M.M]. If it only was not about the Boolean terms!

Well, we can get rid of four of those by just eliminating the diagonal of M:

A = M;
Do[A[[i, i]] = 0.;, {i, 1, Length[M]}];

sum2 = Sum[
   Times[
    A[[i, j]] A[[j, k]] A[[k, l]] A[[l, i]],
    Boole[j =!= l] Boole[k =!= i]
    ]
   , {i, 1, 2^n}, {j, 1, 2^n}, {k, 1, 2^n}, {l, 1, 2^n}];

Only two Boolean factors are left. We can try to ignore one out and then correct it:

sum3 = Plus[
   Sum[
    Times[
     A[[i, j]] A[[j, k]] A[[k, l]] A[[l, i]],
     Boole[k =!= i]
     ]
    , {i, 1, 2^n}, {j, 1, 2^n}, {k, 1, 2^n}, {l, 1, 2^n}],
   -Sum[
     Times[
      A[[i, j]] A[[j, k]] A[[k, l]] A[[l, i]],
      Boole[j == l] Boole[k =!= i]
      ], {i, 1, 2^n}, {j, 1, 2^n}, {k, 1, 2^n}, {l, 1, 2^n}]
   ];

You might think we have not gained much, but the summation over Boole[j == l] can be simplified to

sum3 = Plus[
   Sum[
    Times[
     A[[i, j]] A[[j, k]] A[[k, l]] A[[l, i]],
     Boole[k =!= i]
     ]
    , {i, 1, 2^n}, {j, 1, 2^n}, {k, 1, 2^n}, {l, 1, 2^n}],
   -Sum[
     Times[
      A[[i, j]] A[[j, k]] A[[k, j]] A[[j, i]],
       Boole[k =!= i]
      ], {i, 1, 2^n}, {j, 1, 2^n}, {k, 1, 2^n}]
   ];

And now we can apply the same trick (ignore the Boole expression and then correct it) to both sums. This way I get

sum4 = Plus[
   Sum[A[[i, j]] A[[j, k]] A[[k, l]] A[[l, i]], {i, 1, 2^n}, {j, 1, 2^n}, {k, 1, 2^n}, {l, 1, 2^n}],
   -Sum[A[[i, j]] A[[j, i]] A[[i, l]] A[[l, i]], {i, 1, 2^n}, {j, 1, 2^n}, {l, 1, 2^n}],
   -Sum[A[[i, j]] A[[j, k]] A[[k, j]] A[[j, i]], {i, 1, 2^n}, {j, 1, 2^n}, {k, 1, 2^n}],
   Sum[A[[i, j]] A[[j, i]] A[[i, j]] A[[j, i]], {i, 1, 2^n}, {j, 1, 2^n}]
   ];

A sharpe look at this reveals that this is equal to

sum5 = Plus[
   Tr[A.A.A.A],
   -Diagonal[A.A].Diagonal[A.A],
   -Diagonal[A.A].Total[A Transpose[A]],
   Total[A A Transpose[A] Transpose[A], 2]
   ];

Some of these computations are still a bit redundant. Factoring that in we get the following piece of code:

A = M;
Do[A[[i, i]] = 0.;, {i, 1, Length[M]}];
AT = Transpose[A];
AA = A.A;
B = AT A;
sumFast = Tr[AA.AA] - Diagonal[AA].Diagonal[AA] - Diagonal[AA].Total[B] + Total[B B, 2];

On my machine this takes about 0.01 milliseconds to compute for n = 5, while the original implementation required more than 2 seconds. Let's discuss the complexity.

Obviously, the orginal code had time complexity $O(N^4)$, where $N = 2^n$. The new code needs a matrix-matrix products AA = A.A and AA.AA. Together they cost $O(N^3)$. Total[B] and Total[B B, 2] come at a cost of $O(N^2)$. And the remaining bits have time complexity $O(N)$. So we get this speed-up only partially by complexity reduction. The bulk of the speed-up comes from using standard matrix routines that have been highly optimized.

Btw., the real obstruction for doing this for $n=16$ is that the matrix M already needs 34.4 GB of RAM for storing it. And I also stored a couple of extra matrices of the same size to avoid recomputing things. And each time you increment $n$ by one, you have to quadruple the memory requirements. So this is just not feasible for large n.

Answer 2

Assume we have a dxd matrix. Possible sum indices are 4 non-equal indices out of a Range[d]. Additionally, these may be permuted.

Therefore, to get possible indices, we take all subsets and permute them:

ind=Subsets[Range[d],{4}];
ind= Permutations /@ ind;

Now, to save memory, it is better not to create all permutations before getting the product of matrix elements, but to get the product for every permutation separately. To create a function that calculates the product, given a list of indices, we write:

(mat[[#[[1]], #[[2]]]] mat[[#[[2]], #[[3]]]] mat[[#[[3]], \
#[[4]]]] mat[[#[[4]], #[[1]]]]) &

Finally, we need to flatten the sum over the subsets and add everything together.

Here is an example for d=5. The matrix:

d = 5;
mat = Array[Subscript[m, #1, #2] &, {d, d}];
mat // MatrixForm

Flatten[(mat[[#[[1]], #[[2]] ]] mat[[#[[2]], #[[3]] ]] mat[[#[[3]], \
#[[4]]]] mat[[#[[4]], #[[1]]]]) & /@ Permutations[#] & /@ 
    Subsets[Range[d], {4}] , 1] // Total // Simplify