# Need to find matrix formulation

I have a $$B$$ matrix: $$B = B_{ij}$$

I need to find closed matrix formulation of:

$$\sum_i \sum_j \sum_m \sum_n B_{ij} B_{jm} B_{mn} B_{ni}$$

But I am so confused!

Edit by Henrik:

Originally, it was asked to express

$$\sum_i \sum_j \sum_m \sum_n B_{mi} B_{mj} B_{ni} B_{nj}$$

in terms of matrices.

actually, there is also a condition: i is not equal to j.

• This site is about Mathematica, the software. Should we move your question over to Mathematics?
– Johu
Sep 23, 2018 at 0:38
• Could you also check, if you have not made any mistakes with the indexes.
– Johu
Sep 23, 2018 at 0:44
• @DavidG.Stork I am not sure whether your edit was meaningfull. By swapping some of the indices, you actually changed the question dramatically. For example, Louis' (imho correct) answer became invalid under this modification. Sep 23, 2018 at 7:41
• @LouisB I'd suggest to undelete your post. Imho, it is a perfect answer to the original question. Sep 23, 2018 at 7:53

If this is indeed a Mathematica question, then first note that:

$$\sum _j B_{i,j} B_{j,k}\equiv (B.B)_{i,j}$$

and

$$\sum _i B_{i,i}\equiv \operatorname{Tr}[B]$$

So, the Mathematica equivalent of:

$$\sum _i \sum _j \sum _m \sum _n B_{i,j} B_{j,m} B_{m,n} B_{n,i}$$

is:

Tr[B . B . B . B]


or:

Tr[MatrixPower[B, 4]]


For the original form of the question, note that:

$$B_{i,j}\equiv B^T{}_{j,i}$$

So, the Mathematica equivalent of:

$$\sum _i \sum _j \sum _m \sum _n B_{m,i} B_{m,j} B_{n,i} B_{n,j}$$

is:

Tr[B.Transpose[B].B.Transpose[B]]


The OP added the requirement that terms where $$i=j$$ should not be included. Without explaining why, you can use the following to compute this version:

$$\sum _i \sum _j \sum _m \sum _n B_{m,i} B_{n,i} B_{m,j} B_{n,j} \ \left(1-\delta _{i,j}\right)\equiv \operatorname{Tr}\left[B^T.B.B^T.B\right]-\ \operatorname{Tr}\left[\left(B^T.B\right)^2\right]$$

For your example, $$B=\left( \begin{array}{cc} 1 & 2 \\ 3 & 4 \\ \end{array} \right)$$, we have:

B = {{1, 2}, {3, 4}};
Tr[Transpose[B] . B . Transpose[B] . B] - Tr[ (Transpose[B] . B)^2 ]


392

• thank you so much. firstly, subscripts "k" and "l" should change to "m" and "n", secondly, I checked your answer to my original question with an example, and I didn't find it correct. Sep 24, 2018 at 14:21
• @farhadpti I fixed the typo. How did you check the answer? Sep 24, 2018 at 14:24
• consider B = [1 2; 3 4]. the answer of sigma form and matrix form are 771 and 892 respectively. Sep 24, 2018 at 14:29
• @farhadpti I get 892 for both. Sep 24, 2018 at 14:33
• sorry, Yes you are right! I made a simple mistake in coding!! thank you so much. I updated the question. actually there is a condition that "i" should not be equal to "j". so some elements of B.transpose(B).B.transpose(B) must be removed. is it possible not to calculate them at all? what is the mathematical correct way for that? Sep 24, 2018 at 14:43

Carl's answer is perfect as it is. His use of MatrixPower inspired me to think a little about computational complexity: Matrix-matrix multiplication has complexity $$O(n^3)$$ (in the straight-forward implementations, not these theoretically fancy but practically irrelevant algorithms). If we can get rid of one or two of them, we can even afford a transposition to speed up the computation:

n = 5000;
B = RandomReal[{-1, 1}, {n, n}];
a = Tr[B.B.B.B]; // RepeatedTiming // First
b = Tr[MatrixPower[B, 4]]; // RepeatedTiming // First
c = With[{A = B.B}, Total[A Transpose[A], 2]]; // RepeatedTiming // First
a == b == c


5.7

4.207

2.5

True

For the original question, LouisB's answer (c below) seems to be both correct and efficient:

n = 5000;
B = RandomReal[{-1, 1}, {n, n}];
a = Tr[B.Transpose[B].B.Transpose[B]]; // RepeatedTiming // First
b = With[{u = Flatten[B.Transpose[B]]}, u.u]; // RepeatedTiming // First
c = Total[(Transpose[B].B)^2, 2]; // RepeatedTiming // First
a == b == c


5.8

2.4

2.40

True

This is a closed matrix formulation of $$B^4$$.

• It's the trace of $B^4$, no? Sep 23, 2018 at 2:00