Tensor network Renormalization group I am working on Tensor network Renormalization group. While I am trying to implement TRG on Honeycomb lattice I got stuck on pairing two three rank tensors and represent it as D^2*D^2 matrix. Can you help me to do this one? 
I need to pair two rank three tensors. $A_{ijm} B_{mkl}$. Then I need to find a tensor C and D in the form $C_{lin} D_{jkn}$. How can I do it?
 A: Mathematica implementation. Define the tensor trace function:
tTr[Ts_,s_]:=Activate@TensorContract[Inactive[TensorProduct]@@Ts,s];

The tensor trace can be simply implemented by specifying the list of tensors and the pairs of legs should be contracted. The legs are indexed following the orders of the tensors. For example, the following code takes the 3-leg tensors $A$ and $B$ and contract the legs $(3,4)$ to form the 4-leg tensor $X$.
A = RandomReal[1, {2, 2, 2}];
B = RandomReal[1, {2, 2, 2}];
X = tTr[{A, B}, {{3, 4}}];


To use the tTr function correctly, one needs to understand how to translate the Einstein notation $A_{ijm}\;B_{mkl}$ to the tensor-network language. The rule is quite simple, just attach to every index in the Einstein notation by the leg index 1,2,3... from left to right. For example, one can rewrite $A_{ijm}\;B_{mkl}$ as $A_{(i1)(j2)(m3)}\;\;B_{(m4)(k5)(l6)}\;\;$, then we know $(m3)$ and $(m4)$ are legs to be contracted, so we need to tell tTr to pair up legs $(3,4)$. After the tensor trace, we obtain a 4-leg tensor $X_{ijkl}=A_{ijm}\;B_{mkl}$. 
Then to split the tensor $X_{ijkl}=C_{lin}\;D_{jkn}$, we need to group the legs $li$ and $jk$ first. To help us find the leg indices, we write $X$ as $X_{(i1)(j2)(k3)(l4)}\;\;$, then we know that the legs $li$ are actually legs $(4,1)$ and the legs $jk$ are actually legs $(2,3)$, so we can group the legs and perform the SVD by
{U, S, V} = SingularValueDecomposition@Flatten[X, {{4, 1}, {2, 3}}];

So legs $li$ will go to $U$ and the legs $jk$ will go to $V$ now. We need to absorb the singular value into the left and right unitary matrices, this can be done by
V = V.Sqrt[S];
U = U.Sqrt[S];

For tensor RG, you may want to restrict the dimension of the internal legs at this step. For example, you may do
V = V[[;; d]];
U = U[[;; d]];

where $d$ should be the bond dimension that you specified previously. But we will not do this truncation for now. Finally, we reshape the matrices $U$ and $V$ into 3-leg tensors $C$ and $D$ by
C = ArrayReshape[U, {2, 2, 4}];
D = ArrayReshape[V, {2, 2, 4}];

So we have obtained the new tensors $C$ and $D$. To verify that this is the correct answer, we can contract the tensors $C$ and $D$. By the rewriting $C_{(l1)(i2)(n3)}\;\;D_{(j4)(k5)(n6)}\;\;$, we know this time we need to contract the legs $(3,6)$, so we tell that to tTr
Y = tTr[{C, D}, {{3, 6}}];

Now we have obtained a 4-leg tensor $Y_{lijk}=C_{lin}\;D_{jkn}$. We should have $Y_{lijk}=X_{ijkl}$. But we can not directly ask Mathematica to check if Y == X, because the two tensors now have different leg orderings. To rearrange the legs, we can use the generalized transpose. But we need to tell Mathematica the leg mapping rules. To find the mapping rules from $Y$ to $X$, we can write $Y_{(l1)(i2)(j3)(k4)}\;\;\to X_{(i1)(j2)(k3)(l4)}\;\;$, so we have established the mapping for each leg as $l:(l1)\to(l4)$, $i:(i2)\to(i1)$, $j:(j3)\to(j2)$, $k:(k4)\to(k3)$. We do not need to keep the names $i,j,k,l$, we can just write the mapping as $(1,2,3,4)\to(4,1,2,3)$. We don't even need to tell Mathematica the left-hand-side of the mapping, because by definition, the left-hand-side must always be ordered like 1,2,3..., so we only need to tell Mathematica to transpose the tensor $Y$ by $(4,1,2,3)$, then compare the result with the tensor $X$.
Transpose[Y, {4, 1, 2, 3}] == X

The result is True despite that $A$ and $B$ are randomly generated. So we know that the whole procedure works.
In conclusion, we can realize tensor operations in Mathematica very easily. The following four functions basically covers all our needs:


*

*tTr: perform tensor trace (tensor network contraction),

*Transpose: rearrange tensor legs,

*Flatten: group together tensor legs,

*ArrayReshape: split appart tensor legs (actually it can be used to reshape the tensor to any shape you want).


You can build a tensor RG package on top of these four functions with the SVD solver provided by SingularValueDecomposition.
