# Data-Model-Circuit Tri-Design for Ultra-Light Video Intelligence on Edge
Devices

Yimeng Zhang  Michigan State University, USA ,  Akshay Karkal Kamath  Georgia
Institute of Technology, USA ,  Qiucheng Wu  UC, Santa Barbara, USA ,  Zhiwen
Fan  University of Texas at Austin, USA ,  Wuyang Chen  University of Texas at
Austin, USA ,  Zhangyang Wang  University of Texas at Austin, USA ,  Shiyu
Chang  UC, Santa Barbara, USA ,  Sijia Liu  Michigan State University, USA and
Cong Hao  Georgia Institute of Technology, USA

(2023)

###### Abstract.

In this paper, we propose a data-model-hardware tri-design framework for high-
throughput, low-cost, and high-accuracy multi-object tracking (MOT) on High-
Definition (HD) video stream. First, to enable ultra-light video intelligence,
we propose temporal frame-filtering and spatial saliency-focusing approaches
to reduce the complexity of massive video data. Second, we exploit structure-
aware weight sparsity to design a hardware-friendly model compression method.
Third, assisted with data and model complexity reduction, we propose a
sparsity-aware, scalable, and low-power accelerator design, aiming to deliver
real-time performance with high energy efficiency. Different from existing
works, we make a solid step towards the synergized software/hardware co-
optimization for realistic MOT model implementation. Compared to the state-of-
the-art MOT baseline, our tri-design approach can achieve 12.5$\times$ latency
reduction, 20.9$\times$ effective frame rate improvement, $5.83\times$ lower
power, and $9.78\times$ better energy efficiency, without much accuracy drop.

††journalyear: 2023††copyright: rightsretained††conference: 28th Asia and
South Pacific Design Automation Conference; January 16–19, 2023; Tokyo,
Japan††booktitle: 28th Asia and South Pacific Design Automation Conference
(ASPDAC ’23), January 16–19, 2023, Tokyo, Japan††doi:
10.1145/3566097.3567864††isbn: 978-1-4503-9783-4/23/01

##  1\. introduction

The past decade has witnessed a tremendous success of deep neural networks
(DNNs) (Alom et al., 2018). In this paper, we focus on DNN-based multi-object
tracking (MOT), which has been used as a technological basis for video
intelligence and is drastically improving the quality of human life, e.g.,
autonomous driving (AD) (Yurtsever et al., 2020). Despite a proliferation of
MOT techniques (Kalake et al., 2021), little research effort has been made
towards the data-efficient and co-design solutions for MOT across the full
software/hardware stack. Our work aims to close this gap by developing a high-
throughput, low-cost, and high-accuracy video processing algorithm/hardware
pipeline. We provide a holistic viewpoint on how contemporary MOT methods can
be renovated to deal with massive quantities of data and high-complexity DNNs.

Despite extensive existing works, data and model efficiencies of MOT are
mainly studied either from the algorithm perspective or from the model design
perspective. For example, in (Wang et al., 2020a), the MOT efficiency is
improved by leveraging meta learning (Finn et al., 2017) to solve an instance
detection problem. In addition to algorithms, several works have focused on
developing model compression techniques (such as weight pruning and weight
quantization) to reduce the model complexity of MOT (Stäcker et al., 2021;
Muhawenayo et al., 2021). However, the efficiency of “compressed” models are
evaluated without considering the practical hardware platform, such as low-
power FPGAs. Known FPGA acceleration techniques (Hao et al., 2019; Zhang et
al., 2018) do not explore data sparsity or model sparsity, thereby missing
huge optimization opportunities. Additionally, existing accelerators are
evaluated on the ImageNet dataset with small input image sizes and do not
scale to real-world High-Definition (HD) video frames, limited by the device
and/or design tool capabilities. To the best of our knowledge, there is no
prior work that discusses efficient implementation of MOT on the edge for HD
video processing by fully utilizing data- and model-level sparsity.

Contributions. To close the gap between MOT algorithm design and efficient
implementation, we propose a data-model-hardware tri-design approach. We
summarize our contributions as follows.

$\bullet$ Dynamic frame filtering. We propose a reinforcement learning-based
lightweight algorithm to achieve temporal data reduction.

$\bullet$ Spatial attention focusing. We propose a saliency-guided spatial
data reduction method to eliminate uninformative pixels from both the input
frames and the intermediate feature maps.

$\bullet$ Hardware-aware model compression. We leverage kernel-wise pattern-
aware sparsity of an MOT model to achieve hardware-friendly model compression.

$\bullet$ Tri-design implementation and evaluation. We implement the proposed
tri-design framework on a hardware platform comprising of one Xilinx ZCU104
and two Xilinx Alveo U50 FPGAs and conduct extensive experiments to evaluate
its effectiveness in both accuracy and on-device efficiency metrics. We
demonstrate that our approach can achieve 12.5$\times$ latency reduction,
20.9$\times$ effective frame rate improvement, $5.83\times$ lower power, and
$9.78\times$ better energy efficiency, without much accuracy drop compared to
the state-of-the-art MOT baseline (Pang et al., 2021) on the BDD100K dataset
(Yu et al., 2020).

##  2\. Preliminaries and Motivation

In this paper, we focus on MOT (multiple-object tracking) which deals with
massive video data. Current MOT methods are commonly built upon DNNs and
follow the paradigm of tracking-by-detection (Wang et al., 2020b). QDTrack
(Pang et al., 2021) is the state-of-the-art MOT algorithm introduced on the
BDD100K dataset (Yu et al., 2020), one of the most representative and
challenging self-driving car datasets. QDTrack employs Faster R-CNN (Ren et
al., 2015) with Feature Pyramid Network (FPN) (Lin et al., 2017) as an object
detection backbone network, leverages Contrastive Learning (Wu et al., 2018)
to optimize the backbone network parameters, and utilizes Bi-directional
Softmax (Pang et al., 2021) in the embedding space for object association and
tracking. The focus of this work is not to develop new MOT model algorithms
but to instead build a high-throughput, low-cost, and high-accuracy video
processing implementation pipeline across the full data-processing/DNN
operation/hardware stack.

The computational landscape of MOT has been rapidly evolving on dedicated
computing hardware (e.g., FPGAs and ASICs) (Zhang et al., 2021). Despite the
advancements of hardware accelerators (Li et al., 2019) and software/hardware
co-design techniques (Hao et al., 2019), critical limitations still exist for
real-time video processing. First, there is a sharp increase not only in ML
model size but also computation complexity, especially for HD video frames.
With HD inputs, any two $3\times 3$ convolution layers of ResNet-50 already
have higher complexity than the entire ResNet-50 with ImageNet inputs. This
asserts huge pressure for real-time inference on hardware. Second, both the
hardware devices and design tools exhibit poor scalability. For example, the
most powerful FPGA till date, Xilinx Alveo U55C, has 9,024 DSP slices. It can
provide peak performance of 2,256 GOPs for QDTrack (157.9 GOPs), but the
performance upper bound is 70 ms per frame, which is far from real-time
requirement. Meanwhile, the design tools are strictly limited by the largest
parallelism that can be supported. For example, Xilinx Vitis HLS tool
constrains that the parallelism factor should be lower than 4096 (Aarrestad et
al., 2021), while Intel Design Toolkit does not support more than 256
concurrent kernels. With larger parallelism, hardware designs can easily take
days or even weeks to synthesize, which drastically increases the development
cycle and the time-to-market (TTM). Consequently, there is still a huge gap
between the low-power real-time video processing requirements and the limited
capability of existing hardware acceleration techniques.

As inspired by above, our work advances the design of practical MOT systems by
tackling the following challenges (C1–C3). (C1): The massive video data and
their high frame rates make the implementation of MOT prohibitively expensive
in energy and latency. Thus, our first goal is to develop novel data
complexity reduction techniques that can fully exploit the temporal and
saliency redundancy of video frames over time and space. (C2): SOTA DNN-based
MOT models typically contain gigantic number of weights making them poor
candidates for hardware implementation. To this end, we propose a hardware-
friendly model compression technique for MOT. (C3): A holistic data-model-
hardware full-stack design approach is lacking. It remains elusive how current
MOT solutions can be applied subject to Small Size, Weight, and Power (SWaP)
constraints in practice. To the best of our knowledge, there is no prior work
that addresses the challenges (C1–C3) in a single unified MOT implementation
pipeline.

##  3\. Methodology

Figure 1.  Overview of data-model-hardware tri-design.

###  3.1. Overview of proposed approach.

The proposed data-model-hardware tri-design includes four main components as
shown in Fig. 1.  1 (Data) Temporal frame filtering. It is achieved using
reinforcement learning (RL) to filter unnecessary frames at the front-end data
processing stage and eventually reduce processing power along with latency
costs on hardware.  2 (Data) Spatial saliency focusing. It provides the
second-stage data complexity reduction by peering into the salient region of a
frame to eliminate uninformative pixels.  3 (Model) Hardware-friendly deep
model pruning. It exploits the structured sparsity patterns of weights to
reduce the model size on hardware without much loss of accuracy.  4 (Hardware)
Backend acceleration. Guided by the algorithmic pipeline, it adopts low-power
FPGAs to achieve realistic implementation of MOT with high throughout, low
latency, and high energy efficiency. In a nutshell, components  1 and  2 are
designed for tackling the data complexity challenge (in response to C1);
Component  3 reduces the complexity of back-end model (in response to C2);
Component  4 unveils the practical efficiency of hardware solutions for MOT
(in response to C3).

###  3.2. Data complexity reduction

####  3.2.1. RL-based temporal frame filtering

The embedded MOT solution requires a high throughput in terms of the number of
frames that can be processed in unit time. The frame processing rate should
also consider the practical latency cost. Considering the massive quantities
of video frames, achieving a video processing rate of less than 30ms (for
real-time inference) becomes implausible even on high-end FPGAs. Due to this
hardware constraint, it is necessary to reduce the data complexity through the
lens of temporal video frame dropping. Particularly, in the AD (autonomous
driving) datasets, consecutive video frames may contain largely-overlapped
frame content such as street and road scenes. Thus, temporal frame filter is
not only necessary but also reasonable for real-time MOT. In our preliminary
experiments, we test a random filtering policy with different dropping ratios.
As shown in Fig. 7, dropping frames randomly significantly hampers the object
tracking performance. The incapability of random filtering drives us to design
a lightweight learning-based solution that can self-adjust the frame dropping
scheme based on the task-end tracking performance.

Inspired by above, we propose an RL-based policy network for temporal frame
filtering. Our policy network takes the current frame and the difference
between the current frame and the previous frame as inputs. At the beginning
of the policy network, we adopt a lightweight 3-layer convolutional network to
extract a 2D feature map with the same height and width as the input frames.
In compensation with a shallow network, we adopt 5$\times$5 kernels at each
convolution layer. The numbers of filters of the three convolution layers are
16, 8, and 1. Finally, with a pooling layer and sigmoid layer, we obtain the
probability of dropping the current frame, which indicates the salience of the
current input frame.

We optimize this policy model through a regularized policy gradient method.
Specifically, given a video sequence $\mathbf{v}$ with $n$ frames
$\\{v_{0},v_{1},...,v_{n}\\}$, the policy model $G$ generates a frame
importance score $G(v_{i})\in(0,1)$ for each frame $v_{i}$. $G(v_{i})$ serves
as the probability to retain the current frame and determines the new video
sequence $\mathbf{v}^{\prime}$ with $n^{\prime}$ frames that only consists of
the selected video frames. The training objective function of $G$ is given by

(1) $J=-R(\mathbf{v}^{\prime})\sum_{i=1}^{n}\log
p(G(v_{i}))+\alpha\sum_{i=1}^{n}(G(v_{i})-\mu)^{2},$

where $R(\mathbf{v}^{\prime})$ is the reward of the updated video sequence
$\mathbf{v}^{\prime}$ and $p(G(v_{i}))$ is the probability to sample an
individual frame $v_{i}$. We use the number of ID Switches per frame (i.e.,
IDSw in Sec. 4.1) as the reward function to evaluate the quality of frame
dropping:
$R(\mathbf{v}^{\prime})=-\mathrm{IDSw}(\mathbf{v}^{\prime})/n^{\prime}$, where
$n^{\prime}$ is the number of remaining frames in the video sequence
$\mathbf{v}^{\prime}$. In (1), inspired by (Zhou et al., 2018), we also
regularize the salience scores $G(v_{i})$ with a hyperparameter $\mu$ to avoid
the case of no frame dropping and obtaining high reward for free.

####  3.2.2. Spatial saliency focusing.

Based on the first-stage temporal filtering, we next propose the second-stage
data reduction, which attentively focuses on the salient region of a frame and
its feature map to eliminate uninformative pixels.

Towards the goal of spatial reduction, our MOT model (i.e., QDTrack) achieves
spatial saliency by detecting and characterizing important spatial signals
while monitoring high-rate video streams with low latency. The output is a
localized sub-image for each frame, on which the subsequent backend algorithm
on FPGA can focus. We propose a spatial saliency reduction (SSR) method.
First, it decomposes an input image into multiple patches, and then adopts the
Sobel operator (Jin-Yu et al., 2009) to compute an approximate saliency score
of each patch. This operator exploits two 3 $\times$ 3 convolution kernels on
the given patch to calculate the saliency values (i.e., derivatives). Once the
patch-wise saliency score is obtained, we smoothen it using its four
neighboring patches and create a binary mask indicating the patches of an
input frame that should be discarded.

In addition to input saliency dropping, we further propose feature saliency
dropping, where patches of intermediate features, i.e., activations, can also
be filtered. Our motivations are two-fold. First, even if an input patch has
been dropped, its corresponding feature region might still have responses
“leaked” from neighborhood unfiltered areas, caused by filter stride, padding,
etc. Second, some image regions might be recognized as redundant at the
intermediate feature level rather than at the beginning input level. As shown
in Fig. 2(e) and (g), even if a region is not dropped initially, it leads to
low responses at later stages. Therefore, reducing spatial redundancy beyond
the input level to feature level (Fig. 2(f) and (h)) is valuable and is also
beneficial to hardware implementation. Feature map region dropping is
implemented by interpolating the saliency mask to the spatial resolution of
each layer. The integration of patch dropping and feature dropping further
reduces model GOPs (in a hardware-friendly way) without compromising accuracy.

|  |  |
---|---|---|---
(a) Input frame | (b) Saliency Mask | (c) Masked Input | (d) GT
|  |  |
(e) FM(2,0) | (f) FM(2,0) w/ Mask | (g) FM(3,0) | (h) FM(3,0) w/ Mask
Figure 2. Visualizations of input frame patch dropping (first row) and feature
map patch dropping (second row). (a): An input frame with Ground Truth (GT)
object bounding box shown in (d). (b)-(c): An input saliency mask and the
masked input frame, where the patch dropping ratio is set to 20%, and the
discarded pixels are marked in black. (e)-(f): Feature maps in layer 0 of
block 2 within backbone (ResNet-50), denoted by FM(2,0), without and with
feature patch dropping mask, respectively. (g)-(h): Feature maps in layer 0 of
block 3, denoted by FM(3,0), without and with feature patch dropping mask,
respectively. We can observe that the proposed spatial data reduction approach
could accurately identify the redundant patch regions associated with
background scenes, e.g., sky and street.

###  3.3. Model pruning for computation reduction

The SOTA tracking model is overwhelmingly large. It has been long proven that
model pruning is an efficient way to reduce model complexity. For effective
and efficient hardware implementation, structured pruning is much more
beneficial compared to unstructured pruning, as the latter introduces
additional overhead such as sparse representation and format conversion.
Furthermore, unstructured pruning hinders full utilization of parallel
architectures (Anwar et al., 2017). Therefore, we advocate for and propose
structured pruning to fully accommodate our hardware design.

Table 1. The performance of irregularly pruned model (QDtrack) vs. pruning
ratio (0%, 80%, 85%, 90%). Here the sparse kernel ratio is given by the ratio
of pruned $3\times 3$ kernel weights over the total number of pruned weights.

Metrics | Dense Model | Global pruning
---|---|---
| 0%
---
| 80%
---
| 85%
---
| 90%
---
IDF1 ($\uparrow$) | 0.714 | 0.712 | 0.706 | 0.703
MOTA ($\uparrow$) | 0.637 | 0.631 | 0.627 | 0.624
| Sparse kernel ratio
---
0% | 36.59% | 34.44% | 32.53%

####  3.3.1. Kernel-wise sparsity embedded in irregular weight pruning.

The SOTA pruning method is iterative magnitude pruning (IMP) (Han et al.,
2015), which enables aggressive weight sparsity without loss of model
accuracy, as justified by the lottery ticket hypothesis (Frankle et al.,
2018). In Table 1, we apply IMP to the QDTrack model and indeed observe a
competitive MOT accuracy, measured by IDF1 and MOTA, as illustrated in Sec. 4.
However, the weight sparsity pattern achieved by IMP is unstructured, leading
to the compatibility issue with hardware. Spurred by this fact, we examine if
there exists a structured sparsity pattern hidden in the irregular weight
sparsity achieved by IMP. As shown in Table 1, the sparse kernel ratio, i.e.,
the ratio of pruned $3\times 3$ kernel weights over the total number of pruned
weights takes over 30% in different weight pruning ratios (from $80\%$ to
$90\%$). This implies that kernel-wise sparsity is embedded in irregular
weight pruning on QDTrack.

####  3.3.2. Integration with irregular pattern-aware pruning.

We leverage the irregular pattern-aware pruning method (Ma et al., 2020) to
enforce sparse weights that cannot be covered using kernel-wise sparse
patterns. The rationale behind such irregular pattern-aware pruning

Figure 3.  Pattern pruning.

is that the effective area of a convolution kernel often maintains specific
sparse patterns even if it does not yield kernel-wise sparsity, i.e., kernels
with all zero weights. Thus, we follow (Ma et al., 2020) to pre-define
irregular sparse patterns (Niu et al., 2020) for $3\times 3$ kernels, and
leverage them to conduct irregular but pattern-aware weight pruning. Note that
pruning using a fixed number of pre-defined sparse patterns facilitates
efficient hardware implementation compared with the conventional irregular
weight pruning.

Together with kernel-wise pruning, we propose the following hardware-aware
pruning method:

$\bullet$ Call IMP to achieve the irregular weight sparse mask
$\mathbf{M}_{\mathrm{IMP}}$;

$\bullet$ Extract kernel-wise sparsity mask $\mathbf{M}$ from
$\mathbf{M}_{\mathrm{IMP}}$;

$\bullet$ Conduct irregular pattern-aware pruning on remaining weights
identified by $(\mathbf{1}-\mathbf{M})$, then update $\mathbf{M}$ to obtain a
new mask $\mathbf{M}^{\prime}$;

$\bullet$ Retrain non-zero model weighs under fixed pruning mask
$\mathbf{M}^{\prime}$.

###  3.4. Hardware implementation and acceleration

####  3.4.1. Patch dropping compatible tile-based accelerator design.

Given HD input frames and the large ML model, the intermediate feature map for
even a single layer is too large to be fully stored in the device’s on-chip
memory. Therefore, it is required that the input images and feature maps be
partitioned into tiles and load/compute tile-by-tile, as shown in Fig. 4(a).
Fortunately, incorporating our patch dropping strategy, a tile-based
accelerator design is beneficial since it allows the accelerator to completely
skip the loading and computing of a tile if it is dropped at an early stage;
this not only reduces the computation latency but also significantly reduces
energy consumption, since off-chip data movement is much more power consuming
than on-chip computation (Chen et al., 2016).

####  3.4.2. Pruning-aware and multi-level parallelism.

To achieve extremely high parallelism and to simultaneously exploit structured
pruning, we propose a novel multi-level pruning-aware parallelization
architecture, as illustrated in Fig. 4. Consider one convolution layer as an
example, where the input is a data tile of size $T_{H}\times T_{W}\times
T_{C}$ along with a weight kernel of size $K\times K\times T_{C}$. The
commonly adopted parallel computation is along the $T_{C}$ dimension, i.e.,
$T_{C}$ channels are computed simultaneously (Hao et al., 2019). However, such
a parallelization scheme fails to employ the structured and pattern-based
pruning, since a pruned-away channel in a kernel still needs to be computed
with other channels in parallel. Therefore, in this work, we explore the
parallelism along $T_{H}$ and $T_{W}$ dimensions in two levels, and exploit
the sparsity along $T_{C}$ dimension.

First, row-level parallelism, as shown in Fig. 4(b), allows that at each clock
cycle, one input row is computed with $K$ rows of a kernel, followed by self-
accumulation for partial sum. Second, for column-level parallelism as shown in
Fig. 4(c), within the same clock cycle, $T_{W}$ multiplications are executed,
followed by a $K$-to-$1$ adder tree. Third, since the row-level and column-
level approaches have provided sufficient parallelism, we perform the
convolution channel-by-channel sequentially without sacrificing computation
latency. As shown in Fig. 4(d), if one channel of a filter is pruned, we skip
the data loading and computation for that channel, reducing energy and latency
almost linearly. Specifically, with structured pruning, if $p\%$ channels are
pruned away within a filter, the latency linearly reduces by $p\%$. Similarly,
for pattern-based pruning, our proposed design uses a fixed sparsity pattern
for the entire channel.

Figure 4. Pruning-aware parallelized hardware implementation.

####  3.4.3. Dataflow architecture using multi-FPGA for scalability.

To overcome the device and design tool scalability issue, as mentioned in Sec.
2, we propose a novel dataflow architecture using multiple FPGAs, as
illustrated in Fig. 5. Instead of executing the entire model on a giant FPGA
with extremely high parallelism, which is not feasible under HD frames
(limited by development toolkits), we dissemble the ML model into components
and map them to different FPGAs. As demonstrated in Fig. 5, multiple FPGAs
work in a dataflow manner: 1) The first FPGA executes the backbone ResNet-50
computation for one frame and meanwhile passes the intermediate feature maps
to the second FPGA. Once the first FPGA finishes the backbone computation for
one frame, it would start the computation for the next one. 2) The second FPGA
starts the computation as soon as it receives the first tile of intermediate
feature maps. Such a dataflow architecture not only significantly reduces end-
to-end latency but also increases the throughput by overlapping the
computations of multiple frames. In addition, it resolves the scalability
issue of both the device and the development tool as each FPGA can be
developed independently.

##  4\. Experiments

Figure 5. FPGA cluster for QDTrack implementation. FPGA 1 (Xilinx Alveo U50)
is responsible for ResNet50 backbone; FPGA 2 (U50) is responsible for FPN and
RPN modules; FPGA 3 (low-end ZCU104) processes the object tracking head. Under
each module, we annotate the GOPs before/after pruning.

###  4.1. Experiment setup

We conduct the MOT task under the BDD100K dataset (Yu et al., 2020), one of
the most representative and challenging self-driving car datasets that
incorporates geographic, environmental, weather diversity and intentional
occlusions. It contains 100K videos with 30 frames per second. In our
experiments, the training set has 1,400 videos and the validation set has 200
videos.

To evaluate the accuracy performance of our proposed tri-design approach, we
adopt the standard MOT metrics, ID F1 Score (IDF1) and Multi-Object Tracking
Accuracy (MOTA) . These are two aggregated metrics widely used to evaluate the
overall tracking performance by taking into account object detection accuracy
and tracking consistency. In addition, we use the number of Identity Switches
(IDSw) metric to evaluate the RL reward function in the first stage of data
complexity reduction.

To evaluate the efficiency of hardware implementation of our proposed tri-
design approach, we demonstrate the QDTrack model on a dataflow cluster
composed of three FPGAs: two mid-end Xilinx Alveo U50 boards (5,952 DSPs, 28
MB on-chip memory), and one low-end Xilinx ZCU104 (1,728 DSPs, 1.37 MB on-chip
memory) board. The accelerator is developed using Vitis HLS 2021.1 and
deployed using Vivado 2021.1.

###  4.2. Experiment results

Table 2. Performance comparison of our proposal with its variants and the QDTrack (Pang et al., 2021) baselines (implemented by GPU and FPGA respectively) under various metrics including accuracy (IDF1 and MOTA) and hardware metrics given by on-board latency in the unit of millisecond, effective frame rate (EFR) in the unit of FPS, and power in the unit of Watt. The data reduction is implemented using 40% temporal frame dropping together with 20% spatial patch dropping. The ‘$\times$’ indicates the case of no data/model compression, and ‘—’ means not applicable. The improvement of our proposed tri-design approach over GPU and FPGA baselines are summarized.  Methods | Data/model compression | Metrics |
---|---|---|---
| Data reduction
---
| Pruning
---
|  IDF1 ($\uparrow$)
---
|  MOTA ($\uparrow$)
---
|  Latency ($\downarrow$)
---
|  EFR ($\uparrow$)
---
|  Power ($\downarrow$)
---
|  Energy Efficiency ($\downarrow$)
---
QDTrack (GPU baseline) | $\times$ | $\times$ | 0.714 | 0.637 | 60.9 | 22.5 | 296 W | 13.2 J/frame
QDTrack on FPGA | $\times$ | $\times$ | 0.714 | 0.637 | 554.7 | 1.8 | 50.8 W | 28.2 J/frame
Variant: Frame + patch drop | (40%, 20%) | $\times$ | 0.71 | 0.628 | 443.8 | 2.3 | 50.8 W | 22.0 J/frame
Tri-design (ours) | (40%, 20%) | 90% | 0.704 | 0.617 | 44.4 | 37.6 | 50.8 W | 1.35 J/frame
Improv. over GPU baseline | — | — | -1.40% | -3.14% | 1.37$\times$ | 1.67$\times$ | 5.83$\times$ | 9.78$\times$
Improv. over FPGA baseline | — | — | -1.40% | -3.14% | 12.5$\times$ | 20.9$\times$ | — | 20.9$\times$

####  4.2.1. Results on video frame dropping.

We peer into the effect of the RL-based frame dropping model on our proposed
tri-design pipeline. We compare our frame dropping approach with two
baselines: uniform dropping and random dropping. The former one drops the
frames at a fixed time window for each video, and the latter one randomly
selects and drops frames in a video stream. Fig. 7 shows the frame dropping
performance at different drop rates from $15\%$ to $60\%$. As can be seen, our
RL-based model generally outperforms the random drop model at different
dropping ratios, indicating the effectiveness of the RL-based frame dropping
model. We also observe that the ‘sweet’ dropping ratio is achieved round
$40\%$, which strikes the most graceful balance between spatial reduction and
IDSw. In what follows, unless specified otherwise, we will assume that $40\%$
temporal frame reduction is applied to MOT.

####  4.2.2. Results on saliency-based image patch dropping.

We also studied the effect of spatial patch dropping in our framework. In Fig.
6, we compare the saliency-based method with random dropping versus different
patch drop sizes. As can be observed, the proposed saliency-score based method
consistently outperforms random dropping baseline at different patch size. We
also observe that IDF1 and MOTA have similar performance for sizes smaller
than 175 $\times$ 175, and they start to decrease if the drop patch size
becomes larger than 175 $\times$ 175\. Therefore, we set the drop patch size
as 60 $\times$ 60 in other experiments along with drop ratio as 20% so as to
maintain the IDF1 accuracy around 0.7.

####  4.2.3. Results on model pruning.

We compare our proposed model compression method with two commonly-used
pruning baselines: irregular weight magnitude pruning (IMP) and channel-wise
structured pruning. The results are demonstrated in Fig. 7. As we can see,
channel-wise pruning is quite sensitive to the pruning ratio. By contrast, our
proposed hardware-aware pruning method is capable of achieving MOT performance
similar to IMP across different pruning ratios. Additionally, in Table 2 we
show that when integrating with the spatiotemporal data reduction algorithms,
our proposed hardware-aware pruning method yields back-end acceleration
according to the performance of on-board latency and the effective video frame
rate.

|
---|---
Figure 6. Comparison between random patch dropping and our proposed saliency-
based dropping vs. different drop patch sizes.

####  4.2.4. Tri-design performance.

We next summarize the performance of integrated frame dropping, patch
dropping, and model compression. In Table 2, we demonstrate both MOT accuracy
and hardware implementation efficiency metrics of our proposed tri-design
approach. The baselines include the original QDTrack (Pang et al., 2021)
implementation on GPU Tesla V100 and the Xilinx FPGA Alveo U50 cluster
(Section 3.4, Fig. 5). We present the results of different configurations of
data complexity reduction and model weight pruning; each configuration results
in one variant of the proposed tri-design method. Accuracy. As we can see, the
use of proposed data complexity reduction strategy (i.e., the row ‘Variant’)
does not hamper the MOT accuracies in terms of IDF1 and MOTA scores.
Throughput. Comparing to GPU baseline, our tri-design approach reduces latency
by $1.37\times$ and improves the effective frame rate (EFR) by 1.67$\times$;
comparing to FPGA baseline, our tri-design together with proposed sparsity-
aware multi-level parallelism reduces latency by $12.5\times$ and improves the
EFR by $20.9\times$, demonstrating significant performance boost. Energy
Efficiency. Meanwhile, our approach shows remarkable energy efficiency. The
baseline GPU Tesla V100 has a power consumption of 296W, while our FPGA
cluster has only 50.8 W, leading to $5.83\times$ power reduction. Together
with 1.67$\times$ higher frame rate, the energy efficiency is 9.78$\times$
higher ($5.83\times 1.67$). To summarize, our approach outperforms GPU and
FPGA baselines in all hardware efficiency metrics with negligible accuracy
loss.

Figure 7.  Performance comparison: Left: Temporal filtering; Right: Pruning.

##  5\. Conclusions

We introduce data-model-hardware tri-design for MOT implementation on the edge
devices, which exploits aggressive data reduction, model compression, and
ultra-low-power hardware innovation for a hardware-aware ultra-light algorithm
development. We demonstrate the effectiveness of the proposed tri-design
through extensive experiments. Compared to SOTA MOT baseline, our approach on
Alveo U50 FPGAs achieves 12.5$\times$ latency reduction, 20.9$\times$
effective frame rate improvement, $5.83\times$ lower power, and $9.78\times$
better energy efficiency comparing to Tesla V100 GPU.

## References

  * (1)
  * Aarrestad et al. (2021) Thea Aarrestad et al. 2021\.  Fast convolutional neural networks on FPGAs with hls4ml.  _Machine Learning: Science and Technology_ 2, 4 (2021), 045015\.
  * Alom et al. (2018) Md Zahangir Alom et al. 2018\.  The history began from alexnet: A comprehensive survey on deep learning approaches.  _arXiv preprint arXiv:1803.01164_ (2018).
  * Anwar et al. (2017) Sajid Anwar et al. 2017\.  Structured pruning of deep convolutional neural networks.  _ACM JETC_ 13, 3 (2017), 1–18.
  * Chen et al. (2016) Yu-Hsin Chen et al. 2016\.  Eyeriss: An energy-efficient reconfigurable accelerator for deep convolutional neural networks.  _IEEE journal of solid-state circuits_ 52, 1 (2016), 127–138.
  * Finn et al. (2017) Chelsea Finn et al. 2017\.  Model-agnostic meta-learning for fast adaptation of deep networks. In _ICML_. PMLR, 1126–1135.
  * Frankle et al. (2018) Jonathan Frankle et al. 2018\.  The lottery ticket hypothesis: Finding sparse, trainable neural networks.  _arXiv preprint arXiv:1803.03635_ (2018).
  * Han et al. (2015) Song Han et al. 2015\.  Learning both weights and connections for efficient neural network.  _Advances in neural information processing systems_ 28 (2015).
  * Hao et al. (2019) Cong Hao et al. 2019\.  FPGA/DNN Co-Design: An Efficient Design Methodology for 1oT Intelligence on the Edge. In _2019 56th ACM/IEEE Design Automation Conference (DAC)_. 1–6.
  * Jin-Yu et al. (2009) Zhang Jin-Yu et al. 2009\.  Edge detection of images based on improved Sobel operator and genetic algorithms. In _2009 International Conference on Image Analysis and Signal Processing_. IEEE, 31–35.
  * Kalake et al. (2021) Lesole Kalake et al. 2021\.  Analysis based on recent deep learning approaches applied in real-time multi-object tracking: a review.  _IEEE Access_ 9 (2021).
  * Li et al. (2019) Qin Li et al. 2019\.  Implementing neural machine translation with bi-directional GRU and attention mechanism on FPGAs using HLS. In _Proceedings of the 24th Asia and South Pacific Design Automation Conference_. 693–698.
  * Lin et al. (2017) Tsung-Yi Lin et al. 2017\.  Feature pyramid networks for object detection. In _CVPR_.
  * Ma et al. (2020) Xiaolong Ma et al. 2020\.  An image enhancing pattern-based sparsity for real-time inference on mobile devices. In _ECCV_.
  * Muhawenayo et al. (2021) Gedeon Muhawenayo et al. 2021\.  Compressed object detection.  _arXiv preprint arXiv:2102.02896_ (2021).
  * Niu et al. (2020) Wei Niu et al. 2020\.  Patdnn: Achieving real-time dnn execution on mobile devices with pattern-based weight pruning. In _25th International Conference on Architectural Support for Programming Languages and Operating Systems_.
  * Pang et al. (2021) Jiangmiao Pang et al. 2021\.  Quasi-dense similarity learning for multiple object tracking. In _CVPR_. 164–173.
  * Ren et al. (2015) S. Ren et al. 2015\.  Faster r-cnn: Towards real-time object detection with region proposal networks.  _Advances in neural information processing systems_ 28 (2015).
  * Stäcker et al. (2021) Lukas Stäcker et al. 2021\.  Deployment of Deep Neural Networks for Object Detection on Edge AI Devices with Runtime Optimization. In _ICCV_. 1015–1022.
  * Wang et al. (2020a) Guangting Wang et al. 2020a.  Tracking by instance detection: A meta-learning approach. In _CVPR_. 6288–6297.
  * Wang et al. (2020b) Zhongdao Wang et al. 2020b.  Towards real-time multi-object tracking. In _ECCV_. Springer, 107–122.
  * Wu et al. (2018) Zhirong Wu et al. 2018\.  Unsupervised feature learning via non-parametric instance discrimination. In _CVPR_. 3733–3742.
  * Yu et al. (2020) Fisher Yu et al. 2020\.  Bdd100k: A diverse driving dataset for heterogeneous multitask learning. In _CVPR_. 2636–2645.
  * Yurtsever et al. (2020) Ekim Yurtsever et al. 2020\.  A survey of autonomous driving: Common practices and emerging technologies.  _IEEE access_ 8 (2020), 58443–58469.
  * Zhang et al. (2021) Jinming Zhang et al. 2021\.  A Low-Latency FPGA Implementation for Real-Time Object Detection. In _ISCAS_. 1–5.  https://doi.org/10.1109/ISCAS51556.2021.9401577
  * Zhang et al. (2018) Xiaofan Zhang et al. 2018\.  DNNBuilder: an Automated Tool for Building High-Performance DNN Hardware Accelerators for FPGAs. In _ICCAD_. 1–8.  https://doi.org/10.1145/3240765.3240801
  * Zhou et al. (2018) Kaiyang Zhou et al. 2018\.  Deep reinforcement learning for unsupervised video summarization with diversity-representativeness reward. In _AAAI_ , Vol. 32.


# Multicanonical reweighting for the QCD topological susceptibility

P. Thomas Jahn<EMAIL_ADDRESS>Parikshit M.
Junnarkar<EMAIL_ADDRESS>Guy D. Moore
<EMAIL_ADDRESS>Institut für Kernphysik
(Theoriezentrum), Technische Universität Darmstadt,
Schlossgartenstraße 2, D-64289 Darmstadt, Germany  Daniel Robaina
<EMAIL_ADDRESS>Max-Planck-Institut für Quantenoptik, Hans-
Kopfermann-Straße 1, D-85748 Garching, Germany

###### Abstract

We introduce a reweighting technique which allows for a continuous sampling of
temperatures in a single simulation and employ it to compute the temperature
dependence of the QCD topological susceptibility $\chi_{\mathrm{top}}$ at high
temperatures. The method determines the ratio of $\chi_{\mathrm{top}}$ between
any two temperatures within the explored temperature range. We find that the
results from the method agree with our previous determination and that it is
competitive with but not better than existing methods of determining
$\mathrm{d}\chi_{\mathrm{top}}/\mathrm{d}T$. The method may also be useful in
exploring the temperature dependence of other thermodynamical observables in
QCD in a continuous way.

##  I Introduction

The axion solution to the strong CP problem, proposed more than four decades
ago Peccei and Quinn (1977); Weinberg (1978); Wilczek (1978), solves the fine
tuning problem of the smallness of the $\theta$ parameter in QCD by
introducing a new light ($<\mathrm{meV}$) pseudo Goldstone boson in the
Standard model. Soon after it was proposed, it was also realized that long-
wavelength axions would be produced abundantly early in the Big Bang, and can
serve as a candidate for dark matter Abbott and Sikivie (1983); Preskill  _et
al._ (1983). These proposals have heightened interest both in searching
experimentally for the axion, and in better understanding its phenomenology
and cosmological history; for a review of these topics see for instance Ref.
Irastorza and Redondo (2018). An important input for the cosmological
production of axions and for their contemporary properties is the QCD
topological susceptibility $\chi_{\mathrm{top}}(T)$, which determines the
axion mass:

$m^{2}_{\mathrm{A}}(T)=\frac{\chi_{\mathrm{top}}(T)}{f^{2}_{\mathrm{A}}}.$ (1)

Here $m_{\mathrm{A}}(T)$ is the temperature dependent axion mass and
$f_{\mathrm{A}}$ is the axion decay constant. The low-temperature value of
$\chi_{\mathrm{top}}$ is well determined Grilli di Cortona  _et al._ (2016);
Gorghetto and Villadoro (2019), but the high-temperature regime determines
axion production efficiency; it is particularly important to determine
$\chi_{\mathrm{top}}(T)$ in the temperature regime from 500 MeV to 1200 MeV
Moore (2018). In this temperature range, this can only be achieved by
nonperturbative lattice investigations Borsanyi  _et al._ (2016a). The most
straightforward lattice methods, based on brute-force sampling of the gauge
field configuration space, face difficulties in this temperature range because
topologically nontrivial gauge field configurations become very rare, leading
to a loss of statistical power. New methodologies are needed to overcome this
problem. In recent years one such methodology has been developed Borsanyi  _et
al._ (2016b); Frison  _et al._ (2016) which provides access to high
temperatures. In these calculations, it was shown that the difference of the
expectation of the QCD action in two topological sectors can provide a
determination of $\mathrm{d}(\ln\chi_{\mathrm{top}})/\mathrm{d}(\ln T)$, which
can be integrated to provide the temperature dependence of
$\chi_{\mathrm{top}}$.

Alternatively, one can approach the problem at a fixed, high temperature by
reweighting between topological sectors. A first attempt, based on a fixed
guess for the reweighting function Bonati  _et al._ (2018), explored
temperatures around 500 MeV. Another method, in which the reweighting function
is determined dynamically via an iterative self-consistent technique, was
introduced in Jahn  _et al._ (2018) and further improved in Jahn  _et al._
(2020), where it was shown to be effective in the pure-glue theory up to at
least $7\,T_{\mathrm{c}}$.

Our main motivation for this work is to explore whether a new technique might
improve existing methods for determining $\chi_{\mathrm{top}}$. The method
introduced here is related to but distinct from the technique of Refs Borsanyi
_et al._ (2016b); Frison  _et al._ (2016). We train a single Markov-chain
Monte-Carlo simulation to explore a wide range of temperatures in a detailed-
balance respecting way by replacing the weighting function $\exp[-\beta S]$
with $\exp[-W(S)]$, where $W(S)$ is established by an iterative procedure such
that the resulting ensemble can be reweighted to describe any temperature in a
relatively wide range. We do this separately in the $Q=0$ and $|Q|=1$ (non-
topological and instanton-number $\pm 1$) sectors, and combine with a
determination of $\chi_{\mathrm{top}}(T)$ at one temperature to determine
$\chi_{\mathrm{top}}$ across the full accessible temperature range. We show
that the approach has a similar efficiency to Borsanyi  _et al._ (2016b), with
both (slight) advantages and disadvantages relative to their approach. The
technique can be extended to include fermions if the line of constant physics
(lattice $m(\beta)$) is known. However we find that our approach developed in
Ref. Jahn  _et al._ (2020) seems to afford better numerical efficiency.

In the next section we review the definition of topological susceptibility and
lay the groundwork for our technique. Section III introduces our technique to
carry out a detailed-balance preserving Markov chain over a range of
temperatures. Then Section IV presents our results and Section V closes with a
discussion.

##  II Topological Susceptibility

The topological susceptibility is defined as

$\chi_{\mathrm{top}}(T)=\frac{T}{V}\;\langle
Q^{2}\rangle=\frac{T}{V}\int\mathrm{d}^{4}x\;\mathrm{d}^{4}y\;\langle
q(x)q(y)\rangle,$ (2)

where $Q=\int\mathrm{d}^{4}x\>q(x)$ and $q$ is the topological charge density,
$V$ the spatial volume and $T$ the temperature. In the continuum $Q$ always
takes an integer value while on the lattice one requires a refined definition
of $Q$. Considering the continuum integer topological charge $Q$ one can write
the partition function as a sum over topological sectors:

$Z=\int\mathcal{D}Ae^{-\beta
S}=\sum_{N\in\mathcal{Z}}\int\mathcal{D}Ae^{-\beta
S}\delta(Q-N)\equiv\sum_{N}Z_{N}\,.$ (3)

The susceptibility is then given by :

$\displaystyle\chi_{\mathrm{top}}(T)$ $\displaystyle=$
$\displaystyle\frac{T}{V}\frac{\int\mathcal{D}Ae^{-\beta
S}Q^{2}}{\int\mathcal{D}Ae^{-\beta S}}$ (4) $\displaystyle=$
$\displaystyle\frac{T}{V}\frac{\sum_{N}N^{2}Z_{N}}{\sum_{N}Z_{N}}\,.$

At low temperatures and/or large volumes, this sum will have important
contributions from many $N$ values. But at a sufficiently high temperature,
such that $V/T\;\chi_{\mathrm{top}}(T)\ll 1$, the numerator is dominated by
$N=1$ and $N=-1$ and the denominator is dominated by $N=0$. Renaming
$Z_{1}+Z_{-1}\to Z_{1}$ (the part of the partition function where $Q^{2}=1$,
that is, $Q=\pm 1$) and considering a finite volume we find

$\displaystyle\chi_{\mathrm{top}}(T)\simeq\frac{T}{V}\frac{Z_{1}}{Z_{0}}$
$\displaystyle=\frac{1}{a^{4}N_{x}N_{y}N_{z}N_{\tau}}\frac{Z_{1}}{Z_{0}}$
$\displaystyle\chi_{\mathrm{top}}(\beta)a^{4}(\beta)$
$\displaystyle=\frac{1}{V_{L}}\frac{Z_{1}}{Z_{0}}.$ (5)

Here we have re-expressed the temperature and volume as they would appear in a
lattice calculation, with lattice spacing $a$ and lattice extents
$(N_{\tau},N_{x},N_{y},N_{z})$ in the temporal and the three space directions,
and $V_{L}=N_{\tau}N_{x}N_{y}N_{z}$ the number of lattice sites. The last line
emphasizes that the lattice spacing is a function of the lattice gauge
coupling $\beta=6/g^{2}_{\mathrm{latt}}$. The nontrivial relation between
these two quantities is determined by a scale setting measurement.

In this work, we will present a method which determines the ratio of
susceptibilities at two specified temperatures. With a fixed lattice extent in
the spatial and temporal directions, the two temperatures will correspond to
two gauge couplings $\beta_{\mathrm{h}}$ and $\beta_{\mathrm{c}}$ (h for hot
and c for cold) which will in turn correspond to two different lattice
spacings $a(\beta_{\mathrm{h}})$ and $a(\beta_{\mathrm{c}})$. The ratio of
susceptibilities at two temperatures will then be given as:

$\frac{\chi_{\mathrm{top}}(\beta_{\mathrm{h}})a^{4}(\beta_{\mathrm{h}})}{\chi_{\mathrm{top}}(\beta_{\mathrm{c}})a^{4}(\beta_{\mathrm{c}})}=\frac{Z_{1}(\beta_{\mathrm{h}})}{Z_{1}(\beta_{\mathrm{c}})}\;\frac{Z_{0}(\beta_{\mathrm{c}})}{Z_{0}(\beta_{\mathrm{h}})}.$
(6)

Our method will also allow for the determination of the above ratio for any
pair of temperatures within the prespecified range. We will do so by computing
the ratio of partition functions using two Monte-Carlo simulations: one that
works within the $Q=1$ topological sector and one which works within the $Q=0$
sector, but with each simulation exploring the full range of $\beta$ values,
as we describe in the next section.

##  III Temperature reweighting

Figure 1: A cartoon of the Monte-Carlo sampling of $W[S]$. The green line
represents a $W$ choice which will sample one $\beta$ value, while the blue
line will sample all $\beta$ values between $\beta_{\mathrm{max}}$ and
$\beta_{\mathrm{min}}$.

In a given Monte-Carlo simulation, we intend to sample a range of temperatures
continuously over a prespecified range in pure glue QCD regularized on a
finite space-time lattice. We will accomplish this using a reweighting method
very similar to the old proposal of Berg and Neuhaus Berg and Neuhaus (1991).
In this section we first show how a reweighted Monte-Carlo can be used to
determine the susceptibility; then we show how the reweighted Monte-Carlo can
be carried out; and finally we show how to determine the reweighting function
itself. In practice, the numerical work proceeds in exactly the opposite
order.

###  III.1 Susceptibility using reweighting

In the standard Monte-Carlo simulation one evaluates the partition function
$\mathcal{Z}(\beta)$ shown below at a particular gauge coupling $\beta$ by
generating gauge configurations with the probability distribution
$\mathrm{d}P[U]$ as

$\mathcal{Z}(\beta)=\int\mathcal{D}U\;e^{-\beta
S},\quad\mathrm{d}P[U]=\frac{\mathrm{d}U\;e^{-\beta S}}{\mathcal{Z}(\beta)}.$
(7)

Such a Monte-Carlo simulates a single temperature which is controlled by the
gauge coupling $\beta$ since the lattice dimensions are fixed. In order to
simulate a range of $\beta$ values and therefore a range of temperatures, the
sampling weight in the usual Monte-Carlo in Eq. (7) is replaced with

$\int\mathcal{D}U\;e^{-\beta S}\longrightarrow\int\mathcal{D}U\;e^{-W[S]}.$
(8)

The weight function $W[S]$ is a general function of the action $S$, and its
derivative can be interpreted, approximately, as the gauge coupling $\beta$ to
use at a given action value:

$W^{\prime}[S]\equiv\frac{\mbox{d}W[S]}{\mbox{d}S}=\beta[S].$ (9)

In this sense, the usual Monte-Carlo has simply $W[S]=\beta S$ with the
constant slope simulating the fixed coupling. We show how to cover a range of
temperatures schematically in Figure 1. To get a fixed temperature one picks
$W^{\prime}[S]=\beta$ as indicated by the green line in the figure; to cover a
range of temperatures, one instead uses the blue line in the figure, for which
$W^{\prime}[S]$ varies from $\beta_{\mathrm{max}}$ to $\beta_{\mathrm{min}}$
as $S$ changes from $S_{\mathrm{min}}$ to $S_{\mathrm{max}}$; here
$S_{\mathrm{min}}$ and $S_{\mathrm{max}}$ are the average values of $S$ in
simulations with coupling strengths $\beta_{\mathrm{max}}$ and
$\beta_{\mathrm{min}}$ respectively. For $S_{0}$ an action value between these
limiting choices, $W^{\prime}[S_{0}]$ is chosen to be the $\beta$ value for
which the expectation value of the action would be $S_{0}$. Outside the range
of interest, the weight function $W[S]$ then simulates two different constant
gauge couplings $\beta_{\mathrm{max}}$ and $\beta_{\mathrm{min}}$. With this
choice of $W[S]$, a single Markov chain Monte-Carlo will sample the entire
range of $S$ values between $S_{\mathrm{min}}$ and $S_{\mathrm{max}}$ rather
uniformly. The function $W^{\prime}[S]$ is not known a priori; we will return
to its determination shortly, but will first discuss how a simulation with
such a weight function can be carried out and used.

Once $W[S]$ is known, a Monte-Carlo simulation with $W[S]$ generates an
ensemble which can be reweighted to determine expectation values at a given
$\beta$ value $\beta_{0}$ via

$\displaystyle\mathcal{Z}(\beta_{0})$ $\displaystyle=$
$\displaystyle\int\mathcal{D}Ue^{-W[S]}\;e^{+W[S]-\beta_{0}S}$
$\displaystyle\propto$ $\displaystyle\sum_{i}e^{W[S_{i}]-\beta_{0}S_{i}}\,,$

where $i$ indexes the sampled configurations. In practice this partition
function is used to determine expectation values for operators via

$\langle\mathcal{O}\rangle=\mathcal{Z}^{-1}\int\mathcal{D}Ue^{-\beta_{0}S}\mathcal{O}\simeq\frac{\sum_{i}e^{W[S_{i}]-\beta_{0}S_{i}}\mathcal{O}_{i}}{\sum_{i}e^{W[S_{i}]-\beta_{0}S_{i}}}\,.$
(11)

By establishing two reweighting functions $W[S]$ and $W_{Q}[S]$ for the $Q=0$
and $|Q|=1$ ensembles respectively, we can generate multi-temperature
ensembles, labeled by $i$ and $iQ$, which respectively sample the $Q=0$ and
the $|Q|=1$ ensembles across temperatures. The ratio needed in Eq. (6) is then
given by

$\displaystyle\frac{\chi_{\mathrm{top}}(\beta_{\mathrm{h}})a^{4}(\beta_{\mathrm{h}})}{\chi_{\mathrm{top}}(\beta_{\mathrm{c}})a^{4}(\beta_{\mathrm{c}})}$
(12) $\displaystyle=$
$\displaystyle\mathopen{}\mathclose{{}\left(\frac{\mathopen{}\mathclose{{}\left(\sum_{iQ}e^{W_{Q}[S_{iQ}]}e^{-\beta_{\mathrm{h}}S_{iQ}}}\right)\mathopen{}\mathclose{{}\left(\sum_{i\vphantom{Q}}e^{W[S_{i}]}e^{-\beta_{\mathrm{c}}S_{i}}}\right)}{\mathopen{}\mathclose{{}\left(\sum_{iQ}e^{W_{Q}[S_{iQ}]}e^{-\beta_{\mathrm{c}}S_{iQ}}}\right)\mathopen{}\mathclose{{}\left(\sum_{i\vphantom{Q}}e^{W[S_{i}]}e^{-\beta_{\mathrm{h}}S_{i}}}\right)}}\right)\,.$

There are two subtleties associated with this expression. The first is that
Eq. (III.1) only determines the partition function up to an overall
multiplicative factor. However, this multiplicative factor cancels between the
numerator and denominator expressions computed from the same sample. The
second subtlety is that the partition function $\mathcal{Z}$ also has severe
lattice-spacing dependent renormalizations, which we cannot easily compute.
Fortunately, these cancel because each lattice spacing occurs once in the
numerator and once in the denominator in Eq. (12).

###  III.2 Update algorithm with $W[S]$

In this section we explain the algorithm to perform a Monte-Carlo simulation
with weight function $e^{-W[S]}$. Our aim is to generate a sampling with
probability distribution

$\mbox{d}P[U]=\frac{e^{-W[S[U]]}}{\int\mathcal{D}Ue^{-W[S[U]]}}$ (13)

where $S[U]$ is the standard lattice gauge action and we assume that $W[S[U]]$
is a known differentiable function of the action $S$. Simulating such a weight
requires a slight modification of the standard Hybrid Monte-Carlo (HMC)
algorithm Duane  _et al._ (1987).111It is also straightforward to use a
mixture of heatbath and overrelaxation steps. However this approach does not
generalize to the unquenched theory, so we concentrate on the HMC approach.

As in the standard HMC algorithm, we introduce canonical momenta $\pi_{\mu}$
for the link variables $U_{\mu}$, and define a Hamiltonian for this system as

$\displaystyle\mathcal{H}\big{(}\pi,U\big{)}$ $\displaystyle\equiv$
$\displaystyle\sum_{\mu,x}\frac{1}{2}(\pi_{\mu}(x))^{2}+W[S_{U}]$ (14)
$\displaystyle S_{U}$ $\displaystyle=$
$\displaystyle\sum_{\begin{picture}(4.0,4.0)\put(0.0,0.0){\line(1,0){4.0}}\put(0.0,4.0){\line(1,0){4.0}}
\put(0.0,0.0){\line(0,1){4.0}}\put(4.0,0.0){\line(0,1){4.0}}\end{picture}}\mathopen{}\mathclose{{}\left(1-\frac{1}{3}\mathrm{Tr}\>\begin{picture}(6.0,6.0)\put(0.0,0.0){\line(1,0){6.0}}\put(0.0,6.0){\line(1,0){6.0}}
\put(0.0,0.0){\line(0,1){6.0}}\put(6.0,0.0){\line(0,1){6.0}}\end{picture}}\right)$
(15)

where $S_{U}$ in Eq. (15) is the standard Wilson gauge action written without
the gauge coupling $\beta$ prefactor. The standard HMC would use the same
Hamiltonian but with $W[S_{U}]$ replaced by $\beta S_{U}$, that is, it would
use a strictly linear function for $W[S_{U}]$.

A single HMC update trajectory consists of the standard steps:

  1. 1.

Picking a random canonical momentum $\pi_{\mu}(x)$ from a Gaußian ensemble
independently for each of the elements of the Lie algebra.

  2. 2.

Solving the following Hamilton equations of motion (shown here schematically):

$\displaystyle\frac{\mathrm{d}U}{\mathrm{d}t}$ $\displaystyle=-i\pi U\,,$ (16)
$\displaystyle\frac{\mathrm{d}\pi}{\mathrm{d}t}$
$\displaystyle=iU^{\dagger}\frac{\partial W[S_{U}]}{\partial
U}=iU^{\dagger}\frac{\mathrm{d}W[S_{U}]}{\mathrm{d}S_{U}}\frac{\partial
S_{U}}{\partial U}$ (17)

for a total time $t_{0}$. The derivative with respect to $U_{\mu}$ is a Lie
derivative and in this sense these equations are schematic representation.
Here the time $t$ is a fictitious variable under which the Hamiltonian
$\mathcal{H}$ is conserved.

These Hamiltonian equations are discretized using a time-symmetric solver such
as the leapfrog or Omelyan algorithms. Under these algorithms, one iteratively
solves Eq. (16) for all link variables $U$ at fixed $\pi$, and then solves Eq.
(17) for all $\pi$ variables at fixed $U$. Before applying Eq. (17), we must
compute $S_{U}$ and use the (instantaneous) value of $W^{\prime}[S_{U}]$ in
place of the usual factor of $\beta$ at each time step during the update.

The use of a time-symmetric algorithm is essential, since it ensures the
property that, if the pair $(U_{\mathrm{i}},\pi_{\mathrm{i}})$ (i for initial)
is carried to $(U_{\mathrm{f}},\pi_{\mathrm{f}})$ under the update algorithm,
then the pair $(U_{\mathrm{f}},-\pi_{\mathrm{f}})$ is carried to
$(U_{\mathrm{i}},-\pi_{\mathrm{i}})$ up to roundoff error effects. This is
sufficient to ensure that the algorithm converges, in a Fokker-Planck sense,
to the probability distribution $\exp[-\mathcal{H}]$ provided that we also
include a Metropolis accept/reject step.

  3. 3.

In adding the Metropolis step, the change in the Hamiltonian
$\Delta\mathcal{H}=\mathcal{H}(U_{\mathrm{f}},P_{\mathrm{f}})-\mathcal{H}(U_{\mathrm{i}},P_{\mathrm{i}})$
is compared to a random number drawn uniformly from the interval $[0,1]$:
$\mathcal{R}[0,1]$. Whenever $e^{\Delta\mathcal{H}}<\mathcal{R}[0,1]$, we
accept the change, and proceed with $U_{\mathrm{f}}$ as our new configuration.
Otherwise we revert to $U_{\mathrm{i}}$, that is, we reject the update.

The only differences with respect to the standard HMC algorithm Duane  _et
al._ (1987) are the use of $W^{\prime}[S_{U}]$ the “instantaneous $\beta$
value” in place of $\beta$ in Eq. (17) and the use of $W[S_{U}]$ in place of
$\beta S_{U}$ in the Metropolis accept-reject step. Both modifications are
compatible with the time-symmetry of the update algorithm and therefore
preserve detailed balance. With these modifications, the HMC algorithm now
generates the desired probability distribution.

For the case of a $Q=1$ simulation, an additional accept-reject step is
needed, in which we check to see whether the configuration has fallen down
into the $Q=0$ sector and reject the update if this is the case. In practice
we can buffer every $N$’th configuration and only perform this step after
every $N$ HMC update steps, reverting to the last buffered configuration when
the check fails. We define $Q$ as the lattice sum of an $a^{2}$-improved
topological density definition after $\tau_{\mathrm{F}}=2.4a^{2}$ units of
Wilson flow, as in Jahn  _et al._ (2020). We find that values of $N=5$ or
$N=10$ are adequate to preserve a good acceptance rate.

Lastly we remark on the optimal length of the individual trajectories. The
figure of merit for trajectory length is the mean-squared change in $S_{U}$
per unit numerical effort. The numerical effort is approximately linear in
trajectory length $t_{0}$. For a short trajectory, $\Delta S_{U}$ is linear,
and $(\Delta S_{U})^{2}$ quadratic, in trajectory length; but beyond a certain
(fairly short) trajectory length the action change saturates. Therefore we
started with a study of how $(\Delta S_{U})^{2}$ varies with trajectory length
and chose the value which maximizes $(\Delta S_{U})^{2}/t_{0}$; $t_{0}\simeq
0.75a$. A single trajectory leads to a change of $(\Delta
S_{U}/S_{U})^{2}\simeq 3/N_{\mathrm{dof}}$ where $N_{\mathrm{dof}}=24V_{L}$ is
the number of lattice degrees of freedom (3 polarizations and 8 colors per
site). Therefore, for a well-chosen $W[S_{U}]$ function, since changes to
$S_{U}$ accumulate in a Brownian fashion, the number of updates needed to
explore the full $\beta$ range is of order

$N_{\mathrm{updates}}\sim
N_{\mathrm{dof}}\,\ln^{2}(\beta_{\mathrm{max}}/\beta_{\mathrm{min}}).$ (18)

###  III.3 Choice and determination of $W[S]$

Now we return to the question of how to determine the weight function $W[S]$.
We start by choosing the range of $\beta$ values we want to explore,
$\beta\in[\beta_{\mathrm{min}},\beta_{\mathrm{max}}]$. A short fixed-$\beta$
Markov chain establishes $S_{\mathrm{max}}=\langle
S_{U}\rangle_{\beta_{\mathrm{min}}}$ and $S_{\mathrm{min}}=\langle
S_{U}\rangle_{\beta_{\mathrm{max}}}$. We then follow our procedure in Jahn
_et al._ (2019, 2020) and choose a discrete set of values $S_{i},i=(0,\ldots
N_{i})$ with $S_{0}=S_{\mathrm{min}}$ and $S_{N_{i}}=S_{\mathrm{max}}$;
$W[S_{U}]$ will be determined by its values $W[S_{i}]$. However, because the
update described above works best when both $W[S_{U}]$ and $W^{\prime}[S_{U}]$
are continuous functions, we shall interpolate $W[S_{U}]$ between these points
using a cubic spline function, rather than using a piecewise linear function
as in our previous papers. We extend $W[S]$ above $S_{\mathrm{max}}$ by
choosing $W^{\prime}[S>S_{\mathrm{max}}]=\beta_{\mathrm{min}}$ and similarly
we set $W^{\prime}[S<S_{\mathrm{min}}]=\beta_{\mathrm{max}}$; these are also
the values of the first derivatives used in completing the definition of the
spline function.

We will use the same automated improvement scheme to determine $W[S_{i}]$ as
we introduced in Ref. Jahn  _et al._ (2019) (see also Laine and Rummukainen
(1998)). We generate a Markov chain using the update approach of the previous
subsection, but after each update, we adjust the function $W[S_{U}]$ so as to
make the current $S$ value less likely (implying that the current $S$ is
oversampled). This is done by

  1. 1.

determining $S_{i},S_{i+1}$ such that the current $S$ value lies between them,
$S_{i}<S_{U}<S_{i+1}$

  2. 2.

increasing $W[S_{i}]$ and $W[S_{i+1}]$ by $\Delta
W[S_{i}]=\frac{s_{r}(S_{i+1}-S_{U})}{(S_{i+1}-S_{i})}$ and $\Delta
W[S_{i+1}]=\frac{s_{r}(S_{U}-S_{i})}{(S_{i+1}-S_{i})}$. Here $s_{r}$ is an
update strength which we explain next.

  3. 3.

If $S_{U}$ is out of range, $W[S]$ is not updated.

  4. 4.

If $i=0$ or $i+1=N_{i}$ so we are in the first or last interval, the boundary
value is updated with double strength (since it is only updated half as often
as other values).

The value $s_{r}$ is initially chosen so that, in the time it takes for the
Monte-Carlo to go from the top to the bottom and back (see Eq. (18)), the
$W[S_{i}]$ will change by of order 100. Each time the $S$-value makes its way
from the first interval to the last and back (which we call a “sweep”), we
reduce $s_{r}$; at first we reduce it by a factor of 2, but after the average
$W[S_{i}]$ changes by less than 1 per sweep, we change it by a reduced amount.
The update ends when five sweeps change the average $W[S_{i}]$ by a total of
less than 1.

This approach is inefficient if the initial $W$ function is very far from its
final form. Therefore we improve the initial guess from the simplest approach
(that $W^{\prime}[S]$ is linear). Instead, we choose an intermediate $\beta$
value $\beta_{\mathrm{mid}}$ and perform a fixed-$\beta$ Monte-Carlo
calculation to determine the associated $S_{\mathrm{mid}}$ value. We then fit
$W^{\prime}$ to a quadratic, with values
$(\beta_{\mathrm{max}},\beta_{\mathrm{mid}},\beta_{\mathrm{min}})$ at the
points $(S_{\mathrm{min}},S_{\mathrm{mid}},S_{\mathrm{max}})$, integrate, and
use this to determine starting guesses for the $W[S_{i}]$. With this approach,
the initial and final $W$ functions differ by less than 100, as shown in
Figure 2, which displays the difference between the initial and final $W$
function choice for the specific lattice and $\beta$ range described in the
next section.222We could further improve the initial guess by using more
intermediate values; but using enough values, with precise enough $\langle
S\rangle$ determinations, to determine $W[S]$ to better than $\pm 1$ costs the
same as the $W$ refinement algorithm described here. The complete built
$W^{\prime}$ is shown in Figure 3.

Figure 2:  Difference between the final determined $W[S]$ and the initial
guess. Figure 3: Final choice of $W^{\prime}[S]$ for the lattice study
described in Section IV.

Finally, one must determine $W[S_{U}]$ for the $Q=1$ sector. Here we can take
as an initial guess the $W[S_{U}]$ value determined in the $Q=0$ sector. To
further refine this guess, we shift it by

$W_{Q=1}[S]=W_{Q=0}[S]-11\ln(T/T_{0})$ (19)

where $T$ is the temperature associated with the $\beta$ value described by
the slope $W^{\prime}$ using a scale-setting relation between lattice coupling
$\beta$ and temperature $T$, and $T_{0}$ is a reference temperature which
could for instance be the temperature at $\beta_{\mathrm{min}}$. The factor 11
is the expected temperature dependence of $\chi_{\mathrm{top}}a^{4}$ when the
lattice spacing $a$ varies as $1/T$, at leading perturbative order. Again, Eq.
(19) is only used to refine the initial guess for $W_{Q=1}[S]$; we again
perform an automated $W$-changing Markov chain to improve this guess; however
the initial $s_{r}$ value can be made much smaller.

After the $W[S_{U}]$-setting Markov chains are completed, we freeze the values
of $W_{Q=0,Q=1}[S]$ and use them in detailed-balance respecting Markov chains
which we will use to determine the susceptibility as described in Subsection
III.1.

##  IV Results

Table 1:  Numerical cost of the simulation on a lattice four-volume $V=10\times 32^{2}\times 36$ with $\beta_{\mathrm{min}}=6.9076$ and $\beta_{\mathrm{max}}=8.01951$. Procedure | $Q$ | HMC Trajectories | Sweeps
---|---|---|---
Building | 0 | $4\times 10^{6}$ | 22
Building | 1 | $8\times 10^{6}$ | 17
Measurements | 0 | $6.8\times 10^{6}$ | 45
Measurements | 1 | $4.9\times 10^{6}$ | 35
Figure 4: A histogram of the $S$ distribution of our Markov chain samples in
the $Q=0$ sector (left) and the $Q=1$ sector (right).

In this section we present results of a simulation with aforementioned $W(S)$
sampling. Our goal is only to test the method; we will not study multiple
lattice spacings to attempt a continuum limit. Instead, we choose one lattice
geometry from among the geometries studied in Ref. Jahn  _et al._ (2018,
2020), and which we can therefore use to compare our results for
susceptibility ratios to those obtained via a competing technique. Our goal is
to establish which approach, the one described here or the one described in
Ref. Jahn  _et al._ (2020), is more efficient at establishing the topological
susceptibility at high temperature.

We choose to investigate a lattice with temporal extent $N_{\tau}=10$ and
spatial extent $32^{2}\times 36$, with
$(\beta_{\mathrm{min}},\beta_{\mathrm{max}})=(6.9076,8.01951)$, which
corresponds, according to the scale setting calculation of Ref. Burnier  _et
al._ (2017) which we will use throughout, to
$T_{\beta_{\mathrm{min}}}=2.5\,T_{\mathrm{c}}$ and
$T_{\beta_{\mathrm{max}}}=9.4\,T_{\mathrm{c}}$. This scale-setting relation
involves applying a fit to scale-setting data beyond the range where the
reference has performed simulations, that is, an extrapolation, which means it
may not be absolutely trustworthy. However, by expressing our results in terms
of $\chi_{\mathrm{top}}(\beta)a^{4}(\beta)$, we can remain agnostic about the
relation between $T$ and $\beta$ and just determine how accurately we can
determine $\chi_{\mathrm{top}}a^{4}$ as a function of $\beta$ for this
specific lattice geometry. Getting continuum results over a range of
temperatures will then of course require multiple lattice $N_{\tau}$ values
and a reliable scale setting. However, for the time being our goal is just to
evaluate the precision-to-numerical-cost ratio of the technique, so we leave
this problem for later. The number of updates used in each procedure are
listed in Table 1. Note that an unfortunate choice of $s_{r}$ made the $Q=1$
building unnecessarily inefficient; a more careful choice should have required
a fraction as many trajectories, so that the trajectory count would be
dominated by the measurements.

We begin with a check that our $W[S]$ function correctly generates a rather
uniform sample of configurations across the desired $\beta$ range. We
investigate this by plotting a histogram of the $S$ values measured during the
sampling Markov chain, both for the $Q=0$ and the $Q=1$ ensembles, shown in
Figure 4. The sample is adequately uniform.

$T/T_{\mathrm{c}}$ $\beta$
$\ln(\chi_{\mathrm{top}}a^{4}(T)/\chi_{\mathrm{top}}a^{4}(2.8\,T_{\mathrm{c}}))$
1-$\sigma$ stat. err 3.5 7.1771 -2.25 0.17 4.0 7.2885 -3.68 0.20 5.0 7.4764
-5.79 0.24 7.0 7.7629 -9.39 0.30 9.0 7.9788 -12.20 0.36

Table 2: Several temperatures, the corresponding $\beta$ values using the
scale setting of Ref. Burnier  _et al._ (2017), and our results for the log
susceptibility ratio and its 1-sigma statistical error.

Our main results are presented in Figure 5, which shows how
$\chi_{\mathrm{top}}a^{4}$ changes as a function of $\beta$ across the range
we study, evaluated from our Markov chains using Eq. (12). We have chosen to
use a value near the beginning of the $\beta$ range ($T=2.8\,T_{\mathrm{c}}$)
as the low temperature and to express all other temperatures in relation to
this one. The figure, and further data presented in Table 2, show that the
error bars are smallest between nearby $\beta$ values and grow to around $\pm
0.35$ in $\ln(\chi_{\mathrm{top}})$ for the widest-separated temperatures. Our
results agree within error bars of the results in Ref. Jahn  _et al._ (2020)
for those values where they are directly comparable.

##  V Discussion

Figure 5: Results of the ratio in Eq. 12

We have shown that the method we propose can successfully find the $\beta$
dependence of $\chi_{\mathrm{top}}a^{4}$, and therefore the temperature
dependence of the susceptibility if the line of constant physics (that is,
$a(\beta)$) is known. This determines $\chi_{\mathrm{top}}$ over a range of
temperatures if it is known at the lowest temperature, which is where it is
most easily determined by other approaches.

There are two key questions. Is it more or less effective than the rather
similar approach of Refs Borsanyi  _et al._ (2016b); Frison  _et al._ (2016)?
And how does it compare with the approach of Ref. Jahn  _et al._ (2020)?

The approach of Ref. Borsanyi  _et al._ (2016b); Frison  _et al._ (2016)
computes $\mathrm{d}\ln(\chi_{\mathrm{top}}a^{4})/\mathrm{d}\beta$ at several
$\beta$ values, which it integrates to determine
$\chi_{\mathrm{top}}a^{4}(\beta)$. We review this approach and compare it to
our own in an appendix. To summarize, in a high-statistics determination, the
approaches have essentially the same numerical precision. However, if lower
precision is desired, the numerical cost associated with building the $W(S)$
functions in our approach is essentially “dead weight” which does not
contribute to the statistical power. The other approach does not suffer from
this problem and so it is more efficient for a low-statistics determination.
Our approach has the advantage that it automatically includes all intermediate
temperatures, while the alternative bases the determined
$\mathrm{d}\ln(\chi_{\mathrm{top}}a^{4})/\mathrm{d}\beta$ on a finite set of
values which may leave discrete integration errors. But it is not difficult to
use enough values to render this a minor concern.

Finally, we want to compare the numerical efficiency to the method of Ref.
Jahn  _et al._ (2020). Fortunately, the single lattice we investigated in this
work was also used in that reference, and we can directly compare the error on
$\chi_{\mathrm{top}}a^{4}(\beta_{\mathrm{h}})/\chi_{\mathrm{top}}a^{4}(\beta_{\mathrm{c}})$
found here with the error on the same quantity found there, along with the
number of HMC trajectories needed in each case. In that reference
$\chi_{\mathrm{top}}(T)$ was determined at $\beta=(6.90097,7.30916,7.76294)$,
corresponding to $T=(2.5,4.1,7.0)\,T_{\mathrm{c}}$, a little narrower than the
range considered here. The three determinations required a total of $9.2\times
10^{6}$ trajectories, about half the number which should have been needed
here. The average trajectory length used was also shorter in that reference
than what we used here. The the final errors on $\ln(\chi_{\mathrm{top}})$ in
that study, for this lattice, were $(0.09,0.09,0.08)$ at the three
temperatures. In comparison, in comparing $2.8\,T_{\mathrm{c}}$ to
$7.0\,T_{\mathrm{c}}$ we find statistical errors of $0.30$. To reduce these
errors to the level of the other study would therefore require about 10 times
more statistics in our measurement runs, indicating that the present method is
of order 10 times less numerically efficient. Moreover, Ref. Jahn  _et al._
(2020) finds that the number of trajectories needed for a given statistical
error barely changes as one increases the volume (larger aspect ratio) or
makes the lattice finer (larger $N_{\tau}$ at fixed aspect ratio), whereas we
know that the number of updates needed for the method described in this paper
should scale with the number of lattice sites, see Eq. (18).

We conclude that our method is at least 10 times less efficient than the
single-temperature reweighting approach of Jahn  _et al._ (2020), and will
become still less efficient as one goes closer to the large-volume and
continuum limits. As we understand it, this also implies that it should be
easier in principle to achieve small statistical errors with the approach of
Jahn  _et al._ (2020) than with the approach of Borsanyi  _et al._ (2016b);
Frison  _et al._ (2016).

###### Acknowledgements.

The authors acknowledge support by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) through the CRC-TR 211 “Strong-interaction matter
under extreme conditions” – project number 315477589 – TRR 211. We also thank
the GSI Helmholtzzentrum and the TU Darmstadt and its Institut für Kernphysik
for supporting this research. Calculations were conducted on the Lichtenberg
high performance computer of the TU Darmstadt. This work was performed using
the framework of the publicly available openQCD-1.6 package ope .

##  Appendix A Comparison with the slope method

Our approach is closely related to the approach of Refs. Borsanyi  _et al._
(2016b); Frison  _et al._ (2016), and as we understand it, the errors per
numerical effort are nearly the same. To explain this conclusion, we start
with a quick review of their approach, and then look at the issue of
statistical power in each approach.

Their approach also seeks to compute Eq. (6) and then use a determination of
$\chi_{\mathrm{top}}(\beta_{\mathrm{c}})$ to determine $\chi_{\mathrm{top}}$
at other temperatures. Taking the log of Eq. (6) we find

$\displaystyle\ln\frac{\chi_{\mathrm{top}}(\beta_{\mathrm{h}})a^{4}(\beta_{\mathrm{h}})}{\chi_{\mathrm{top}}(\beta_{\mathrm{c}})a^{4}(\beta_{\mathrm{c}})}$
$\displaystyle=\Big{(}\ln(Z_{1}(\beta_{\mathrm{h}}))-\ln(Z_{1}(\beta_{\mathrm{c}})\Big{)}$
$\displaystyle\phantom{=}{}-\Big{(}\ln(Z_{0}(\beta_{\mathrm{h}}))-\ln(Z_{0}(\beta_{\mathrm{c}})\Big{)}\,.$
(20)

Now note that

$\displaystyle Z(\beta)$ $\displaystyle=\int\mathcal{D}A_{\mu}e^{-\beta
S}\quad\Rightarrow$ $\displaystyle-\frac{\partial\ln Z}{\partial\ln\beta}$
$\displaystyle=\frac{1}{Z}\int\mathcal{D}A_{\mu}e^{-\beta S}\beta
S=\langle\beta S\rangle$ (21)

the $\beta$ dependence of $\ln Z$ is set by the expectation value of the
action. Therefore

$\ln Z_{1}(\beta_{\mathrm{h}})-\ln
Z_{1}(\beta_{\mathrm{c}})=-\int_{\ln\beta_{\mathrm{c}}}^{\ln\beta_{\mathrm{h}}}\langle\beta
S\rangle_{1}\>\mathrm{d}(\ln\beta)$ (22)

and the ratio we want is

$\ln\frac{\chi_{\mathrm{top}}(\beta_{\mathrm{h}})a^{4}(\beta_{\mathrm{h}})}{\chi_{\mathrm{top}}(\beta_{\mathrm{c}})a^{4}(\beta_{\mathrm{c}})}=-\int_{\ln\beta_{\mathrm{c}}}^{\ln\beta_{\mathrm{h}}}\mathopen{}\mathclose{{}\left(\langle\beta
S\rangle_{1}-\langle\beta S\rangle_{0}}\right)\mathrm{d}(\ln\beta)\,.$ (23)

Refs. Borsanyi  _et al._ (2016b); Frison  _et al._ (2016) evaluate
$\langle\beta S\rangle$ in both $Q=0$ and $Q=1$ ensembles at a number of
$\beta$ values, which are then used to estimate this integral by, eg, the
trapezoid rule.

To compute Eq. (23) using our approach, first write

$\displaystyle\hskip
6.00006pt\int_{\ln\beta_{\mathrm{c}}}^{\ln\beta_{\mathrm{h}}}-\langle\beta
S\rangle\;\mathrm{d}(\ln\beta)$
$\displaystyle=\int_{\ln\beta_{\mathrm{c}}}^{\ln\beta_{\mathrm{h}}}\frac{-\sum_{i}e^{W[S_{i}]-\beta
S_{i}}\beta S_{i}}{\sum_{i}e^{W[S_{i}]-\beta S_{i}}}\mathrm{d}(\ln\beta)$
$\displaystyle=\int_{\ln\beta_{\mathrm{c}}}^{\ln\beta_{\mathrm{h}}}\frac{\mathrm{d}}{\mathrm{d}\ln\beta}\ln\mathopen{}\mathclose{{}\left(\sum_{i}e^{W[S_{i}]-\beta
S_{i}}}\right)\mathrm{d}(\ln\beta)$
$\displaystyle=\ln\frac{\sum_{i}e^{W[S_{i}]-\beta_{\mathrm{h}}S_{i}}}{\sum_{i}e^{W[S_{i}]-\beta_{\mathrm{c}}S_{i}}}$
(24)

so applying Eq. (11) to Eq. (23) leads directly to Eq. (12). Therefore the
real difference between the approaches is whether Eq. (23) is estimated based
on interpolating results for several temperatures, or using a single Markov
chain which spans all temperatures.

Now consider the statistical power of each approach. The accuracy of a Monte-
Carlo evaluation of $\langle\beta S\rangle$ is set by the variance of $\beta
S$ and the number of independent configurations used. The variance should be
reasonably approximated as that for $N_{\mathrm{dof}}$ Gaußian random
variables: $\sigma_{\beta S}^{2}\simeq N_{\mathrm{dof}}/2$. Therefore order-1
errors in $\langle\beta S\rangle$ require $N_{\mathrm{dof}}/2$ evaluations.
Since the expectation value determines an integrand, this is multiplied by the
integration range, so $N_{\mathrm{dof}}/2$ evaluations return an error in the
ratio of partition functions which is of order
$\ln(\beta_{\mathrm{max}})-\ln(\beta_{\mathrm{min}})$. Evaluating
$\langle\beta S\rangle$ at multiple $\beta$ values leads to a larger error at
each evaluation, but because each is responsible for a narrower $\Delta\beta$
range and the errors are uncorrelated, the final statistical uncertainty is
independent of the number of $\beta$ values used in the evaluation and depends
only on the total number of Markov steps and the width of the $\beta$ range
considered. The final error estimate is
$\Delta\ln\chi=\ln(\beta_{\mathrm{max}}/\beta_{\mathrm{min}})\,\sqrt{N_{\mathrm{dof}}/2N_{\mathrm{updates}}}$.
The error rises by $\sqrt{2}$ and $N_{\mathrm{updates}}$ is doubled when we
recall that separate simulations are needed in the $Q=0$ and $Q=1$ sectors.

In comparison, we see in Eq. (18) that our approach can explore the full
$\beta$ range, leading to order-1 errors in $\ln\chi$, in
$N_{\mathrm{updates}}\sim
N_{\mathrm{dof}}\,\ln^{2}(\beta_{\mathrm{max}}/\beta_{\mathrm{min}})$.
Therefore the two approaches produce errors per unit numerical effort which
are the same up to an order-1 factor. In a numerical experiment on a toy
problem ($N$ independent Gaußian random variables $x$ with action $S=\sum
x^{2}/2$) we find that the order-1 factor is in fact 1, so the two approaches
have the same statistical power per compute time, provided that $W[S]$ is well
determined and neglecting the computational effort expended in evaluating it.

We should also remark on how each approach is extended to full (unquenched)
QCD. In each case the main challenge is dealing with the way quark masses must
be varied with the lattice spacing and therefore with $\beta$: $m=m(\beta)$
(which must also be determined as part of the scale setting). This added
$\beta$ dependence changes Eq. (A), replacing $\langle\beta
S\rangle\to\langle\beta
S+\beta\,\mathrm{d}m/\mathrm{d}\beta\,\bar{\psi}\psi\rangle$. In our approach
one must replace $W[S]\to W[S]+\bar{\psi}(\not{D}+m(W^{\prime}[S]))\psi$ where
we use $W^{\prime}$ in place of $\beta$ for the scale dependence of $m$. This
amends Eq. (17) by the addition of the standard fermionic force term and by
the replacement
$\mathrm{d}W/\mathrm{d}S_{U}\to(\mathrm{d}W/\mathrm{d}S_{U})+(\mathrm{d}m/\mathrm{d}\beta)(\mathrm{d}W^{\prime}/\mathrm{d}S_{U})\langle\bar{\psi}\psi\rangle$
where $\langle\bar{\psi}\psi\rangle$ is the sum of the $\bar{\psi}\psi$ value
over all sites in the current configuration. Finally, in Eq. (12), the
$W-\beta S$ reweighting must be complemented by a determinant-ratio from the
$S$-dependent mass to the physical mass for the desired $\beta$ value.

## References

  * Peccei and Quinn (1977) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977).
  * Weinberg (1978) S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek (1978) F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
  * Abbott and Sikivie (1983) L. F. Abbott and P. Sikivie, Phys. Lett. B120, 133 (1983).
  * Preskill  _et al._ (1983) J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B 120, 127 (1983).
  * Irastorza and Redondo (2018) I. G. Irastorza and J. Redondo, Prog. Part. Nucl. Phys. 102, 89 (2018).
  * Grilli di Cortona  _et al._ (2016) G. Grilli di Cortona, E. Hardy, J. P. Vega, and G. Villadoro, JHEP 01, 034 (2016).
  * Gorghetto and Villadoro (2019) M. Gorghetto and G. Villadoro, JHEP 03, 033 (2019), arXiv:1812.01008 [hep-ph] .
  * Moore (2018) G. D. Moore,  _Proceedings, 35th International Symposium on Lattice Field Theory (Lattice 2017): Granada, Spain, June 18-24, 2017_ , EPJ Web Conf. 175, 01009 (2018).
  * Borsanyi  _et al._ (2016a) S. Borsanyi, M. Dierigl, Z. Fodor, S. D. Katz, S. W. Mages, D. Nogradi, J. Redondo, A. Ringwald, and K. K. Szabo, Phys. Lett. B 752, 175 (2016a).
  * Borsanyi  _et al._ (2016b) S. Borsanyi _et al._ , Nature 539, 69 (2016b).
  * Frison  _et al._ (2016) J. Frison, R. Kitano, H. Matsufuru, S. Mori, and N. Yamada, JHEP 09, 021 (2016).
  * Bonati  _et al._ (2018) C. Bonati, M. D’Elia, G. Martinelli, F. Negro, F. Sanfilippo, and A. Todaro, JHEP 11, 170 (2018), [JHEP18,170(2020)], arXiv:1807.07954 [hep-lat] .
  * Jahn  _et al._ (2018) P. T. Jahn, G. D. Moore, and D. Robaina, Phys. Rev. D98, 054512 (2018), arXiv:1806.01162 [hep-lat] .
  * Jahn  _et al._ (2020) P. T. Jahn, G. D. Moore, and D. Robaina, (2020), arXiv:2002.01153 [hep-lat] .
  * Laio  _et al._ (2016) A. Laio, G. Martinelli, and F. Sanfilippo, JHEP 07, 089 (2016).
  * Wang and Landau (2001) F. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001).
  * Berg and Neuhaus (1991) B. A. Berg and T. Neuhaus, Phys. Lett. B267, 249 (1991).
  * Duane  _et al._ (1987) S. Duane, A. D. Kennedy, B. J. Pendleton, and D. Roweth, Phys. Lett. B 195, 216 (1987).
  * Jahn  _et al._ (2019) P. T. Jahn, G. D. Moore, and D. Robaina, Eur. Phys. J. C79, 510 (2019), arXiv:1805.11511 [hep-lat] .
  * Laine and Rummukainen (1998) M. Laine and K. Rummukainen, Nucl. Phys. B535, 423 (1998).
  * Burnier  _et al._ (2017) Y. Burnier, H. T. Ding, O. Kaczmarek, A. L. Kruse, M. Laine, H. Ohno, and H. Sandmeyer, JHEP 11, 206 (2017), arXiv:1709.07612 [hep-lat] .
  * (23) http://luscher.web.cern.ch/luscher/openQCD/index.html.