Need to find matrix formulation

Question

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.

Answer 1

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]]

Addendum

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

Answer 2

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

Answer 3

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