# Combining Bayesian reconstruction entropy with maximum entropy method for
analytic continuations of matrix-valued Green’s functions

Songlin Yang  Liang Du<EMAIL_ADDRESS>Li Huang<EMAIL_ADDRESS>College of Physics and Technology, Guangxi Normal University, Guilin, Guangxi
541004, China

###### Abstract

The Bayesian reconstruction entropy is considered an alternative to the
Shannon-Jaynes entropy, as it does not exhibit the asymptotic flatness
characteristic of the Shannon-Jaynes entropy and obeys the scale invariance.
It is commonly utilized in conjunction with the maximum entropy method to
derive spectral functions from Euclidean time correlators produced by lattice
QCD simulations. This study expands the application of the Bayesian
reconstruction entropy to the reconstruction of spectral functions for
Matsubara or imaginary-time Green’s functions in quantum many-body physics.
Furthermore, it extends the Bayesian reconstruction entropy to implement the
positive-negative entropy algorithm, enabling the analytic continuations of
matrix-valued Green’s functions on an element-wise manner. Both the diagonal
and off-diagonal components of the matrix-valued Green’s functions are treated
equally. Benchmark results for the analytic continuations of synthetic Green’s
functions indicate that the Bayesian reconstruction entropy, when combined
with the preblur trick, demonstrates comparable performance to the Shannon-
Jaynes entropy. Notably, it exhibits greater resilience to noises in the input
data, particularly when the noise level is moderate.

##  I Introduction

In quantum many-body physics, single-particle Green’s function holds
significant importance Negele and Orland (1998); Coleman (2016). It is
typically computed using numerical methods such as quantum Monte Carlo Gull
_et al._ (2011); Foulkes  _et al._ (2001); Gubernatis  _et al._ (2016) and
functional renormalization group Metzner  _et al._ (2012); C. Platt and
Thomale (2013), or many-body perturbative techniques like random phase
approximation Harl and Kresse (2009); Chen  _et al._ (2017) and GW
approximation Hedin (1965); Aryasetiawan and Gunnarsson (1998); Onida  _et
al._ (2002); Gukelberger  _et al._ (2015), within the framework of imaginary
time or Matsubara frequency formulation. This quantity provides valuable
insights into the dynamic response of the system. However, it cannot be
directly compared to real-frequency experimental observables. Therefore, it is
necessary to convert the calculated values to real frequencies, which presents
the fundamental challenge of the analytic continuation problem. In theory, the
retarded Green’s function $G^{R}(\omega)$ and the spectral function
$A(\omega)$ can be derived from the imaginary-time Green’s function $G(\tau)$
or Matsubara Green’s function $G(i\omega_{n})$ by performing an inverse
Laplace transformation Negele and Orland (1998); Coleman (2016). Due to the
ill-conditioned nature of the transformation kernel, solving it directly is
nearly impossible. In fact, the output for this inverse problem is highly
sensitive to the input. Even minor fluctuations or noises in $G(\tau)$ or
$G(i\omega_{n})$ can result in nonsensical changes for $A(\omega)$ or
$G^{R}(\omega)$, making analytic continuations extremely unstable.

In recent decades, numerous techniques have been developed to address the
challenges associated with analytic continuation problems. These methods
include Padé approximation (PA) Beach  _et al._ (2000); Schött  _et al._
(2016a); Vidberg and Serene (1977); Markó  _et al._ (2017), maximum entropy
method (MaxEnt) Jarrell and Gubernatis (1996); Gubernatis  _et al._ (1991);
Silver  _et al._ (1990), stochastic analytic continuation (SAC) Sandvik (1998,
2016); Shao and Sandvik (2023); Beach (2004); Fuchs  _et al._ (2010); Vafayi
and Gunnarsson (2007); Reichman and Rabani (2009); Ghanem and Koch (2020a, b);
Syljuåsen (2008), stochastic optimization method (SOM) Mishchenko  _et al._
(2000); Bao  _et al._ (2016); Goulko  _et al._ (2017); Krivenko and Harland
(2019); Krivenko and Mishchenko (2022), stochastic pole expansion (SPX) Huang
and Liang (2023a, b), sparse modeling (SpM) Otsuki  _et al._ (2017); Motoyama
_et al._ (2022), Nevanlinna analytical continuation (NAC) Fei  _et al._
(2021a, b), causal projections Huang  _et al._ (2023), Prony fits Zhang and
Gull (2023); Ying (2022a, b), and machine learning assistant approaches
Fournier  _et al._ (2020); Yoon  _et al._ (2018); Huang and Yang (2022); Zhang
_et al._ (2022); Yao  _et al._ (2022); Arsenault  _et al._ (2017); Nichols
_et al._ (2022), and so on. Each method has its own advantages and
disadvantages. There is no doubt that the MaxEnt method is particularly
popular due to its ability to maintain a balance between computational
efficiency and accuracy Jarrell and Gubernatis (1996); Gubernatis  _et al._
(1991); Silver  _et al._ (1990). Several open source toolkits, such as
$\Omega$MAXENT Bergeron and Tremblay (2016), QMEM Sim and Han (2018), MAXENT
by Kraberger _et al._ Kraberger  _et al._ (2017), MAXENT by Levy _et al._ Levy
_et al._ (2017), ALF Assaad  _et al._ (2022), ACFlow Huang (2023), eDMFT Haule
_et al._ (2010), and ana_cont Kaufmann and Held (2023), have been developed to
implement this method.

The MaxEnt method is rooted in Bayesian inference Papoulis (1989). Initially,
the spectral function $A(\omega)$ is constrained to be non-negative and is
interpreted as a probability distribution. Subsequently, the MaxEnt method
endeavors to identify the most probable spectrum by minimizing a functional
$Q[A]$, which comprises the misfit functional $L[A]$ and the entropic term
$S[A]$ [see Eq. (10)]. Here, $L[A]$ quantifies the difference between the
input and reconstructed Green’s function, while $S[A]$ is employed to regulate
the spectrum Jarrell and Gubernatis (1996). The MaxEnt method is well-
established for the diagonal components of matrix-valued Green’s functions, as
the corresponding spectral functions are non-negative (for fermionic systems).
However, the non-negativity of the spectral functions cannot be assured for
the off-diagonal components of matrix-valued Green’s functions. As a result,
conventional MaxEnt method fails in this situation.

Currently, there is a growing emphasis on obtaining the full matrix
representation of the Green’s function (or the corresponding self-energy
function) on the real axis De Raychaudhury  _et al._ (2007); Sim and Han
(2018). This is necessary for calculating lattice quantities of interest in
Dyson’s equation Haule  _et al._ (2010). Therefore, there is a strong demand
for a reliable method to perform analytic continuation of the entire Green’s
function matrix Huang  _et al._ (2023); Zhang and Gull (2023); Fei  _et al._
(2021b). To meet this requirement, several novel approaches have been
suggested in recent years as enhancements to the conventional MaxEnt method.
The three most significant advancements are as follows: (1) _Auxiliary Green’s
function algorithm_. It is possible to create an auxiliary Green’s function by
combining the off-diagonal and diagonal elements of the matrix-valued Green’s
function to ensure the positivity of the auxiliary spectrum Tomczak and
Biermann (2007); Gull and Millis (2014); Reymbaut  _et al._ (2015, 2017); Gull
_et al._ (2013); Yue and Werner (2023). The analytic continuation of the
auxiliary Green’s function can then be carried out using the standard MaxEnt
method. (2) _Positive-negative entropy algorithm_. Kraberger _et al._ proposed
that the off-diagonal spectral functions can be seen as a subtraction of two
artificial positive functions. They extended the entropic term to relax the
non-negativity constraint presented in the standard MaxEnt method Kraberger
_et al._ (2017); Kaufmann and Held (2023). This enables simultaneous treatment
of both the diagonal and off-diagonal elements, thereby restoring crucial
constraints on the mathematical properties of the resulting spectral
functions, including positive semi-definiteness and Hermiticity. (3) _Maximum
quantum entropy algorithm_. Recently, Sim and Han generalized the MaxEnt
method by reformulating it with quantum relative entropy, maintaining the
Bayesian probabilistic interpretation Papoulis (1989); Jarrell and Gubernatis
(1996). The matrix-valued Green’s function is directly continued as a single
object without any further approximation or ad-hoc treatment Sim and Han
(2018).

Although various algorithms have been proposed to enhance the usefulness of
the MaxEnt method Kraberger  _et al._ (2017); Bergeron and Tremblay (2016);
Sim and Han (2018), further improvements in algorithms and implementations are
always beneficial. As mentioned above, the essence of the MaxEnt method lies
in the definition of the entropic term $S[A]$. Typically, it takes the form of
Shannon-Jaynes entropy (dubbed $S_{\text{SJ}}$) in the realm of quantum many-
body physics Jaynes (1957). However, alternative forms such as the Tikhonov
regulator Tikhonov  _et al._ (1995), the positive-negative entropy Kraberger
_et al._ (2017), and the quantum relative entropy Sim and Han (2018) are also
acceptable. It is worth noting that in the context of lattice QCD simulations,
the extraction of spectral functions from Euclidean time correlators is of
particular important. The MaxEnt method is also the primary tool for the
analytic continuation of lattice QCD data Asakawa  _et al._ (2001). Burnier
and Rothkopf introduced an enhanced dimensionless prior distribution, known as
the Bayesian reconstruction entropy (dubbed $S_{\text{BR}}$) Burnier and
Rothkopf (2013); Rothkopf (2022), which demonstrates superior asymptotic
behavior compared to $S_{\text{SJ}}$ Jaynes (1957). They discovered that by
using $S_{\text{BR}}$ in conjunction with the MaxEnt method, a significant
improvement in the reconstruction of spectral functions for Euclidean time
correlators can be achieved. They latter extended $S_{\text{BR}}$ to support
Bayesian inference of non-positive spectral functions in quantum field theory
Rothkopf (2017). But to the best of our knowledge, the application of
$S_{\text{BR}}$ has been limited to post-processing for lattice QCD data thus
far. Hence, some questions naturally arise. How effective is $S_{\text{BR}}$
for solving analytic continuation problems in quantum many-body simulations?
Is it truly superior to $S_{\text{SJ}}$? Keeping these questions in mind, the
primary objective of this study is to broaden the application scope of
$S_{\text{BR}}$. We at first verify whether $S_{\text{BR}}$ can handle
Matsubara (or imaginary-time) Green’s functions. Then, we generalize
$S_{\text{BR}}$ to implement the positive-negative entropy algorithm Kraberger
_et al._ (2017) for analytic continuations of matrix-valued Green’s functions.
The simulated results indicate that $S_{\text{BR}}$ works quite well in
transforming imaginary-time or Matsubara Green’s function data to real
frequency, no matter whether the Green’s function is a matrix or not. We find
that $S_{\text{BR}}$ has a tendency to sharpen the peaks in spectral
functions, but this shortcoming can be largely overcame by the preblur trick
Kaufmann and Held (2023); Kraberger  _et al._ (2017). Overall, the performance
of $S_{\text{BR}}$ is comparable to that of $S_{\text{SJ}}$ Jaynes (1957) for
the examples involved.

The rest of this paper is organized as follows. Section II introduces the
spectral representation of the single-particle Green’s function. Furthermore,
it explains the fundamental formalisms of the MaxEnt method, the Bayesian
reconstruction entropy, and the principle of the positive-negative entropy
algorithm. Section III is dedicated to various benchmark examples, including
analytic continuation results for synthetic single-band Green’s functions and
multi-orbital matrix-valued Green’s functions. The spectral functions obtained
by $S_{\text{SJ}}$ and $S_{\text{BR}}$ are compared with the exact spectra. In
Section IV, we further examine the preblur trick and the auxiliary Green’s
function algorithm for $S_{\text{SJ}}$ and $S_{\text{BR}}$. The robustness of
$S_{\text{BR}}$ with respect to noisy Matsubara data is discussed and compared
with that of $S_{\text{SJ}}$. A concise summary is presented in Section V.
Finally, Appendix A introduces the goodness-of-fit functional. The detailed
mathematical derivations for $S_{\text{SJ}}$ and $S_{\text{BR}}$ are provided
in Appendices B and C, respectively.

##  II Method

###  II.1 Spectral representation of single-particle Green’s function

It is well known that the single-particle imaginary-time Green’s function
$G(\tau)$ can be defined as follows:

$G(\tau)=\langle\mathcal{T}_{\tau}d(\tau)d^{\dagger}(0)\rangle,$ (1)

where $\tau$ is imaginary time, $\mathcal{T}_{\tau}$ means the time-ordered
operator, $d$ ($d^{\dagger}$) denotes the annihilation (creation) operator
Negele and Orland (1998); Coleman (2016). Given $G(\tau)$, Matsubara Green’s
function $G(i\omega_{n})$ can be constructed via direct Fourier
transformation:

$G(i\omega_{n})=\int^{\beta}_{0}d\tau~{}e^{-i\omega_{n}\tau}G(\tau),$ (2)

where $\beta$ is the inverse temperature of the system ($\beta\equiv 1/T$),
and $\omega_{n}$ is the Matsubara frequency. Note that $\omega_{n}$ is equal
to $(2n+1)\pi/\beta$ for fermions and $2n\pi/\beta$ for bosons ($n$ is an
integer).

Let us assume that the spectral function of the single-particle Green’s
function is $A(\omega)$, then the spectral representations of $G(\tau)$ and
$G(i\omega_{n})$ read:

$G(\tau)=\int^{+\infty}_{-\infty}K(\tau,\omega)A(\omega)d\omega,$ (3)

and

$G(i\omega_{n})=\int^{+\infty}_{-\infty}K(i\omega_{n},\omega)A(\omega)d\omega,$
(4)

respectively. Here, $K(\tau,\omega)$ and $K(i\omega_{n},\omega)$ are called
the kernel functions. Their definitions are as follows:

$K(\tau,\omega)=\frac{e^{-\tau\omega}}{1\pm e^{-\beta\omega}},$ (5)

and

$K(i\omega_{n},\omega)=\frac{1}{i\omega_{n}-\omega}.$ (6)

In the right-hand side of Eq. (5), + is for fermionic correlators and - is for
bosonic correlators Jarrell and Gubernatis (1996). In the subsequent
discussion, our focus will be on fermionic correlators.

We observe that Eqs. (3) and (4) are classified as Fredholm integral equations
of the first kind. When $A(\omega)$ is given, it is relatively straightforward
to compute the corresponding $G(\tau)$ and $G(i\omega_{n})$ by numerical
integration methods. However, the inverse problem of deducing $A(\omega)$ from
$G(\tau)$ or $G(i\omega_{n})$ is challenging due to the ill-conditioned nature
of the kernel function $K$. To be more specific, the condition number of $K$
is very large because of the exponential decay of $K$ with $\omega$ and
$\tau$. Consequently, direct inversion of $K$ becomes impractical from a
numerical standpoint Gubernatis  _et al._ (1991). Furthermore, the $G(\tau)$
or $G(i\omega_{n})$ data obtained from finite temperature quantum Monte Carlo
simulations is not free from errors Gull  _et al._ (2011); Foulkes  _et al._
(2001); Gubernatis  _et al._ (2016), further complicating the solution of Eqs.
(3) and (4).

###  II.2 Maximum entropy method

The cornerstone of the MaxEnt method is Bayes’ theorem. Let us treat the input
Green’s function $G$ and the spectrum $A$ as two events. According to Bayes’
theorem, the posterior probability $P[A|G]$ can be calculated by:

$P[A|G]=\frac{P[G|A]P[A]}{P[G]},$ (7)

where $P[G|A]$ is the likelihood function, $P[A]$ is the prior probability,
and $P[G]$ is the evidence Papoulis (1989). $P[G|A]$ is assumed to be in
direct proportion to $e^{-L[A]}$, where the misfit functional $L[A]$ reads:

$L[A]=\frac{1}{2}\chi^{2}[A]=\frac{1}{2}\sum^{N}_{i=1}\left(\frac{G_{i}-\tilde{G}_{i}[A]}{\sigma_{i}}\right)^{2}.$
(8)

Here, $N$ is the number of input data points, $\sigma$ denotes the standard
deviation of $G$, $\tilde{G}[A]$ is the reconstructed Green’s function via
Eqs. (3) and (4), and $\chi^{2}[A]$ is called the goodness-of-fit functional
(see Appendix A for more details). On the other hand, $P[A]$ is supposed to be
in direct proportion to $e^{\alpha S[A]}$, where $\alpha$ is a regulation
parameter and $S$ is entropy. Since the evidence $P[G]$ can be viewed as a
normalization constant, it is ignored in what follows. Thus,

$P[A|G]\propto e^{Q},$ (9)

where $Q$ is a functional of $A$ Jarrell and Gubernatis (1996). It is defined
as follows:

$Q[A]=\alpha S[A]-L[A]=\alpha S[A]-\frac{1}{2}\chi^{2}[A].$ (10)

The basic idea of the MaxEnt method is to identify the optimal spectrum
$\hat{A}$ that maximizes the posterior probability $P[A|G]$ (or equivalently
$Q[A]$). In essence, our goal is to determine the most favorable $\hat{A}$
that satisfies the following equation:

$\frac{\partial Q}{\partial A}\Big{|}_{A=\hat{A}}=0.$ (11)

Eq. (11) can be easily solved by the standard Newton’s algorithm Nocedal and
Wright (2006). In Appendices A, B and C, all terms in the right-hand side of
Eq. (10) are elaborated in detail. We also explain how to solve Eq. (11)
efficiently.

###  II.3 Bayesian reconstruction entropy

The entropy $S$ is also known as the Kullback-Leibler distance sometimes. In
principle, there are multiple choices for it. Perhaps $S_{\text{SJ}}$ is the
most frequently used in quantum many-body physics and condensed matter physics
Jaynes (1957); Jarrell and Gubernatis (1996). It reads:

$S_{\text{SJ}}[A]=\int
d\omega\left[A(\omega)-D(\omega)-A(\omega)\log\left(\frac{A(\omega)}{D(\omega)}\right)\right].$
(12)

Here, $D(\omega)$ is called the default model, providing the essential
features of spectra. Usually $D(\omega)$ is a constant or a Gaussian function,
but it can be determined by making use of high-frequency behavior of input
data as well Bergeron and Tremblay (2016). If both $A(\omega)$ and $D(\omega)$
have the same normalization, Eq. (12) reduces to:

$S_{\text{SJ}}[A]=-\int
d\omega~{}A(\omega)\log\left(\frac{A(\omega)}{D(\omega)}\right).$ (13)

The $S_{\text{BR}}$ introduced by Burnier and Rothkopf Burnier and Rothkopf
(2013); Rothkopf (2022) is dominant in high-energy physics and particle
physics. It reads:

$S_{\text{BR}}[A]=\int
d\omega\left[1-\frac{A(\omega)}{D(\omega)}+\log\left(\frac{A(\omega)}{D(\omega)}\right)\right].$
(14)

Note that $S_{\text{SJ}}$ is constructed from four axioms Jarrell and
Gubernatis (1996). While for $S_{\text{BR}}$, two of these axioms are replaced
Burnier and Rothkopf (2013). First of all, the scale invariance is
incorporated in $S_{\text{BR}}$. It means that $S_{\text{BR}}$ only depends on
the ratio between $A$ and $D$. The integrand in Eq. (14) is dimensionless,
such that the choice of units for $A$ and $D$ will not change the result of
the spectral reconstruction. Second, $S_{\text{BR}}$ favors choosing smooth
spectrum. The spectra that deviate between two adjacent frequencies,
$\omega_{1}$ and $\omega_{2}$, should be penalized.

###  II.4 Positive-negative entropy

The formulas discussed above are only correct for positive definite spectra
with a finite norm, i.e.,

$\int d\omega~{}A(\omega)>0.$ (15)

However, the norm could be zero for off-diagonal elements of matrix-valued
Green’s functions due to the fermionic anti-commutation relation Negele and
Orland (1998); Huang  _et al._ (2023):

$\int d\omega~{}A(\omega)=0.$ (16)

This suggests that the spectral function is not positive definite any more.
The equations for $S_{\text{SJ}}$ and $S_{\text{BR}}$ [i.e., Eqs. (12) and
(14)] should be adapted to this new circumstance. The positive-negative
entropy algorithm proposed by Kraberger _et al._ is a graceful solution to
this problem Kraberger  _et al._ (2017); Kaufmann and Held (2023). They
rewrite the off-diagonal spectral function $A(\omega)$ as the subtraction of
two positive definite spectra:

$A(\omega)=A^{+}(\omega)-A^{-}(\omega).$ (17)

Here $A^{+}(\omega)$ and $A^{-}(\omega)$ are independent, but have the same
norm. Then the resulting entropy can be split into two parts:

$S^{\pm}[A]=S[A^{+},A^{-}]=S[A^{+}]+S[A^{-}].$ (18)

$S^{\pm}[A]$ (or $S[A^{+},A^{-}]$) is called the positive-negative entropy,
which was first used for the analysis of NMR spectra Laue  _et al._ (1985);
Sibisi (1990). The expressions of positive-negative entropy for
$S_{\text{SJ}}$ and $S_{\text{BR}}$ are as follows:

$S^{\pm}_{\text{SJ}}[A]=\int
d\omega~{}\left[\sqrt{A^{2}+4D^{2}}-2D-A\log\left(\frac{\sqrt{A^{2}+4D^{2}}+A}{2D}\right)\right],$
(19)

and

$S^{\pm}_{\text{BR}}[A]=\int
d\omega~{}\left[2-\frac{\sqrt{A^{2}+D^{2}}+D}{D}+\log\left(\frac{\sqrt{A^{2}+D^{2}}+D}{2D}\right)\right].$
(20)

The default models for $A^{+}$ and $A^{-}$ are $D^{+}$ and $D^{-}$,
respectively. To derive Eqs. (19) and (20), $D^{+}=D^{-}=D$ is assumed. For a
detailed derivation, please see Appendices B and C. Similar to
$S_{\text{BR}}[A]$, $S^{\pm}_{\text{BR}}[A]$ also exhibits scale invariance.
In other words, it depends on $A/D$ only.

###  II.5 Entropy density

Figure 1: Comparison of entropy densities for different entropic terms. The
entropy $S$ is related to the entropy density $s$ via Eq. (21). For the
positive-negative entropy $S^{\pm}$, $D^{+}=D^{-}=D$ is assumed.

The integrand in the expression for entropy $S$ is referred as entropy density
$s$, i.e.,

$S[A]=\int d\omega~{}s[A(\omega)].$ (21)

Next, we would like to make a detailed comparison about the entropy densities
of different entropic terms. Supposed that $x=A/D$, the expressions for
various entropy densities are collected as follows:

$s_{\text{SJ}}(x)=x-1-x\log{x},$ (22) $s_{\text{BR}}(x)=1-x+\log{x},$ (23)
$s^{\pm}_{\text{SJ}}(x)=\sqrt{x^{2}+4}-2-x\log{\frac{\sqrt{x^{2}+4}+x}{2}},$
(24)
$s^{\pm}_{\text{BR}}(x)=2-(\sqrt{x^{2}+1}+1)+\log{\frac{\sqrt{x^{2}+1}+1}{2}}.$
(25)

They are visualized in Figure 1.

Clearly, all the entropy densities are strictly convex and non-positive (i.e.,
$s\leq 0$ and $s^{\pm}\leq 0$). The ordinary entropy density $s(x)$ is just
defined for positive $x$. $s_{\text{BR}}(x)$ becomes maximal at $x=1$, and
exhibits a similar quadratic behavior around its maximum as
$s_{\text{SJ}}(x)$. In the case $x\to 0$ ($A\ll D$), $s_{\text{BR}}(x)$ is not
suppressed, while $s_{\text{SJ}}(x)\to-1$. Thus, $s_{\text{BR}}(x)$ avoids the
asymptotic flatness inherent in $s_{\text{SJ}}(x)$ Burnier and Rothkopf
(2013). The positive-negative entropy density $s^{\pm}(x)$ is valid for any
$x$ ($x\in\mathbb{R}$). Both $s^{\pm}_{\text{BR}}(x)$ and
$s^{\pm}_{\text{SJ}}(x)$ are even functions. They also exhibit quadratic
behaviors around $x=0$, and $s^{\pm}_{\text{BR}}(x)\geq
s^{\pm}_{\text{SJ}}(x)$. In the limit of $\alpha\to\infty$, the goodness-of-
fit functional $\chi^{2}[A]$ has negligible weight, and the entropic term
$\alpha S[A]$ becomes dominant Bergeron and Tremblay (2016). Thus, the
maximization of $Q[A]$ (equivalently, the maximization of $S[A]$) yields
$\hat{A}=D$ (at $x=1$) for conventional entropy, in contrast to $\hat{A}=0$
(at $x=0$) for the positive-negative entropy Kraberger  _et al._ (2017).

##  III Results

Figure 2: Analytic continuations of single-band Green’s functions. (a)-(c) For
single off-centered peak case. (d)-(f) For two Gaussian peaks with a gap case.
In panels (a) and (d), the calculated and exact spectra are compared. In
panels (b) and (e), the $\log_{10}(\chi^{2})-\log_{10}(\alpha)$ curves are
plotted. The inverted triangle symbols indicate the optimal $\alpha$
parameters, which are determined by the $\chi^{2}$-kink algorithm Bergeron and
Tremblay (2016). The deviations of the reproduced Green’s function from the
raw ones are shown in panels (c) and (f). Note that in panel (c), $R$ is a
function of Matsubara frequency. However, $R$ is $\tau$-dependent in panel
(c).  Figure 3: Analytic continuations of two-band matrix-valued Green’s
functions. For analytic continuations for the diagonal elements of the Green’s
function, $S_{\text{SJ}}$ and $S_{\text{BR}}$ are used in the calculations.
While for the off-diagonal elements of the Green’s functions,
$S^{\pm}_{\text{SJ}}$ and $S^{\pm}_{\text{BR}}$ are employed. In the MaxEnt-BR
simulations, the preblur trick is always applied, and the blur parameter $b$
is fixed to 0.2. The default models for off-diagonal calculations are
constructed by $A_{11}$ and $A_{22}$. See main text for more details. (a)-(c)
$\theta=0.1$; (d)-(f) $\theta=0.5$; (g)-(i) $\theta=0.9$.

###  III.1 Computational setups

When the conventional MaxEnt method is combined with $S_{\text{SJ}}$ (or
$S^{\pm}_{\text{SJ}}$) Jaynes (1957), it is called the MaxEnt-SJ method. On
the other hand, a combination of the MaxEnt method with $S_{\text{BR}}$ (or
$S^{\pm}_{\text{BR}}$) Burnier and Rothkopf (2013); Rothkopf (2017, 2022) is
called the MaxEnt-BR method. We implemented both methods in the open source
analytic continuation toolkit, namely ACFlow Huang (2023). All the benchmark
calculations involved in this paper were done by this toolkit. The simulated
results for the MaxEnt-SJ and MaxEnt-BR methods will be presented in this
section. Next, we will explain the computational details.

We always start from an exact spectrum $A(\omega)$, which consists of one or
more Gaussian-like peaks. Then $A(\omega)$ is used to construct the imaginary-
time Green’s function $G(\tau)$ via Eq. (3), or the Matsubara Green’s function
$G(i\omega_{n})$ via Eq. (4). Here, for the sake of simplicity, the kernel
function $K$ is assumed to be fermionic. Since the realistic Green’s function
data from finite temperature quantum Monte Carlo simulations is usually noisy,
multiplicative Gaussian noise will be manually added to the clean $G(\tau)$ or
$G(i\omega_{n})$ to imitate this situation. The following formula is adopted:

$G_{\text{noisy}}=G_{\text{clean}}[1+\delta N_{\mathbb{C}}(0,1)],$ (26)

where $\delta$ denotes the noise level of the input data and $N_{\mathbb{C}}$
is the complex-valued normal Gaussian noise Huang  _et al._ (2023). Unless
stated otherwise, the standard deviation of $G$ is fixed to $10^{-4}$. The
numbers of input data for $G(\tau)$ and $G(i\omega_{n})$ are 1000 and 50,
respectively. During the MaxEnt simulations, the $\chi^{2}$-kink algorithm
Bergeron and Tremblay (2016) is used to determine the regulation parameter
$\alpha$. The Bryan algorithm Bryan (1990) is also tested. It always gives
similar results, which will not be presented in this paper. Once
$S_{\text{BR}}$ and $S^{\pm}_{\text{BR}}$ are chosen, the preblur trick
Kaufmann and Held (2023); Kraberger  _et al._ (2017) is used to smooth the
spectra. The blur parameter is adjusted case by case. Finally, the
reconstructed spectrum is compared with the true solution. We also adopt the
following quantity to quantify the deviation of the reconstructed Green’s
function from the input one:

$R_{i}=100\left|\frac{\Delta
G_{i}}{G_{i}}\right|=100\left|\frac{G_{i}-\tilde{G}_{i}[\hat{A}]}{G_{i}}\right|,$
(27)

where $\tilde{G}$ is evaluated by $\hat{A}$ via Eq. (3) or Eq. (4).

###  III.2 Single-band Green’s functions

At first, the analytic continuations of single-band Green’s functions are
examined. The exact spectral functions are constructed by a superposition of
some Gaussian peaks. We consider two representative spectra. (i) Single off-
centered peak. The spectral function is:

$A(\omega)=\exp\left[-\frac{(\omega-\epsilon)^{2}}{2\Gamma^{2}}\right],$ (28)

where $\epsilon=0.5$, $\Gamma=1.0$, $\delta=10^{-4}$, $\sigma=10^{-4}$, and
$\beta=20.0$. The input data is Matsubara Green’s function, which contains 50
Matsubara frequencies ($N=50$). The blur parameter for the MaxEnt-BR method is
0.45. (ii) Two Gaussian peaks with a gap. The spectral function is:

$\displaystyle A(\omega)$ $\displaystyle=$
$\displaystyle\frac{f_{1}}{\sqrt{2\pi}\Gamma_{1}}\exp\left[-\frac{(\omega-\epsilon_{1})^{2}}{2\Gamma_{1}^{2}}\right]$
(29) $\displaystyle+$
$\displaystyle\frac{f_{2}}{\sqrt{2\pi}\Gamma_{2}}\exp\left[-\frac{(\omega-\epsilon_{2})^{2}}{2\Gamma_{2}^{2}}\right],$

where $f_{1}=f_{2}=1.0$, $\epsilon_{1}=-\epsilon_{2}=2.0$,
$\Gamma_{1}=\Gamma_{2}=0.5$, $\delta=10^{-4}$, $\sigma=10^{-3}$, and
$\beta=5.0$. The synthetic data is imaginary-time Green’s function, which
contains 1000 imaginary time slices ($N=1000$). The blur parameter for the
MaxEnt-BR method is 0.30.

The analytic continuation results are shown in Fig. 2. For case (i), it is
clear that the exact spectrum and Matsubara Green’s function are well
reproduced by both the MaxEnt-SJ method and the MaxEnt-BR method [see Fig.
2(a) and (c)]. For case (ii), the MaxEnt-SJ method works quite well. The
MaxEnt-BR method can resolve the gap and the positions of the two peaks
accurately. However, it slightly overestimates their heights. We also observe
that the difference between $G(\tau)$ and $\tilde{G}(\tau)$ becomes relatively
apparent in the vicinity of $\tau=\beta/2$ [see Fig. 2(d) and (f)]. Figure
2(b) and (e) plot $\log_{10}(\chi^{2})$ as a function of $\log_{10}(\alpha)$.
The plots can be split into three distinct regions Bergeron and Tremblay
(2016): (i) _Default model region_ ($\alpha\to\infty$). $\chi^{2}$ is
relatively flat and goes to its global maximum. The entropic term $\alpha S$
is much larger than $\chi^{2}$. The obtained spectrum $A(\omega)$ resembles
the default model $D(\omega)$. (ii) _Noise-fitting region_ ($\alpha\to 0$).
$\chi^{2}$ approaches its global minimum, but it is larger than $\alpha S$. In
this region, the calculated spectrum $A(\omega)$ tends to fit the noise in
$G(\tau)$ or $G(i\omega_{n})$. (iii) _Information-fitting region_. $\chi^{2}$
increases with the increment of $\alpha$. It is comparable with $\alpha S$.
The optimal $\alpha$ parameter is located at the crossover between the noise-
fitting region and the information-fitting region. For $S_{\text{BR}}$ and
$S_{\text{SJ}}$, their $\chi^{2}(\alpha)$ curves almost overlap at the default
model region and the noise-fitting region, but differ at the information-
fitting region. It indicates that the optimal $\alpha$ parameters for
$S_{\text{SJ}}$ are larger than those for $S_{\text{BR}}$.

###  III.3 Matrix-valued Green’s functions

Next, the analytic continuations of matrix-valued Green’s functions are
examined. We consider a two-band model. The initial spectral function is a
diagonal matrix:

$\mathbf{A}^{\prime}(\omega)=\left[\begin{array}[]{cc}A^{\prime}_{11}(\omega)&0\\\
0&A^{\prime}_{22}(\omega)\\\ \end{array}\right].$ (30)

Here $A^{\prime}_{11}(\omega)$ and $A^{\prime}_{22}(\omega)$ are constructed
by using Eq. (29). For $A^{\prime}_{11}(\omega)$, the parameters are:
$f_{1}=f_{2}=0.5$, $\epsilon_{1}=1.0$, $\epsilon_{2}=2.0$, $\Gamma_{1}=0.20$,
and $\Gamma_{2}=0.70$. For $A^{\prime}_{22}(\omega)$, the parameters are:
$f_{1}=f_{2}=0.5$, $\epsilon_{1}=-1.0$, $\epsilon_{2}=-2.1$,
$\Gamma_{1}=0.25$, and $\Gamma_{2}=0.60$. The true spectral function is a
general matrix with non-zero off-diagonal components. It reads:

$\mathbf{A}(\omega)=\left[\begin{array}[]{cc}A_{11}(\omega)&A_{12}(\omega)\\\
A_{21}(\omega)&A_{22}(\omega)\\\ \end{array}\right].$ (31)

The diagonal matrix $\mathbf{A}^{\prime}(\omega)$ is rotated by a rotation
matrix $\mathbf{R}$ to generate $\mathbf{A}(\omega)$. The rotation matrix
$\mathbf{R}$ is defined as follows:

$\mathbf{R}=\left[\begin{array}[]{cc}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\\\ \end{array}\right],$ (32)

where $\theta$ denotes the rotation angle. In the present work, we consider
three different rotation angles. They are: (i) $\theta=0.1$, (ii)
$\theta=0.5$, and (iii) $\theta=0.9$. Then $\mathbf{A}(\omega)$ is used to
construct the matrix-valued Green’s function $\mathbf{G}(i\omega_{n})$ by
using Eq. (4). The essential parameters are: $\delta=0.0$, $\sigma=10^{-4}$,
and $\beta=40.0$. The input data contains 50 Matsubara frequencies. For
analytic continuations of the diagonal elements of $\mathbf{G}(i\omega_{n})$,
the default models are Gaussian-like. However, for the off-diagonal elements,
the default models are quite different. They are evaluated by
$D_{12}(\omega)=\sqrt{A_{11}(\omega)A_{22}(\omega)}$, where $A_{11}(\omega)$
and $A_{22}(\omega)$ are the calculated spectral functions for diagonal
elements Kraberger  _et al._ (2017); Kaufmann and Held (2023). This means that
we have to carry out analytic continuations for the diagonal elements at
first, and then use the diagonal spectra to prepare the default models for the
off-diagonal elements. For the MaxEnt-BR method, the preblur algorithm is used
to smooth the spectra. The blur parameter is fixed to 0.2.

The simulated results are illustrated in Fig. 3. For the diagonal spectral
spectral functions ($A_{11}$ and $A_{22}$), both the MaxEnt-SJ and MaxEnt-BR
methods can accurately resolve the peaks that are close to the Fermi level.
However, for the high-energy peaks, the MaxEnt-BR method tends to overestimate
their heights and yield sharper peaks [see the peaks around $\omega=2.0$ in
Fig. 3(a), (d), and (g)]. For the off-diagonal spectral functions, only
$A_{12}$ is shown, since $A_{21}$ is equivalent to $A_{12}$. We observe that
the major features of $A_{12}$ are well captured by both the MaxEnt-SJ and
MaxEnt-BR methods. Undoubtedly, the MaxEnt-SJ method works quite well in all
cases. The MaxEnt-BR method exhibits good performance when the rotation angle
is small [$\theta=0.1$, see Fig. 3(c)]. But some wiggles emerge around
$\omega=\pm 2.0$ when the rotation angle is large [$\theta=0.5$ and
$\theta=0.9$, see Fig. 3(f) and (i)]. We test more rotation angles
($0.0<\theta<2.0$). It seems that these wiggles won’t be enhanced unless the
blur parameter $b$ is decreased.

##  IV Discussions

Figure 4: Test of the preblur trick. A single-band Green’s function is
treated. The exact spectral function is generated by using Eq. (29). The model
parameters are the same as those used in Sec. III.2. Only the analytic
continuation results obtained by the MaxEnt-BR method are shown.  Figure 5:
Test of the preblur trick. A matrix-valued Green’s function is treated via the
positive-negative entropy approach. The initial spectral function is generated
by using Eq. (29). The model parameters are the same as those used in Sec.
III.3. The rotation angle $\theta$ is 0.1. Only the off-diagonal spectral
functions ($A_{12}$) obtained by the MaxEnt-BR method are shown.  Figure 6:
Analytic continuations of matrix-valued Green’s functions via the auxiliary
Green’s function approach. The model parameters are the same as those used in
Sec. III.3. For the MaxEnt-BR simulations, the preblur trick is always
applied, and the blur parameter $b$ is fixed to 0.2. Only the results for the
off-diagonal elements are presented ($A_{12}$). The upper panels show the
spectral functions for the auxiliary Green’s functions [see Eq. (35)]. The
lower panels show the off-diagonal spectral functions evaluated by Eq. (36).
(a)-(b) $\theta=0.1$; (c)-(d) $\theta=0.5$; (e)-(f) $\theta=0.9$.  Figure 7:
Robustness of the MaxEnt method with respect to the noisy Matsubara data. Here
we treat the off-diagonal elements ($G_{12}$) of the two-band matrix-valued
Green’s functions only via the positive-negative entropy method. The model
parameters are the same as those used in Sec. III.3. The rotation angle
$\theta$ is 0.5. The blur parameter for the MaxEnt-BR method is fixed to 0.2.
(a)-(g) Off-diagonal spectral functions ($A_{12}$) extracted from noisy
Matsubara data. $\delta$ denotes the noise level of the input data. It varies
from $10^{-8}$ (tiny noise) to $10^{-2}$ (large noise). (h) Integrated real
axis error as a function of $\delta$. It measures the deviations from the
exact spectral function, see Eq. (37).

In previous section, analytic continuations for single-band Green’s functions
and matrix-valued Green’s functions by using the MaxEnt-BR method have been
demonstrated. However, there are still some important issues that need to be
clarified. In this section, we would like to further discuss the preblur
algorithm, the auxiliary Green’s algorithm, and the noise tolerance for the
MaxEnt-BR method.

###  IV.1 Effect of preblur

The benchmark results shown in Section III imply that the MaxEnt-BR method has
a tendency to generate sharp peaks or small fluctuations in the high-energy
regions of the spectra. The preblur algorithm is helpful to alleviate this
phenomenon. In the context of analytic continuation, the preblur algorithm was
introduced by Kraberger _et al._ Kraberger  _et al._ (2017). The kernel
function $K(i\omega_{n},\omega)$ [see Eq. (6)] is “blurred” by using the
following expression:

$K_{b}(i\omega_{n},\omega)=\int^{\infty}_{-\infty}d\omega^{\prime}~{}K(i\omega_{n},\omega^{\prime})g_{b}(\omega-\omega^{\prime}).$
(33)

Here, $K_{b}(i\omega_{n},\omega)$ is the blurred kernel, which is then used in
Eq. (8) to evaluate the $\chi^{2}$-term. $g_{b}(\omega)$ is a Gaussian
function:

$g_{b}(\omega)=\frac{\exp{(-\omega^{2}/2b^{2})}}{\sqrt{2\pi}b},$ (34)

where $b$ is the blur parameter.

In Figures 4 and 5, we analyze the effects of the preblur trick for two
typical scenarios: (i) Positive spectral function with two separate Gaussian
peaks. (ii) Complicated spectral function for off-diagonal Green’s function.
Note that the model parameters for generating the exact spectra are taken from
Sections III.2 and III.3, respectively. It is evident that the MaxEnt-BR
method without the preblur trick ($b=0$) has trouble resolving the spectra
accurately. It usually favors sharp peaks (see Fig. 4). Sometimes it may lead
to undesirable artifacts, such as the side peaks around $\omega=\pm 1.5$ in
Fig. 5. The preblur algorithm can remedy this problem to some extent. The
major peaks are smoothed, and the artificial side peaks are suppressed upon
increasing $b$. But it’s not the case that bigger is always better. There is
an optimal $b$. For case (i), the optimal $b$ is 0.3. A larger $b$ ($b>0.3$)
will destroy the gap, inducing a metallic state. For case (ii), the optimal
$b$ is 0.2. If $b$ is further increased, it is difficult to obtain a stable
solution. Here, we should emphasize that these unphysical features seen in the
spectra are not unique for the MaxEnt-BR method. Actually, a similar tendency
was already observed in the early applications of the MaxEnt-SJ method in
image processing tasks. In order to address this problem, Skilling suggested
that the spectrum can be expressed as a Gaussian convolution: $A=g_{b}\star
h$, where $h$ is a “hidden” function Skilling (1991). Then the entropy is
evaluated from $h(\omega)$, instead of $A(\omega)$. In fact, Skilling’s
approach is equivalent to the preblur trick.

###  IV.2 Auxiliary Green’s functions

As stated before, in addition to the positive-negative entropy method, the
auxiliary Green’s algorithm can be used to continue the off-diagonal elements
of the Green’s functions Tomczak and Biermann (2007); Gull and Millis (2014);
Reymbaut  _et al._ (2015, 2017); Gull  _et al._ (2013); Yue and Werner (2023).
Its idea is quite simple. At first, an auxiliary Green’s function is
constructed as follows:

$G_{\text{aux}}(i\omega_{n})=G_{11}(i\omega_{n})+G_{22}(i\omega_{n})+2G_{12}(i\omega_{n}).$
(35)

Supposed that $A_{\text{aux}}(\omega)$ is the corresponding spectral function
for $G_{\text{aux}}(i\omega_{n})$. It is easy to prove that
$A_{\text{aux}}(\omega)$ is positive semi-definiteness. So, it is safe to
perform analytic continuation for $G_{\text{aux}}(i\omega_{n})$ by using the
traditional MaxEnt-SJ or MaxEnt-BR method. Next, the off-diagonal spectral
function $A_{12}(\omega)$ can be evaluated as follows:

$A_{12}(\omega)=\frac{A_{\text{aux}}(\omega)-A_{11}(\omega)-A_{22}(\omega)}{2}.$
(36)

This algorithm is not restricted to the particle-hole symmetric case, as
assumed in Ref. Gull and Millis (2014). Here, we employ the two-band’s example
presented in Section III.3 again to test the combination of the auxiliary
Green’s function algorithm and the MaxEnt-BR method (and the MaxEnt-SJ
method). The model and computational parameters are kept, and the analytic
continuation results are plotted in Fig. 6.

Three rotation angles are considered in the simulations. We find the auxiliary
Green’s function algorithm is numerically unstable. (i) $\theta=0.1$. When the
rotation angle is small, the amplitude (absolute value) of $A_{12}(\omega)$ is
small. Both the MaxEnt-SJ method and the MaxEnt-BR method can resolve the
major peaks around $\omega=\pm 1.0$. But they will produce apparent
oscillations at higher energies. Increasing the $b$ parameter further could
lead to wrong peaks. It seems that these unphysical features are likely due to
the superposition of errors in $A_{\text{aux}}(\omega)$, $A_{11}(\omega)$, and
$A_{22}(\omega)$. The performance of the MaxEnt-BR method is worse than that
of the MaxEnt-SJ method. (ii) $\theta=0.5$ and $\theta=0.9$. When the rotation
angle is moderate or large, the spectra obtained by the MaxEnt-SJ method are
well consistent with the exact solutions. By using the MaxEnt-BR method,
though the fluctuations in high-energy regions are greatly suppressed, small
deviations from the exact spectra still exist. Overall, the auxiliary Green’s
function algorithm is inferior to the positive-negative entropy method in the
examples studied. Especially when the rotation angle is small, the auxiliary
Green’s function algorithm usually fails, irrespective of which form of
entropy is adopted.

###  IV.3 Robustness with respect to noisy input data

Analytic continuation is commonly used on noisy Monte Carlo data. The
precision of the Monte Carlo data, which strongly depends on the sampling
algorithm and the estimator used, is rarely better than $10^{-5}$ Gull  _et
al._ (2011); Foulkes  _et al._ (2001); Gubernatis  _et al._ (2016). The
distributions of inherent errors in the Monte Carlo data are often Gaussian.
This poses a strict requirement for the noise tolerance of the analytic
continuation methods. Previous studies have demonstrated that the MaxEnt-SJ
method is robust to noise. Here, we would like to examine the robustness of
the MaxEnt-BR method with respect to noisy Matsubara data. We reconsider the
analytic continuation of the off-diagonal elements of matrix-valued Green’s
functions. As mentioned above, the synthetic Matsubara data for this case is
assumed to be noiseless ($\delta=0.0$). Now the noise is manually added by Eq.
(26) and the noise level $\delta$ is changed from $10^{-8}$ to $10^{-2}$. The
other computational parameters are the same as those used in Section III.3. In
order to estimate the sensitivity of the analytic continuation method to the
noisy data, a new quantity, namely integrated real axis error, is introduced:

$err(\delta)=\int d\omega~{}\left|A(\omega)-\hat{A}_{\delta}(\omega)\right|.$
(37)

Here, $\hat{A}_{\delta}(\omega)$ means the reconstructed (optimal) spectral
function under the given noise level $\delta$.

In Figure 7(a)-(g), the convergence of the spectral functions by the MaxEnt-BR
method for simulated Gaussian errors with varying magnitude is shown. When
$\delta=10^{-2}$, the main peaks at $\omega\approx\pm 1.0$ are roughly
reproduced, while the side peaks at $\omega\approx 2.0$ are completely smeared
out. When $\delta=10^{-3}$, the major characteristics of the off-diagonal
spectrum are successfully captured. As the simulated errors decrease further
($\delta<10^{-3}$), the MaxEnt-BR method rapidly converges to the exact
solution. For comparison, the results for the MaxEnt-SJ method are also
presented. It is evident that the MaxEnt-SJ method is less robust than the
MaxEnt-BR method when the noise level is moderate. It fails to resolve the
satellite peaks around $\omega=\pm 2.0$. Only when $\delta\leq 10^{-6}$, the
MaxEnt-SJ method can recover the exact spectrum. Figure 7(h) exhibits the
integrated error $err(\delta)$. When $10^{-6}<\delta<10^{-3}$, the MaxEnt-BR
method exhibits better robustness with respect to noise than the MaxEnt-SJ
method.

##  V Concluding remarks

In summary, we extend the application scope of $S_{\text{BR}}$ to analytic
continuations of imaginary-time Green’s functions and Matsubara Green’s
functions in quantum many-body physics and condense matter physics. It is
further generalized to the form of positive-negative entropy to support
analytic continuation of matrix-valued Green’s function, in which the positive
semi-definiteness of the spectral function is broken. We demonstrate that the
MaxEnt-BR method, in conjunction with the preblur algorithm and the positive-
negative entropy algorithm, is capable of capturing the primary features of
the diagonal and off-diagonal spectral functions, even in the presence of
moderate levels of noise. Overall, its performance is on par with that of the
MaxEnt-SJ method in the examples studied. Possible applications of the MaxEnt-
BR method in the future include analytic continuations of anomalous Green’s
functions and self-energy functions Yue and Werner (2023); Gull and Millis
(2014), bosonic response functions (such as optical conductivity and spin
susceptibility) Gunnarsson  _et al._ (2010); Schött  _et al._ (2016b), and
frequency-dependent transport coefficients with non-positive spectral weight
Reymbaut  _et al._ (2017, 2015), etc. Further investigations are highly
desirable.

Finally, the MaxEnt-BR method, together with the MaxEnt-SJ method, has been
integrated into the open source software package ACFlow Huang (2023), which
may be useful for the analytic continuation community.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
(No. 12274380 and No. 11934020), and the Central Government Guidance Funds for
Local Scientific and Technological Development (No. GUIKE ZY22096024).

##  Appendix A Goodness-of-fit functional

Given the spectral function $A(\omega)$, the Matsubara Green’s function
$\tilde{G}[A]$ could be reconstructed by using the following equation (the
$\tilde{~{}}$ symbol is used to distinguish the reconstructed Green’s function
from the input Green’s function):

$\tilde{G}_{n}[A]=\int d\omega~{}K(i\omega_{n},\omega)A(\omega),$ (38)

where $K(i\omega_{n},\omega)$ denotes the kernel function, and $n$ is the
index for Matsubara frequency. For the sake of simplicity, Eq. (38) is
reformulated into its discretization form:

$\tilde{G}_{n}[A]=\sum_{m}K_{nm}A_{m}\Delta_{m},$ (39)

where $K_{nm}\equiv K(i\omega_{n},\omega_{m})$, $A_{i}\equiv A(\omega_{i})$,
and $\Delta_{i}$ means the weight of the mesh at real axis.

The goodness-of-fit functional $\chi^{2}[A]$ measures the distance between the
input Green’s function $G$ and the reconstructed Green’s function
$\tilde{G}[A]$. Its expression is as follows:

$\chi^{2}[A]=\sum^{N}_{n=1}\frac{(G_{n}-\tilde{G}_{n}[A])^{2}}{\sigma^{2}_{n}}$
(40)

Here we just assume that there are $N$ data points, and $\sigma_{n}$ denotes
the standard derivative (i.e., the error bar) of $G_{n}$. Substituting Eq.
(39) into Eq. (40), then we arrive:

$\chi^{2}[A]=\sum^{N}_{n=1}\frac{(G_{n}-\sum_{m}K_{nm}A_{m}\Delta_{m})^{2}}{\sigma^{2}_{n}}.$
(41)

The first derivative of $\chi^{2}[A]$ with respect to $A_{i}$ reads:

$\frac{\partial\chi^{2}[A]}{\partial
A_{i}}=2\sum^{N}_{n=1}\frac{K_{ni}\Delta_{i}}{\sigma^{2}_{n}}\left(\sum_{m}K_{nm}A_{m}\Delta_{m}-G_{n}\right).$
(42)

##  Appendix B Shannon-Jaynes entropy

In this appendix, the technical details for $S_{\text{SJ}}$,
$S^{\pm}_{\text{SJ}}$, and the MaxEnt-SJ method are reviewed. Most of the
equations presented in this appendix have been derived in Refs. Jarrell and
Gubernatis (1996), Kraberger  _et al._ (2017), and Kaufmann and Held (2023).
We repeat the mathematical derivation here so that this paper is self-
contained.

###  B.1 Entropy

The Shannon-Jaynes entropy is defined as follows:

$S_{\text{SJ}}[A]=\int
d\omega\left[A(\omega)-D(\omega)-A(\omega)\log\left(\frac{A(\omega)}{D(\omega)}\right)\right],$
(43)

The discretization form of Eq. (43) is:

$S_{\text{SJ}}[A]=\sum_{i}\Delta_{i}\left[A_{i}-D_{i}-A_{i}\log\left(\frac{A_{i}}{D_{i}}\right)\right].$
(44)

The first derivative of $S_{\text{SJ}}[A]$ with respect to $A_{i}$ reads:

$\frac{\partial S_{\text{SJ}}[A]}{\partial
A_{i}}=-\Delta_{i}\log\left(\frac{A_{i}}{D_{i}}\right).$ (45)

###  B.2 Parameterization of spectral function

A singular value decomposition can be readily performed for the kernel
function $K$:

$K=U\xi V^{T},$ (46)

where $U$ and $V$ are column-orthogonal matrices and $\xi$ is the vector of
singular values. So, the matrix element of $K$ reads:

$K_{ni}=\sum_{m}U_{nm}\xi_{m}V_{im}.$ (47)

The columns of $V$ can be understood as basis functions for the spectral
function. That is to say they span a so-called singular value space. And in
the singular value space the spectrum can be parameterized through:

$A_{l}=D_{l}\exp\left(\sum_{m}V_{lm}u_{m}\right).$ (48)

It is easy to prove that:

$\frac{\partial A_{l}}{\partial
u_{i}}=D_{l}\exp\left(\sum_{m}V_{lm}u_{m}\right)V_{li}=A_{l}V_{li}.$ (49)

###  B.3 Newton’s method

Next, we will derive the equations for $\\{u_{m}\\}$. Let us recall the
stationary condition [see also Eq. (11)]:

$\frac{\partial Q[A]}{\partial A_{i}}=0.$ (50)

It is actually:

$\alpha\frac{\partial S_{\text{SJ}}[A]}{\partial
A_{i}}-\frac{1}{2}\frac{\partial\chi^{2}[A]}{\partial A_{i}}=0.$ (51)

Substituting Eqs. (45) and (42) into the above equation:

$\alpha\Delta_{i}\log\left(\frac{A_{i}}{D_{i}}\right)+\sum^{N}_{n=1}\frac{K_{ni}\Delta_{i}}{\sigma^{2}_{n}}\left(\sum_{m}K_{nm}A_{m}\Delta_{m}-G_{n}\right)=0.$
(52)

Eliminating $\Delta_{i}$:

$\alpha\log\left(\frac{A_{i}}{D_{i}}\right)+\sum^{N}_{n=1}\frac{K_{ni}}{\sigma^{2}_{n}}\left(\sum_{m}K_{nm}A_{m}\Delta_{m}-G_{n}\right)=0.$
(53)

Substituting Eqs. (48) and (47) into Eq. (53):

$\displaystyle\alpha\sum_{m}V_{im}u_{m}+$ (54)
$\displaystyle~{}~{}\sum^{N}_{n=1}\frac{1}{\sigma^{2}_{n}}\sum_{m}U_{nm}\xi_{m}V_{im}\left(\sum_{l}K_{nl}A_{l}\Delta_{l}-G_{n}\right)=0.$

Now $\sum_{m}V_{im}$ can be removed from the two terms in the left-hand side.
Finally we get:

$\alpha
u_{m}+\sum^{N}_{n=1}\frac{1}{\sigma^{2}_{n}}\xi_{m}U_{nm}\left(\sum_{l}K_{nl}A_{l}\Delta_{l}-G_{n}\right)=0.$
(55)

Since $A_{l}$ depends on $\\{u_{m}\\}$ as well, Eq. (55) is actually a non-
linear equation about $\\{u_{m}\\}$. In the ACFlow package Huang (2023), the
Newton’s method is adopted to solve it. So we have to evaluate the following
two variables and pass them to the Newton’s algorithm Nocedal and Wright
(2006):

$f_{m}=\alpha
u_{m}+\xi_{m}\sum^{N}_{n=1}\frac{1}{\sigma^{2}_{n}}U_{nm}\left(\sum_{l}K_{nl}A_{l}\Delta_{l}-G_{n}\right),$
(56) $J_{mi}=\frac{\partial f_{m}}{\partial
u_{i}}=\alpha\delta_{mi}+\xi_{m}\sum^{N}_{n=1}\frac{1}{\sigma^{2}_{n}}U_{nm}\sum_{l}K_{nl}A_{l}\Delta_{l}V_{li},$
(57)

where $J_{mi}$ can be considered as a Jacobian matrix. Note that Eq. (49) is
applied to derive Eq. (57).

Since the calculations for $f_{m}$ and $J_{mi}$ are quite complicated and
time-consuming, the Einstein summation technology is used to improve the
computational efficiency. The following three variables should be precomputed
and stored during the initialization stage:

$B_{m}=\sum^{N}_{n=1}\frac{1}{\sigma^{2}_{n}}\xi_{m}U_{nm}G_{n},$ (58)
$W_{ml}=\sum_{pn}\frac{1}{\sigma^{2}_{n}}U_{nm}\xi_{m}U_{np}\xi_{p}V_{lp}\Delta_{l}D_{l},$
(59) $W_{mli}=W_{ml}V_{li}.$ (60)

Clearly, the three variables only depend on the input Green’s function $G$,
the default model $D$, and the singular value decomposition of the kernel
function $K$. Now $f_{m}$ and $J_{mi}$ can be reformulated in terms of them:

$f_{m}=\alpha u_{m}+\sum_{l}W_{ml}w_{l}-B_{m},$ (61)
$J_{mi}=\alpha\delta_{mi}+\sum_{l}W_{mli}w_{l},$ (62)

where

$w_{l}=\exp\left(\sum_{m}V_{lm}u_{m}\right).$ (63)

###  B.4 Positive-negative entropy

For matrix-valued Green’s function, the spectral functions of the off-diagonal
components could exhibit negative weight. However, the ordinary MaxEnt method
is only rigorous for non-negative spectrum Jarrell and Gubernatis (1996). In
order to remedy this problem, one could imagine that the off-diagonal spectral
functions originate from a subtraction of two artificial positive functions
Kraberger  _et al._ (2017),

$A=A^{+}-A^{-},$ (64) $A^{-}=A^{+}-A.$ (65)

Assuming the independence of $A^{+}$ and $A^{-}$, the resulting entropy
$S_{\text{SJ}}[A^{+},A^{-}]$ is the sum of the respective entropies:

$\displaystyle S_{\text{SJ}}[A^{+},A^{-}]=$ $\displaystyle\int
d\omega\left[A^{+}-D-A^{+}\log\left(\frac{A^{+}}{D}\right)\right]+$ (66)
$\displaystyle\int
d\omega\left[A^{-}-D-A^{-}\log\left(\frac{A^{-}}{D}\right)\right].$

Thus, $S_{\text{SJ}}[A^{+},A^{-}]$ is called the positive-negative entropy in
the literature Kraberger  _et al._ (2017); Kaufmann and Held (2023). We at
first eliminate $A^{-}$

$S_{\text{SJ}}[A,A^{+}]=\int
d\omega\left[A^{+}-D-A^{+}\log\left(\frac{A^{+}}{D}\right)\right]+\int
d\omega\left[(A^{+}-A)-D-(A^{+}-A)\log\left(\frac{A^{+}-A}{D}\right)\right],$
(67)

and write:

$Q[A,A^{+}]=\alpha S_{\text{SJ}}[A,A^{+}]-\frac{1}{2}\chi^{2}[A].$ (68)

Since we are searching for a maximum of $Q$ with respect to $A$ and $A^{+}$,
we also apply:

$\frac{\partial Q[A,A^{+}]}{\partial A^{+}}=\frac{\partial
S_{\text{SJ}}[A,A^{+}]}{\partial A^{+}}=0,$ (69)

to eliminate $A^{+}$. The following equations show some intermediate steps:

$\int
d\omega\left[\log\left(\frac{A^{+}}{D}\right)+\log\left(\frac{A^{+}-A}{D}\right)\right]=0,$
(70) $\log\left(\frac{A^{+}}{D}\right)+\log\left(\frac{A^{+}-A}{D}\right)=0,$
(71) $\frac{A^{+}(A^{+}-A)}{D^{2}}=1.$ (72)

Finally, we obtain:

$A^{+}=\frac{\sqrt{A^{2}+4D^{2}}+A}{2},$ (73)
$A^{-}=\frac{\sqrt{A^{2}+4D^{2}}-A}{2}.$ (74)

Substituting Eqs. (73) and (74) into Eq. (66) to eliminate $A^{+}$ and
$A^{-}$, we obtain:

$\displaystyle S_{\text{SJ}}[A^{+},A^{-}]$ $\displaystyle=$ $\displaystyle\int
d\omega\Bigg{[}\sqrt{A^{2}+4D^{2}}-2D$ (75) $\displaystyle-$ $\displaystyle
A\log\left(\frac{\sqrt{A^{2}+4D^{2}}+A}{2D}\right)\Bigg{]}.$

Its discretization form becomes:

$\displaystyle S_{\text{SJ}}[A^{+},A^{-}]$ $\displaystyle=$
$\displaystyle\sum_{m}\Bigg{[}\sqrt{A_{m}^{2}+4D_{m}^{2}}-2D_{m}$ (76)
$\displaystyle-$ $\displaystyle
A_{m}\log\left(\frac{\sqrt{A_{m}^{2}+4D_{m}^{2}}+A_{m}}{2D_{m}}\right)\Bigg{]}\Delta_{m}.$

Note that Eq. (75) is exactly the same with Eq. (19).
$S_{\text{SJ}}[A^{+},A^{-}]$ is just $S^{\pm}_{\text{SJ}}[A]$. The first
derivative of $S_{\text{SJ}}[A^{+},A^{-}]$ with respect to $A_{i}$ is:

$\frac{\partial S_{\text{SJ}}[A^{+},A^{-}]}{\partial
A_{i}}=-\Delta_{i}\log\left(\frac{A^{+}_{i}}{D_{i}}\right).$ (77)

###  B.5 Parameterization for positive-negative spectral functions

According to Eq. (73) and Eq. (74), we find:

$A^{+}A^{-}=D^{2}.$ (78)

Inspired by Eq. (48), it is naturally to parameterize $A^{+}$ and $A^{-}$ in
the singular value space as well:

$A^{+}_{i}=D_{i}\exp\left(\sum_{m}V_{im}u_{m}\right),$ (79)
$A^{-}_{i}=D_{i}\exp\left(-\sum_{m}V_{im}u_{m}\right).$ (80)

Hence,

$A_{i}=D_{i}\exp\left(\sum_{m}V_{im}u_{m}\right)-D_{i}\exp\left(-\sum_{m}V_{im}u_{m}\right)$
(81) $A_{i}=D_{i}\left(w_{i}-\frac{1}{w_{i}}\right).$ (82) $\frac{\partial
A^{+}_{l}}{\partial u_{i}}=A^{+}_{l}V_{li}.$ (83) $\frac{\partial
A^{-}_{l}}{\partial u_{i}}=-A^{-}_{l}V_{li}.$ (84) $\frac{\partial A}{\partial
u_{i}}=(A^{+}_{l}+A^{-}_{l})V_{li}=D_{l}\left(w_{l}+\frac{1}{w_{l}}\right)V_{li}.$
(85)

By using the Einstein summation notations as defined in Section B.3, we
immediately obtain:

$f_{m}=\alpha u_{m}+\sum_{l}W_{ml}\left(w_{l}-\frac{1}{w_{l}}\right)-B_{m},$
(86)
$J_{mi}=\alpha\delta_{mi}+\sum_{l}W_{mli}\left(w_{l}+\frac{1}{w_{l}}\right).$
(87)

##  Appendix C Bayesian reconstruction entropy

In this appendix, the technical details for $S_{\text{BR}}$,
$S^{\pm}_{\text{BR}}$, and the MaxEnt-BR method are discussed.

###  C.1 Entropy

The Bayesian reconstruction entropy reads:

$S_{\text{BR}}[A]=\int
d\omega\left[1-\frac{A(\omega)}{D(\omega)}+\log\left(\frac{A(\omega)}{D(\omega)}\right)\right].$
(88)

Its discretization form is:

$S_{\text{BR}}[A]=\sum_{i}\Delta_{i}\left[1-\frac{A_{i}}{D_{i}}+\log\left(\frac{A_{i}}{D_{i}}\right)\right].$
(89)

Its first derivative with respect to $A_{i}$ is:

$\frac{\partial S_{\text{BR}}[A]}{\partial
A_{i}}=-\Delta_{i}\left(\frac{1}{D_{i}}-\frac{1}{A_{i}}\right).$ (90)

###  C.2 Parameterization of spectral function

The original parameterization of $A$ [see Eq. (48)] can not be used here. A
new parameterization scheme for $A$ is necessary. Assumed that

$\frac{1}{D_{i}}-\frac{1}{A_{i}}=\sum_{m}V_{im}u_{m},$ (91)

then we have:

$A_{l}=\frac{D_{l}}{1-D_{l}\sum_{m}V_{lm}u_{m}}.$ (92)

So, the first derivative of $A_{l}$ with respect to $u_{i}$ reads:

$\frac{\partial A_{l}}{\partial u_{i}}=A_{l}A_{l}V_{li}.$ (93)

###  C.3 Newton’s method

Next, we would like to derive $f_{m}$ and $J_{mi}$ for the Bayesian
reconstruction entropy. Substituting Eq. (90) into Eq. (51):

$\alpha\Delta_{i}\left(\frac{1}{D_{i}}-\frac{1}{A_{i}}\right)+\sum^{N}_{n=1}\frac{K_{ni}\Delta_{i}}{\sigma^{2}_{n}}\left(\sum_{m}K_{nm}A_{m}\Delta_{m}-G_{n}\right)=0.$
(94)

Eliminating $\Delta_{i}$:

$\alpha\left(\frac{1}{D_{i}}-\frac{1}{A_{i}}\right)+\sum^{N}_{n=1}\frac{K_{ni}}{\sigma^{2}_{n}}\left(\sum_{m}K_{nm}A_{m}\Delta_{m}-G_{n}\right)=0.$
(95)

Substituting Eq. (91) into Eq. (95):

$\displaystyle\alpha\sum_{m}V_{im}u_{m}+$ (96)
$\displaystyle~{}~{}\sum^{N}_{n=1}\frac{1}{\sigma^{2}_{n}}\sum_{m}U_{nm}\xi_{m}V_{im}\left(\sum_{l}K_{nl}A_{l}\Delta_{l}-G_{n}\right)=0.$

Eliminating $\sum_{m}V_{im}$ again:

$\alpha
u_{m}+\sum^{N}_{n=1}\frac{1}{\sigma^{2}_{n}}\xi_{m}U_{nm}\left(\sum_{l}K_{nl}A_{l}\Delta_{l}-G_{n}\right)=0.$
(97)

Clearly, Eq. (97) is the same with Eq. (55). Thus, the equations for $f_{m}$
and $J_{mi}$ are as follows:

$f_{m}=\alpha
u_{m}+\xi_{m}\sum^{N}_{n=1}\frac{1}{\sigma^{2}_{n}}U_{nm}\left(\sum_{l}K_{nl}A_{l}\Delta_{l}-G_{n}\right),$
(98) $J_{mi}=\frac{\partial f_{m}}{\partial
u_{i}}=\alpha\delta_{mi}+\xi_{m}\sum^{N}_{n=1}\frac{1}{\sigma^{2}_{n}}U_{nm}\sum_{l}K_{nl}A_{l}A_{l}\Delta_{l}V_{li}.$
(99)

Note that Eq. (98) is the same with Eq. (56). However, Eq. (99) differs from
Eq. (57) by a factor $A_{l}$ in the second term of the right-hand side. To
derive Eq. (99), Eq. (93) is used. Now the Einstein summation notation is
adopted to simplify the calculations of $f_{m}$ and $J_{mi}$ again:

$f_{m}=\alpha u_{m}+\sum_{l}W_{ml}w_{l}-B_{m},$ (100)
$J_{mi}=\alpha\delta_{mi}+\sum_{l}W_{mli}D_{l}w_{l}w_{l},$ (101)

where

$w_{l}=\frac{1}{1-D_{l}\sum_{m}V_{lm}u_{m}}.$ (102)

Here, the definitions of $W_{ml}$, $W_{mli}$, and $B_{m}$ could be found in
Section B.3.

###  C.4 Positive-negative entropy

Next, we would like to generalize the Bayesian reconstruction entropy to
realize the positive-negative entropy algorithm to support the analytic
continuation of matrix-valued Green’s function Kraberger  _et al._ (2017). The
positive-negative entropy for the Bayesian reconstruction entropy is defined
as follows:

$\displaystyle S_{\text{BR}}[A^{+},A^{-}]$ $\displaystyle=$ $\displaystyle\int
d\omega\left[1-\frac{A^{+}}{D}+\log\left(\frac{A^{+}}{D}\right)\right]$ (103)
$\displaystyle+$ $\displaystyle\int
d\omega\left[1-\frac{A^{-}}{D}+\log\left(\frac{A^{-}}{D}\right)\right].$

Since $A=A^{+}-A^{-}$, Eq. (103) can be transformed into:

$\displaystyle S_{\text{BR}}[A,A^{+}]$ $\displaystyle=$ $\displaystyle\int
d\omega\left[1-\frac{A^{+}}{D}+\log\left(\frac{A^{+}}{D}\right)\right]$ (104)
$\displaystyle+$ $\displaystyle\int
d\omega\left[1-\frac{A^{+}-A}{D}+\log\left(\frac{A^{+}-A}{D}\right)\right].$

Then we should eliminate $A^{+}$ in Eq. (104). Because

$\frac{\partial S_{\text{BR}}[A,A^{+}]}{\partial A^{+}}=0=\int
d\omega\left(\frac{1}{A^{+}}+\frac{1}{A^{+}-A}-\frac{2}{D}\right),$ (105)

we immediately get:

$A^{+}=\frac{\sqrt{A^{2}+D^{2}}+D+A}{2},$ (106)

and

$A^{-}=\frac{\sqrt{A^{2}+D^{2}}+D-A}{2}.$ (107)

From the definitions of $A^{+}$ and $A^{-}$, it is easily to prove:

$2A^{+}A^{-}=(A^{+}+A^{-})D,$ (108)

and

$A^{+}+A^{-}=\sqrt{A^{2}+D^{2}}+D.$ (109)

We would like to express $S_{\text{BR}}[A^{+},A^{-}]$ in terms of $A$:

$S_{\text{BR}}[A^{+},A^{-}]=\int
d\omega\left[2-\frac{A^{+}}{D}-\frac{A^{-}}{D}+\log\left(\frac{A^{+}A^{-}}{D^{2}}\right)\right].$
(110)

Substituting Eq. (108) into Eq. (110):

$S_{\text{BR}}[A^{+},A^{-}]=\int
d\omega\left[2-\frac{A^{+}+A^{-}}{D}+\log\left(\frac{A^{+}+A^{-}}{2D}\right)\right].$
(111)

Substituting Eq. (109) into Eq. (111):

$S_{\text{BR}}[A^{+},A^{-}]=\int
d\omega\left[2-\frac{\sqrt{A^{2}+D^{2}}+D}{D}+\log\left(\frac{\sqrt{A^{2}+D^{2}}+D}{2D}\right)\right].$
(112)

Its discretization form reads:

$S_{\text{BR}}[A^{+},A^{-}]=\sum_{m}\left[2-\frac{\sqrt{A_{m}^{2}+D_{m}^{2}}+D_{m}}{D_{m}}+\log\left(\frac{\sqrt{A_{m}^{2}+D_{m}^{2}}+D_{m}}{2D_{m}}\right)\right]\Delta_{m}.$
(113)

Note that Eq. (112) is exactly the same with Eq. (20).
$S_{\text{BR}}[A^{+},A^{-}]$ is just $S^{\pm}_{\text{BR}}[A]$. Then the first
derivative of $S_{\text{BR}}[A^{+},A^{-}]$ with respect to $A_{i}$ is:

$\frac{\partial S[A^{+},A^{-}]}{\partial
A_{i}}=-\Delta_{i}\left(\frac{1}{D_{i}}-\frac{1}{A_{i}^{+}}\right)$ (114)

This equation is very similar to Eq. (90).

###  C.5 Parameterization for positive-negative spectral functions

Now we assume that:

$A^{+}_{l}=\frac{D_{l}}{1-D_{l}\sum_{m}V_{lm}u_{m}},$ (115)

and

$A^{-}_{l}=\frac{D_{l}}{1+D_{l}\sum_{m}V_{lm}u_{m}}.$ (116)

We also introduce:

$w^{+}_{l}=\frac{1}{1-D_{l}\sum_{m}V_{lm}u_{m}},$ (117)

and

$w^{-}_{l}=\frac{1}{1+D_{l}\sum_{m}V_{lm}u_{m}}.$ (118)

So $A^{+}_{l}=D_{l}w^{+}_{l}$, $A^{-}_{l}=D_{l}w^{-}_{l}$, and
$A_{l}=D_{l}(w^{+}_{l}-w^{-}_{l})$. It is easy to verify that the definitions
of $A^{+}_{l}$ and $A^{-}_{l}$ obey Eq. (108). The first derivatives of
$A^{+}_{l}$ and $A^{-}_{l}$ with respect to $u_{i}$ are:

$\frac{\partial A^{+}_{l}}{\partial u_{i}}=A^{+}_{l}A^{+}_{l}V_{li},$ (119)

and

$\frac{\partial A^{-}_{l}}{\partial u_{i}}=-A^{-}_{l}A^{-}_{l}V_{li}.$ (120)

The first derivative of $A_{l}$ with respect to $u_{i}$ is:

$\frac{\partial A_{l}}{\partial
u_{i}}=A^{+}_{l}A^{+}_{l}V_{li}+A^{-}_{l}A^{-}_{l}V_{li}.$ (121)

After some simple algebra, we easily get:

$f_{m}=\alpha u_{m}+\sum_{l}W_{ml}(w^{+}_{l}-w^{-}_{l})-B_{m},$ (122)
$J_{mi}=\alpha\delta_{mi}+\sum_{l}W_{mli}D_{l}(w^{+}_{l}w^{+}_{l}+w^{-}_{l}w^{-}_{l}).$
(123)

## References

  * Negele and Orland (1998) John W. Negele and Henri Orland,  _Quantum Many-Particle Systems_ (Perseus Books, 1998).
  * Coleman (2016) Piers Coleman,  _Introduction to Many-Body Physics_ (Cambridge University Press, 2016).
  * Gull  _et al._ (2011) Emanuel Gull, Andrew J. Millis, Alexander I. Lichtenstein, Alexey N. Rubtsov, Matthias Troyer, and Philipp Werner, “Continuous-time Monte Carlo methods for quantum impurity models,” Rev. Mod. Phys. 83, 349–404 (2011).
  * Foulkes  _et al._ (2001) W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, “Quantum Monte Carlo simulations of solids,” Rev. Mod. Phys. 73, 33–83 (2001).
  * Gubernatis  _et al._ (2016) James Gubernatis, Naoki Kawashima, and Philipp Werner,  _Quantum Monte Carlo Methods: Algorithms for Lattice Models_ (Cambridge University Press, 2016).
  * Metzner  _et al._ (2012) Walter Metzner, Manfred Salmhofer, Carsten Honerkamp, Volker Meden, and Kurt Schönhammer, “Functional renormalization group approach to correlated fermion systems,” Rev. Mod. Phys. 84, 299–352 (2012).
  * C. Platt and Thomale (2013) W. Hanke C. Platt and R. Thomale, “Functional renormalization group for multi-orbital Fermi surface instabilities,” Adv. Phys. 62, 453–562 (2013).
  * Harl and Kresse (2009) Judith Harl and Georg Kresse, “Accurate Bulk Properties from Approximate Many-Body Techniques,” Phys. Rev. Lett. 103, 056401 (2009).
  * Chen  _et al._ (2017) Guo P. Chen, Vamsee K. Voora, Matthew M. Agee, Sree Ganesh Balasubramani, and Filipp Furche, “Random-Phase Approximation Methods,” Annu. Rev. Phys. Chem. 68, 421–445 (2017).
  * Hedin (1965) Lars Hedin, “New Method for Calculating the One-Particle Green’s Function with Application to the Electron-Gas Problem,” Phys. Rev. 139, A796–A823 (1965).
  * Aryasetiawan and Gunnarsson (1998) F Aryasetiawan and O Gunnarsson, “The GW method,” Rep. Prog. Phys. 61, 237 (1998).
  * Onida  _et al._ (2002) Giovanni Onida, Lucia Reining, and Angel Rubio, “Electronic excitations: density-functional versus many-body Green’s-function approaches,” Rev. Mod. Phys. 74, 601–659 (2002).
  * Gukelberger  _et al._ (2015) Jan Gukelberger, Li Huang, and Philipp Werner, “On the dangers of partial diagrammatic summations: Benchmarks for the two-dimensional Hubbard model in the weak-coupling regime,” Phys. Rev. B 91, 235114 (2015).
  * Beach  _et al._ (2000) K. S. D. Beach, R. J. Gooding, and F. Marsiglio, “Reliable Padé analytical continuation method based on a high-accuracy symbolic computation algorithm,” Phys. Rev. B 61, 5147–5157 (2000).
  * Schött  _et al._ (2016a) Johan Schött, Inka L. M. Locht, Elin Lundin, Oscar Grånäs, Olle Eriksson, and Igor Di Marco, “Analytic continuation by averaging Padé approximants,” Phys. Rev. B 93, 075104 (2016a).
  * Vidberg and Serene (1977) H. J. Vidberg and J. W. Serene, “Solving the Eliashberg equations by means of N-point Padé approximants,” J. Low Temp. Phys. 29, 179–192 (1977).
  * Markó  _et al._ (2017) Gergely Markó, Urko Reinosa, and Zsolt Szép, “Padé approximants and analytic continuation of Euclidean $\mathrm{\Phi}$-derivable approximations,” Phys. Rev. D 96, 036002 (2017).
  * Jarrell and Gubernatis (1996) Mark Jarrell and J.E. Gubernatis, “Bayesian inference and the analytic continuation of imaginary-time quantum Monte Carlo data,” Phys. Rep. 269, 133–195 (1996).
  * Gubernatis  _et al._ (1991) J. E. Gubernatis, Mark Jarrell, R. N. Silver, and D. S. Sivia, “Quantum Monte Carlo simulations and maximum entropy: Dynamics from imaginary-time data,” Phys. Rev. B 44, 6011–6029 (1991).
  * Silver  _et al._ (1990) R. N. Silver, D. S. Sivia, and J. E. Gubernatis, “Maximum-entropy method for analytic continuation of quantum Monte Carlo data,” Phys. Rev. B 41, 2380–2389 (1990).
  * Sandvik (1998) Anders W. Sandvik, “Stochastic method for analytic continuation of quantum Monte Carlo data,” Phys. Rev. B 57, 10287–10290 (1998).
  * Sandvik (2016) Anders W. Sandvik, “Constrained sampling method for analytic continuation,” Phys. Rev. E 94, 063308 (2016).
  * Shao and Sandvik (2023) Hui Shao and Anders W. Sandvik, “Progress on stochastic analytic continuation of quantum Monte Carlo data,” Phys. Rep. 1003, 1–88 (2023).
  * Beach (2004) K. S. D. Beach, “Identifying the maximum entropy method as a special limit of stochastic analytic continuation,” (2004), arXiv:0403055 [cond-mat.str-el] .
  * Fuchs  _et al._ (2010) Sebastian Fuchs, Thomas Pruschke, and Mark Jarrell, “Analytic continuation of quantum Monte Carlo data by stochastic analytical inference,” Phys. Rev. E 81, 056701 (2010).
  * Vafayi and Gunnarsson (2007) K. Vafayi and O. Gunnarsson, “Analytical continuation of spectral data from imaginary time axis to real frequency axis using statistical sampling,” Phys. Rev. B 76, 035115 (2007).
  * Reichman and Rabani (2009) David R. Reichman and Eran Rabani, “Analytic continuation average spectrum method for quantum liquids,” J. Chem. Phys. 131, 054502 (2009).
  * Ghanem and Koch (2020a) Khaldoon Ghanem and Erik Koch, “Extending the average spectrum method: Grid point sampling and density averaging,” Phys. Rev. B 102, 035114 (2020a).
  * Ghanem and Koch (2020b) Khaldoon Ghanem and Erik Koch, “Average spectrum method for analytic continuation: Efficient blocked-mode sampling and dependence on the discretization grid,” Phys. Rev. B 101, 085111 (2020b).
  * Syljuåsen (2008) Olav F. Syljuåsen, “Using the average spectrum method to extract dynamics from quantum Monte Carlo simulations,” Phys. Rev. B 78, 174429 (2008).
  * Mishchenko  _et al._ (2000) A. S. Mishchenko, N. V. Prokof’ev, A. Sakamoto, and B. V. Svistunov, “Diagrammatic quantum Monte Carlo study of the Fröhlich polaron,” Phys. Rev. B 62, 6317–6336 (2000).
  * Bao  _et al._ (2016) F. Bao, Y. Tang, M. Summers, G. Zhang, C. Webster, V. Scarola, and T. A. Maier, “Fast and efficient stochastic optimization for analytic continuation,” Phys. Rev. B 94, 125149 (2016).
  * Goulko  _et al._ (2017) Olga Goulko, Andrey S. Mishchenko, Lode Pollet, Nikolay Prokof’ev, and Boris Svistunov, “Numerical analytic continuation: Answers to well-posed questions,” Phys. Rev. B 95, 014102 (2017).
  * Krivenko and Harland (2019) Igor Krivenko and Malte Harland, “TRIQS/SOM: Implementation of the stochastic optimization method for analytic continuation,” Comput. Phys. Commun. 239, 166–183 (2019).
  * Krivenko and Mishchenko (2022) Igor Krivenko and Andrey S. Mishchenko, “TRIQS/SOM 2.0: Implementation of the stochastic optimization with consistent constraints for analytic continuation,” Comput. Phys. Commun. 280, 108491 (2022).
  * Huang and Liang (2023a) Li Huang and Shuang Liang, “Stochastic pole expansion method for analytic continuation of the Green’s function,” Phys. Rev. B 108, 235143 (2023a).
  * Huang and Liang (2023b) Li Huang and Shuang Liang, “Reconstructing lattice QCD spectral functions with stochastic pole expansion and Nevanlinna analytic continuation,” (2023b), arXiv:2309.11114 [hep-lat] .
  * Otsuki  _et al._ (2017) Junya Otsuki, Masayuki Ohzeki, Hiroshi Shinaoka, and Kazuyoshi Yoshimi, “Sparse modeling approach to analytical continuation of imaginary-time quantum monte carlo data,” Phys. Rev. E 95, 061302 (2017).
  * Motoyama  _et al._ (2022) Yuichi Motoyama, Kazuyoshi Yoshimi, and Junya Otsuki, “Robust analytic continuation combining the advantages of the sparse modeling approach and the padé approximation,” Phys. Rev. B 105, 035139 (2022).
  * Fei  _et al._ (2021a) Jiani Fei, Chia-Nan Yeh, and Emanuel Gull, “Nevanlinna Analytical Continuation,” Phys. Rev. Lett. 126, 056402 (2021a).
  * Fei  _et al._ (2021b) Jiani Fei, Chia-Nan Yeh, Dominika Zgid, and Emanuel Gull, “Analytical continuation of matrix-valued functions: Carathéodory formalism,” Phys. Rev. B 104, 165111 (2021b).
  * Huang  _et al._ (2023) Zhen Huang, Emanuel Gull, and Lin Lin, “Robust analytic continuation of Green’s functions via projection, pole estimation, and semidefinite relaxation,” Phys. Rev. B 107, 075151 (2023).
  * Zhang and Gull (2023) Lei Zhang and Emanuel Gull, “Minimal Pole Representation and Controlled Analytic Continuation of Matsubara Response Functions,” (2023), arXiv:2312.10576 [cond-mat.str-el] .
  * Ying (2022a) Lexing Ying, “Pole Recovery From Noisy Data on Imaginary Axis,” J. Sci. Comput. 92, 107 (2022a).
  * Ying (2022b) Lexing Ying, “Analytic continuation from limited noisy Matsubara data,” J. Comput. Phys. 469, 111549 (2022b).
  * Fournier  _et al._ (2020) Romain Fournier, Lei Wang, Oleg V. Yazyev, and QuanSheng Wu, “Artificial Neural Network Approach to the Analytic Continuation Problem,” Phys. Rev. Lett. 124, 056401 (2020).
  * Yoon  _et al._ (2018) Hongkee Yoon, Jae-Hoon Sim, and Myung Joon Han, “Analytic continuation via domain knowledge free machine learning,” Phys. Rev. B 98, 245101 (2018).
  * Huang and Yang (2022) Dongchen Huang and Yi-feng Yang, “Learned optimizers for analytic continuation,” Phys. Rev. B 105, 075112 (2022).
  * Zhang  _et al._ (2022) Rong Zhang, Maximilian E. Merkel, Sophie Beck, and Claude Ederer, “Training biases in machine learning for the analytic continuation of quantum many-body Green’s functions,” Phys. Rev. Res. 4, 043082 (2022).
  * Yao  _et al._ (2022) Juan Yao, Ce Wang, Zhiyuan Yao, and Hui Zhai, “Noise enhanced neural networks for analytic continuation,” Mach. Learn.: Sci. Technol. 3, 025010 (2022).
  * Arsenault  _et al._ (2017) Louis-François Arsenault, Richard Neuberg, Lauren A Hannah, and Andrew J Millis, “Projected regression method for solving Fredholm integral equations arising in the analytic continuation problem of quantum physics,” Inv. Prob. 33, 115007 (2017).
  * Nichols  _et al._ (2022) Nathan S. Nichols, Paul Sokol, and Adrian Del Maestro, “Parameter-free differential evolution algorithm for the analytic continuation of imaginary time correlation functions,” Phys. Rev. E 106, 025312 (2022).
  * Bergeron and Tremblay (2016) Dominic Bergeron and A.-M. S. Tremblay, “Algorithms for optimized maximum entropy and diagnostic tools for analytic continuation,” Phys. Rev. E 94, 023303 (2016).
  * Sim and Han (2018) Jae-Hoon Sim and Myung Joon Han, “Maximum quantum entropy method,” Phys. Rev. B 98, 205102 (2018).
  * Kraberger  _et al._ (2017) Gernot J. Kraberger, Robert Triebl, Manuel Zingl, and Markus Aichhorn, “Maximum entropy formalism for the analytic continuation of matrix-valued Green’s functions,” Phys. Rev. B 96, 155128 (2017).
  * Levy  _et al._ (2017) Ryan Levy, J.P.F. LeBlanc, and Emanuel Gull, “Implementation of the maximum entropy method for analytic continuation,” Comput. Phys. Commun. 215, 149–155 (2017).
  * Assaad  _et al._ (2022) F. F. Assaad, M. Bercx, F. Goth, A. Götz, J. S. Hofmann, E. Huffman, Z. Liu, F. Parisen Toldin, J. S. E. Portela, and J. Schwab, “The ALF (Algorithms for Lattice Fermions) project release 2.0. Documentation for the auxiliary-field quantum Monte Carlo code,” SciPost Phys. Codebases , 1 (2022).
  * Huang (2023) Li Huang, “ACFlow: An open source toolkit for analytic continuation of quantum Monte Carlo data,” Comput. Phys. Commun. 292, 108863 (2023).
  * Haule  _et al._ (2010) Kristjan Haule, Chuck-Hou Yee, and Kyoo Kim, “Dynamical mean-field theory within the full-potential methods: Electronic structure of ${\text{CeIrIn}}_{5}$, ${\text{CeCoIn}}_{5}$, and ${\text{CeRhIn}}_{5}$,” Phys. Rev. B 81, 195107 (2010).
  * Kaufmann and Held (2023) Josef Kaufmann and Karsten Held, “ana_cont: Python package for analytic continuation,” Comput. Phys. Commun. 282, 108519 (2023).
  * Papoulis (1989) Athanasios Papoulis,  _Probability and Statistics_ (Prentice Hall, 1989).
  * De Raychaudhury  _et al._ (2007) M. De Raychaudhury, E. Pavarini, and O. K. Andersen, “Orbital Fluctuations in the Different Phases of ${\mathrm{LaVO}}_{3}$ and ${\mathrm{YVO}}_{3}$,” Phys. Rev. Lett. 99, 126402 (2007).
  * Tomczak and Biermann (2007) Jan M. Tomczak and Silke Biermann, “Effective band structure of correlated materials: the case of VO2,” J. Phys.: Condens. Matter 19, 365206 (2007).
  * Gull and Millis (2014) E. Gull and A. J. Millis, “Pairing glue in the two-dimensional Hubbard model,” Phys. Rev. B 90, 041110 (2014).
  * Reymbaut  _et al._ (2015) A. Reymbaut, D. Bergeron, and A.-M. S. Tremblay, “Maximum entropy analytic continuation for spectral functions with nonpositive spectral weight,” Phys. Rev. B 92, 060509 (2015).
  * Reymbaut  _et al._ (2017) A. Reymbaut, A.-M. Gagnon, D. Bergeron, and A.-M. S. Tremblay, “Maximum entropy analytic continuation for frequency-dependent transport coefficients with nonpositive spectral weight,” Phys. Rev. B 95, 121104 (2017).
  * Gull  _et al._ (2013) Emanuel Gull, Olivier Parcollet, and Andrew J. Millis, “Superconductivity and the Pseudogap in the Two-Dimensional Hubbard Model,” Phys. Rev. Lett. 110, 216405 (2013).
  * Yue and Werner (2023) Changming Yue and Philipp Werner, “Maximum entropy analytic continuation of anomalous self-energies,” (2023), arXiv:2303.16888 [cond-mat.supr-con] .
  * Jaynes (1957) E. T. Jaynes, “Information Theory and Statistical Mechanics,” Phys. Rev. 106, 620–630 (1957).
  * Tikhonov  _et al._ (1995) A. N. Tikhonov, A. V. Goncharsky, V. V. Stepanov, and A. G. Yagola,  _Numerical Methods for the Solution of Ill-Posed Problems_ (Springer Netherlands, 1995).
  * Asakawa  _et al._ (2001) M. Asakawa, Y. Nakahara, and T. Hatsuda, “Maximum entropy analysis of the spectral functions in lattice QCD,” Prog. Part. Nucl. Phys. 46, 459–508 (2001).
  * Burnier and Rothkopf (2013) Yannis Burnier and Alexander Rothkopf, “Bayesian Approach to Spectral Function Reconstruction for Euclidean Quantum Field Theories,” Phys. Rev. Lett. 111, 182003 (2013).
  * Rothkopf (2022) Alexander Rothkopf, “Bayesian inference of real-time dynamics from lattice QCD,” Front. Phys. 10, 1028995 (2022).
  * Rothkopf (2017) Alexander Rothkopf, “Bayesian inference of nonpositive spectral functions in quantum field theory,” Phys. Rev. D 95, 056016 (2017).
  * Nocedal and Wright (2006) Jorge Nocedal and Stephen J. Wright,  _Numerical Optimization_ (Springer New York, NY, 2006).
  * Laue  _et al._ (1985) E. D. Laue, J. Skilling, and J. Staunton, “Maximum entropy reconstruction of spectra containing antiphase peaks,” J. Magn. Reson. 63, 418–424 (1985).
  * Sibisi (1990) Sibusiso Sibisi, “Quantified Maxent: An NMR Application,” in  _Maximum Entropy and Bayesian Methods_, edited by Paul F. Fougère (Springer Netherlands, Dordrecht, 1990) pp. 351–358.
  * Bryan (1990) R. K. Bryan, “Maximum entropy analysis of oversampled data problems,” Eur. Biophys. J. 18, 165–174 (1990).
  * Skilling (1991) J Skilling, “Fundamentals of MaxEnt in data analysis,” in  _Maximum Entropy in Action: A Collection Of Expository Essays_ (Oxford University Press, 1991).
  * Gunnarsson  _et al._ (2010) O. Gunnarsson, M. W. Haverkort, and G. Sangiovanni, “Analytical continuation of imaginary axis data for optical conductivity,” Phys. Rev. B 82, 165125 (2010).
  * Schött  _et al._ (2016b) J. Schött, E. G. C. P. van Loon, I. L. M. Locht, M. I. Katsnelson, and I. Di Marco, “Comparison between methods of analytical continuation for bosonic functions,” Phys. Rev. B 94, 245140 (2016b).


# A Real-Time Wrong-Way Vehicle Detection Based on YOLO and Centroid Tracking
††thanks: This is funded by ICT incubator project in Chittagong University of
Engineering and Technology under the ICT Division of Ministry of Information
and Communication Technology, Bangladesh.

Zillur Rahman1, Amit Mazumder Ami, Muhammad Ahsan Ullah  Department of
Electrical and Electronic Engineering
Chittagong University of Engineering and Technology
Chittagong-4349, Bangladesh
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Wrong-way driving is one of the main causes of road accidents and traffic jam
all over the world. By detecting wrong-way vehicles, the number of accidents
can be minimized and traffic jam can be reduced. With the increasing
popularity of real-time traffic management systems and due to the availability
of cheaper cameras, the surveillance video has become a big source of data. In
this paper, we propose an automatic wrong-way vehicle detection system from
on-road surveillance camera footage. Our system works in three stages: the
detection of vehicles from the video frame by using the You Only Look Once
(YOLO) algorithm, track each vehicle in a specified region of interest using
centroid tracking algorithm and detect the wrong-way driving vehicles. YOLO is
very accurate in object detection and the centroid tracking algorithm can
track any moving object efficiently. Experiment with some traffic videos shows
that our proposed system can detect and identify any wrong-way vehicle in
different light and weather conditions. The system is very simple and easy to
implement.

###### Index Terms:

vehicle, YOLO, centroid, wrong-way, computer vision

††publicationid: pubid: 978-1-7281-7366-5/20/$31.00 ©2020 IEEE

##  I Introduction

Road accident is a very common issue in a densely populated country like
Bangladesh. In 2019, there were 4702 road accidents resulting in 5227 deaths
including many children and students[1]. The capacity of the roadway is not
sufficient for the growing number of vehicles and thus imbalance is created.
The drivers do not follow traffic rules and take advantage of driving in the
wrong-side in times of red traffic signals. It increases traffic on one side
and hampers traffic flow greatly [2]. It also increases the possibility of
head-on collision several times. About 355 people die every year due to the
crashes in wrong-way driving in the United States [3]. Fig. 1 shows the
traffic jam caused by wrong-way driving. This is a common scenario in most of
the cities of Bangladesh. So, it is essential to stop drivers from driving on
the wrong side. To ensure it, those who don’t follow traffic rules need to
find out and strict law should be applied.

Figure 1: Heavy traffic jam due to wrong-way driving[4]

Computer vision is used in Intelligent Transportation Systems (ITS) since the
beginning of this technology. It is very useful in the traffic surveillance
system. The video footage from the on-road cameras can be analyzed to find out
the cause of road accidents, traffic jam and also traffic rules violations,
counting vehicle and speed. In the past times, it was very difficult to use
this technology due to the expensive computation and unavailability of data.
But, the recent improvement of different machine learning and deep learning
algorithm and availability of powerful Graphics Processing Unit (GPU) and the
cheaper camera has made the whole system more efficient.

The system we present here has three stages. In the first stage, every vehicle
in the video frame is detected using the YOLO object detection algorithm and a
bounding box is generated for each detected vehicle. Then, the bounding boxes
are fed to the centroid based moving object tracking algorithm. The algorithm
tracks each vehicle independently in a specified region of interest (ROI).
Finally, the direction of the vehicle is determined by calculating its
centroid’s height in each frame and detect whether it moves in the wrong
direction or not. If the vehicle is on the wrong side, then the system will
capture an image of the vehicle.

In the next section, the previous works in this field are described. In
section III, the details of the YOLO algorithm are described. In section IV,
the details of our proposed system are described. The results are presented in
section V and section VI describes the conclusion.

##  II RELATED WORKS

Computer vision technology is being used in many intelligent traffic
monitoring systems. So far, many systems have been developed for wrong-way
vehicle detection. The existing methods can be classified into sensor-based
and motion pattern-based. Sensor-based methods [5] use a magnetic sensor to
detect the direction of vehicles. The earth’s magnetic field gets disturbed in
the presence of a vehicle and the sensor gives the signal according to the
variation of the magnetic field. From the signal, the direction of the vehicle
can be detected. But in this system [5], the magnetic field can be varied by
other reasons too. In that case, it won’t give a satisfactory result. The
motion pattern-based method uses optical flow measurement. The method proposed
in [6] uses optical flow calculation to detect the direction and compare it
with the modeled lane direction. But this method faces a problem due to
occlusion. Improved optical flow estimation [7] was used and the authors used
background subtraction and Lucas-Kanade method to detect the wrong-way
vehicles. But this system is not going to provide a satisfactory result in
different light conditions because the background subtraction method is highly
affected by the shadow of the vehicle itself.

##  III YOLO

There are a lot of methods available to detect the vehicles with bounding
boxes. Among those methods, for real-time applications, Mask-RCNN [8] and
YOLO[9] are most popular. The Mask-RCNN method creates bounding boxes using a
process called selective search. But for real-time application, YOLO provides
more speed than Mask-RCNN and it is the state of the art. There are several
versions of YOLO. In this paper, we used the most recent YOLOv3 [10]. The main
thing about YOLO is that it looks at the image once and derives each object
with class probability and bounding box and that is why it is very fast
comparing with other object detection algorithm. First, it divides the images
into the SxS grid cell. Each cell can have B bounding boxes [9]. If there are
C number of classes then, the total bounding box will be SxSxB and each box
will have (5+C) attributes. Four attributes for the box coordinate and one is
for confidence score that this box contains one of the C classes. There is a
threshold confidence score. For our system, we used 0.5 as the threshold. So,
it will only show the box with a confidence score of more than 0.5 with the
name of the class. All of these operations are completed in a single pass
through the network. There is a problem when multiple detections of the same
object happen. To overcome this, a process called non-max suppression is used.
In this process, first, a box with the highest confidence score is selected.
Then, the boxes which have overlap more than a threshold value with the
previously selected box will be removed. This process ensures only one box per
object.

Figure 2: YOLO divides the image into the SxS grid and predicts B bounding box
for each cell [9]

In YOLOv3, the authors present a new deeper convolutional neural network (CNN)
to extract the features. It has 53 convolutional layers. Batch norm and Leaky
ReLU activation function were applied to each convolutional layer. The
features of the input image are extracted by these convolutional layers. No
softmax layer was used because it was not necessary. For the class prediction,
binary cross-entropy was used during the training.

YOLOv3 is trained on the COCO dataset which has 80 different classes. For our
system, we only used different classes of vehicles including motorbike, bus,
truck, and cars. The previous versions of YOLO did not give good results for
small objects but YOLOv3 shows very good results in that case.

##  IV PROPOSED METHOD

This section describes the details of the whole proposed system. The details
of each stage will be described consecutively in the following subsections.

###  IV-A Virtual Region of Interest

We want to track the vehicles in a certain region of the video frame. That is
why we selected a rectangular area in the frame which is shown in Fig. 3 in
red color. When the vehicle is in the region, it will be tracked. When the
vehicle is out of this region, it will be removed from the tracked vehicle
list. As this region depends on the camera view, it must be done manually. The
region should have a reasonable physical distance or there will be a problem
to track the large vehicles.

Figure 3: The region of interest for vehicle tracking

###  IV-B Vehicle Tracking

To track each vehicle, we use the centroid tracking algorithm. This algorithm
takes the bounding box as the input. So first, the bounding boxes are
generated using YOLO. Then, those boxes are fed to the centroid tracker. When
the center of each vehicle that means the center of the corresponding bounding
box enters the region of interest, it is given a unique identification number
which is shown in Fig. 4(a). In the next frame, the center of all the objects
move in another place or maybe not have any movement which is shown in Fig.
4(b). The centroid tracking algorithm is based on an assumption which is that
each object will move very little in between the subsequent frame. So, if we
can relate any new centroid which has the minimum distance with an old
centroid, we can say that this object is previously identified and the new
centroid of that object will be updated. This is shown in Fig. 5(a). To do
this, all possible Euclidean distance between each pair of the new centroids
(yellow color) and the old centroids (red color) are computed. The Euclidean
distance d between two pixels $(x_{i},y_{i})and(x_{j},y_{j})$ is as follow:

$d=\sqrt{(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}}$ (1)

(a)

(b)

Figure 4: (a) The first frame (b) second frame

(a)

(b)

Figure 5: (a) updated frame (b) when a new frame has more object

If the current frame has fewer objects than the previous frame, which is shown
in Fig. 4(b), that means one or more objects are disappeared from the region
of interest. As their centroids are not in the region of interest, the ID of
those vehicles will be removed from the tracked vehicle list.

But, when the number of objects in the current frame becomes more than the
previous frame which is shown in Fig. 5(b) then there must be a new object.
The algorithm will be updated and old objects will be assigned new centroids.
The remaining objects will be identified as new objects and given a new
identity number.

In this way, our proposed system tracks each vehicle independently which is
shown in Fig. 06. The figure shows that a CNG vehicle is in the region of
interest. The algorithm successfully tracked that vehicle and gave an identity
number. On the right side of the image, it is seen that a pedestrian is
detected by the YOLO but is not tracked. Because, we feed only the vehicle
classes to the centroid tracker and remove others.

Figure 6: A CNG vehicle is tracked

###  IV-C Wrong-way vehicle detection

The final step of our proposed system is to detect the wrong-way vehicle. Our
system already can track each vehicle that is in our specified region of
interest. The centroid of every tracked vehicle has a height from the top of
the frame. When the vehicle is registered first and given an identity number,
the height of the centroid $H_{1}$ is computed and stored in a file
corresponding to the identity number.

In the next frame, the height of the centroid $H_{2}$ is computed and stored
in another file along with its identity number. This $H_{2}$ will be updated
in each consecutive frame. If the vehicles move, the $H_{1}$ and $H_{2}$ of a
vehicle will not be equal. By comparing these two heights, our system will
predict the direction of the vehicle. In our system, we defined that if the
vehicle moves away from the camera, it will be detected as a wrong-way
vehicle. So, if $H_{1}<H_{2}$ then, the vehicle is coming towards the camera
and is in the right way. If otherwise, our system will detect it as a wrong-
way vehicle. The opposite can also be defined just by changing the condition.
After the detection of such a vehicle, an image of the frame will be captured
automatically for further inspection. The flowchart of the whole system is
shown in Fig. 7.

Figure 7: The flowchart of the whole system

## RESULTS

To check our system, we captured three videos from Chittagong city of
Bangladesh each including only one wrong-way vehicle. The videos were captured
from roadside and resolution was 1280 x 720 pixels. In Fig. 8, it is seen that
our system successfully detected a wrong-way vehicle. Each wrong-way vehicle
from three videos is successfully detected. So, the accuracy of the system is
almost 100%. As the other related systems [5] [6] [7] are not based on deep
learning algorithm, it is not possible to examine them with our used dataset
and experimental setup and also the authors of those work did not mention the
accuracy of their system specifically, so it would not be appropriate to
compare them with our system.

We designed the system to detect the wrong-way vehicle from only one side of
the road. The camera will focus on only one side. In our future work, we will
extend this and implement the system on both sides of the road using a single
camera. In that case, the camera will be placed in the middle of the road
focusing on both sides.

The system was performed in Python language with OpenCV deep learning library.
As the YOLO algorithm is computationally expensive, the frame per second is
2.5 in AMD Ryzen 7 CPU in the windows platform. As it is not necessary to take
into account every frame, we took one frame and skip the next four. The GPU
version is much faster and by using an Nvidia GPU in the Linux platform, the
speed can be increased by several times.

Figure 8: A motorbike is driving on the wrong way

## CONCLUSIONS

This paper presented a system that can detect wrong-way vehicles and mark them
from on-road video footage. In the first stage, every vehicle in the video
frame is detected using the YOLO object detector as it is very accurate and
faster than any object detector algorithm. Then, the bounding boxes generated
by the YOLO are fed to the centroid based tracking algorithm. The tracking
algorithm tracks each vehicle in the specified region of interest. Then, by
computing the centroid height of each vehicle in consecutive frames, the
direction of vehicles can be determined. The system shows a very promising
result as we tested it using traffic footage we collected from the roads. The
whole system’s accuracy depends on the detection of the YOLO. If YOLO detects
correctly, our proposed method will have an accuracy of almost 100%.

The limitation of this system is with the centroid tracking algorithm. The
centroids of the object must lie close together between subsequent frames or
the ID number might be switched due to overlap of one object to another.
Though this problem doesn’t affect our system much.

## References

  * [1] “Bangladesh road accidents in 2019 statistics by Nirapad Sarak Chai,”The Daily Star,05-Jan-2020. [Online]. Available: https://www.thedailystar.net/frontpage/bangladesh-road-accidents-in-2019-stats-by-nirapad-sarak-chai-1849588.[Accessed Feb. 20, 2020].
  * [2] “Traffic Jam in Dhaka: ‘Driving on wrong side a disease’, says Obaidul.”The Daily Star,29-May-2018.[Online].Available:https://www.thedailystar.net/country/dhaka-traffic-jam-driving-wrong-side-disease-obaidul-quader-1583323. [Accessed Feb. 21, 2020].
  * [3] M. Pour-rouholamin and J. Shaw, “Current Practices of Safety Countermeasures for Wrong-Way Driving,” TRB 94th Annu. Meet. Compend. Pap., no. January, pp. 1–14, 2015.
  * [4] “ways to avoid traffic congestion Archives - GoMetro Transport App.” [Online]. Available: http://www.getgometro.com/tag/ways-to-avoid-traffic-congestion/. [Accessed: 25-Feb-2020].
  * [5] M. J. Caruso and L. S. Withanawasam, “Vehicle Detection and Compass Applications using AMR Magnetic Sensors,” Sensors Expo, vol. 477, p. 39, 1999.
  * [6] “WRONG WAY DRIVERS DETECTION BASED ON OPTICAL FLOW Gonc Institute for System and Robotics Dep . of Electrical Engineering and Computers University of Coimbra - Portugal,” no. 1, pp. 141–144, 2007.
  * [7] S. V.-U. Ha, H.-H. Nguyen, H. M. Tran, and P. Ho-Thanh, “Improved optical flow estimation in wrong way vehicle detection,” J. Inf. Assur. Secur., vol. 9, no. 5, pp. 165–169, 2014.
  * [8] K. He, G. Gkioxari, P. Dollar, and R. Girshick, “Mask r-CNN,” in 2017 IEEE International Conference on Computer Vision (ICCV), 2017, pp. 2980–2988.
  * [9] J. Redmon, S. Divvala, R. Girshick, and A. Farhadi, “You only look once: Unified, real-time object detection,” in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 779–788.
  * [10] J. Redmon and A. Farhadi, “Yolov3: An incremental improvement,” arXiv preprint arXiv:1804.02767,2018.