# NUTA: Non-uniform Temporal Aggregation for Action Recognition

Xinyu Li equally contribute Chunhui Liu11footnotemark: 1 Bing Shuai  Yi Zhu
Hao Chen  Joseph Tighe
Amazon Web Services
{xxnl, chunhliu, bshuai, yzaws, hxen<EMAIL_ADDRESS>
###### Abstract

In the world of action recognition research, one primary focus has been on how
to construct and train networks to model the spatial-temporal volume of an
input video. These methods typically uniformly sample a segment of an input
clip (along the temporal dimension). However, not all parts of a video are
equally important to determine the action in the clip. In this work, we focus
instead on learning where to extract features, so as to focus on the most
informative parts of the video. We propose a method called the non-uniform
temporal aggregation (NUTA), which aggregates features only from informative
temporal segments. We also introduce a synchronization method that allows our
NUTA features to be temporally aligned with traditional uniformly sampled
video features, so that both local and clip-level features can be combined.
Our model has achieved state-of-the-art performance on four widely used large-
scale action-recognition datasets (Kinetics400, Kinetics700, Something-
something V2 and Charades). In addition, we have created a visualization to
illustrate how the proposed NUTA method selects only the most relevant parts
of a video clip.

##  1 Introduction

A key challenge in action recognition is how to learn a feature that captures
the relevant spatial and motion queues in an efficient and compact
representation. This problem has been well studied, primarily using
convolution neural networks (CNNs) [12], from frame-based methods [16] to
segment based temporal information aggregation [35] to the I3D based methods
[2]. These methods have primarily focused on how to perform feature extraction
that captures a complete spatial-temporal description of the video. The
majority of these methods treat each frame, or point in time, with equal
weight, but not all parts of the video are equally important and thus it is
also key that we develop feature extraction methods that can determine where
to extract features from. Some recent works have started to look at this
problem with non-local modules [36], directional convolutions [18] or
ensembles of features at different sample rates [8, 42], but these methods
still focus more on how to comprehensively extract features rather than where
to extract the features. In this work we leverage all the advancements that
have been made in the past years on how to extract motion features and focus
our attention on the question of where to extract these features from.

There are two major challenges to answer the “where to” question. First, non-
uniformly extracting features efficiently from only informative temporal
episodes is challenging as it is required to look at the whole video to
determine which parts are informative. Some recent work [17, 40] have proved
that a recognition system can benefit from selecting the informative frames
rather than simply taking the uniformly sampled frames as inputs. However,
these systems treat the frame selection and feature extraction as two separate
stages and thus the frame selection can not benefit from the later feature
extraction thus reducing the descriptive power of the network and adding
redundancy in the two stages.

Figure 1:  A comparison of different temporal modeling methods. The temporal
pooling/striding and temporal convolution aggregate features uniformly along
time. Our non-uniform temporal aggregation (NUTA) extracts features non-
uniformly from informative parts of the video clip.

Second, exchanging information between different feature sampling strategies
efficiently and comprehensively remains a challenge. Other works have tackled
this using kernels [20, 24, 36], LSTMs [6] or most recently feature pyramids
[42] but each of these methods only partially deals with this information
exchange problem.

In order to address the first challenge (how to extract feature from only
representative temporal instances), we introduce a non-uniform temporal
aggregation (NUTA) method. In contrast to the average pooling operator or
temporal stride convolution kernels (Figure 1), the NUTA module generates a
temporal map that non-uniformly projects $T$ temporal instances to
$T^{\prime}$ ($T^{\prime}<T)$. The NUTA has receptive field in temporal
dimension as large as the entire clip so the NUTA module is able to generate
clip-level features. Differing from previous two-stage systems that perform
frame selection first and then do recognition on the selected frames [17, 40],
our NUTA module performs temporal information selection at the feature level
in an end-to-end manner. Differing from previous work that focus on “stacking”
features of the same kind [8, 36] to construct the comprehensive feature
representation, NUTA works by focusing on a sub-set of the frames that are the
most informative for the task. To address the second challenge (information
exchange between network branches), we propose a two-branch network with a
novel temporal synchronization strategy. The two-branch design maintains the
uniformly aggregated local feature branch and the non-uniformly aggregated
clip-level feature branch separately. Note that the purpose of our two-branch
design is to separate the feature descriptors instead of using multiple input
modalities to encourage different feature extraction behaviors as in [8, 42].
The temporal synchronization strategy temporally aligns the features from the
two branches to allow the information from each branch to flow to the other.
Such design not only improves the action recognition performance by removing
noise from both the uniform and non-uniform branches, but also reduces the
computational requirements of our network, allowing our two-branch design to
maintain similar computational cost to the original single branch network.

We tested our model on four of the most popular action recognition datasets:
Kinetics 400, Kinetics 700, Something Something V2 and Charades. Our NUTA
network achieves the state-of-the-art (SOTA) or comparable performance on all
four datasets with less computation cost. We perform a detailed ablation study
that shows our non-uniform aggregation method is effective across a wide range
of configurations. We also visualized our model to better demonstrate how our
NUTA behaves and how it helps answer the “where to” question. To summarize,
our contributions are:

  1. 1.

We propose a NUTA module to answer the “where to” question, which is able to
identify and sample only the informative parts of the video clips;

  2. 2.

We propose a two-branch NUTA network for action recognition, which learns
local and clip-level features by fusing the feature extracted from the best in
bread video modeling methods and our novel NUTA module;

  3. 3.

We propose information fusion and temporal synchronization methods that allow
both branches to specialize on their own aspects of video modelling while
sharing their feaures after each stage of the network;

  4. 4.

We present experimental results on four common, large scale action recognition
datasets, which shows state-of-the-art or better performance on each.

##  2 Related Work

The early video action recognition work used 2D CNNs [16, 29] to extract and
aggregate frame-wise features with no effective temporal modeling. Therefore,
the question of how to perform temporal modeling became a key area of
interest. RNNs and LSTMs [14] were used to model temporal associations among
features extracted by 2D networks [19, 43, 5, 21, 15, 23]. However, LSTM based
methods often significantly increased the computation with little benefit to
accuracy. Segment level predictions methods such as TSN [35] and rank pooling
[10, 1] based methods were proposed to better capture the temporal feature
evolution at the video level. Although quite successful, these methods still
followed the idea of aggregating 2D features instead of performing pixel-level
3D feature learning. More recently, 3D CNNs [2, 32, 34, 41, 33, 27, 25, 34]
have gained significant popularity due to their ability to perform spatio-
temporal feature modeling. The 3D CNNs effectively address the “how to”
question and are now used as basic building blocks for most temporal modeling
tasks. Most of recent works have continued to focus on “how to” perform
temporal modeling more comprehensively: the non-local network [36] leverages
non-local connections to establish pixel level long-term spatio-temporal
association; the TEA [44] and TAM [24] methods use a local and global kernel
for long-shot term feature learning; the SlowFast [8] and TPN [42] methods
ensemble feature at different frame-rate for better performance. There have
been a handful of papers that have tried to answer the “where to” question.
Early research tried to generate long-term feature with LSTM and temporal
attention [31, 6] but none of these methods ended up achieving competitive
performance. Two recent methods, SCSampler [17] and multi-agent sampler [40]
select only informative clips from a video for action recognition and have
been able to demonstrate better performance than uniformly sampling the frames
from a video.

Different from these previous work, we answer the “where to” question by
proposing the non-uniform temporal aggregation module which learn the feature
non-uniformly from informative temporal episodes. Instead of adding features,
e.g. from non-local connection or from multiple branches that takes different
types of input [36, 8], our network reduce redundant information for better
recognition. Different from two-stage system that first select representative
frames and run classifier on selected frames [17], our proposed method runs
end-to-end and performs feature level information selection.

Many two-branch architectures have been proposed to fuse information from
different sources, e.g. input frames at different resolution [16], RGB and
optical flow images [35, 29, 9], and video clips at different sample rate [8,
44]. Differing from most previous two-branch systems that fuse features from
different input modalities or representations, our two-branch design models
both different and complementary features from the same input.

##  3 Methodology

Figure 2: An overview of proposed NUTA network: a two-branch design (top). A
detailed view of our NUTA module and its interaction with uniform branch is
show on the bottom.

###  3.1 Overview

We propose the Non-uniform Temporal Aggregation (NUTA) network that follows a
two-branch design with one uniform branch that maintains the relative temporal
spacing of the input clip and captures temporally local information, and one
non-uniform branch for clip-level temporal feature selection (Figure 2). We
adopt a 3D convolutional style network for our uniform branch to learn local
spatial-temporal features. We propose a novel NUTA module and build our non-
uniform branch by stringing multiple NUTA modules together. Each NUTA module
is able to model the full clip-level temporal feature by sampling non-
uniformly along the temporal dimension. Our NUTA module does not perform any
spatial modeling, thus to enable the non-uniform branch to model the clip at
different spatial scales we transfer information between the two branches
after each stage. This information exchange between our uniform and non-
uniform branches is non trivial because each operates on different samples
from the input clip. To overcome this issue we introduce our temporal
synchronization strategy, which synchronizes the features from the uniform and
non-uniform branches bidirectionally. The rest of this section introduces the
components of our NUTA module including the non-uniform down-sampling, feature
synchronization and feature fusion strategies.

###  3.2 Non-uniform Temporal Down-sampling

A core feature of our NUTA network is our non-uniform sampling of the video
clip feature along the temporal dimension. This component is to only sample
the most informative parts to feature across time. We achieve this via our
down-sampling module that is parameterized as a learned projection map, which
maps a feature of temporal dimension $T$ to $T^{\prime}$
($T^{\prime}=\frac{T}{2}$ in this paper). Given a feature tensor
$\mathbf{F}\in\mathbb{R}^{C\times T\times W\times H}$, we apply a self-
attention-like module to learn the projection map $\mathbf{M}$ as:

$\mathbf{M}=\operatorname{softmax}\left(\Gamma[\phi(\mathbf{\operatorname{pool}(F}))]\times\Gamma[\theta(\mathbf{F})]^{T}\right)$
(1)

where $\operatorname{pool}(*)$ denotes temporal max-pooling with kernel size
of 2, functions $\phi(*)$ and $\theta(*)$ are convolution operators, and
$\Gamma[*]$ denotes a reshape and permutation operation. Then, the clip-level,
non-uniformly sampled feature $\mathbf{F_{NUTA}}$ is generated by applying the
temporal projection map as:

$\mathbf{F_{NUTA}}=\operatorname{Conv}\left(\Gamma^{-1}[\mathbf{M}\times\Gamma[\delta(\mathbf{F})]]\right)$
(2)

where the $\delta(*)$ denotes a convolution, the $\Gamma^{-1}$ is the inverse
reshape of $\Gamma$, and $\operatorname{Conv}$ denotes a convolution to
compress the channel dimension for information fusion with the next stage. In
practice, we use a $(3,1,1)$ convolution kernel for $\phi(*)$, $\theta(*)$ and
$\delta(*)$. We find such convolution gives a better performance than commonly
used $(1,1,1)$ kernel, which may due to it gives larger capacity on temporal
modeling.

Depending on the complexity of the scene and action classification needed,
there may be some set of frames that are most informative for one type of
action and a completely different set for another type of action. To allow our
network to have the flexibility of picking different sets of temporal periods
for each channel we extend our projection map by making it multi-head. Each
head is then free to select its own set of frames and extract features from
those frames. We implement this using our $\Gamma[*]$ operation that split
$\phi(\mathbf{F})\in\mathbb{R}^{T\times C\times W\times H}$ to $n$ heads along
the channels, such that $\Gamma[\phi(\mathbf{F})]\in\mathbb{R}^{n\times
T\times\frac{C}{n}WH}$. To make our design computationally efficient, we use
the group convolution and empirically set groups as 64 (see section ablation
study for details).

Although it may appear similar, our proposed temporal modeling is
significantly different from non-local network [36] as our generated
projection map $\mathbf{M}\in\mathbb{R}^{n\times\frac{T}{2}\times T}$
represents a temporal association based on spatial features while the non-
local kernel $\mathbf{M^{\prime}}\in\mathbb{R}^{TWH\times TWH}$ represents a
pixel level spatial temporal association learned from channels.

###  3.3 Temporal Synchronization

The output feature from our NUTA module
$\mathbf{F_{NUTA}}_{n}\in\mathbb{R}^{C\times\frac{T}{2}\times W\times H}$ is
not compatible with the uniform branch feature
($\mathbf{F_{res}}_{n}\in\mathbb{R}^{C\times T\times W\times H}$) for two
reasons: 1) the temporal dimensions do not match ($\frac{T}{2}\neq T$) and 2)
they are not sampled from the same set of frames. To handle this
incompatibility we propose to synchronize the features from the non-uniform
branch to the uniform branch (Figure 2 yellow arrow) as:

$\mathbf{\mathbf{F}_{res}}=\operatorname{Conv}\left(\Gamma^{-1}[\mathbf{M}\times\Gamma[\zeta(\mathbf{F_{res}})]]\right)+\operatorname{pool}(\mathbf{F_{res}})$
(3)

where $\zeta$ denotes a 3D convolution with kernel $(3,1,1)$ similar to
$\theta$ and $\phi$. Operation
$\operatorname{Conv}\left(\Gamma^{-1}[\mathbf{M}\times\Gamma[\zeta(\mathbf{F_{res}})]]\right)$
can be intuitively understood as applying the same temporal projection to the
uniform branch feature $\mathbf{F_{res}}$, so that feature is synchronized
with clip-level feature from the non-uniform branch. Note that although we
apply the same temporal projection map to the feature from uniform branch, the
following 3D resnet layers will still generate the local features. The
temporal synchronization removes the non-informative features from uniform
branch, which helps improve the action recognition performance and reduces the
computation cost (see ablation study section).

###  3.4 Information Fusion

The temporal synchronization handles the synchronization and combination of
the non-uniform features with the uniform branch, however it doesn’t encode
the spatial information learned in the uniform branch into the non-uniform
branch. To address this, we fuse the uniform branch feature from res-stage $n$
and non-uniform feature from NUTA-stage $n-1$ and feed them to NUTA-stage $n$
as follows:

$\mathbf{F}_{n}=f[\mathbf{F_{NUTA}}_{n-1},\mathbf{F_{res}}_{n}]$ (4)

where $\mathbf{F}_{n}$ stands for the fused feature at $n_{th}$ stage.
$\mathbf{F_{NUTA}}_{n-1}$ denotes the feature from non-uniform branch (NUTA
module) and $\mathbf{F_{res}}_{n}$ is the feature from the uniform branch
(ResNet feature backbone). We consider the following three fusion strategies
$f$: (1). Non-local fusion ($f_{n}$): where we use $\mathbf{F_{NUTA}}_{n-1}$
as query to pull information from $\mathbf{F_{res}}_{n}$. (2). Feature sum-up
($f_{s}$): where we apply the 1D convolution to $\mathbf{F_{NUTA}}_{n-1}$ and
apply the pixel wise addition to $\mathbf{F_{res}}_{n}$. (3). Channel
concatenate ($f_{c}$): where we concatenate $\mathbf{F_{res}}_{n}$ with
$\mathbf{F_{NUTA}}$ along channel dimension. Our results show that $f_{c}$
gives best performance at lowest computational cost. Note that the initial
$\mathbf{F_{NUTA}}$ is generated by spatially down-sampling the uniform branch
feature as shows in Figure 2.

###  3.5 Implementation Details

We initialize the uniform branch with Res2D or I3D backbone networks using
ImageNet pre-trained weights. Our experiments show that the training process
converges faster if the uniform branch is initialized with Kinetics 400 pre-
trained weights. For the other components, we used kaiming initialization
[12]. The model is trained on 8 machines with 8 GPUs on each machine (64 Tesla
V100 GPUs in total). The batch size is set to 8 per GPU without using the
synchronized batch normalization. If the batch size is smaller than 8 per GPU,
we train the model with synchronized batch normalization.

All models are trained for 180 epochs. We use an initial learning rate of 0.4
(equivalent to 0.05 on a single instance) and weights decay of 1e-4. Learning
rate is dropped by scale of 10 at epoch 60, 120 and 150. To avoid over-
fitting, we apply the commonly used data augmentation techniques including:
random resized crop (i.e., with short side randomly resized to [256, 320] and
random crops of 224 following [30, 36]), random temporal sampling (i.e., we
keep the original video frame rate, when a video has less frames than needed,
we pad the missing frames with last frame) and random horizontal flipping.
Besides data augmentation, we also use dropout [13] before the final linear
layer with a dropout ratio of $0.6$.

##  4 Experimental Results

###  4.1 Dataset

We test our method on four of the most widely used datasets:
Kinetics 400 [4] consists of approximately 240k training and 20k validation
videos trimmed to 10 seconds from 400 human action categories. Similarly to
other works, we report top-1 classification accuracy on the validation set.
Kinetics 700 [3] contains approximately 650k video clips that covers 700 human
action classes. Each action class has at least 600 video clips. Each clip is
human annotated with a single action class and lasts around 10s. We report
top-1 classification accuracy on the validation set.
Something-Something V2 [11] dataset consists of 174 actions and contains
220,847 videos. Following other works, we report top-1 classification accuracy
on the validation set. Something-Something dataset is fairly different from
the other three datasets as it does not depict people performing common
activities but instead manipulations of one or two objects such as opening
something, covering something with something, moving something behind
something. This dataset in particular requires strong temporal modeling as
many activities cannot be simply inferred based on spatial features.
Charades [28] has about 9.8k training videos and 1.8k validation videos in 157
classes in a multi-label classification setting with longer activities
spanning about 30 seconds on average. Performance is measured in mean Average
Precision (mAP).

Table 1: Results on Kinetics-400 dataset. We report top 1(%) on the validation set. I3D* stands for the I3D network with temporal down sampling [36]. The ‘Input’ column indicates what frames of the 64 frame clip are actually sent to the network. $n\times\tau$ input indicates we feed $n$ frames to the network sampled every $\tau$ frames. Model | Input | GFLOPs | Top1
---|---|---|---
R2D50 [36] | $32\times 2$ | - | 69.9
R2D101 [36] | $32\times 2$ | - | 71.3
R2D50-NL [36] | $32\times 2$ | - | 73.5
TSM [22] | 8 | 33 | 74.1
NUTA R2D-50 | $32\times 2$ | 54 | 74.2
NUTA R2D-101 | $32\times 2$ | 66 | 75.5
I3D50* [2] | $32\times 2$ | 33 | 74.0
I3D101* [2] | $32\times 2$ | 58 | 77.4
I3D50 [42] | $32\times 2$ | 168 | 75.7
NL50 [36] | $32\times 2$ | 282 | 76.5
NL101 [36] | $32\times 2$ | 544 | 77.7
TEA50 [20] | $16\times 2$ | 70 | 76.1
CIDC [18] | $32\times 2$ | 101 | 75.5
FG 152 [26] | $32\times 2$ | - | 78.8
SF50 [8] | $32\times 2$ | 66 | 77.0
SF101 [8] | $32\times 2$ | 106 | 77.5
TPN-50 [42] | $32\times 2$ | 256 | 77.7
TPN-101 [42] | $32\times 2$ | 458 | 78.9
X3D-L [7] | $16\times 5$ | 24.8 | 77.5
NUTA network I3D-50* | $32\times 2$ | 50 | 76.2
NUTA network I3D-50 | $16\times 4$ | 98 | 76.8
NUTA network I3D-50 | $32\times 2$ | 197 | 77.2
NUTA network I3D-101* | $32\times 2$ | 111 | 77.6
NUTA network I3D-101 | $16\times 4$ | 186 | 78.3
NUTA network I3D-101 | $32\times 2$ | 372 | 78.9

###  4.2 Comparison to State-of-the-art

####  4.2.1 Kinetics 400

Following previous works, at inference we run 30-view (10 temporal crops and 3
spatial crops) through our network and average their predictions to procure
the final classification score. We report results on the validation set of
Kinetics 400 in Table 1. We report the top 1 accuracy and GFLOPs (Giga
Floating-Point Operations) required to compute results on one view.

Our method achieves comparable performance to state-of-the-art methods with
fewer FLOPs and notably outperforms other methods that have a similar goal of
answering where to focus the feature extraction across the video clip. In
comparison to the non-local model that emphasizes long-term temporal modeling
[36], our network achieves 1.2% higher accuracy with roughly half of the
computation. Looking at our I3D-50 variant, we outperform the most recent TEA
network [20] by 0.7% at roughly same FLOPs using 16 frames input. Compared to
state-of-the-art networks [8, 42], our model achieves similar performance with
30% less FLOPs compared with TPN [42]. The results demonstrate that dropping
non-informative frames is an efficient way for temporal feature aggregation in
comparison to fusing features at multiple different scales. The proposed NUTA
network boosts the performance of an I3D network by 3.2% but only increases
computation by about 8%, which demonstrates the proposed network architecture
is computationally efficient.

NUTA is generalizable to different backbones. From the table we show that NUTA
network works well with both the I3D variant that includes temporal down-
sampling [36] and the one that has no temporal down-sampling [42]. Of the
methods in Tabel 1, X3D [7], which is based on neural architecture search,
arguably has the best trade-off between performance and cost. Our NUTA network
is fully compatible with it and other 3D convolution based backbones as well.

The NUTA module also works with 2D backbones. The results in the top rows of
Table 1 show that our model is able to outperform the recently proposed 2D TSM
[22] network with the similar computational cost, which demonstrates that the
non-uniform temporal information aggregation is effective on both 2D and 3D
feature extraction pipelines.

####  4.2.2 Kinetics 700

Kinetics 700 [3] is the latest Kinetics dataset for video classification. We
follow the same 30-view evaluation protocol as used in Kinetics 400 and report
the FLOPs and top1 accuracy.

Table 2: Results on Kinetics-700 dataset. We report top 1 accuracy (%) on validation set. K4&K6 stand for Kinetics 400 and Kinetics 600. SF-NL 101 stands for slowfast network with non-local block [8]. $n\times\tau$ indicates we feed $n$ frames to the network sampled every $\tau$ frames.  Model | pre-train | FLOPs | Top1
---|---|---|---
I3D50 [3] | N/A | N/A | 58.7
SF101-NL $8\times 8$ | K4&K6 | 115G | 70.6
SF101-NL $16\times 8$ | K4&K6 | 234G | 71.0
NUTA network I3D-50 ($32\times 2$) | K400 | 197G | 68.9
NUTA network I3D-101 ($16\times 4$) | K400 | 186G | 69.5
NUTA network I3D-101 ($32\times 2$) | K400 | 372G | 70.3

Our experiments show a consistent performance trend on Kinetics 700 dataset.
As listed in Table 2, our model achieves comparable performance to slowfast
[8], and significantly outperforms the baseline I3D by 10%. Note that the
slowfast model is pre-trained on both Kinetics 400 and 600 while our model is
only pre-trained on Kinetics 400, which may explain the small performance gap
between our model and slowfast with non-local [8].

####  4.2.3 Something-Something V2

Something-something dataset is unique in that it shows people manipulating
objects, rather than performing common actions. We believe that distinguishing
those manipulation actions needs more motion context than other datasets so
it’s a good dataset to validate the effectiveness of temporal modelling.
Following previous works [22]. during inference, we use the center crop of
size $224\times 224$ from 8 segments to compute the classification score on
the validation set.

Table 3: Results on Something-something V2 validation set. The results are generated by taking the center crop of 1 clip/video [22] as input to the network. Model | Init. | FLOPs | Top1
---|---|---|---
I3D50 | K400 | 33G | 50.0
TRN [44] | ImgNet | 42G | 55.5
CIDC [18] | K400 | 92G | 56.1
TSM [22] | K400 | 33G | 59.1
SF [39] | K400 | 66G | 60.9
SF (multigrid) [39] | K400 | 66G | 61.2
TPN [42] | K400 | 56G | 62.0
NUTA network I3D-50 | K400 | 49G | 61.5
NUTA network I3D-101 | K400 | 98G | 62.1
NUTA network I3D-101 | K700 | 98G | 63.0

As the results in Table 3 show, our model achieves state-of-the-art
recognition performance (single model with RGB input) on something-something
v2 dataset with 63.0% top1 accuracy. Our NUTA network is able to outperform
I3D and I3D based approaches [39, 18] with lower computational cost. It is
worth mentioning that our NUTA network generalizes well to different datasets.
Previous works have demonstrated that I3D is not as effective in modeling the
temporal relations [22, 18]. The state-of-the-art methods on Kinetics e.g. TPN
[42] and slowfast [42] must either use the TSM [22] backbone or the multigrid
training trick to achieve high performance on the Something-something dataset,
showing that for Something-something dataset a specialize adaptation is
usually required to transfer a high performing Kinetics model over to the
task. We apply our NUTA method directly to the unmodified I3D backbone and are
able to get a 11.5% performance boost over the I3D baseline. We believe this
results demonstrate that our selection of informative features is not only
useful in improving performance, it generalizes more easily to different
domains of video understanding.

####  4.2.4 Charades

Charades is a dataset [28] of longer duration video annotated for a multi-
label action classification problem. Table 4 shows our model achieves very
competitive recognition accuracy that is comparable to state-of-the-art
methods while at significantly lower computational cost. Compared to both
slowfast-101 and slowfast-101 with non-local [8], we achieve slightly lower
performance (41.2 mAP) but with 18% fewer FLOPs. We achieve state-of-the-art
performance when we initialize the model with K700 pre-trained weights, but
still at much lower computational cost comparing with previous methods [38,
8].

Table 4: Results on Charades dataset. The results are generated from an input of 32 frames sampled 1 out of 4 frames. Model | Init. | FLOPs | mAP
---|---|---|---
Non-local [36] | K400 | 544G | 37.5
STRG [37] | K400 | N/A | 39.7
LFB-NL [38] | K400 | 529G | 42.5
SF [8] | K400 | 213G | 42.1
SF-NL [8] | K400 | 234G | 42.5
X3D-XL [7] | K400 | 48.4G | 43.4
NUTA network 101 | K400 | 186G | 41.2
NUTA network 101 | K700 | 186G | 43.1
Table 5:  Ablation studies on Kinetics 400 dataset. We use an I3D-50 backbone.
I3D and I3D* stand for the I3D network without and with temporal down sampling
respectively. The evaluation is performed on 30 views with 32 frame input
unless specified.

Model | GFLOPs | Top1
---|---|---
I3D* [36] | 33 | 74.0
I3D [42] | 168 | 75.7
NUTA network I3D | 197 | 77.2
(a) Comparison between I3D with and without temporal down-sampling, and our
NUTA network

Model | Top1
---|---
I3D (backbone) | 75.7
\+ NUTA (only I3D features) | 76.3
\+ I3D and NUTA features | 77.2
(b) NUTA network performance with different feature extraction strategies

Input | GFLOPs | Top1
---|---|---
$8\times 8$ | 49 | 75.5
$16\times 4$ | 98 | 76.8
$32\times 2$ | 197 | 77.2
(c) Model performance with different input clip lengths

Input | GFLOPs | Top1
---|---|---
$128^{2}$ | 52 | 76.3
$256^{2}$ | 197 | 77.2
$312^{2}$ | 395 | 77.6
(d) Model performance with different input frame resolution

Backbone | Number of groups | Top1
---|---|---
I3D* | 16 | 75.6
I3D* | 32 | 76.0
I3D* | 64 | 76.2
(e) Model performance with different NUTA group head configurations

Backbone | Add NTUA at | Top1
---|---|---
I3D | NUTA3+NUTA4+NUTA5 | 76.7
I3D | NUTA4+NUTA5 | 77.2
I3D | NUTA5 | 75.7
(f) Performance comparison of including different NUTA stages

###  4.3 Ablation Study

####  4.3.1 NUTA vs. uniform temporal aggregation

We compare the performance of an I3D model with temporal down-sampling
(temporal pooling and striding ) [36], an I3D without temporal down-sampling
[42] and an I3D with our proposed NUTA network (Table 5f(a)). The results show
that the I3D without temporal down-sampling outperforms the I3D with temporal
downsampling. We believe this is because it is able to maintain more temporal
information. Our NUTA network is able to further outperform I3D without
temporal down-sampling by an additional 1.5%, which shows that non-uniformly
aggregating the informative temporal information helps with the action
understanding. Note that the NUTA models only add roughly 10% additional GLOPs
since the NUTA modules perform their own temporal down-sampling, significantly
reducing its computational overhead.

####  4.3.2 NUTA network branch analysis

To understand the feature from both the uniform and non-uniform branches we
look at different configurations of our network using the I3D backbone as a
starting point. First we take the baseline by using the I3D backbone as the
uniform branch without adopting an non-uniform branch. Then we construct a
NUTA network by attaching our NUTA modules to form a non-uniform branch but
still perform classification from only the features coming out of the uniform
branch (I3D features). Such setup outperforms the baseline by 0.7% (Table
5f(b)), and we conjecture this is because the features from uniform branch
utilize the temporal mapping learned by the non-uniform branch from NUTA
module to drop non-informative frames. Adding the features from non-uniform
branch classification further improves the Top1 accuracy by 0.9%, showing that
those two branches learn features that are complementary for better action
recognition.

####  4.3.3 Frame rate

We test the NUTA network at different frame-rates following previous papers
[8, 36, 42]. Our results (Table 5f(c)) shows that dense sampling gives a
better performance but requires significantly higher computational cost. Our
method particularly benefits from more densely-sampled frames as this prevents
information loss at the early stages of the network, thus giving our NUTA
units sufficient information for temporal feature aggregation.

####  4.3.4 Frame resolution

Table 5f(d) compares our NUTA network performance across different input
resolutions. Generally a larger input resolution leads to a better
performance, but will significantly increase the computational cost. We used
$256^{2}$ for all other experiments and GFLOPs estimation.

####  4.3.5 NUTA configuration

We study the impact of number of groups used in the 3D convolutions (equation
1) in our NUTA module. Our results (Table 5f(e)) show that 64 groups give the
best performance.

####  4.3.6 Number of NUTA units

Table 5f(f) compares including (or not including) different NUTA stages. The
results show that adding more NUTA stages does not neccessarily lead to better
performance. Adding NUTA stage 4 and 5 gives the best performance (Table
5f(f)). A possible explanation is that the NUTA module aggregates features by
focusing on important temporal instances but stacking too many NUTA modules
will compress feature too aggressively (since each NUTA will down-sample
temporal dimension by 2).

##  5 Visualization

Figure 3: Visualization of the temporal projection map learned by our NUTA
module. Our NUTA module performs uniform aggregation when information is
smoothly distributed across clip (top clip) but is able to focus on
informative temporal segments when there is irrelevant segments in the video
(middle and bottom clips).

To understand how the proposed NUTA unit works to aggregate temporal features,
we visualize the temporal projection matrix (equation 1). We notice that: when
the input clip has smooth transition, which indicates the information is
roughly uniformly distributed over time, the NUTA unit performs uniform
sampling (Figure 3 top row). Most previous works based on 3D convolutions also
perform temporal information aggregation in this way. When the frames are
highly repetitive or contain only background without useful information for
action recognition, the NUTA unit is able to skip the non-informative frames
and focus more on the representative features (e.g. Figure 3 middle row,
information from first half of the clip is skipped). When the input features
have scene changes or contain noises (e.g. transition frames), the temporal
mapping generated by NUTA stays to focus on representative information by
skipping noises (e.g. Figure 3 bottom row, the NUTA is able to focus on the
start and end of the video while skipping the middle frames).

##  6 Conclusion

In this paper, we have addressed the “where to extract feature” question by
proposing a non-uniform temporal aggregation (NUTA) module, which is able to
select the informative features. We have further proposed the two-branch
network with a uniform branch that learns local feature and non-uniform branch
that learns clip-level features. Our experimental results have showed that the
proposed NUTA network achieves state-of-the-art accuracy on four public action
recognition datasets. Besides, we have demonstrated how our NUTA network works
by visualizing the intermediate temporal projection matrices. One future
direction is to apply the proposed NUTA to other tasks (e.g. action
detection).

## References

  * [1] Hakan Bilen, Basura Fernando, Efstratios Gavves, Andrea Vedaldi, and Stephen Gould.  Dynamic image networks for action recognition.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 3034–3042, 2016.
  * [2] Joao Carreira.  Quo vadis, action recognition.  A new model and the kinetics dataset. CoRR, abs/1705.07750, 2:3, 2017.
  * [3] Joao Carreira, Eric Noland, Chloe Hillier, and Andrew Zisserman.  A short note on the kinetics-700 human action dataset.  arXiv preprint arXiv:1907.06987, 2019.
  * [4] J. Carreira and A. Zisserman.  Quo vadis, action recognition? a new model and the kinetics dataset.  In 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 4724–4733, July 2017.
  * [5] Jeffrey Donahue, Lisa Anne Hendricks, Sergio Guadarrama, Marcus Rohrbach, Subhashini Venugopalan, Kate Saenko, and Trevor Darrell.  Long-term recurrent convolutional networks for visual recognition and description.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2625–2634, 2015.
  * [6] Wenbin Du, Yali Wang, and Yu Qiao.  Recurrent spatial-temporal attention network for action recognition in videos.  IEEE Transactions on Image Processing, 27(3):1347–1360, 2017.
  * [7] Christoph Feichtenhofer.  X3d: Expanding architectures for efficient video recognition.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 203–213, 2020.
  * [8] Christoph Feichtenhofer, Haoqi Fan, Jitendra Malik, and Kaiming He.  Slowfast networks for video recognition.  In Proceedings of the IEEE International Conference on Computer Vision, pages 6202–6211, 2019.
  * [9] Christoph Feichtenhofer, Axel Pinz, and Richard P Wildes.  Spatiotemporal multiplier networks for video action recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4768–4777, 2017.
  * [10] Basura Fernando, Efstratios Gavves, José Oramas, Amir Ghodrati, and Tinne Tuytelaars.  Rank pooling for action recognition.  IEEE transactions on pattern analysis and machine intelligence, 39(4):773–787, 2016.
  * [11] Raghav Goyal, Samira Ebrahimi Kahou, Vincent Michalski, Joanna Materzynska, Susanne Westphal, Heuna Kim, Valentin Haenel, Ingo Fruend, Peter Yianilos, Moritz Mueller-Freitag, et al.  The” something something” video database for learning and evaluating visual common sense.  In ICCV, volume 1, page 3, 2017.
  * [12] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [13] Geoffrey E Hinton, Nitish Srivastava, Alex Krizhevsky, Ilya Sutskever, and Ruslan R Salakhutdinov.  Improving neural networks by preventing co-adaptation of feature detectors.  arXiv preprint arXiv:1207.0580, 2012.
  * [14] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural computation, 9(8):1735–1780, 1997.
  * [15] Yueyu Hu, Chunhui Liu, Yanghao Li, and Jiaying Liu.  Temporal perceptive network for skeleton-based action recognition.  In BMVC.
  * [16] Andrej Karpathy, George Toderici, Sanketh Shetty, Thomas Leung, Rahul Sukthankar, and Li Fei-Fei.  Large-scale video classification with convolutional neural networks.  In Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, pages 1725–1732, 2014.
  * [17] Bruno Korbar, Du Tran, and Lorenzo Torresani.  Scsampler: Sampling salient clips from video for efficient action recognition.  In Proceedings of the IEEE International Conference on Computer Vision, pages 6232–6242, 2019.
  * [18] Xinyu Li, Bing Shuai, and Joseph Tighe.  Directional temporal modeling for action recognition.  arXiv preprint arXiv:2007.11040, 2020.
  * [19] Xinyu Li, Yanyi Zhang, Jianyu Zhang, Moliang Zhou, Shuhong Chen, Yue Gu, Yueyang Chen, Ivan Marsic, Richard A Farneth, and Randall S Burd.  Progress estimation and phase detection for sequential processes.  Proceedings of the ACM on interactive, mobile, wearable and ubiquitous technologies, 1(3):73, 2017.
  * [20] Yan Li, Bin Ji, Xintian Shi, Jianguo Zhang, Bin Kang, and Limin Wang.  Tea: Temporal excitation and aggregation for action recognition.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 909–918, 2020.
  * [21] Zhenyang Li, Kirill Gavrilyuk, Efstratios Gavves, Mihir Jain, and Cees GM Snoek.  Videolstm convolves, attends and flows for action recognition.  Computer Vision and Image Understanding, 166:41–50, 2018.
  * [22] Ji Lin, Chuang Gan, and Song Han.  Tsm: Temporal shift module for efficient video understanding.  In Proceedings of the IEEE International Conference on Computer Vision, 2019.
  * [23] Chunhui Liu, Yanghao Li, Yueyu Hu, and Jiaying Liu.  Online action detection and forecast via multitask deep recurrent neural networks.  In 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 1702–1706. IEEE, 2017.
  * [24] Zhaoyang Liu, Limin Wang, Wayne Wu, Chen Qian, and Tong Lu.  Tam: Temporal adaptive module for video recognition.  arXiv preprint arXiv:2005.06803, 2020.
  * [25] Brais Martinez, Davide Modolo, Yuanjun Xiong, and Joseph Tighe.  Action recognition with spatial-temporal discriminative filter banks.  In Proceedings of the IEEE International Conference on Computer Vision, pages 5482–5491, 2019.
  * [26] Brais Martinez, Davide Modolo, Yuanjun Xiong, and Joseph Tighe.  Action recognition with spatial-temporal discriminative filter banks.  In Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV), October 2019.
  * [27] Zhaofan Qiu, Ting Yao, and Tao Mei.  Learning spatio-temporal representation with pseudo-3d residual networks.  In proceedings of the IEEE International Conference on Computer Vision, pages 5533–5541, 2017.
  * [28] Gunnar A. Sigurdsson, Gül Varol, Xiaolong Wang, Ivan Laptev, Ali Farhadi, and Abhinav Gupta.  Hollywood in homes: Crowdsourcing data collection for activity understanding.  ArXiv e-prints, 2016.
  * [29] Karen Simonyan and Andrew Zisserman.  Two-stream convolutional networks for action recognition in videos.  In Advances in neural information processing systems, pages 568–576, 2014.
  * [30] Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  In International Conference on Learning Representations, 2015.
  * [31] Sijie Song, Cuiling Lan, Junliang Xing, Wenjun Zeng, and Jiaying Liu.  Spatio-temporal attention-based lstm networks for 3d action recognition and detection.  IEEE Transactions on image processing, 27(7):3459–3471, 2018.
  * [32] Graham W Taylor, Rob Fergus, Yann LeCun, and Christoph Bregler.  Convolutional learning of spatio-temporal features.  In European conference on computer vision, pages 140–153. Springer, 2010.
  * [33] Du Tran, Lubomir Bourdev, Rob Fergus, Lorenzo Torresani, and Manohar Paluri.  Learning spatiotemporal features with 3d convolutional networks.  In Proceedings of the IEEE international conference on computer vision, pages 4489–4497, 2015.
  * [34] Du Tran, Heng Wang, Lorenzo Torresani, Jamie Ray, Yann LeCun, and Manohar Paluri.  A closer look at spatiotemporal convolutions for action recognition.  In Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, pages 6450–6459, 2018.
  * [35] Limin Wang, Yuanjun Xiong, Zhe Wang, Yu Qiao, Dahua Lin, Xiaoou Tang, and Luc Van Gool.  Temporal segment networks: Towards good practices for deep action recognition.  In European conference on computer vision, pages 20–36. Springer, 2016.
  * [36] Xiaolong Wang, Ross Girshick, Abhinav Gupta, and Kaiming He.  Non-local neural networks.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 7794–7803, 2018.
  * [37] Xiaolong Wang and Abhinav Gupta.  Videos as space-time region graphs.  In Proceedings of the European conference on computer vision (ECCV), pages 399–417, 2018.
  * [38] Chao-Yuan Wu, Christoph Feichtenhofer, Haoqi Fan, Kaiming He, Philipp Krahenbuhl, and Ross Girshick.  Long-term feature banks for detailed video understanding.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 284–293, 2019.
  * [39] Chao-Yuan Wu, Ross Girshick, Kaiming He, Christoph Feichtenhofer, and Philipp Krahenbuhl.  A multigrid method for efficiently training video models.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 153–162, 2020.
  * [40] Wenhao Wu, Dongliang He, Xiao Tan, Shifeng Chen, and Shilei Wen.  Multi-agent reinforcement learning based frame sampling for effective untrimmed video recognition.  In Proceedings of the IEEE International Conference on Computer Vision, pages 6222–6231, 2019.
  * [41] Saining Xie, Chen Sun, Jonathan Huang, Zhuowen Tu, and Kevin Murphy.  Rethinking spatiotemporal feature learning: Speed-accuracy trade-offs in video classification.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 305–321, 2018.
  * [42] Ceyuan Yang, Yinghao Xu, Jianping Shi, Bo Dai, and Bolei Zhou.  Temporal pyramid network for action recognition.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 591–600, 2020.
  * [43] Joe Yue-Hei Ng, Matthew Hausknecht, Sudheendra Vijayanarasimhan, Oriol Vinyals, Rajat Monga, and George Toderici.  Beyond short snippets: Deep networks for video classification.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4694–4702, 2015.
  * [44] Bolei Zhou, Alex Andonian, Aude Oliva, and Antonio Torralba.  Temporal relational reasoning in videos.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 803–818, 2018.


# Finite temperature tensor network algorithm for frustrated two-dimensional
quantum materials

Philipp Schmoll  Dahlem Center for Complex Quantum Systems and Institut für
Theoretische Physik, Freie Universität Berlin, Arnimallee 14, 14195 Berlin,
Germany  Christian Balz  ISIS Neutron and Muon Source, Rutherford Appleton
Laboratory, Didcot OX11 0QX, UK  Bella Lake  Helmholtz-Zentrum Berlin für
Materialien und Energie, Hahn-Meitner-Platz 1, 14109 Berlin, Germany  Institut
für Festkörperphysik, Technische Universität Berlin, Hardenbergstraße 36,
D-10623 Berlin, Germany Jens Eisert  Dahlem Center for Complex Quantum Systems
and Institut für Theoretische Physik, Freie Universität Berlin, Arnimallee 14,
14195 Berlin, Germany  Helmholtz-Zentrum Berlin für Materialien und Energie,
Hahn-Meitner-Platz 1, 14109 Berlin, Germany  Augustine Kshetrimayum  Dahlem
Center for Complex Quantum Systems and Institut für Theoretische Physik, Freie
Universität Berlin, Arnimallee 14, 14195 Berlin, Germany  Helmholtz-Zentrum
Berlin für Materialien und Energie, Hahn-Meitner-Platz 1, 14109 Berlin,
Germany  Theory Division, Saha Institute of Nuclear Physics, 1/AF Bidhannagar,
Kolkata 700 064, India

###### Abstract

Aimed at a more realistic classical description of natural quantum systems, we
present a two-dimensional tensor network algorithm to study finite temperature
properties of frustrated model quantum systems and real quantum materials. For
this purpose, we introduce the infinite projected entangled simplex operator
ansatz to study thermodynamic properties. To obtain state-of-the-art
benchmarking results, we explore the highly challenging spin-1/2 Heisenberg
anti-ferromagnet on the Kagome lattice, a system for which we investigate the
melting of the magnetization plateaus at finite magnetic field and
temperature. Making close connection to actual experimental data of real
quantum materials, we go on to studying the finite temperature properties of
Ca10Cr7O28. We compare the magnetization curve of this material in the
presence of an external magnetic field at finite temperature with classically
simulated data. As a first theoretical tool that incorporates both thermal
fluctuations as well as quantum correlations in the study of this material,
our work contributes to settling the existing controversy between the
experimental data and previous theoretical works on the magnetization process.

##  I Introduction

Simulating complex quantum materials is considered to be one of the hardest
problems in modern physics. Density functional theory is arguably the most
popular approach to date for calculating the electronic structure of molecules
and extended materials Kohn and Sham (1965); Tkatchenko and Scheffler (2009).
In situations in which strong correlations are expected to be dominant,
however, its applicability can be limited. Ultimately, the core computational
challenge in the numerical simulation of strongly correlated quantum materials
arises from the exponential scaling of the size of the Hilbert space with the
system size. Thus, it comes as no surprise that the _exact diagonalization_
(ED) technique can only study small sizes and therefore may fail to capture
the important physics of emergent many-body phenomena. Mean-field techniques
are also unsuitable in the study of quantum materials as they neglect the most
crucial ingredient in describing these systems: quantum entanglement. While
_quantum Monte Carlo_ constitutes a versatile tool for simulating unfrustrated
strongly correlated systems Becca and Sorella (2017), they suffer from severe
limitations for frustrated quantum systems due to the sign problem. In this
respect, _tensor network_ techniques have emerged as a powerful alternative
for studying challenging many-body problems which does not suffer from any of
those limitations Orús (2014); Verstraete  _et al._ (2008); Bridgeman and
Chubb (2017); Eisert  _et al._ (2010).

The success of one-dimensional tensor networks, also known as _matrix product
states_ (MPS) White (1992); Schollwöck (2005, 2011); Perez-Garcia  _et al._
(2007), in describing one-dimensional phases of matter have provided much
impetus to the development of two-dimensional tensor network algorithms. While
the situation is much more intricate and challenging in two spatial
dimensions, such tensor network algorithms, also known as _projected entangled
pair states_ (PEPS) or _iPEPS_ Verstraete and Cirac (2004); Jordan  _et al._
(2008); Verstraete  _et al._ (2008) in its infinite instance tackling directly
the thermodynamic limit, have recently matured and have been employed
successfully to study various challenging problems in two dimensions. This
includes finding ground states of frustrated systems and real quantum
materials Liao  _et al._ (2017); Picot  _et al._ (2016a); Picot and Poilblanc
(2015); Kshetrimayum  _et al._ (2016); Boos  _et al._ (2019); Kshetrimayum
_et al._ (2020a); Astrakhantsev  _et al._ (2021) and non-equilibrium systems
Kshetrimayum  _et al._ (2017); Czarnik  _et al._ (2019a); Hubig and Cirac
(2019); Kshetrimayum  _et al._ (2020b, 2021); Dziarmaga (2021, 2022). While
most of the efforts has been dedicated towards identifying ground states of
closed quantum systems, in order to accurately capture the physics of quantum
materials in realistic conditions in the lab, one needs to include the effects
of temperature. With this aim, there have been several recent works on two-
dimensional _finite temperature tensor network algorithms_ Czarnik  _et al._
(2012); Czarnik and Dziarmaga (2015); Kshetrimayum  _et al._ (2019); Czarnik
_et al._ (2019b, 2021); Mondal  _et al._ (2020). Most of these works have,
however, focused on paradigmatic, theoretical models such as the Ising, Kitaev
or Heisenberg models, and mostly models that are defined on the square
lattice.

In this work, we develop a two-dimensional tensor network algorithm for
studying finite temperature properties of existing quantum materials, thus
mimicking experimental studies as closely as possible. We start by describing
our method and then present results on two important instances of strongly
correlated systems: (i) the paradigmatic spin-$1/2$ Kagome Heisenberg anti-
ferromagnet both in the absence and presence of an external magnetic field and
(ii) the real quantum material Ca10Cr7O28 that features a bilayer Kagome
structure.

##  II Method

Our method substantially advances the algorithm proposed in Ref. Kshetrimayum
_et al._ (2019) by extending it to the more challenging realm of frustrated
quantum systems and real quantum materials. This step renders it possible to
directly compare experimental data and theoretical tensor network simulations,
as we do here. We will now review the underlying annealing algorithm and
highlight the improvements. In order to simulate a quantum system at finite
temperature $\beta\coloneqq 1/T>0$, we describe it by an (unnormalized)
thermal quantum state

$\displaystyle\rho(\beta)=\mathrm{e}^{-\beta H},$ (1)

where $H$ is the full local many-body Hamiltonian. To obtain such a Gibbs
state, we start from an infinite temperature state, i.e., $\rho(\beta=0)$ and
cool down the system to the desired temperature $\beta^{-1}>0$. The initial
state is simply a tensor product of identities, the (unnormalized) single-
particle thermal state in the limit $T\rightarrow\infty$. The evolution to the
desired temperature can be generated by suitably many small temperature steps
$\delta\beta$, so that the full quantum state is obtained for $N\in\mathbb{N}$
by

$\displaystyle\rho(\beta)=\rho(\delta\beta)^{N}=\left(\mathrm{e}^{-\delta\beta
H}\right)^{N}$ (2)

with $\delta\beta\coloneqq\beta/N$ and $\rho(\delta\beta)$ as what we call the
_infinitesimal thermal density matrix_ (ITDM). This cooling is implemented by
a simple update technique Jiang  _et al._ (2008a); Kshetrimayum  _et al._
(2019). The simple update is adopted here for its numerical stability and
efficiency Liao  _et al._ (2017); Picot and Poilblanc (2015); Kshetrimayum
_et al._ (2019); Gauthé and Mila (2022), particularly relevant while working
on systems with large physical dimensions, which seems a necessity for the
demanding task considered here (see the Appendix).

Instead of directly cooling down to $\beta>0$, it is advantageous to cool down
to $\beta/2$ and evaluate the Gibbs state as

$\displaystyle\rho(\beta)=\rho(\beta/2)^{\dagger}\rho(\beta/2)\ .$ (3)

This ensures that the resulting operator is _positive semi-definite_ and hence
reflects a valid quantum state, which is otherwise not guaranteed in tensor
network implementations due to truncation effects Chen  _et al._ (2018a);
Werner  _et al._ (2016); Chen  _et al._ (2018b). Eq. (3) is _the_ main
difference to the underlying algorithm presented in Ref. Kshetrimayum  _et
al._ (2019) and is the crucial improvement which enables the simulation of
frustrated systems (along with using the correct tensor network structure of
the underlying lattice of the model as we discuss in the next paragraph).
Thus, we have the freedom of evolving up to only $N/2$ steps thereby saving a
factor of two in the number of annealing steps or evolving up to $N$ steps
with each step size being $\delta\beta/2$. The latter choice is adopted in our
simulations and decreases the Trotter error from
$\mathcal{O}({\delta\beta}^{2})$ to $\mathcal{O}({\delta\beta}^{2}/4)$.

We will now introduce the tensor network representation of the Gibbs state:
The _infinite projected entangled simplex operator_ (iPESO) shown in Fig. 1.
It is the operator version of the _infinite projected entangled simplex state_
(iPESS) proposed in Ref. Xie  _et al._ (2014), applied to the simulation of
thermal density matrices. In both tensor networks, the quantum correlations
inherently present on the Kagome triangles are efficiently and accurately
captured by exploiting the structure of its dual, the honeycomb lattice.

Figure 1: iPESO ansatz for the simulation of Gibbs states on the Kagome
lattice (shown in light blue).

Green tensors represent the lattice sites of the Kagome lattice, with two
physical indices for the density matrix (as opposed to a quantum state, for
which tensors only have a single physical index). They are connected by purely
virtual simplex tensors shown in grey. The accuracy with which the iPESO
approximates the thermal density matrix is controlled by the bond dimension of
the virtual bulk indices, denoted as $\chi_{B}$. It is important to note that
$\chi_{B}$ needs to be chosen sufficiently large to prevent truncation effects
in the ITDM. This leads to a minimal bond dimension of $p^{2}$, where $p$ is
the dimension of the Hilbert space of the local physical degrees of freedom (a
detailed explanation is given in the Appendix). For the final simulations of
the targeted real material we choose the bond dimension such that the total
truncation error is below $\sim 10^{-5}$, see Fig. 11. Expectation values are
then directly computed in the tensor network representation of the thermal
state according to

$\displaystyle\langle\hat{O}\rangle={N_{\rho}}^{-1}\text{Tr}\bigl{[}\rho(\beta/2)^{\dagger}\,\hat{O}\,\rho(\beta/2)\bigr{]},$
(4)

with a normalization factor
$N_{\rho}\coloneqq\textrm{Tr}[\rho(\beta/2)^{\dagger}\,\rho(\beta/2)]$.
Expectation values can be computed by either using the simple update mean-
field environment, or by a full _corner transfer matrix renormalization group
(CTMRG)_ procedure Nishino and Okunishi (1996); Orús and Vidal (2009); Orús
(2012), which captures quantum correlations more faithfully. For the latter,
the environment bond dimension $\chi_{E}$ controls the approximations in the
contraction of the infinite two-dimensional lattice. Details for both the
simple update cooling and calculations of expectation values are presented in
detail in the Appendix. The smallest possible unit cell of the iPESO consists
of three lattice site tensors and two simplex tensors, as presented in Fig. 1.
Besides this structure, we also employ a nine-site unit cell in our numerical
simulation. This is required to capture thermal states with larger structures
that are not commensurate with three-site translational invariance.

##  III Models and results

###  III.1 Kagome Heisenberg anti-ferromagnet

The first application of the developed iPESO method is the finite temperature
study of the frustrated spin-$1/2$ Heisenberg anti-ferromagnet on the Kagome
lattice, a paradigmatic model that has been a topic of intense study in the
community Sachdev (1992); Ran  _et al._ (2007); Jiang  _et al._ (2008b, 2012);
Liao  _et al._ (2017); Götze  _et al._ (2011); Iqbal  _et al._ (2013). Its
Hamiltonian is given by

$\displaystyle H=J\sum_{\langle
i,j\rangle}\vec{S}_{i}\cdot\vec{S}_{j}-h_{z}\sum_{i}S_{i}^{z},$ (5)

where $\vec{S}_{i}$ are spin-$1/2$ operators on site $i$ and $\langle
i,j\rangle$ denotes nearest-neighbours in the underlying lattice, $h_{z}$ is a
magnetic field applied along the $z$-axis. In the following, we employ the
iPESO method to study the model at $J=1.0$ over a large temperature range,
choosing an infinitesimal temperature step $\delta\beta=10^{-3}$. In the main
panel of Fig. 2 we show the thermal state energy for a three-site iPESO ansatz
at bulk bond dimensions up to $\chi_{B}=10$. These results are computed using
a CTMRG procedure with individual environment bond dimensions $\chi_{E}$ such
that expectation values are well converged.

Figure 2: Thermal state energy for the spin-$1/2$ Heisenberg model on the
Kagome lattice at various bulk bond dimensions $\chi_{B}$, using CTMRG
environments. The dashed line corresponds to the $T=0$ ground state energy of
Ref. Läuchli  _et al._ (2019). (Top inset) Convergence of the thermal state
energy with $\chi_{B}$ at $T=0.01$. (Bottom inset) Specific heat $C$ at
$\chi_{B}=16$, using mean-field environments.

The energy of the thermal state approaches the ground state energy at $T=0$
for low temperatures, as shown in the top inset. This, along with the
vanishing magnetization when approaching the ground state (not shown here),
indicates that the annealing procedure does not get stuck in local minima and
flows towards the correct ground state as would have been obtained using
ground state optimization. While the thermal state energies in Fig. 2 have
been computed with CTMRG environments, we note that the accuracy is not
affected while using the mean-field environment of the simple update.
Therefore, we compute the _heat capacity_ $C\coloneqq\partial U/\partial T$ at
a higher bulk bond dimension $\chi_{B}=16$, using these environments. The
result is shown in the bottom inset and matches previous finite temperature
studies of the model Chen  _et al._ (2018b). In the analysis of the model and
for the targeted real material we therefore restrict to mean-field
environments, since CTMRG calculations are limited to inexpressively small
environment bond dimensions in those cases.

We further use our method to study the effect of temperature on the
magnetization behaviour of the Heisenberg model in Eq. (5). It is known that
four different magnetization plateaus at values
$m_{z}/m_{S}=[1/9,1/3,5/9,7/9]$ of the saturation magnetization $m_{S}=1/2$
appear at $T=0$ upon tuning the magnetic field $h_{z}$ Picot  _et al._
(2016b); Kshetrimayum  _et al._ (2016). We focus our study on the most
prominent $m_{z}/m_{S}=1/3$ magnetic plateau and classically simulate the
Kagome Heisenberg model in a field for various temperatures in the range of
$T\,\in\,[0.001,100]$ using a nine-site iPESO at bond dimension $\chi_{B}=12$.
The results in Fig. 3 are complemented with regular iPESS simulations at $T=0$
and $\chi_{B}=12$, where the magnetization plateaus appear most prominently.

Figure 3: Magnetization of the Kagome Heisenberg anti-ferromagnet for various
temperatures $T$ and magnetic fields $h_{z}$. Additionally, we show the $T=0$
magnetization curve obtained via iPESS simulations in black. (Inset) Melting
of the $m_{z}/m_{S}=1/3$ plateau.

Of all the $T=0$ magnetization plateaus, we find that only the $1/3$ plateau
survives at finite temperature up to $T\sim 2\cdot 10^{-2}$, where $T$ is in
units of the interaction strength $J$. Above this temperature, the plateau
starts melting and disappears. A cross section of Fig. 3 along constant
temperature slices of $1/T=[\infty,100,80,60,40]$ reveals the melting of the
plateau, shown as the inset on the face of the cube. We observe that the
melting is stronger at the low-field end of the plateau, which is in agreement
with a recent exact diagonalization study of the melting of the magnetization
plateaus Schlüter  _et al._ (2022). Our results serve as an important guide to
the experimental study of the magnetization process of closely related real
materials such as Herbertsmithite ZnCu3(OH)6Cl2 and its relatives.

###  III.2 Real material Ca10Cr7O28

The material Ca10Cr7O28 has a breathing bilayer Kagome structure, with
alternating ferro- and anti-ferromagnetic Heisenberg interactions on
neighbouring triangles, defined by

$\displaystyle\begin{split}H&=\sum_{k=1}^{2}\left[J_{\bigtriangledown
k}\sum_{\langle i,j\rangle}\vec{S}_{i}\cdot\vec{S}_{j}+J_{\bigtriangleup
k}\sum_{\langle i,j\rangle}\vec{S}_{i}\cdot\vec{S}_{j}\right]\\\
&+J_{\text{inter}}\sum_{\langle
i,j\rangle}\vec{S}_{i}\cdot\vec{S}_{j}-h_{z}\sum_{i}S_{i}^{z}\end{split},$ (6)

where $J_{\bigtriangledown k}$ and $J_{\bigtriangleup k}$ are the intra-Kagome
couplings in the two layers $k=1$ and $k=2$ while $J_{\text{inter}}$ denotes
the coupling between the two layers. The lattice structure is shown in Fig. 4.

Figure 4: Breathing bilayer Kagome lattice with different couplings for the
compound material Ca10Cr7O28. The double bonds denote inter-layer coupling
with coupling strength $J_{\text{inter}}$.

The different coupling parameters have been determined from neutron scattering
experiments of the real material in Ref. Balz  _et al._ (2017). These values
(in meV) are $J_{\bigtriangledown 1}=+0.09(2)$, $J_{\bigtriangleup
1}=-0.27(3)$, $J_{\bigtriangledown 2}=-0.76(5)$, $J_{\bigtriangleup
2}=+0.11(3)$ and $J_{\text{inter}}=-0.08(4)$, where the numbers in round
brackets indicate the uncertainties.

The bilayer Kagome structure can be mapped to a single layer by combining the
two spins in the different layers to a single physical site, so that a regular
iPESO ansatz with an enlarged local physical dimension of $d=4$ can be used.
The large physical dimension of the system also means that the minimum bond
dimension required for accurate simulations of the material needs to be
sufficiently high ($\chi_{B}\geq 16$). This becomes a bottleneck in computing
expectation values using CTMRG routines. For this reason, we have adopted the
mean-field environment calculation which takes into account the quantum
correlations within a cluster. We find that even with such approximations in
computing the magnetization and heat capacity, our results are compatible with
the experimental data and provide novel insights into the earlier discrepancy
between theory and experiment data as we show below.

####  III.2.1 Magnetization behaviour

The magnetic properties of this real material have been investigated
previously at zero temperature using tensor networks and compared to
experimental measurements by some of the current authors Kshetrimayum  _et
al._ (2020a). While results at small magnetic fields were in good agreement,
the magnetization curve has shown a significant discrepancy between simulation
and experiment at large values of the field. Such a discrepancy has also been
observed when comparing the experimental data to theoretical mean-field
calculations Balz  _et al._ (2017). Experimental data show that the
magnetization of this material increases rapidly for small external magnetic
fields up to $1$ or $2\text{\,}\mathrm{T}$ above which the slope flattens and
saturation is achieved for a field value of $\approx$12\text{\,}\mathrm{T}$$.
In contrast, our previous tensor network simulation predicted saturation at a
much smaller value of the external magnetic field of
$\approx$1\text{\,}\mathrm{T}$$. This theoretical investigation, while
quantum, has ignored thermal fluctuations. For comparison, molecular dynamics
simulation which take temperature effects into account result in a better
agreement Pohle  _et al._ (2021), however the quantum fluctuations are ignored
by this method. Our present technique encompasses both the quantum properties
of the material as well as the effect of finite temperature.

We now investigate the quantum material Ca10Cr7O28 in the presence of magnetic
fields between $0\text{\,}\mathrm{T}$ and $12\text{\,}\mathrm{T}$ at various
temperatures. The TN simulations are done with a nine-site iPESO at bond
dimension $\chi_{B}=30$ and with $\delta\beta=10^{-2}$, results are shown in
Fig. 5.

Figure 5: Magnetization curve of Ca10Cr7O28 for the Hamiltonian given in Eq.
(6) and magnetic fields between $0\text{\,}\mathrm{T}$ and
$12\text{\,}\mathrm{T}$ at different temperatures. With decreasing
temperature, the slope increases and the magnetization saturates earlier. The
inset shows the magnetization at fixed field of
$h_{z}=$0.1\text{\,}\mathrm{T}$$ over temperature.

Our theoretical results are plotted for temperatures
$T=$1.81\text{\,}\mathrm{K}$$ and $T=$2.98\text{\,}\mathrm{K}$$. This is then
compared against the experimental data measured at
$T=$1.8\text{\,}\mathrm{K}$$ and $T=$3.0\text{\,}\mathrm{K}$$. The conversion
factors between the theoretical calculations and experiment are shown in the
Appendix. From the plots in Fig. 5, we see that as we increase the
temperature, the field value at which the magnetization saturates becomes
larger and approaches the experimental findings. Overall, we find that the
magnetization curve shows a strong dependence on the temperature and the
saturation sets in quicker for low $T$. This seems to indicate that the
earlier discrepancy between theory and experimental data has been largely due
to neglecting finite temperature effects in the theory simulations. For
comparison, we show the $T=$0\text{\,}\mathrm{K}$$ magnetization curve
obtained with a nine-site iPESS at $\chi_{B}=24$ with mean-field environments.

Inelastic neutron scattering has previously revealed that the spin liquid
ground state of this material is destroyed by a magnetic field of
$1\text{\,}\mathrm{T}$ Balz  _et al._ (2017). However, heat capacity
measurements could show that magnetic fields of up to $0.5\text{\,}\mathrm{T}$
leave the spin liquid ground state intact as indicated by featureless $C/T$
curves Balz  _et al._ (2017). To check that the effect of small fields is
correctly captured by our model, we have computed the magnetization as a
function of temperature at fixed field strength
$h_{z}=$0.1\text{\,}\mathrm{T}$$ and contrasted it with measured experimental
data. The comparison is shown in the inset of Fig. 5, and the model indeed
reproduces the experimental curve without any anomalies that would indicate a
phase transition into a magnetically ordered ground state.

####  III.2.2 Heat capacity and entropy

Finally, we compute the magnetic heat capacity from the thermal state energy
$U$ according to $C\coloneqq\partial U/\partial T$ for two different values of
the magnetic field $h_{z}=$2.0\text{\,}\mathrm{T}$$ and
$h_{z}=$3.0\text{\,}\mathrm{T}$$ and compare it with the experimental data.
Results are shown in Fig. 6. The conversion factors between the theoretical
calculations and experiment are again shown in the Appendix.

Figure 6: Magnetic heat capacity $C/T$ as a function of temperature for
different strengths of the magnetic fields $h_{z}$. Integrating this quantity
over $T$ yields the entropy, which accurately approaches the theoretical value
of $S_{\textrm{th}}=6R\ln(2)$, with $R$ the ideal gas constant.

As well as the heat capacity data, we have also computed the thermodynamic
entropy by integrating the heat capacity. This is shown in the inset of Fig.
6. We find good agreement between our theoretical predictions and experimental
data.

The heat capacity of Ca10Cr7O28 at intermediate fields is characterized by a
broad and smooth peak of a Schottky anomaly due to the excitations that become
gapped by the magnetic field. The position of this peak shifts to higher
temperature with increasing field both in the model and the experimental data.
This is in qualitative agreement with an increasing gap due to the Zeeman term
in the Hamiltonian, see Eq. 6. Integrating $C/T$ to obtain the magnetic
entropy shows that the model does well in capturing the total possible entropy
for spin-1/2 over the temperature range up to $15\text{\,}\mathrm{K}$.

##  IV Conclusions and outlook

In this work, we have presented a two-dimensional tensor network algorithm for
studying finite temperature properties for the highly challenging realm of
frustrated systems and two-dimensional quantum materials. We achieve this by
introducing the infinite projected entangled simplex operator algorithm. Our
algorithm explicitly preserves the positive semi-definiteness of the Gibbs
state represented by the iPESO. We use our technique to benchmark against
finite temperature properties of the well-known, paradigmatic model of the
spin-1/2 Kagome Heisenberg anti-ferromagnet and obtain very competitive state-
of-the art results for the thermal state energy and heat capacity. We also
study the melting of the magnetization plateaus of this model in the presence
of external magnetic field at finite temperature. By focusing on the most
prominent $1/3$ plateau, we find that it starts melting and disappears at
temperature $T\sim 2\cdot 10^{-2}$. Moreover, the plateau starts melting from
the lower end of the field, an observation that was also made recently in an
independent exact diagonalization study Schlüter  _et al._ (2022).

Finally, we have investigated the finite temperature properties of the quantum
material Ca10Cr7O28 using our tensor network technique. This is particularly
important due to a recent discrepancy in the magnetization process predicted
by theoretical simulations compared to experimental findings. As a first
theoretical study of this real material that includes both quantum
correlations and finite temperature effects, we find a strong temperature
dependence of the magnetization curve of this material in the presence of an
external magnetic field. We find that on systematically increasing the
temperature, our theoretical simulations approach the experimental data which
was collected at finite temperature. We provide a direct comparison of the
theoretical magnetization data with the experimental data at
$T=$1.8\text{\,}\mathrm{K}$$ and $T=$3.0\text{\,}\mathrm{K}$$ and find them to
be in surprisingly good but not quite perfect agreement. We also computed the
magnetic heat capacity ($C/T$) as a function of temperature at different field
strengths $h_{z}=$2.0\text{\,}\mathrm{T}$$ and
$h_{z}=$3.0\text{\,}\mathrm{T}$$ as well as the entropy $S$. For all these
quantities, we find good agreement with the experimental data.

One can argue that the agreement is striking, given that the Hamiltonian in
Eq. (6) has only been recovered by neutron scattering techniques to finite
precision considering five Heisenberg interactions, while Dzyaloshinskii-
Moriya interactions have been excluded. Furthermore, there are truncation
errors in the classical simulation. One can argue that the present analysis
allows to cross-benchmark quantum experiments with classical simulations. The
findings can also be seen as an invitation, however, to use high-precision
tools of _Hamiltonian learning_ to better identify the actual underlying
microscopic Hamiltonian, given data from Gibbs states Haah  _et al._ (2021);
Anshu  _et al._ (2021), possibly even based on tensor networks akin the
approach taken in Ref. Wilde  _et al._ (2022). These steps would further
contribute to an engineering perspective of studying realistic strongly
correlated quantum materials with tensor networks.

We believe our work to be an important step towards bridging the gap between
theoretical simulations and experimental studies of quantum materials. It
would be straightforward to extend our algorithm to other lattices and
geometries that may suit other quantum materials. By incorporating both
quantum correlations and finite temperature effects, we have now made direct
comparison between experimental and theoretical data possible.

##  V Acknowledgements

The authors are thankful for discussions with Ji-Yao Chen, Dante Kennes,
Corinna Kollath, David Luitz, Jan Naumann, Román Orús, Matteo Rizzi and Anne-
Maria Visuri. The authors would like to thank the HPC Service of ZEDAT, Freie
Universität Berlin, for computing time Bennett  _et al._ (2020). The FUB team
acknowledges funding by the Deutsche Forschungsgemeinschaft (CRC 183 on
‘Entangled states of matter’ and FOR 2724 on ‘Thermal machines in the quantum
world’), the Helmholtz Association, and the BMBF (MUNIQC-ATOMS), for which
this work constitutes method development. B. L. acknowledges the support of
Deutsche Forschungsgemeinschaft through project B06 of SFB 1143 on ‘Correlated
magnetism: From frustration to topology’ (ID 247310070).

## References

  * Kohn and Sham (1965) W. Kohn and L. J. Sham, “Self-consistent equations including exchange and correlation effects,” Phys. Rev. 140, A1133–A1138 (1965).
  * Tkatchenko and Scheffler (2009) A. Tkatchenko and M. Scheffler, “Accurate molecular van der waals interactions from ground-state electron density and free-atom reference data,” Phys. Rev. Lett. 102, 073005 (2009).
  * Becca and Sorella (2017) F. Becca and S. Sorella,  _Quantum Monte Carlo approaches for correlated systems_ (Cambridge University Press, Cambridge, 2017).
  * Orús (2014) R. Orús, “A practical introduction to tensor networks: Matrix product states and projected entangled pair states,” Ann. Phys. 349, 117–158 (2014).
  * Verstraete  _et al._ (2008) F. Verstraete, J. I. Cirac, and V. Murg, “Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems,” Adv. Phys. 57, 143 (2008).
  * Bridgeman and Chubb (2017) J. C. Bridgeman and C. T. Chubb, “Hand-waving and interpretive dance: An introductory course on tensor networks,” J. Phys. A 50, 223001 (2017).
  * Eisert  _et al._ (2010) J. Eisert, M. Cramer, and M. B. Plenio, “Area laws for the entanglement entropy,” Rev. Mod. Phys. 82, 277 (2010).
  * White (1992) S. R. White, “Density matrix formulation for quantum renormalization groups,” Phys. Rev. Lett. 69, 2863–2866 (1992).
  * Schollwöck (2005) U. Schollwöck, “The density-matrix renormalization group,” Rev. Mod. Phys. 77, 259–315 (2005).
  * Schollwöck (2011) U. Schollwöck, “The density-matrix renormalization group in the age of matrix product states,” Ann. Phys. 326, 96–192 (2011).
  * Perez-Garcia  _et al._ (2007) D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac, “Matrix product state representations,” Quantum Inf. Comput. 7, 401 (2007).
  * Verstraete and Cirac (2004) F. Verstraete and I. Cirac, “Renormalization algorithms for quantum-many body systems in two and higher dimensions,” arXiv , cond–mat/0407066 (2004).
  * Jordan  _et al._ (2008) J. Jordan, R. Orus, G. Vidal, F. Verstraete, and J. I. Cirac, “Classical simulation of infinite-size quantum lattice systems in two spatial dimensions,” Phys. Rev. Lett. 101, 250602 (2008).
  * Liao  _et al._ (2017) H. J. Liao, Z. Y. Xie, J. Chen, Z. Y. Liu, H. D. Xie, R. Z. Huang, B. Normand, and T. Xiang, “Gapless spin-liquid ground state in the $S=1/2$ Kagome anti-ferromagnet,” Phys. Rev. Lett. 118, 137202 (2017).
  * Picot  _et al._ (2016a) T. Picot, M. Ziegler, R. Orús, and D. Poilblanc, “Spin-$S$ Kagome quantum anti-ferromagnets in a field with tensor networks,” Phys. Rev. B 93, 060407 (2016a).
  * Picot and Poilblanc (2015) T. Picot and D. Poilblanc, “Nematic and supernematic phases in Kagome quantum anti-ferromagnets under the influence of a magnetic field,” Phys. Rev. B 91, 064415 (2015).
  * Kshetrimayum  _et al._ (2016) A. Kshetrimayum, T. Picot, R. Orús, and D. Poilblanc, “Spin-$\frac{1}{2}$ Kagome XXZ model in a field: Competition between lattice nematic and solid orders,” Phys. Rev. B 94, 235146 (2016).
  * Boos  _et al._ (2019) C. Boos, S. P. G. Crone, I. A. Niesen, P. Corboz, K. P. Schmidt, and F. Mila, “Competition between intermediate plaquette phases in SrCu2(BO3)2 under pressure,” Phys. Rev. B 100, 140413 (2019).
  * Kshetrimayum  _et al._ (2020a) A. Kshetrimayum, C. Balz, B. Lake, and J. Eisert, “Tensor network investigation of the double layer Kagome compound Ca10Cr7O28,” Ann. Phys. (N.Y.) 421, 168292 (2020a).
  * Astrakhantsev  _et al._ (2021) N. Astrakhantsev, F. Ferrari, N. Niggemann, T. Müller, A. Chauhan, A. Kshetrimayum, P. Ghosh, N. Regnault, R. Thomale, J. Reuther, T. Neupert, and Y. Iqbal, “Pinwheel valence bond crystal ground state of the spin-$\frac{1}{2}$ Heisenberg anti-ferromagnet on the Shuriken lattice,” Phys. Rev. B 104, L220408 (2021).
  * Kshetrimayum  _et al._ (2017) A. Kshetrimayum, H. Weimer, and R. Orús, “A simple tensor network algorithm for two-dimensional steady states,” Nature Comm. 8, 1291 (2017).
  * Czarnik  _et al._ (2019a) P. Czarnik, J. Dziarmaga, and P. Corboz, “Time evolution of an infinite projected entangled pair state: An efficient algorithm,” Phys. Rev. B 99, 035115 (2019a).
  * Hubig and Cirac (2019) C. Hubig and J. I. Cirac, “Time-dependent study of disordered models with infinite projected entangled pair states,” SciPost Phys. 6, 31 (2019).
  * Kshetrimayum  _et al._ (2020b) A. Kshetrimayum, M. Goihl, and J. Eisert, “Time evolution of many-body localized systems in two spatial dimensions,” Phys. Rev. B 102, 235132 (2020b).
  * Kshetrimayum  _et al._ (2021) A. Kshetrimayum, M. Goihl, D. M. Kennes, and J. Eisert, “Quantum time crystals with programmable disorder in higher dimensions,” Phys. Rev. B 103, 224205 (2021).
  * Dziarmaga (2021) J. Dziarmaga, “Time evolution of an infinite projected entangled pair state: Neighborhood tensor update,” Phys. Rev. B 104, 094411 (2021).
  * Dziarmaga (2022) J. Dziarmaga, “Time evolution of an infinite projected entangled pair state: A gradient tensor update in the tangent space,” Phys. Rev. B 106, 014304 (2022).
  * Czarnik  _et al._ (2012) P. Czarnik, L. Cincio, and J. Dziarmaga, “Projected entangled pair states at finite temperature: Imaginary time evolution with ancillas,” Phys. Rev. B 86, 245101 (2012).
  * Czarnik and Dziarmaga (2015) P. Czarnik and J. Dziarmaga, “Variational approach to projected entangled pair states at finite temperature,” Phys. Rev. B 92, 035152 (2015).
  * Kshetrimayum  _et al._ (2019) A. Kshetrimayum, M. Rizzi, J. Eisert, and R. Orús, “Tensor network annealing algorithm for two-dimensional thermal states,” Phys. Rev. Lett. 122, 070502 (2019).
  * Czarnik  _et al._ (2019b) P. Czarnik, A. Francuz, and J. Dziarmaga, “Tensor network simulation of the Kitaev-Heisenberg model at finite temperature,” Phys. Rev. B 100, 165147 (2019b).
  * Czarnik  _et al._ (2021) P. Czarnik, M. M. Rams, P. Corboz, and J. Dziarmaga, “Tensor network study of the $m=1/2$ magnetization plateau in the Shastry-Sutherland model at finite temperature,” Phys. Rev. B 103, 075113 (2021).
  * Mondal  _et al._ (2020) S. Mondal, A. Kshetrimayum, and T. Mishra, “Two-body repulsive bound pairs in a multibody interacting Bose-Hubbard model,” Phys. Rev. A 102, 023312 (2020).
  * Jiang  _et al._ (2008a) H. C. Jiang, Z. Y. Weng, and T. Xiang, “Accurate determination of tensor network state of quantum lattice models in two dimensions,” Phys. Rev. Lett. 101, 090603 (2008a).
  * Gauthé and Mila (2022) O. Gauthé and F. Mila, “Thermal Ising transition in the spin-$1/2$ ${J}_{1}$-${J}_{2}$ Heisenberg model,” Phys. Rev. Lett. 128, 227202 (2022).
  * Chen  _et al._ (2018a) B.-B. Chen, L. Chen, Z. Chen, W. Li, and A. Weichselbaum, “Exponential thermal tensor network approach for quantum lattice models,” Phys. Rev. X 8, 031082 (2018a).
  * Werner  _et al._ (2016) A. H. Werner, D. Jaschke, P. Silvi, M. Kliesch, T. Calarco, J. Eisert, and S. Montangero, “Positive tensor network approach for simulating open quantum many-body systems,” Phys. Rev. Lett. 116, 237201 (2016).
  * Chen  _et al._ (2018b) X. Chen, S.-J. Ran, T. Liu, C. Peng, Y.-Z. Huang, and G. Su, “Thermodynamics of spin-1/2 Kagomé Heisenberg anti-ferromagnet: algebraic paramagnetic liquid and finite-temperature phase diagram,” Science Bulletin 63, 1545–1550 (2018b).
  * Xie  _et al._ (2014) Z. Y. Xie, J. Chen, J. F. Yu, X. Kong, B. Normand, and T. Xiang, “Tensor renormalization of quantum many-body systems using projected entangled simplex states,” Phys. Rev. X 4, 011025 (2014).
  * Nishino and Okunishi (1996) T. Nishino and K. Okunishi, “Corner transfer matrix renormalization group method,” J. Phys. Soc. Jap. 65, 891–894 (1996).
  * Orús and Vidal (2009) R. Orús and G. Vidal, “Simulation of two-dimensional quantum systems on an infinite lattice revisited: Corner transfer matrix for tensor contraction,” Phys. Rev. B 80, 094403 (2009).
  * Orús (2012) R. Orús, “Exploring corner transfer matrices and corner tensors for the classical simulation of quantum lattice systems,” Phys. Rev. B 85, 205117 (2012).
  * Sachdev (1992) S. Sachdev, “Kagome- and triangular-lattice Heisenberg anti-ferromagnets: Ordering from quantum fluctuations and quantum-disordered ground states with unconfined bosonic spinons,” Phys. Rev. B 45, 12377–12396 (1992).
  * Ran  _et al._ (2007) Y. Ran, M. Hermele, P. A. Lee, and X.-G. Wen, “Projected-wave-function study of the spin-$1/2$ Heisenberg model on the Kagomé lattice,” Phys. Rev. Lett. 98, 117205 (2007).
  * Jiang  _et al._ (2008b) H. C. Jiang, Z. Y. Weng, and D. N. Sheng, “Density matrix renormalization group numerical study of the Kagome anti-ferromagnet,” Phys. Rev. Lett. 101, 117203 (2008b).
  * Jiang  _et al._ (2012) H.-C. Jiang, Z. Wang, and L. Balents, “Identifying topological order by entanglement entropy,” Nature Phys. 8, 902–905 (2012).
  * Götze  _et al._ (2011) O. Götze, D. J. J. Farnell, R. F. Bishop, P. H. Y. Li, and J. Richter, “Heisenberg anti-ferromagnet on the kagome lattice with arbitrary spin: A higher-order coupled cluster treatment,” Phys. Rev. B 84, 224428 (2011).
  * Iqbal  _et al._ (2013) Y. Iqbal, F. Becca, S. Sorella, and D. Poilblanc, “Gapless spin-liquid phase in the Kagome spin-$\frac{1}{2}$ Heisenberg anti-ferromagnet,” Phys. Rev. B 87, 060405 (2013).
  * Läuchli  _et al._ (2019) A. M. Läuchli, J. Sudan, and R. Moessner, “$S=\frac{1}{2}$ Kagome Heisenberg anti-ferromagnet revisited,” Phys. Rev. B 100, 155142 (2019).
  * Picot  _et al._ (2016b) T. Picot, M. Ziegler, R. Orús, and D. Poilblanc, “Spin-$S$ Kagome quantum anti-ferromagnets in a field with tensor networks,” Phys. Rev. B 93, 060407 (2016b).
  * Schlüter  _et al._ (2022) H. Schlüter, J. Richter, and J. Schnack, “Melting of magnetization plateaus for Kagomé and square-Kagomé lattice anti-ferromagnets,” J. Phys. Soc. Jap. 91, 094711 (2022).
  * Balz  _et al._ (2017) C. Balz, B. Lake, A. T. M. Nazmul Islam, Y. Singh, J. A. Rodriguez-Rivera, T. Guidi, E. M. Wheeler, G. G. Simeoni, and H. Ryll, “Magnetic Hamiltonian and phase diagram of the quantum spin liquid ${\mathrm{Ca}}_{10}{\mathrm{Cr}}_{7}{\mathrm{O}}_{28}$,” Phys. Rev. B 95, 174414 (2017).
  * Pohle  _et al._ (2021) R. Pohle, H. Yan, and N. Shannon, “Theory of ${\mathrm{Ca}}_{10}{\mathrm{Cr}}_{7}{\mathrm{O}}_{28}$ as a bilayer breathing-Kagome magnet: Classical thermodynamics and semiclassical dynamics,” Phys. Rev. B 104, 024426 (2021).
  * Haah  _et al._ (2021) J. Haah, R. Kothari, and E. Tang, “Optimal learning of quantum Hamiltonians from high-temperature Gibbs states,” (2021), arXiv:2108.04842.
  * Anshu  _et al._ (2021) A. Anshu, S. Arunachalam, T. Kuwahara, and M. Soleimanifar, “Sample-efficient learning of quantum many-body systems,” Nature Phys. 17, 931–935 (2021).
  * Wilde  _et al._ (2022) F. Wilde, A. Kshetrimayum, I. Roth, D. Hangleiter, R. Sweke, and J. Eisert, “Scalably learning quantum many-body Hamiltonians from dynamical data,” (2022), arXiv:2209.14328.
  * Bennett  _et al._ (2020) L. Bennett, B. Melchers, and B. Proppe, “Curta: A general-purpose high-performance computer at ZEDAT, Freie Universität Berlin,” http://dx.doi.org/10.17169/refubium-26754 (2020).
  * Singh  _et al._ (2010) S. Singh, R. N. C. Pfeifer, and G. Vidal, “Tensor network decompositions in the presence of a global symmetry,” Phys. Rev. A 82, 050301 (2010).
  * Schmoll  _et al._ (2020) P. Schmoll, S. Singh, M. Rizzi, and R. Orús, “A programming guide for tensor networks with global SU(2) symmetry,” Ann. Phys. 419, 168232 (2020).
  * Schuch  _et al._ (2007) N. Schuch, M. M. Wolf, F. Verstraete, and J. I. Cirac, “Computational complexity of projected entangled pair states,” Phys. Rev. Lett. 98, 140506 (2007).
  * Haferkamp  _et al._ (2020) J. Haferkamp, D. Hangleiter, J. Eisert, and M. Gluza, “Contracting projected entangled pair states is average-case hard,” Phys. Rev. Research 2, 013010 (2020).
  * Phien  _et al._ (2015) H. N. Phien, J. A. Bengua, H. D. Tuan, P. Corboz, and R. Orús, “Infinite projected entangled pair states algorithm improved: Fast full update and gauge fixing,” Phys. Rev. B 92, 035142 (2015).
  * Czarnik  _et al._ (2019c) P. Czarnik, J. Dziarmaga, and P. Corboz, “Time evolution of an infinite projected entangled pair state: An efficient algorithm,” Phys. Rev. B 99, 035115 (2019c).

##  Appendix A Details of the tensor network algorithm

###  A.1 Simple update

The simple update describes an efficient, yet approximate scheme to do the
annealing respectively the imaginary time evolution of the initial density
matrix. It essentially implements the evolution in Eq. (8) together with a
local truncation to keep the bulk bond dimension fixed. Without the loss of
generality, we consider a Hamiltonian with local interactions in the form of

$\displaystyle H=H_{\bigtriangledown}+H_{\bigtriangleup},$ (7)

where $H_{\bigtriangledown}$ and $H_{\bigtriangleup}$ are three-spin
interactions on the two types of triangles of the Kagome lattice,
respectively. Making use of a first-order _Suzuki-Trotter decomposition_ , the
imaginary time evolution to evolve the thermal density matrix
$\rho(\beta)\rightarrow\rho(\beta+\delta\beta)$ can be approximated by
applying the operator

$\displaystyle\begin{split}U(\delta\beta)&=\mathrm{e}^{-\delta\beta
H_{\bigtriangledown}}\mathrm{e}^{-\delta\beta
H_{\bigtriangleup}}+O(\delta\beta^{2})\\\ &\approx
U_{\bigtriangledown}(\delta\beta)U_{\bigtriangleup}(\delta\beta)\end{split}$
(8)

to both three-site configurations in the tensor network. In Fig. 7 we
illustrate the evolution of the iPESO with the three-body gate
$U_{\bigtriangledown}(\delta\beta)$. This step involves three lattice tensors,
as well as the simplex tensor $\bigtriangledown$.

Figure 7: Simple update step for the (imaginary) time evolution of the iPESO
ansatz. After applying the three-body gate to the physical indices, a
truncated higher-order SVD is used to decompose the nine-index tensor and
restore the lattice and simplex tensors.

After the gate has been contracted with the tensors, a _higher-order singular
value decomposition_ (HOSVD) with subsequent truncation is used to separate
the network back into simplex and lattice tensors. Since the truncation is
based only on the singular values for the three indices, it is purely local.
In a similar fashion, the simplex $\bigtriangleup$ is updated alongside the
three lattice site tensors by applying the three-body gate
$U_{\bigtriangleup}(\delta\beta)$. After both steps have been performed, we
obtain the thermal density matrix $\rho(\beta+\delta\beta)$, represented by a
three-site iPESO. This process is repeated for a fixed number of steps, such
that the final thermal density matrix represents the quantum system at the
desired (inverse) temperature. Naturally, this can be extended to Hamiltonians
with less or more than three-site interactions. Additionally, physical
symmetries of the Hamiltonian (like $U(1)$ or $SU(2)$) can be readily directly
incorporated, exploiting symmetry-preserving tensors Singh  _et al._ (2010);
Schmoll  _et al._ (2020).

###  A.2 Bond dimension considerations

As it is common in tensor network applications, the bond dimension controls
the precision of the simulations. Here, we aim at presenting a discussion of
the minimal bond dimensions required in order to obtain meaningful results.
The infinitesimal thermal density matrix

$\displaystyle\rho(\delta\beta)=\prod_{\langle
i,j,k\rangle}\mathrm{e}^{-\delta\beta
H_{i,j,k}}+\mathcal{O}(\delta\beta^{2}),$ (9)

here to first order in the Suzuki-Trotter decomposition, can be constructed by
applying the Trotterized Hamiltonian gates $\exp(-\delta\beta H_{i,j,k})$ onto
the infinite temperature Gibbs state $\rho(\beta=0)$, as shown in Fig. 8 (the
infinite temperature state is simply a tensor product of identity matrices).

Figure 8: Construction of the infinitesimal thermal density operator
$\rho(\delta\beta)$ by a decomposition of the Trotter gates (see Eq. (8)). An
exact representation (apart from the inevitable Trotter error) can only be
achieved without truncating the virtual bulk bond dimension, which is
therefore at least $p^{2}$.

An accurate representation (within the inevitable Trotter error) of this state
is only possible if the resulting iPESO tensors are not truncated. Since the
infinite temperature state has a bond dimension of unity, the infinitesimal
thermal density matrix necessarily has bond dimension $p^{2}$, where $p$ is
the physical dimension of the system. Naturally, cooling the state down to
lower temperatures can only produce meaningful results, if the bond dimension
is larger than the minimally necessary one.

###  A.3 Effective environments and expectation values

In order to evaluate physical observables and compute expectation values
accurately, the infinite two-dimensional iPESO tensor network needs to be
contracted. It is known that this task cannot be performed exactly classically
efficiently both in worst case and average case complexity Schuch  _et al._
(2007); Haferkamp  _et al._ (2020), without an exponential increase in
computation time, so that approximate methods must be employed. Here we
utilize the so-called _corner transfer matrix renormalization group_ (CTMRG)
to compute the effective environment tensors for every lattice site. To this
end, we coarse-grain the iPESO network to an iPEPO network, the operator form
of the famous _infinite projected entangled pair state_ (iPEPS), as visualized
in Fig. 9.

Figure 9: Coarse-graining of a three-site iPESO into a single-site iPEPO
tensor. For larger iPESO unit cells the resulting iPEPO network will have a
larger unit cell, too.

After coarse-graining, the environment surrounding each local thermal density
matrix can be conveniently computed using a standard CTMRG procedure. To this
end, the contraction of the infinite square lattice is approximated by a set
of fix-point environment tensors, as shown in Fig. 10.

Figure 10: A directional CTMRG routine is used to approximate the contraction
of the infinite square lattice by a set of fixed-point environment tensors,
denoted in grey.

This enables both accurate calculations of expectation values and would be
essential in devising a sophisticated update procedure that includes all
quantum correlations in the system — the so-called _full update_ Phien  _et
al._ (2015); Czarnik  _et al._ (2019c). Since the two physical indices are
traced over, the procedure is a straightforward extension of a regular CTMRG
routine for a two-dimensional iPEPS wave function. In order to ensure that the
thermal density matrix is reflected by a positive semi-definite operator, a
double-layer approach is taken in contrast to the original proposal in Ref.
Kshetrimayum  _et al._ (2019).

###  A.4 Truncation effects in the simple update

The annealing scheme adopted in this study is based on the simple update,
which requires truncations in order to keep the bulk bond dimension $\chi_{B}$
constant. Moreover, the choice of the infinitesimal cooling step $\delta\beta$
controls the unavoidable error in the Trotterization, and the number of
annealing steps which include a truncation. For the iPESO simulations of
Ca10Cr7O28 we choose a step size of $\delta\beta=10^{-2}$, which leads to the
accumulated truncation errors shown in Fig. 11, for several temperatures.

The accumulated truncation error is given by the sum of the discarded weights
of all simple update annealing steps. The discarded weight is the sum of the
discarded squared singular values in the _singular value decomposition_ (SVD)
Schollwöck (2011). Naturally, the truncation error decreases with increasing
magnetic field, since the thermal states become closer to a reduced density
matrix of a product state. Moreover, it increases with decreasing temperature,
because more cooling steps are necessary to reach lower temperatures. In
general, the large bond dimension of $\chi_{B}=30$ keeps the accumulated
truncation error low enough for our results to be meaningful down to the
lowest temperatures we consider.

Figure 11: Accumulated truncation error in the simple update annealing for the
bilayer Kagome compound Ca10Cr7O28 at $\chi_{B}=30$ and for
$\delta\beta=10^{-2}$.

##  Appendix B Conversion between experiment and simulation

In order to compare the simulated tensor network data with measured
experimental data, we need to apply the correct conversion factors. Since the
coupling constants in the Hamiltonian are given in units of $\mathrm{meV}$ and
we set $k_{B}=1$, the temperature $T_{\text{sim}}$ is in $\mathrm{meV}$, too.
The proper conversion to $\mathrm{K}$ is, therefore, given by

$\displaystyle\frac{T_{\text{exp}}}{T_{\text{sim}}}=\frac{$1\text{\,}\mathrm{meV}$}{k_{B}}=\frac{$1.602\text{\times}{10}^{-22}\text{\,}\mathrm{J}$}{k_{B}}\approx$11.6\text{\,}\mathrm{K}$.$
(10)

Furthermore, we need to convert the heat capacity $C$ between simulated and
measured data. The tensor network data is given per spin in units of
$\mathrm{meV}\text{\,}{\mathrm{K}}^{-1}$. In order to convert it, a factor of

$\displaystyle\frac{C_{\text{exp}}}{C_{\text{sim}}}=6\cdot$1.602\text{\times}{10}^{-22}\text{\,}\mathrm{J}$\cdot
N_{A}\approx$578.8\text{\,}\mathrm{J}\text{\,}{\mathrm{mol}}^{-1}$$ (11)

with $N_{A}=$6.022\text{\times}{10}^{23}\text{\,}{\mathrm{mol}}^{-1}$$ the
Avogadro constant, is required. The additional factor of six stems from the
fact, that one formula unit of Ca10Cr7O28 has six spin-$1/2$ chromium ions.