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.
 A: 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

A: 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

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