# Stack operation of tensor networks

Tianning Zhang<EMAIL_ADDRESS>Science, Mathematics and
Technology, Singapore University of Technology and Design, 8 Somapah Road,
487372 Singapore  Tianqi Chen  School of Physical and Mathematical Sciences,
Nanyang Technological University, 639798 Singapore  Erping Li  College of
Information Science and Electronic Engineering, Zhejiang University, Hangzhou
310027, China  Bo Yang  School of Physical and Mathematical Sciences, Nanyang
Technological University, 639798 Singapore  L. K. Ang  Science, Mathematics
and Technology, Singapore University of Technology and Design, 8 Somapah Road,
487372 Singapore

###### Abstract

The tensor network, as a facterization of tensors, aims at performing the
operations that are common for normal tensors, such as addition, contraction
and stacking. However, due to its non-unique network structure, only the
tensor network contraction is so far well defined. In this paper, we propose a
mathematically rigorous definition for the tensor network stack approach, that
compress a large amount of tensor networks into a single one without changing
their structures and configurations. We illustrate the main ideas with the
matrix product states based machine learning as an example. Our results are
compared with the _for_ loop and the _efficient coding_ method on both CPU and
GPU.

##  I Introduction

Tensors are multi-dimensional arrays describing multi-linear relationships in
a vector space. The linear nature of the tensors implies that two tensors may
be added, multiplied, or stacked. The stack operation of tensors place all the
tensors with the same shape row by row, resulting in a higher-dimensional
tensor with an extra stack dimension. For example, stacking scalars will lead
to a vector and stacking vectors gives a matrix. The tensor operation is
designed for the fundamental calculation in modern chips such as graphics
processing unit (GPU) and tensor processing unit (TPU) [1]. Thus, computations
with a combination of tensor operations would be more efficient as compared to
the conventional ways. For example, there are two ways to calculate the dot
product between a hundred vectors $v^{i}_{n\times 1}(i=1,2,\cdots,100)$ and a
core vector $c_{n\times 1}$: either by applying the dot operation $100$ times
via a _for_ loop or align all the vectors into a big matrix $M_{100\times n}$
and apply a matrix-vector dot product. The later method is much faster in GPU
or TPU due to their intrinsic design.

The tensor network [2, 3, 4, 5] is a decomposition of a very large tensor into
a network structure of smaller tensors, which has seen recent applications in
machine learning[6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. As it is merely the
decomposition of the tensors, the tensor networks were expected to be able to
maintain the same operations as mentioned previously. However, due to its
arbitrary decomposition into a network structure, only the contraction
operation is well-defined. Other operations such as ‘stack’ and the batched
tensor networks have not been studied thoroughly so far. In fact, the stack
operation is quite an important part of the modern machine learning tasks, by
treating hundreds of input tensors as a single higher-dimensional batched
tensor at training and inference steps. Recent tensor network machine learning
approaches includes the matrix product states (MPS) [6, 11], the string bond
states (SBS) [12], projected entangled pair states (PEPS) [16, 17], the tree
tensor network (TTN) [9, 13, 18] and others[19, 20] that either use a naïve
iteration loop, or skip the batched tensor network via a programming trick (we
note it as the efficient coding). Several studies have explored the concept of
stack operation of tensor network to a certain extent. For instance, Ref. [21]
puts forward the ‘A-tensor’ to represent the new local unit of MPS after
adding the extra local state space. It represents the local tensor of the
symmetric tensor network, which is regarded as the stack of different
‘multiplets’. This approach is developed for the scenario that the system has
certain symmetries so that the computational cost can be well controlled. It
did not explicitly reveal the concept of stacking hundreds of tensor networks
into one network. A more recent work [22] has put forward the ‘add’ operation
of MPS for solving the many-body Schödinger’s equation. It represents the
coefficients of the many-body wave function as tensor networks and ‘add’ them
into one compact form. The ‘add’ operation is similar to our definition of the
stacked tensor network with an additional reshape. Therefore, it is useful to
discuss about what is the stack of tensor network and how it performs in
machine learning tasks. It is expected that the batched tensor network method
may be an alternative option for a faster machine learning program when
efficient coding is unavailable. It can also provide the theoretical
background for future chips designs based on tensor network contractions.

In this paper, a regular tensor network stack operation method is put forward
to compress two or more tensor networks with the same structure into one
tensor network with one more stack bond dimension. The remaining part of the
stacked tensor network is the same as that of the original tensors. We
illustrate the main procedures with the realization of the stack operation on
a one-dimensional tensor network system, i.e., the MPS. However, our result
can be extended to arbitrarily shaped tensor network such as the PEPS. The
rest of the paper is organized as follows. In Sec. II, we present the concept
of the stack operation for tensor networks, as well as the proof using MPS as
an example. In Sec. III, we benchmark the speed and compare the machine
learning performance for three batch computing methods for the MPS machine
learning task: the naive loop method in Sec. III.1, the batch tensor network
in Sec. III.2, and the efficient coding method in Sec. III.3. In Sec. III.4,
we further discuss about the relationship between the batch tensor network and
the efficient coding method. We draw our conclusion and provide future
perspective in Sec. IV.

##  II Batch contraction of tensor networks

In this section, we provide the definition of the regular stack operation of
tensor networks used in the paper. We will first discuss some basic concepts
of tensor in Sec II.1 and tensor networks in Sec II.2. We then present the
details of the stack operation for MPS in Sec. II.3 and its generalized
formation about higher dimensional tensor network in Sec. II.4.

###  II.1 The stack operation

A tensor $\mathcal{T}$ is a specifically organized high-dimensional
collections of real or complex numbers. For example, a rank-$k$ tensor has $k$
indices. Therefore, a rank-$1$ tensor has only one index, and it is just a
vector; a rank-$3$ tensor is a higher-order tensor with three independent
indices. In this paper, we will use a tuple with parentheses
$(L_{1},L_{2},\dots,L_{K})$ to represent the rank-$K$ tensor

$\mathcal{T}\equiv\mathcal{T}_{p_{1},p_{2},\dots,p_{K}}=(L_{1},L_{2},\dots,L_{K}).$

where $p_{i}$ is the $i^{\text{th}}$ local index which takes the value from
$1$, $2$ to $L_{i}$, i.e. for each $i$, the local dimension is $L_{i}$.
Therefore, a vector which has $N$ dimension could be represented as a tuple
like $(N,)$ and a matrix whose size is $N\times M$ is $(N,M)$.

‘Stack’ is a pervasive operation in tensor analysis. It could merge a series
of tensors having the same shape into a higher rank tensor. For example, if we
stack $N$ vectors $(M,)$, then we get a $(N,M)$. If we stack $K$ matrix
$(N,M)$, then we get a rank-$3$ tensor $(K,N,M)$ [Fig. 1].

Figure 1: Illustration of stacking $4$ vectors of $(6,)$ into a single matrix
$(4,6)$.

In computer science, the stack operation generally converts the list of
tensors having the same shape into a very compact form that is memory
efficient. For example in Matlab, for a certain problem which requires going
through all the elements, transforming it into matrices production is much
more efficient than using a for-loop operation for the multiplication of row
vectors.

For modern machine learning or deep learning, the feeding input data is always
encoded as tensors. In order to deal with thousands of inputs, those tensors
are stacked together and reformulated into a higher-rank tensor called _batch_
before passing it into the model. This type of higher-rank tensor is also
referred to as stack operation in machine learning.

###  II.2 Basics of tensor networks

Tensor network $\\{T^{i}\\}$ is the graph factorization of a very large
tensors by some decomposition processes like SVD or QR. The correspondingly
tensor results from contracting all the auxiliary bonds of the tensor network.

$\mathcal{T}\xrightleftharpoons[\text{contraction}]{\text{decomposition}}\\{T^{i}\\}$

where $T_{i}$ represent the small local tensor in tensor network.

A graphical representation of tensors is shown in Fig. 2(a). In general, a
rank-$k$ tensor has $k$ indices. Therefore, a vector as the rank-$1$ has one
dangling leg, and a rank-$3$ tensor has three dangling legs, as shown in Fig.
2(a). The contraction between two tensors can thus be defined by contracting
the same index between two arbitrary tensors. In Fig. 2(a), a rank-$1$ tensor
$M_{k}$ is contracted with another rank-$3$ tensor $N_{i,j,k}$ as $M_{k}\cdot
N_{i,j,k}$, and it ends up with a rank-$2$ tensor, i.e. a matrix.

In the context of quantum physics, a many-body wave function can be
represented as a one-dimensional rank-$N$ tensor network where each leg
represents the physical index for each site, as shown graphically in Fig.
2(b). This is usually referred to as the matrix product state (MPS). The MPS
can therefore be seen as a connected array of local rank-$3$ tensors. Each
tensor $M_{j}$ has three indices: the physical index $\sigma_{j}$, and two
auxiliary indices $D_{j}$ and $D_{j+1}$ that are contracted and therefore
implicit 111Note that for the first and the last site, one of the auxiliary
indices is a dummy index and therefore can be neglected graphically.. An MPS
with open boundary condition can be written as

$\displaystyle|{\psi}\rangle=\sum_{\\{\bm{\sigma}\\}}\left[M_{0}^{\sigma_{0}}M_{1}^{\sigma_{1}}\dots
M_{n-1}^{\sigma_{n-1}}\right]|\sigma_{0},\dots,\sigma_{n-1}\rangle$ (1)

Figure 2:  Diagrammatic representations for tensor networks: (a) Diagrammatic
representations for tensors $M_{i}$ and $N_{i,j,k}$, and the contraction
between them $M_{i}\cdot N_{i,j,k}$; (b) Diagrammatic representation for one-
dimensional tensor network MPS: $M_{j}$ is the local rank-$3$ tensor and each
is contracted via auxiliary bond dimensions (vertical solid lines);
$\sigma_{j}$ is the index for the physical dimension.

###  II.3 Stack operation for MPS

A tensor $\mathcal{T}$ can be converted to a tensor network $\\{T^{i}\\}$
which consists of contracted local tensors $T_{i}$ via decomposition such as
the tensor train method [24]. Mathematically, one could perform operations
such as product, trace, contraction, splitting and grouping, etc. on both of
them.

However, the stack operation can only be directly applied to one single tensor
network rather than to multiple ones separately. More precisely, we want to
deal with such a problem as shown in Fig. 3: given a series of tensor networks
$\\{T^{i}\\}_{\alpha}$ with the same configuration, we want to get their
stacked tensor network representation with the same physical structure and
shape.

Figure 3:  Schematics of stack operations for MPS: Illustration for the $n$
identical MPS with auxiliary bond dimension $D$ which ends up with a new MPS
with same physics dimension. One of sites should has an extra bond dimension
$n$ to fulfill the stack dimension requirement. (Here we set the last site)

One direct method to realize this is to convert those tensor network
$\\{T^{i}\\}$ firstly to the corresponding tensor $(k,k,\dots,k)$, and then
decompose the stacked tensor $(n,k,k,\dots,k)$. However, this deviates from
the original intention of representing a tensor as a tensor network. Note that
we use a tensor network to avoid storing or expressing the full tensor in
computing, which is memory-consuming and may be prohibited in large systems.
Thus we would like to explore other approach that can help us efficiently to
establish the representation of the stacked tensor network without accessing
any contraction.

We show in Fig. 4: each rank-$3$ tensor in MPS is $(k,nD,nD)$. Each matrix
$(nD,nD)$ along physical dimension $k$, is a block diagonal matrix. The $i$-th
block $(D,D)$ is exactly the $i$-th rank-$3$ tensor in stack list at the same
physical dimension $k$ (so it is a matrix). There are minor differences
between the start and the end of the MPS. The start tensor unit is a rank-$2$
matrix with a row stack for the first unit in stack list. The end tensor is a
diagonal rank-$3$ tensor $(k,nD,n)$ with each diagonal $D\times 1$ block
filled by the end unit in stack list. A full mathematical proof is given in
Appendix. A.

Figure 4:  The diagram result shows the result of the ‘stack’ operation for
MPS. It will allocate a block diagonal local rank-$3$ tensor: for each
physical dimension, it gives a block diagonal matrix where the size of each
block is a matrix $(D,D)$. Note the first tensor unit’s block is a vector
$(D,)$ and the last unit’s block is a matrix $(D,1)$

###  II.4 Batch operation for general tensor networks

For other types of tensor network like PEPS, the unit tensor is a $4+1$ rank
tensor. The formula shown in II.3 would be more complicated. We post the
results here directly: For any tensor network consisting of $M$ tensor unit

$\\{T^{i}=(k,L_{1},L_{2},\dots,L_{m})|i=1,2,\dots,M\\}$

its $n$-batched tensor network version is another tensor network

$\\{T^{i}_{\text{batch}}=(k,nL_{1},nL_{2},\dots,nL_{m})|i=1,2,\dots,M\\}$

with one extra The extra index represents the stack dimension(corresponding to
the tensor stack) is free to add since each local tensor is block-diagonal. We
only need to unsqueeze/reshape the diagonal block from $(k,L_{1},L_{2},\dots)$
to $(k,1,L_{1},L_{2},\dots)$. For the MPS case, we could put it at the right
end, as shown in Fig. 4.

The sub-tensor $(nL_{1},nL_{2},\dots,nL_{m})$ along the physical indices for
each tensor is a diagonal block tensor. For example, the $j$-th diagonal block
is

$(T^{i}_{\text{batch}})_{k,p_{1},p_{2},\dots,p_{m}}=T^{i}_{j}$

where $p_{i}$ is the local index for index $i$ which takes the value from
$p_{i}\in[j*L_{i},(j+1)*L_{i}]$ and $T^{i}_{j}$ is the $T_{i}$ of the $j$-th
tensor network in tensor network list $\\{T_{i}\\}_{\alpha}$.

Meanwhile, one unit (for example, take the $N$-th local tensor) should be
‘unsqueezed’ to generate the dangling leg for the stack number indices

${T}^{N}_{batch}=(k,nL_{1},nL_{2},\dots,nL_{m},n).$

For example, the MPS case in II.3, will generate a row vector which is just
the ‘diagonal block’ form for the rank-$1$ tensor. And the last MPS unit would
get an external ‘leg’ with dimension $N$ is the requirement of the batch.

##  III Application

In this section, we study the potential applications of the stack of tensor
network. We apply this stack operation method to the tensor network machine
learning. For a tensor network based machine learning task, both the model and
input data can be collectively represented as a tensor or tensor network. In
fact, here, the model actually refers to the trainable parameters, and we
would note it as machine learning core in the following context. The forward
processing to calculate a response signal $y$ corresponding to the input $x$
is an inner product in a linear or nonlinear function space, as the
interaction between the machine learning core and the input data. The lost
function is designed based on the type of task. For example, we would measure
the distance between the real signal $\hat{y}$ and the calculated signal $y$;
for unsupervised learning, we might design the structure loss of hundreds of
signals $y$ to maximize the distance matrix. In this work, we take supervised
machine learning as an example. However, such a method is flexible for any
machine learning framework since it is a modeling technology. Here, by using
the language from quantum physics, we denote the machine learning core as
$|\mathcal{C}\rangle$ and one of the input data as
$|I_{j}=I_{j}(\psi)\rangle$. As both the core $|\mathcal{C}\rangle$ and the
input $|I_{j}\rangle$ are represented as tensors, this processing can be batch
executed in parallel. This machine learning task can then be efficiently
represented as

$\displaystyle\langle\mathcal{C}|(B,I_{j})\rangle=(B,\langle\mathcal{C}|I_{j}\rangle)=(B,\alpha_{j})$
(2)

where $\alpha_{j}$ is the resulting output data, and $|{(B,I_{j})}\rangle$ is
the batch representation of $B$ stacked tensors or a tensor network, where
$|{(B,I_{j})}\rangle$ has an explicit format.

In the following, we consider the one-dimensional case of the MPS machine
learning task as an example. As shown in Fig. 5, both the core
$|\mathcal{C}\rangle$ and the inputs $|I_{j}\rangle$ are the MPS, and the
output responsive signal $\alpha_{j}$ is the result obtained by contracting
both physical bonds from the two MPS.

Figure 5:  Schematics of inputting individual data as MPS
($|{I_{j}\left(\psi\right)}\rangle$) onto the core model
($\langle{\mathcal{C}}|$): each operation follows the arrow with the same
color, and results in an outcome $\alpha_{j}$.

Below we show that there are three ways to accomplish this machine learning
task: the loop (LP) method, the batched tensor network (BTN) method and the
efficient coding (EC) method. The core MPS $|\mathcal{C}\rangle$ consist of
$L$ local tensor units with $k$ being the dimension for all physical bonds and
$V$ being the dimension for all auxiliary bonds. Each input MPS
$|I_{j}\rangle$ consists of $L$ local tensor units with $k$ being the
dimension for all physical bonds and $D$ being the dimension for all auxiliary
bonds.

###  III.1 Naïve loop method

As shown in Fig. 5, for $B$ input samples, one needs to perform $B$
contractions to obtain $B$ outcomes. The most obvious way to accomplish this
is to perform contractions one by one. When we utilize the auto differential
feature from the libraries PyTorch [25] and Tensorflow [26], the gradient
information would be accumulated in every step. At the end of the loop, we can
acquire the batch-gradient by averaging over the sum, and apply it to the
gradient descent method for weight updating.

###  III.2 Batched Tensor Network Method

As discussed in Sec. II, we can stack all the inputs MPS to a batched tensor
network which behaves like another MPS. $|{(B,I_{j})}\rangle=|{I(\psi_{\rm
stack})}\rangle$. Assume the number of input is $B$, although the auxiliary
bond dimension of the batched MPS $|{(B,I_{j})}\rangle$ gets $B$ times larger
[see Fig. 4], the action of the inner product

$\displaystyle\langle\mathcal{C}|(B,I_{j})\rangle=\langle\mathcal{C}|{I(\psi_{\rm
stack})}\rangle$ (3)

can now be written as one single tensor contraction, as shown in Fig. 6

Figure 6:  Schematics of contraction between the core MPS
$\langle{\mathcal{C}}|$ (upper), and the batched MPS $|{(B,I_{j})}\rangle$
(down).

. With the help of the format of the batched tensor network, we can now
compress $B$ times contraction into one single contraction. Our advantage is
obvious: software libraries such as PyTorch and computational units such as
GPU are intrinsically optimized for tensor operation, so contracting a large
tensor is usually much faster than the matrix multiplication using a for loop.
As you can see in the benchmark shown in Fig. 8, when using the GPU
environment, the BTN method can get constant time complexity while the LP
method is the linear function of batch size.

Notice that we need to calculate an optimized contraction path to efficiently
contract the tensor network as shown in Fig. 6 while the LP method can simply
contract physical bond first and then sweep from left to right. When $B$ gets
larger, contracting physical bond first in order requires storing an MPS with
auxiliary bond equal to $nVB$, which is quite memory-consuming. Such a method
would cost extra memory allocation for the block diagonal tensor. The memory
increment is $O(B^{m})$, where $m$ is the number of the auxiliary bond. For
example, the $rank$-3 tensor of the batched MPS will allocate a $k\times
nB\times nB$ memory to store only $B\times k\times n\times n$ valid parameter.
If an autodifferential action is needed in the following process, the
intermediate tensor needed for gradient calculation will also become $B$ times
larger. Thus, although the speed is faster when the batch goes larger, the
memory will gradually blow up. The BTN method is a typical ‘memory swap speed’
method.

One possible solution is to treat the diagonal block tensor as a sparse
tensor, so we only need to store the valid indices and values. Meanwhile, if
we contract the physical bonds first, it will return a matrix chain whose unit
is also a sparse diagonal matrix of the shape $(nBV,nBV)$. This implies that
we can block-wise compute the matrix chain, so that each sparse unit performs
like a compact tensor $(B,nV,nV)$. At this juncture, it turns out that we can
then use the popular coding method called efficient coding, which is widely
used in modern tensor network machine learning task [11, 16, 18, 6, 14].

###  III.3 Efficient coding

The efficient coding (EC) does the ‘local batch contraction’ for every tensor
unit calculation in the tensor network contraction.

We first batch contract each physical index and get a tensor train consist of
batched units as shown in Fig. 7, then we perform the batch contraction on all
the auxiliary bonds for this batched tensor train.

Figure 7:  Illustration of efficient coding for MPS. Each unit from the core
MPS is batch contracted with input MPS (upper), and they form a batched unit
as a local tensor in the resulting MPS (down).

The batch contraction is realized by the built-in function einsim which is
provided in many modern scientific software libraries such as Numpy [27],
PyTorch [25] and Tensorflow [26]:

$\displaystyle\text{einsum('kij,Bknm}\rightarrow\text{Binjm',core,batch)}$ (4)
$\displaystyle\text{einsum('Bij,Bjk}\rightarrow\text{Bik',core,batch)}$ (5)

For those without einsim, one alternative solution is to implement a highly
efficient ‘batch matrix multiply’ function. The fundamental realization of
PyTorch.einsim is in fact the ‘batch matrix multiply’ called PyTorch.bmm.
Notice such a method is not new for tensor network machine learning, and Refs.
[11, 14, 16] has already realized successful learning algorithms based on
this.

The memory allocated for EC is much less than that for the BTN, which is the
major reason why the BTN method is slower than EC.

Also, the dense BTN may require an optimized contraction path for a large
batch case, while EC for MPS would only allow contraction of the physics bond
first, followed by the contraction for the auxiliary bonds as the tensor train
contraction. The standard tensor train contraction method sweeps from the
left-hand side of the MPS to the right. In particular, when we require all MPS
in the uniform shape, i.e., all the tensor units share the same auxiliary bond
dimension, we can then stack those units together($(B,L,V,V)$) and split it by
odd indices ($(B,L/2,V,V)$) or even indices ($(B,L/2,V,V)$). In doing so, we
can get the half-length tensor units ($(B,L/2,V,V)$) for the next iteration by
contracting the odd part and even part.

###  III.4 Comparison and discussion

Figure 8:  Speed-up benchmark results for different methods as a function of
batch size with: (a) Numpy on CPU; (b) PyTorch on CPU; (c) PyTorch on GPU; (d)
Tensorflow on GPU. The red, blue and green solid lines correspond to the
resulting obtained using the loop method (LP), efficient coding (EC), and the
batch tensor network method (BTN) respectively. All the panels are plotted in
$\log-\log$ format. All the result is test on a Intel(R) Core(TM) i7-8700K CPU
and 1080Ti GPU.

In Fig. 8, we show the speed-up benchmark results for all three different
methods introduced from Sec. III.1 to Sec. III.3 with respect to different
batch size in a $\log-\log$ plot. We compare the performance using different
libraries in Python. We mainly test three batch contraction methods: the Loop
method (LP), the batch tensor network method (BTN), and the efficient coding
method (EC) on the MPS batch contraction task in Fig. 5. The system size is a
$L=20$ units MPS with each unit assigned $V=6$, $k=3$ and $D=1$ [see Fig. 7].
Batch contraction task size from $1$ to $1000$. Overall, EC has the best
speed-up performance over other methods for almost all libraries, especially
at large batch size. Due to the pre-processing requirement, the BTN can be
slower than others in small-batch cases. The CPU is not efficient in storing
and processing larger matrices or tensors, while the GPU is deeply optimized
to deal with tensor data structure. Thus, the BTN gets similarly linear time
complexity like LP in CPU but maintain a constant time complexity in GPU. Note
the BTN will allocate tensors with $(nB,nB)$. Thus the larger the batch size,
the larger the intermediate tensor. We can see that the PyTorch(CPU) gets
better optimization with large tensor than NumPy. In Fig. 8(c) and (d), both
BTN and EC have constant time complexity in GPU since the GPU would compress
the time complexity of matrix contraction to a constant level. It also imply
that there is a close relationship between EC and BTN: As we discussed at the
end of the Sec. III.2, the EC is the another realization of the sparse matrix
BTN. The redundancy memory requirement for the diagonal block tensor slows
down the BTN method, but such redundancy requirement is eliminated in EC.

For EC and LP method, we contract the physical bond first, and obtain a tensor
chain-like structure

$(B,V)-(B,V,V)-(B,V,V)-\dots-(B,V,V)-(B,V,O)$

where $B=1$ for the LP method and $B>1$ for the EC method. Notice that we
assume $D=1$ here for the traditional MPS machine learning task. The memory
complexity costs $(L-1)\times BV^{2}+BVO+BV$ as the $O(B)$. Since the output
bond is at the right-hand side, the optimal path is contracting from left end
to the right end. Here, we use the auto differential engine to automatically
calculate the gradient for each unit, then the computer records all the
intermediates during the forward calculation. We therefore obtain the $(B,V)$
batched vector calculated by a batched vector-matrix dot for the first two
units

$\displaystyle(B,V)\times(B,V,V)\rightarrow(B,V)\quad\text{for the rest}$
$\displaystyle(B,V)\times(B,V,O)\rightarrow(B,O)\quad\text{for the last}$

So the memory cost for EC and LP method is straightforward as saving $L-1$
times batched vector-matrix dot result, which is $(L-2)\times BV+BO$. If we
apply the same contraction strategy to the BTN method, the resulted chain is

$(BV,)-(BV,BV)-\dots-(BV,BV)-(BV,BO).$

The memory complexity then costs $(L-1)\times B^{2}V^{2}+B^{2}DVO+BV$ memory
as the $O(B^{2})$. Surprisingly, the intermediate cost is the same:
$(L-2)\times BV+BO$. However, the left-right sweep path is not always the
optimal path for the BTN method, as we discussed in Sec. III.2. In our
experiment, the memory complexity for the optimal path consumes much less than
$O(B^{2})$. In Fig. 9, we give an example of the memory cost for increasing
batch size which takes a $L=21$ units of MPS with each unit assigned with
$V=50,k=3,D=1$ as the system parameters. The output $O=10$ classes leg is set
at the right end. The batch size increases from $B=10$ to $B=1500$. We also
plot the $B$ times larger curve of the $EC$, and we simply labeled it as
$EC*B$. We could see that although the BTN method takes more memory than the
EC method, the memory complexity is much less than $O(B^{2})$, which makes it
an available option for machine learning. Also, for other problems which
require small batches but larger auxiliary bond dimensions, the BTN will
perform better than the EC method as shown in Fig. 10 with a
$B=10,k=3,D=1,O=10$ and $V=10\rightarrow 3000$ MPS learning system.
Contracting the auxiliary bonds first and keeping comparably smaller
intermediates is a better choice when $B<<V$. This is the reason why we see
the BTN method can allocate much less memory than the EC method for large
auxiliary dimensions in Fig. 10

Figure 9: The memory needed for the increasing batch size MPS machine learning
problem. The system is an MPS of $L=21$ units, with each unit assigned with
$V=50,k=3,D=1,O=10$. The green curve($BTN$) is the required memory for BTN
method; The blue curve($EC$) is the required memory for EC method; The red
curve($EC*B$) is the $B$ times larger curve of $EC$.  Figure 10: The memory
needed for the increasing auxiliary bond dimension MPS machine learning
problem. The system is an MPS with $L=21$ units, with each unit assigned with
$B=10,k=3,D=1,O=10$. The green curve($BTN$) is the required memory for BTN
method; The blue curve($EC$) is the required memory for EC method.  Figure 11:
Evaluation of the train accuracy of MPS machine learning tasks on MNIST
dataset using three batch contraction method: (red solid curve) the batch
tensor network method (BTN); (green dashed curve) the efficient coding; (blue
dotted curve) the loop method (LP). Figure 12:  Evaluation of the valid
accuracy of MPS machine learning tasks on MNIST dataset using three batch
contraction method: (red solid curve) the batch tensor network method (BTN);
(green dashed curve) the efficient coding; (blue dotted curve) the loop method
(LP).

We also demonstrate a practice MPS machine learning task in Fig. 11 and Fig.
12. The dataset we used here is the MNIST handwritten digit database [28]. The
machine learning task requires predicting precisely the number from $0$ to $9$
among the dataset. The core tensor network is a ring MPS consist of $24\times
24$ unit with each unit is $(V,k,V)$ except the last one $(C,V,k,V)$. Here,
$V=20$ is the bond dimension, $k=2$ is the physics bond dimension, and $C=10$
is the number of classes. The bond dimension of the input tensor network is
$D=1$. We use the ring MPS as it is easier to train than the open boundary
MPS. Fig. 12 shows the valid accuracy v.s. time. The valid accuracy is
measured after each epoch and only plots $5$ epochs since the model get over-
fitting later. Each epoch will swap the whole dataset once. The best valid
accuracy for this setup is around 92%, while the state-of-the-art MPS machine
learning task [16] can reach 99% with a large bond 222The question on how to
perform training in MPS machine learning, for example, how to initial the
weight and pre-process data, is beyond the scope of this work, and we would
refer to Ref. [11, 16, 18, 6, 14] for reference..

The result is tested on an Intel(R) Core(TM) $i7-8700K$ CPU and $1080Ti$ GPU
and is shown in Fig. 11. The optimizer we used here is Adadelta[30] with the
learning rate $\text{\rm lr=}0.001$. The batch size is $100$. The random seed
is $1$. The test environment is PyTorch. The $x$-axis in Fig. 11 is the time
cost. Notice these three methods are the only different approaches to realize
the batch contraction so that they won’t influence the contraction result.
That is why these three curves in Fig. 11 are exactly the same. The EC method
is the best choice for the tensor network machine learning. The LP method is
the most inefficient method, which would cost one hour to traverse the dataset
while the EC only needs $60$ seconds and BTN takes $300$ seconds.

##  IV Conclusions

We have investigated the stack realization of the tensor network as the spread
concept of the _stack_ operation for tensors. The resulting stacked/batched
tensor network consists of block-diagonal local tensors larger bonds. We
connected this method to the efficient coding batch contraction technique
which is widely used in tensor network machine learning. An MPS machine
learning task has been used as an example to validate its function and
benchmark the performance for different batch sizes, different numerical
libraries, and different chips. Our algorithm provides an alternative way of
realizing fast tensor network machine learning in GPU and other advanced
chips. All the code used in this work is available at [31].

Several possible future directions may be of significant interest. Our work
reveals an intrinsic connection between the block-diagonal tensor network
units and the batched contraction scheme, which is potentially useful for
helping researchers realize faster tensor network implementation with the
block-diagonal design like [17]. Secondly, our algorithm comes without
performing SVD to the stacked tensor network; It would be great to define a
batch SVD and check its performance and error when truncating the singular
values for larger data size. Thirdly, the efficient coding method requires the
contraction of the physical dimensions first. However, the batch method
provides the possibility to start contraction in other dimensions first, which
may be useful for the contraction on a two-dimensional PEPS or through tensor
renormalization group method [32, 33, 34]. Finally, the BTN method takes
advantage of compressing thousands of contractions into one contracting graph.
From the hardware design perspective, with a specially designed platform
optimized for tensor network contraction, the BTN method may intrinsically
accelerate the tensor network machine learning or any task involving a stacked
tensor network. Meanwhile, the tensor network states in many-body quantum
physics has a non-trivial relationship with the quantum circuits. So the idea
of representing several tensor contracting into one single contraction also
provides the possibility of efficiently computing on NISQ-era quantum
computers.

###### Acknowledgements.

This work was supported by USA Office of Naval Research Global
(N62909-19-1-2047) and SUTD-ZJU Visiting Professor (VP 201303). T.Z
acknowledges the support of Singapore Ministry of Education PhD Research
Scholarship.

##  Appendix A Mathematical proof for stack operation of two MPS

For simplicity, we consider performing the stack operation of two MPS:
$|{\psi_{\alpha}}\rangle,|{\psi_{\beta}}\rangle$, which can be written as

$\displaystyle\begin{bmatrix}|{\psi_{\alpha}}\rangle,&|{\psi_{\beta}}\rangle\end{bmatrix}=|{\psi_{\alpha}}\rangle\otimes|{0}\rangle+|{\psi_{\beta}}\rangle\otimes|{1}\rangle$
(6)

and

$\displaystyle|\psi_{\alpha}\rangle=\sum_{\\{\bm{\sigma}\\}}\mathrm{Tr}\left[M_{\alpha}^{\sigma_{0}}M_{\alpha}^{\sigma_{1}}\dots
M_{\alpha}^{\sigma_{n-1}}\right]|\sigma_{0},\dots,\sigma_{n-1}\rangle$ (7)
$\displaystyle|\psi_{\beta}\rangle=\sum_{\\{\bm{\sigma}\\}}\mathrm{Tr}\left[M_{\beta}^{\sigma_{0}}M_{\beta}^{\sigma_{1}}\dots
M_{\beta}^{\sigma_{n-1}}\right]|\sigma_{0},\dots,\sigma_{n-1}\rangle$

are the two MPS’s. We seek to find a new MPS $|{\phi}\rangle$ which exactly
represents the stacked operation from Eq. (6). We denote the new MPS as

$\displaystyle|\phi\rangle=\sum_{\\{\bm{\sigma}\\}}\mathrm{Tr}\left[\tilde{M}^{\sigma_{0}}\tilde{M}^{\sigma_{1}}\dots\tilde{M}^{\sigma_{n-1}}\right]|\sigma_{0},\dots,\sigma_{n-1}\rangle$
(8)

To see the detailed structure of each $\tilde{M}^{\sigma_{i}}$ in the above
expression, we expand both terms on the right-hand side of Eq. (6)

$\displaystyle|\psi_{\alpha}\rangle|0\rangle$ (9)
$\displaystyle=\sum_{\sigma_{0},\dots,\sigma_{n-1}}\mathrm{Tr}\left[M_{\alpha}^{\sigma_{0}}M_{\alpha}^{\sigma_{1}}\dots
M_{\alpha}^{\sigma_{n-1}}\right]|\sigma_{0},\dots,\sigma_{n-1}\rangle|0\rangle$
$\displaystyle=\sum_{\sigma_{0},\dots,{\sigma_{n}}}\mathrm{Tr}\left[M_{\alpha}^{\sigma_{0}}M_{\alpha}^{\sigma_{1}}\dots
M_{\alpha}^{\sigma_{n-1}}\right]\delta_{0}^{\sigma_{n}}|\sigma_{0},\dots,\sigma_{n-1}\rangle|\sigma_{n}\rangle$
$\displaystyle=\sum_{\sigma_{0},\dots,{\sigma_{n}}}\mathrm{Tr}\left[M_{\alpha}^{\sigma_{0}}M_{\alpha}^{\sigma_{1}}\dots
M_{\alpha}^{\sigma_{n-1}}\Delta_{0}^{\sigma_{n}}\right]|\sigma_{0},\dots,\sigma_{n-1},\sigma_{n}\rangle$

where $\Delta_{0}^{\sigma_{n}}=\delta_{0}^{\sigma_{n}}\mathbb{I}$, and
$\mathbb{I}$ is the identity matrix of which the row size is the same as the
column size of the previous site matrix $M_{\sigma}^{\sigma_{n-1}}$.

Similarily, we have

$\displaystyle|\psi_{\beta}\rangle|1\rangle$ (10)
$\displaystyle=\sum_{\sigma_{0},\dots,{\sigma_{n}}}\mathrm{Tr}\left[M_{\beta}^{\sigma_{0}}M_{\beta}^{\sigma_{1}}\dots
M_{\beta}^{\sigma_{n-1}}\Delta_{1}^{\sigma_{n}}\right]|\sigma_{0},\dots,\sigma_{n-1},\sigma_{n}\rangle$

By summing up Eq. (9) and (10), we obtain

$\displaystyle|{\psi_{\alpha}}\rangle\otimes|{0}\rangle+|{\psi_{\beta}}\rangle\otimes|{1}\rangle$
(11)
$\displaystyle=\sum_{\\{\bm{\sigma}\\}}\mathrm{Tr}\left[M^{\sigma_{0}}_{\alpha}M^{\sigma_{1}}_{\alpha}\dots
M^{\sigma_{n-1}}_{\alpha}\Delta_{0}^{\sigma_{n}}\right]|\sigma_{0},\dots,\sigma_{n-1},\sigma_{n}\rangle$
$\displaystyle+\sum_{\\{\bm{\sigma}\\}}\mathrm{Tr}\left[M^{\sigma_{0}}_{\beta}M^{\sigma_{1}}_{\beta}\dots
M^{\sigma_{n-1}}_{\beta}\Delta_{1}^{\sigma_{n}}\right]|\sigma_{0},\dots,\sigma_{n-1},\sigma_{n}\rangle$
$\displaystyle=\sum_{\\{\bm{\sigma}\\}}\mathrm{Tr}\left[M^{\sigma_{0}}_{\alpha}M^{\sigma_{1}}_{\alpha}\dots
M^{\sigma_{n-1}}_{\alpha}\Delta_{0}^{\sigma_{n}}\right.$
$\displaystyle\left.+M^{\sigma_{0}}_{\beta}M^{\sigma_{1}}_{\beta}\dots
M^{\sigma_{n-1}}_{\beta}\Delta_{1}^{\sigma_{n}}\right]|{\sigma_{0},\sigma_{1},\dots,\sigma_{n-1},\sigma_{n}}\rangle$
$\displaystyle=\sum_{\\{\bm{\sigma}\\}}\mathrm{Tr}\left\\{\left[\begin{array}[]{c}M_{\alpha}^{\sigma_{0}}\\\
M_{\beta}^{\sigma_{0}}\end{array}\right]^{T}\left[\begin{array}[]{cc}M_{\alpha}^{\sigma_{0}}&0\\\
0&M_{\beta}^{\sigma_{0}}\\\ \end{array}\right]\cdots\right.$ (16)
$\displaystyle\left.\left[\begin{array}[]{cc}M_{\alpha}^{\sigma_{n-1}}&0\\\
0&M_{\beta}^{\sigma_{n-1}}\\\
\end{array}\right]\left[\begin{array}[]{cc}\Delta_{0}^{\sigma_{n}}&0\\\
0&\Delta_{1}^{\sigma_{n}}\\\
\end{array}\right]\right\\}|{\sigma_{0},\sigma_{1},\dots,\sigma_{n-1},\sigma_{n}}\rangle$
(21)
$\displaystyle=\sum_{\\{\bm{\sigma}\\}}\mathrm{Tr}\left[\tilde{M}^{\sigma_{0}}\tilde{M}^{\sigma_{1}}\dots\tilde{M}^{\sigma_{n-1}}\tilde{\Delta}^{\sigma_{n}}\right]|{\sigma_{0},\sigma_{1},\dots,\sigma_{n-1},\sigma_{n}}\rangle$
$\displaystyle=\sum_{\\{\bm{\sigma}\\}}\mathrm{Tr}\left[\tilde{M}^{\sigma_{0}}\tilde{M}^{\sigma_{1}}\dots\tilde{M}^{\sigma_{n-1}}\right]|{\sigma_{0},\sigma_{1},\dots,\sigma_{n-1},\sigma_{n}}\rangle$

Note that for the first site $\sigma_{0}$, the resulting matrix
$\tilde{M}^{\sigma_{0}}$ is stacked in row geometry [see the first site of the
resulting stacked MPS in Fig. 4], all the other matrices
$\tilde{M}^{\sigma_{j}}(j=1,\cdots,n-1)$ are block diagonal [Fig. 4], as the
stack operation is essentially equivalent to a direct sum of local matrices.
In the final step of the above equation, to maintain the form of the stacked
MPS as shown in Eq. (8), we have absorbed $\tilde{\Delta}^{\sigma_{n}}$ into
$\tilde{M}^{\sigma_{n-1}}$ and reformulate it as

$\displaystyle\tilde{M}^{\sigma_{n-1}}\rightarrow\tilde{M}^{\sigma_{n-1}}\tilde{\Delta}^{\sigma_{n}}$
(22) $\displaystyle=\left[\begin{array}[]{cc}M_{\alpha}^{\sigma_{n-1}}&0\\\
0&M_{\beta}^{\sigma_{n-1}}\\\
\end{array}\right]\left[\begin{array}[]{cc}\Delta_{0}^{\sigma_{n}}&0\\\
0&\Delta_{1}^{\sigma_{n}}\\\ \end{array}\right]$ (27)
$\displaystyle=\left[\begin{array}[]{cc}M_{\alpha}^{\sigma_{n-1}}\Delta_{0}^{\sigma_{n}}&0\\\
&M_{\beta}^{\sigma_{n-1}}\Delta_{1}^{\sigma_{n}}\\\ \end{array}\right]$ (30)

which is still block diagonal, but each block with size $D\times 1$.

The above process is the stack operation for two MPS, and it is easy to
extrapolate to the stack operation of $n$ MPS. For simplicity, we assign MPS
with the same auxiliary dimension $D$ for each bond. It is not necessary as
they can be different. One only needs to change the diagonal block indices
from $(D,D)$ to $(D_{i},D_{j})$.

## References

  * Jouppi  _et al._ [2017] N. P. Jouppi, C. Young, N. Patil, D. Patterson, G. Agrawal, R. Bajwa, S. Bates, S. Bhatia, N. Boden, A. Borchers, R. Boyle, P.-l. Cantin, C. Chao, C. Clark, J. Coriell, M. Daley, M. Dau, J. Dean, B. Gelb, T. V. Ghaemmaghami, R. Gottipati, W. Gulland, R. Hagmann, C. R. Ho, D. Hogberg, J. Hu, R. Hundt, D. Hurt, J. Ibarz, A. Jaffey, A. Jaworski, A. Kaplan, H. Khaitan, D. Killebrew, A. Koch, N. Kumar, S. Lacy, J. Laudon, J. Law, D. Le, C. Leary, Z. Liu, K. Lucke, A. Lundin, G. MacKean, A. Maggiore, M. Mahony, K. Miller, R. Nagarajan, R. Narayanaswami, R. Ni, K. Nix, T. Norrie, M. Omernick, N. Penukonda, A. Phelps, J. Ross, M. Ross, A. Salek, E. Samadiani, C. Severn, G. Sizikov, M. Snelham, J. Souter, D. Steinberg, A. Swing, M. Tan, G. Thorson, B. Tian, H. Toma, E. Tuttle, V. Vasudevan, R. Walter, W. Wang, E. Wilcox, and D. H. Yoon, In-datacenter performance analysis of a tensor processing unit, SIGARCH Comput. Archit. News 45, 1–12 (2017).
  * Verstraete  _et al._ [2008] F. Verstraete, V. Murg, and J. Cirac, Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems, Advances in Physics 57, 143 (2008).
  * Orús [2014] R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Annals of Physics 349, 117 (2014).
  * Orús [2019] R. Orús, Tensor networks for complex quantum systems, Nature Reviews Physics 1, 538 (2019).
  * Cirac  _et al._ [2021] J. I. Cirac, D. Pérez-García, N. Schuch, and F. Verstraete, Matrix product states and projected entangled pair states: Concepts, symmetries, theorems, Rev. Mod. Phys. 93, 045003 (2021).
  * Stoudenmire and Schwab [2016] E. Stoudenmire and D. J. Schwab, Supervised learning with tensor networks, in  _Advances in Neural Information Processing Systems_, Vol. 29, edited by D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett (Curran Associates, Inc., 2016).
  * Gao and Duan [2017] X. Gao and L.-M. Duan, Efficient representation of quantum many-body states with deep neural networks, Nature Communications 8, 662 (2017).
  * Deng  _et al._ [2017] D.-L. Deng, X. Li, and S. Das Sarma, Quantum entanglement in neural network states, Phys. Rev. X 7, 021021 (2017).
  * Stoudenmire [2018] E. M. Stoudenmire, Learning relevant features of data with multi-scale tensor networks, Quantum Science and Technology 3, 034003 (2018).
  * Chen  _et al._ [2018] J. Chen, S. Cheng, H. Xie, L. Wang, and T. Xiang, Equivalence of restricted boltzmann machines and tensor network states, Phys. Rev. B 97, 085104 (2018).
  * Han  _et al._ [2018] Z.-Y. Han, J. Wang, H. Fan, L. Wang, and P. Zhang, Unsupervised generative modeling using matrix product states, Phys. Rev. X 8, 031012 (2018).
  * Glasser  _et al._ [2019] I. Glasser, N. Pancotti, and J. I. Cirac, From probabilistic graphical models to generalized tensor networks for supervised learning (2019), arXiv:1806.05964 [quant-ph] .
  * Liu  _et al._ [2019] D. Liu, S.-J. Ran, P. Wittek, C. Peng, R. B. García, G. Su, and M. Lewenstein, Machine learning by unitary tensor network of hierarchical tree structure, New Journal of Physics 21, 073059 (2019).
  * Efthymiou  _et al._ [2019] S. Efthymiou, J. Hidary, and S. Leichenauer, Tensornetwork for machine learning (2019), arXiv:1906.06329 [cs.LG] .
  * Roberts  _et al._ [2019] C. Roberts, A. Milsted, M. Ganahl, A. Zalcman, B. Fontaine, Y. Zou, J. Hidary, G. Vidal, and S. Leichenauer, Tensornetwork: A library for physics and machine learning (2019), arXiv:1905.01330 [physics.comp-ph] .
  * Cheng  _et al._ [2021] S. Cheng, L. Wang, and P. Zhang, Supervised learning with projected entangled pair states, Phys. Rev. B 103, 125117 (2021).
  * Vieijra  _et al._ [2022] T. Vieijra, L. Vanderstraeten, and F. Verstraete, Generative modeling with projected entangled-pair states (2022), arXiv:2202.08177 [quant-ph] .
  * Cheng  _et al._ [2019] S. Cheng, L. Wang, T. Xiang, and P. Zhang, Tree tensor networks for generative modeling, Phys. Rev. B 99, 155131 (2019).
  * Novikov  _et al._ [2015] A. Novikov, D. Podoprikhin, A. Osokin, and D. P. Vetrov, Tensorizing neural networks, in  _Advances in Neural Information Processing Systems_, Vol. 28, edited by C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett (Curran Associates, Inc., 2015).
  * Guo  _et al._ [2018] C. Guo, Z. Jie, W. Lu, and D. Poletti, Matrix product operators for sequence-to-sequence learning, Phys. Rev. E 98, 042114 (2018).
  * Weichselbaum [2012] A. Weichselbaum, Non-abelian symmetries in tensor networks: A quantum symmetry space approach, Annals of Physics 327, 2972 (2012).
  * Hong  _et al._ [2022] R. Hong, Y.-X. Xiao, J. Hu, A.-C. Ji, and S.-J. Ran, Functional tensor network solving many-body schrödinger equation (2022).
  * Note [1] Note that for the first and the last site, one of the auxiliary indices is a dummy index and therefore can be neglected graphically.
  * Oseledets [2011] I. V. Oseledets, Tensor-train decomposition, SIAM Journal on Scientific Computing 33, 2295 (2011).
  * Paszke  _et al._ [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, Pytorch: An imperative style, high-performance deep learning library, in  _Advances in Neural Information Processing Systems 32_, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019) pp. 8024–8035.
  * Abadi  _et al._ [2015] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, TensorFlow: Large-scale machine learning on heterogeneous systems (2015), software available from tensorflow.org.
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Array programming with NumPy, Nature 585, 357 (2020).
  * LeCun  _et al._ [1998] Y. LeCun, C. Cortes, and C. J. Burges, Mnist handwritten digit database (1998).
  * Note [2] The question on how to perform training in MPS machine learning, for example, how to initial the weight and pre-process data, is beyond the scope of this work, and we would refer to Ref. [11, 16, 18, 6, 14] for reference.
  * Zeiler [2012] M. D. Zeiler, Adadelta: An adaptive learning rate method (2012), arXiv:1212.5701 [cs.LG] .
  * [31] https://github.com/veya2ztn/Stack_of_Tensor_Network.
  * Levin and Nave [2007] M. Levin and C. P. Nave, Tensor renormalization group approach to two-dimensional classical lattice models, Phys. Rev. Lett. 99, 120601 (2007).
  * Gu  _et al._ [2008] Z.-C. Gu, M. Levin, and X.-G. Wen, Tensor-entanglement renormalization group approach as a unified method for symmetry breaking and topological phase transitions, Phys. Rev. B 78, 205116 (2008).
  * Evenbly and Vidal [2015] G. Evenbly and G. Vidal, Tensor network renormalization, Phys. Rev. Lett. 115, 180405 (2015).