# Källén–Lehmann Spectral Representation of the Scalar SU(2) Glueball

David Dudala,b<EMAIL_ADDRESS>a KU Leuven Campus Kortrijk–Kulak,
Department of Physics, Etienne Sabbelaan 53 bus 7657, 8500 Kortrijk, Belgium
b Ghent University, Department of Physics and Astronomy, Krijgslaan 281-S9,
9000 Gent, Belgium
c CFisUC, Department of Physics, University of Coimbra, 3004-516 Coimbra,
Portugal Orlando Oliveirac<EMAIL_ADDRESS>a KU Leuven Campus Kortrijk–Kulak,
Department of Physics, Etienne Sabbelaan 53 bus 7657, 8500 Kortrijk, Belgium
b Ghent University, Department of Physics and Astronomy, Krijgslaan 281-S9,
9000 Gent, Belgium
c CFisUC, Department of Physics, University of Coimbra, 3004-516 Coimbra,
Portugal Martin Roelfsa<EMAIL_ADDRESS>Corresponding author a KU
Leuven Campus Kortrijk–Kulak, Department of Physics, Etienne Sabbelaan 53 bus
7657, 8500 Kortrijk, Belgium
b Ghent University, Department of Physics and Astronomy, Krijgslaan 281-S9,
9000 Gent, Belgium
c CFisUC, Department of Physics, University of Coimbra, 3004-516 Coimbra,
Portugal

###### Abstract

The estimation of the Källén-Lehmann spectral density from gauge invariant
lattice QCD two point correlation functions is proposed, and explored via an
appropriate inversion method. As proof of concept the SU(2) glueball spectrum
for the quantum numbers $J^{PC}=0^{++}$ is investigated for various values of
the lattice spacing. The spectral density and the glueball spectrum are
estimated using the published data of [1]. Our estimates for the ground state
mass are in good agreement with the traditional approach published therein,
which is based on the large time exponential behaviour of the correlation
functions. Furthermore, the spectral density also contains hints of excites
states in the spectrum and, despite using non-optimized gauge configurations,
the method is able to estimate their mass values.

##  1 Introduction

A major effort to understand the dynamics of strong interactions is the
computation of the hadronic spectra. In QCD, hadrons, namely mesons and
baryons, are seen as bound states of quarks and gluons. Besides the
conventional hadronic states, multiple-quark and pure glue bound states are
also predicted by theory. Of these exotic states, so far only multi-quark
states have been identified experimentally [2]. For glueballs there are a
number of candidate states [3, 4, 5], but no established observation.
Glueballs continue to attract a great deal of theoretical [6, 7, 8, 9, 10, 11,
12] and experimental attention, with an ongoing effort to identify these pure
glue states unambiguously [13, 14, 15]. Additionally, glueball states have
also been discussed as possible candidates for dark matter [1, 16].

A first principles approach to access the hadron spectra, which has been
advancing since the beginning of the eighties from last century [17], is
lattice QCD. In a typical lattice QCD computation of the bound state mass, an
appropriate two point correlation function is evaluated and, from its large
time decay behaviour and the associated slope, the ground state mass is
estimated [18]. It is common practise in the computation of bound masses to
rely on techniques that improve the signal to noise ratio of the Monte Carlo
simulation. These techniques allow access not only to the ground state mass
with the chosen quantum numbers, but, sometimes, also allow the first excited
state to be extracted. In practice however, the computation of the excited
states masses has proven to be a difficult task.

Herein, we aim to discuss an alternative way of accessing the particle masses:
via the computation of the Källén-Lehmann spectral representation associated
with the momentum space particle propagator [19]. A possible advantage of
using spectral representations compared to a conventional lattice calculation
is that it does not necessarily require the use of smearing, or other
techniques, to improve the Monte Carlo signal to noise ratio, see e.g. [20,
21] and references therein. As seen below, the computation of the spectral
function allows simultaneous extraction of the ground state and, to some
extent, also of the 1st excited state.

The interest in the Källén-Lehmann spectral representation goes beyond the
determination of particle spectra. Besides accessing the spectra at zero and
finite temperature, the spectral representation is linked with the analytical
structure of the associated propagator [22], and it allows the computation of
thermodynamical and transport properties [23]. Furthermore, in QCD, or any
other confining theory, it can help to understand the confinement mechanism
[24, 25, 26, 27, 28, 29, 30].

In the current work we focus on the SU(2) glueball states with quantum numbers
$J^{PC}=0^{++}$, i.e. the scalar glueball. However, in principle, the
procedure can be extended to other quantum numbers and other gauge groups. As
discussed below, the masses of the glueball $0^{++}$ obtained from the
spectral function are in good agreement with the estimates of a more
conventional mass calculation.

Let us describe our procedure to access the Källén-Lehmann spectral
representation $\rho(\omega)$ from a two point correlation function
$G(p^{2})$. The relation between these two functions reads

$G(p^{2})=\int_{0}^{\infty}\frac{2\omega\rho(\omega)\differential{\omega}}{\omega^{2}+p^{2}}=\int_{-\infty}^{\infty}\frac{\rho(\omega)\differential{\omega}}{\omega-
ip}\ .$

This definition assumes that the integrations are well defined [31] and,
therefore, $\rho(\omega)$ approaches zero sufficiently fast as
$|\omega|\rightarrow\infty$. Certain correlators contain polynomial terms that
diverge for large $p$ as happens e.g. for the glueball operator, something
which follows immediately from power counting and a dimensional analysis.
However in these cases it is still possible to write down a sensible
Källén–Lehmann spectral representation, if one first subtracts the polynomial
part [31]. This corresponds to adding appropriate contact counterterms. So, if
$G(p)$ has a polynomial part, the propagator can be written as

$\displaystyle G(p^{2})$
$\displaystyle=\sum_{k=0}^{n-1}a_{k}(p^{2}-\bar{p}^{2})^{k}+(-p^{2}+\bar{p}^{2})^{n}\int_{0}^{\infty}\frac{2\omega\tilde{\rho}(\omega)\differential{\omega}}{\omega^{2}+p^{2}}$
(1)
$\displaystyle=\sum_{k=0}^{n-1}a_{k}(p^{2}-\bar{p}^{2})^{k}+(-p^{2}+\bar{p}^{2})^{n}\int_{-\infty}^{\infty}\frac{\tilde{\rho}(\omega)\differential{\omega}}{\omega-
ip}$

with

$\displaystyle a_{n}$
$\displaystyle=\frac{1}{n!}\partialderivative[n]{G(p^{2})}{(p^{2})}\Bigg{|}_{p^{2}=\bar{p}^{2}},$
(2) $\displaystyle\tilde{\rho}(\omega)$
$\displaystyle=\frac{\rho(\omega)}{(\omega^{2}+\bar{p}^{2})^{n}},$ (3)

and $\bar{p}^{2}$ is a reference momentum scale at which the subtraction is
done. A derivation of the above relations can be found in appendix A. For the
scalar glueball a dimensional analysis [10] shows that $n=3$. The spectral
function $\rho(\omega)$ can be obtained, if the subtraction of the polynomial
part can be performed. However, performing subtractions on numerical data is
extremely sensitive to the choice of $\bar{p}$ and results in relative rapid
variation of $\\{a_{k}\\}$ with the subtraction point [32].

An elegant way to perform the subtractions, without actually having to do
these, is by considering the Fourier transform of $G(p^{2})$ and look at the
Schwinger function defined as

$C(\tau)=\mathcal{F}\left\\{G(p)\right\\}(\tau)=\int_{-\infty}^{\infty}G(p^{2})\bigg{|}_{\vec{p}=0,p_{4}\neq
0}~{}e^{-ip_{4}\tau}\differential{p_{4}}.$ (4)

Setting the subtraction reference momentum $\bar{p}=0$, yields

$\displaystyle C(\tau)$ $\displaystyle=\mathcal{F}\left\\{G(p)\right\\}(\tau)$
(5)
$\displaystyle=\mathcal{L}\left\\{\rho(\omega)\right\\}(|\tau|)+2\pi\sum_{k=0}^{n-1}a_{k}(-1)^{k}\delta^{(2k)}(\tau)$
(6)
$\displaystyle\quad+4\pi(-1)^{n+1}\sum_{k=2}^{n}\delta^{(2(n-k)+2)}(\tau)\int_{0}^{\infty}\differential{\omega}\omega^{2k-3}\tilde{\rho}(\omega),$

see appendix A. The important observation however, is that
$C(\tau)=\mathcal{L}\left\\{\rho(\omega)\right\\}(|\tau|)$ when
$\tau~{}\neq~{}0$, and equal to a sum of (derivatives of) Dirac delta
functions when $\tau=0$. Therefore, $\rho(\omega)$ can be recovered by taking
$C(\tau)$ for $\tau>0$, and inverting the Laplace transformation.

Because the inverse Laplace transform is an ill-defined numerical problem,
regularization is necessary in order to perform the inversion. We shall use
Tikhonov regularization, similar to our previously published method [25, 33].
However, because glueballs are observable particles, their spectral density
$\rho(\omega)$ is non-negative. Therefore we shall implement Tikhonov
regularization using Non-Negative Least Squares (NNLS) [34], to ensure a
positive spectral function $\rho(\omega)\geq 0$. Tikhonov regularization is
not the only possible way to regularize the inversion, and a number of
different regularization strategies have been explored in literature to access
the spectral function for particle correlators by various authors [35, 36, 37,
38, 39, 40, 30, 28, 41, 42]. An advantage of the Tikhonov regularization being
that it keeps the optimization function a quadratic function, which translates
into solving a modified linear system of equations. The most widely used
approach is based on the maximum entropy method [35]. We will therefore use
the same litmus test as in [35] to benchmark our approach, before turning to
the actual glueball case.

##  2 The Numerical Method

In a lattice simulation the propagator111Note the change in notation from
$G(p^{2})$ to $G(p)$, where $p$ now stands for the four component of the
momentum. See Eq. (4) for the definition of the Schwinger function. $G(p_{n})$
is computed on a finite set of evenly spaced momenta $p_{n}$. Given a data set
$\\{G(p_{n})\\}:=\\{G(p_{0}),\ldots,G(p_{N-1})\\}$, the Schwinger function is
computed using DFFT, resulting in a data set $\\{C(\tau_{k})\\}$, where

$C(\tau_{k})=\sum_{n=0}^{N-1}G(p_{n})e^{-i2\pi kn/N}.$

The Laplace transformation to access the spectral function

$\displaystyle
C(\tau_{k})=\mathcal{L}\left\\{\rho(\omega)\right\\}(\tau_{k})=\int_{0}^{\infty}e^{-\omega\tau_{k}}\rho(\omega)\differential{\omega},$
(7)

can be written as a matrix equation

$\boldsymbol{C}=\boldsymbol{K}\boldsymbol{\rho},$

with the elements of $\boldsymbol{K}$ defined as

$\boldsymbol{K}_{k\ell}:=e^{-\omega_{\ell}\tau_{k}}\Delta\omega.$

Since $\boldsymbol{\rho}$ needs to be obtained, and a direct solution is
impossible due to the near zero singular values of $\mathbf{K}$, the original
problem is replaced by the minimisation of the Tikhonov regularizing
functional

$J_{\alpha}=\norm{\boldsymbol{K}\boldsymbol{\rho}-\boldsymbol{C}}_{2}^{2}+\alpha^{2}\norm{\boldsymbol{\rho}-\boldsymbol{\rho}^{*}}_{2}^{2}$
(8)

where $\alpha^{2}>0$ is the Tikhonov parameter and $\boldsymbol{\rho}^{*}$ is
a prior estimate for $\boldsymbol{\rho}$. In order to impose the constraint
$\boldsymbol{\rho}\geq 0$, define

$\displaystyle\boldsymbol{A}=\matrixquantity(\boldsymbol{K}\\\
\alpha\mathbb{1}),\quad\boldsymbol{b}=\matrixquantity(\boldsymbol{C}\\\
\alpha\boldsymbol{\rho}^{*}),$ (9)

and rewrite $J_{\alpha}$ as

$J_{\alpha}=\norm{\boldsymbol{A}\boldsymbol{\rho}-\boldsymbol{b}}_{2}^{2}.$
(10)

It is straightforward to show that (8) and (10) are equivalent. However, the
formulation of the problem using Eq. (10) can be solved with a non-negative
least squares (NNLS) solver such that $\boldsymbol{\rho}\geq 0$ is guaranteed
[34].

#### Determination of $\alpha$

The regularization parameter $\alpha$, in essence, provides a soft threshold
to the singular values of $\boldsymbol{K}$, such that the smallest singular
values no longer cause numerical issues. Choosing $\alpha$ is a delicate
affair, since setting it too small means the problem remains ill-defined,
whereas setting it too large destroys a lot of the information contained in
the data.

Our preferred criterion for $\alpha^{2}$ relies on the Morozov discrepancy
principle [33], which states that $\alpha^{2}$ should be chosen such that

$\norm{\boldsymbol{K}\boldsymbol{\rho}-\boldsymbol{C}}_{2}^{2}=\sum_{i}\sigma_{i}^{2},$
(11)

where $\sum_{i}\sigma_{i}^{2}$ is the total variance in the data. The
$\alpha^{2}$ obeying eq. 11 is guaranteed to be unique [43], and means that
the quality of the reconstruction is identical to the quality of the data.

#### Construction of $\boldsymbol{K}$

The matrix $\boldsymbol{K}$ should perform the Laplace transform as truthfully
as possible, and therefore $\omega$ should range from
$[0,\omega_{\text{max}})$, where $\omega_{\text{max}}$ is sufficiently large
compared to any features that might appear in $\rho(\omega)$. In order to
ensure this we choose to sample $\omega$ evenly in logarithmic space from
$[10^{-5},10^{5}]$ $\text{\,}\mathrm{GeV}$ in $N_{\omega}$ steps. However, it
should be remarked that the peak positions are very consistent, in both linear
and logarithmic space, provided $\omega_{\text{max}}$ is large enough. For the
problems under consideration, the peak positions do not change significantly
after $\omega>100$ $\text{\,}\mathrm{GeV}$. A numerical study is provided in
appendix B.

##  3 Results and Discussion

In this section we detail the results of applying the method of section 2 to
various mock and real data. In section 3.1 the method is applied to a toy
model based on a vector-meson spectral density to establish its reliability.
Then, in section 3.2 the method is applied to recent lattice SU(2) propagator
data for the $0^{++}$ glueball [1].

###  3.1 Meson toy-model

Figure 1: Reconstruction the toy-model spectral density function for various
$d$ and $N$. The dashed black line is the original spectral function, while
the blue full curve is the reconstructed spectral function as given by the
Tikhonov regularized NNLS method.

In order to investigate the reliability of the method we consider a realistic
toy-model, based on a vector-meson model decay into hadrons, as used before in
[35]. This particular model needs a single subtraction, and therefore provides
an excellent test of the method. To allow comparison of our results to those
of [35], we shall use the exact same process to generate the mock data. The
meson spectral density function is given by

$\rho(\omega)=\frac{2}{\pi}\bqty{F^{2}_{\rho}\frac{\Gamma_{\rho}m_{\rho}}{(\omega^{2}-m_{\rho}^{2})^{2}+\Gamma_{\rho}^{2}m_{\rho}^{2}}+\frac{1}{8\pi}\pqty{1+\frac{\alpha_{s}}{\pi}}\frac{1}{1+e^{(\omega_{0}-\omega)/\delta}}},$
(12)

with an energy-dependent width

$\Gamma_{\rho}(\omega)=\frac{g_{\rho\pi\pi}^{2}}{48\pi}m_{\rho}\pqty{1-\frac{4m_{\pi}^{2}}{\omega^{2}}}^{3/2}\theta\pqty{\omega-2m_{\pi}}.$
(13)

The empirical values of the parameters are

$\displaystyle m_{\rho}$ $\displaystyle=0.77\text{ GeV}\quad$ $\displaystyle
m_{\pi}=0.14\text{ GeV}$ $\displaystyle g_{\rho\pi\pi}$ $\displaystyle=5.45$
$\displaystyle F_{\rho}=\frac{m_{\rho}}{g_{\rho\pi\pi}}$
$\displaystyle\omega_{0}$ $\displaystyle=1.3\text{ GeV}$
$\displaystyle\delta=0.2\text{ GeV}$

As $\omega\to\infty$, this model behaves like
$\rho(\omega\to\infty)=(1/4\pi^{2})(1+\alpha_{s}/\pi)$. Therefore the integral
section 1 does not converge, and a single subtraction has to be performed,
that means $n=1$ in the notation used in eq. 2.

Assuming $\alpha_{s}=0.3$, the value $(1/4\pi^{2})(1+\alpha_{s}/\pi)=0.0277$
can be used as the prior, but identical to [35] we shall use the slightly
smaller value $\rho_{\text{prior}}~{}=~{}0.0257$.

In order to generate mock data, we compute $C_{\text{orig}}(\tau_{k})$ as the
Laplace transform of $\rho(\omega)$, on $N$ points $\tau_{k}$ spaced by
$\Delta\tau=0.085\text{ fm}=0.43078\text{ GeV}^{-1}$. The standard deviation
of the noise is chosen as

$\sigma(\tau_{k})=d\,C_{\text{orig}}(\tau_{k})\frac{\tau_{k}}{\Delta\tau},$

where $d$ is a parameter which controls the noise level, identical to that of
[35]. The mock data set is then generated as

$C(\tau_{k})=\mathcal{N}(\mu=C_{\text{orig}}(\tau_{k}),\sigma^{2}=\sigma\pqty{\tau_{k}}^{2}).$

These mock data sets were then inverted, ignoring $C(\tau_{0})$, for various
values of $N$ and $d$, using the method of section 2, to test the robustness
of the method. No inversions without positivity constraints were performed, as
significant positivity violations were observed in initial trials. Here we
used $N_{\omega}=1000$ in the construction of $\boldsymbol{K}$. The results
are shown in fig. 1. Both more data-points, or less noise, are found to
improve the quality of the reconstruction. Interestingly, more data-points do
lead to a greater variance in the reconstructions, though their quality is
improved with the number of data-points. A direct comparison with [35, Figure
4] is complicated by the absence of uncertainties on that Figure, but the
performance of the methods seems comparable using the naked eye.

###  3.2 SU(2) Glueball data

Yamanaka et al. [1], see also [44], have provided us with the Schwinger
functions for the SU(2) pure Yang-Mills glueballs, using lattice simulations
for $\beta=$ 2.1, 2.2, 2.3, 2.4, 2.5. The Schwinger functions were computed
using the raw data, that is, without any smearing applied. The lattice
volumes, the number of configurations, the number of Schwinger function time
slices $N$, and the uncertainties in each data set, are shown in table 1.

The simulations of Yamanaka et al. rely on large ensembles of configurations,
resulting in data sets with very low statistical uncertainties. This makes the
results of the Tikhonov regularized inversion reliable, as seen in section
3.1. Additionally, the number of Schwinger function time slices $N$ for these
data sets is $7$, $9$, or $13$, meaning the variance in the inversion will be
small, as was observed in section 3.1. For further details about the data
sets, we refer the reader to [1].

$\beta$ | Volume | Configurations | $N$ | $\max(C(\tau))/\expectationvalue{\sigma_{C}(\tau)}$ | $a\sqrt{\sigma}$ | $a[\Lambda^{-1}]$
---|---|---|---|---|---|---
$2.1$ | $10^{3}\times 12$ | 1,000,000 | 7 | $3.28\times 10^{-4}$ | $0.608(16)$ | $0.356(27)$
$2.2$ | $12^{4}$ | 9,999,990 | 7 | $1.04\times 10^{-4}$ | $0.467(10)$ | $0.273(20)$
$2.3$ | $14^{3}\times 16$ | 4.100,000 | 9 | $1.38\times 10^{-4}$ | $0.3687(22)$ | $0.216(15)$
$2.4$ | $16^{3}\times 24$ | 2,030,000 | 13 | $1.55\times 10^{-4}$ | $0.2660(21)$ | $0.156(11)$
$2.5$ | $20^{3}\times 24$ | 520,000 | 13 | $3.04\times 10^{-4}$ | $0.1881(28)$ | $0.110(8)$

Table 1: The glueball data sets of [1]. Selection of information taken from
[1, Table I, Table II].

figure[.45][][b]  figure[.45][][b]

Figure 2: $C(\tau)/C(\tau_{1})$ for all data sets.  Figure 3: $\rho(\omega)$
for all data sets, using the $ip$-method. The $y$-axis has been normalized to
give the last peak an intensity of 1. Clearly positivity constraints should be
imposed.

The Schwinger functions for the various simulations can be seen in fig. 3,
whereas fig. 3 shows the spectral density functions as obtained without any
constraints on $\rho(\omega)$, using the $ip$-method published in [33].

Here we used $N_{\omega}=10000$ in the construction of $\boldsymbol{K}$ to
create high quality reconstructions. However, this is unnecessarily large, and
was done purely to improve the aesthetics of fig. 5. See appendix B for the
dependence of the ground state mass on $N_{\omega}$. As fig. 3 illustrates,
there are significant positivity violations and rapid oscillations when
positivity is not enforced, especially for small $\omega$, a tell-tale sign of
overfitting. On the other hand, when positivity is imposed, the resulting
spectral functions in fig. 5 all display a clear mass gap that correspond to a
ground state mass of $\sim 1.4\text{\,}\mathrm{GeV}$. Moreover, the infrared
oscillations are gone.

The exact mass values are given in table 2, which lists the $0^{++}$ ground
state mass as obtained from the spectral density function, as well as those
calculated by Yamanaka et al. [1].

figure[][][t]  table[][][t]  $\beta$ $\omega\pm\text{HWHM}\;/\text{GeV}$ 2.1
$1.382+0.049-0.041$ $3.712+1.048-0.486$ 2.2 $1.444\pm 0.003$
$3.354+0.023-0.031$ 2.3 $1.401+0.241-0.138$ $3.941+1.929-1.106$ 2.4
$1.505+0.085-0.064$ $2.670+0.431-0.307$ $6.118+2.644-1.520$ 2.5
$1.392+0.103-0.079$ $3.798+1.833-0.950$

Figure 4: $\rho(\omega)$ for all datasets, subject to $\rho(\omega)~{}\geq~{}0$. The $y$-axis has been normalized to give the maximum an intensity of 1. Figure 5: Maxima of fig. 5 in order of increasing $\omega$$\text{\,}\mathrm{/}\mathrm{G}\mathrm{eV}$, with corresponding left and right Half Width at Half Maximum (HWHM). |  Traditional [1] | Spectral representation
---|---|---
$\beta$ | $am_{\phi}$ | $m_{\phi}[\Lambda]$ |  $m_{\phi}$ /GeV | $am_{\phi}$ | $m_{\phi}[\Lambda]$ |  $m_{\phi}$ /GeV
2.1 | 1.853(13) | 5.21(39) | 1.34(56) | 1.910(52) | 5.364(15) | 1.382(38)
2.2 | 1.517(10) | 5.55(41) | 1.43(34) | 1.532(32) | 5.613(12) | 1.444(30)
2.3 | 1.241(6) | 5.75(54) | 1.48(54) | 1.174(8) | 5.436(38) | 1.401(10)
2.4 | 0.924(8) | 5.93(43) | 1.53(16) | 0.910(8) | 5.832(53) | 1.505(14)
2.5 | 0.696(6) | 6.32(46) | 1.63(23) | 0.595(10) | 5.408(88) | 1.392(23)
Table 2: The $0^{++}$ glueball masses as presented in Table 3 of [1] compared
with the spectral representation method. The physical units were calculated
assuming $\sqrt{\sigma}=0.44$ GeV.

The values of $a\sqrt{\sigma}$ and $a[\Lambda^{-1}]$ are given in table 1. A
string tension $\sqrt{\sigma}$ of $0.44\text{\,}\mathrm{GeV}$ was assumed to
convert to physical units.

Since the toy model study from section 3.1 suggests that a small number of
data points, combined with high quality data, leads to a low variance in the
inversion, we assume that the uncertainties in $a\sqrt{\sigma}$ or
$a[\Lambda^{-1}]$ are the dominant sources of uncertainty in the calculated
masses. The uncertainties in the glueball masses as given in table 2 were
therefore calculated by propagating the uncertainties in $a\sqrt{\sigma}$ or
$a[\Lambda^{-1}]$ respectively.

The scalar glueball masses as obtained by [1], and those extracted using the
spectral density method, agree within one standard deviation. This is actually
quite remarkable, since we based ourselves on the unsmeared data, whilst the
mass esimates from [1] were obtained after smearing, which usually improves
the ground state signal.

In addition to the ground state, the spectra presented in fig. 5 also show a
clear indication of an excited state, though its position shows much greater
variance with $\beta$ than that of the ground state. Figure 5 lists the
$\omega/\text{GeV}$ values of all the local maxima in fig. 5, and the left and
right Half Width at Half Maximum (HWHM) values of each local maximum.

Except for $\beta=2.4$, all conditions indicate a single excited state at
$3.3-3.9$ $\text{\,}\mathrm{GeV}$. The only exception is $\beta=2.4$, which
has a second maximum at $2.670(25)\text{\,}\mathrm{GeV}$ and a third one at
$6.118(57)\text{\,}\mathrm{GeV}$. The relatively wide peaks for the excited
states indicate that we might need a better signal, as there are in principle
no decay channels open for these 1st excited states in pure gluodynamics, not
even to the $0^{++}$ ground state.

We do not have a clear-cut understanding of why the $\beta=2.2$ (very sharp
1st excited state) and the $\beta=2.4$ (2nd excited state and relative
dislocation of the 1st excited state peak) are a bit different, but we dare to
speculate it has either to do with the much larger configuration number for
$\beta=2.2$, or with the $\beta=2.4$ data set encompassing a considerably
larger lattice volume, combined with still a quite large statistical ensemble,
relative to the other $\beta$-values considered. We should also keep in mind
that we are aiming at reliable continuum estimates, which implies $\beta$
should be sufficiently large, so probably values of $\beta<2.4$ are suboptimal
from this viewpoint. The size of the (inverse) lattice spacing also determines
the maximally accessible momentum scale, forcing our spectral function
estimates to vanish in the UV, in contradistinction with the continuum tail
going as $\omega^{4}$, as follows from power counting or explicit computation
[45].

Interestingly, re-expressing the $\beta=2.4$ groundstate and 1st excited state
mass in units of string tension, we get $3.420(32)\sqrt{\sigma}$ and
$6.069(56)\sqrt{\sigma}$ respectively, values that compare quite well with an
independent SU(2) estimate giving $3.781(23)\sqrt{\sigma}$ for the ground
state and $6.126(38)\sqrt{\sigma}$ for the 1st excited state222We thank A.
Athenodorou and M. Teper for sharing these yet unpublished results with us..

##  4 Summary and Conclusion

In this work we discuss the computation of the Källén-Lehmann spectral density
function $\rho$ from two point correlation functions $G(p^{2})$ that can have
a (divergent) polynomial part. The extraction of the finite part of $G(p^{2})$
that is associated with $\rho$ requires the subtraction of that polynomial
part. To avoid doing the subtraction directly on the numerical data, an
operation that can be cumbersome and a source of large uncertainties, we
consider the Schwinger function $C(\tau)$ instead of $G(p^{2})$. Once the
Schwinger function has been calculated from $G(p^{2})$, subtractions are
performed simply by removing the data point $C(0)$. The spectral density is
then obtained by taking the inverse Laplace transform of $C(\tau)$ for
$\tau>0$. Further, the inversion method as described in [33], that does not
constrain $\rho$, is compared with the reformulation of the problem imposing a
positive $\rho$, that relies on the non-negative least square method to
perform the inversion.

The numerical method was first tested on mock data from a given spectral
density. The test followed the same rules as used in the original maximum
entropy method inversion paper [35], with similar results.

The inversion method was then applied to extract the Källén-Lehmann spectral
density function from the lattice SU(2) glueball propagator data with quantum
numbers $J^{PC}=0^{++}$, based on the unsmeared data of [1]. We found that the
inversion method is robust and allows estimation of the ground state and
excited states masses, albeit that the latter must be considered with some
more care. The comparison of the mass estimates between a conventional lattice
approach [1] and the method described herein shows that the latter achieves a
better statistically accurate result for the ground state, despite the lack of
smearing. Further, the spectral method is able to give information in the
$0^{++}$ glueball spectra not appearing from the conventional approach invoked
in [1].

The new spectral method is sufficiently general to be applied to any
spectroscopic calculation. In future research, we will scrutinize the usage of
smeared data to boost the contribution of the lowest lying states, which
should allow us to access at least the 1st excited state with a better
accuracy. Afterwards, the analysis performed in the current manuscript might
also be extended to the analysis of the spectra of other glueball states with
different quantum numbers.

##  5 Acknowledgements

The authors are indebted to N. Yamanaka, H. Iida, A. Nakamura and M. Wakayama
for sharing the Schwinger function estimates based on their data of [1]. We
also thank A. Athenodorou and M. Teper for useful communication.

O.O. was partly supported by the FCT – Fundação para a Ciência e a Tecnologia,
I.P., under project numbers UIDB/04564/2020 and UIDP/04564/2020, while the
work of D.D. and M.R. was supported by KU Leuven IF project C14/16/067.

##  Appendix A Subtracted correlator

A correlator $G(p^{2})$ can be expressed in terms of a Källén–Lehmann spectral
representation as

$G(p^{2})=\int_{0}^{\infty}\frac{2\omega\rho(\omega)\differential{\omega}}{\omega^{2}+p^{2}},$
(14)

if this integral converges. However, if the correlator contains divergences
polynomial in $p^{2}$, it is still possible to find a Källén–Lehmann spectral
representation, if one first subtracts the polynomial part. In order to do
this, consider the Taylor expansion of $G(p)$ around $\bar{p}$:

$\displaystyle
G(p^{2})=\sum_{k=0}^{n-1}(p^{2}-\bar{p}^{2})^{k}\partialderivative[k]{G(p^{2})}{(p^{2})}\Bigg{|}_{p^{2}=\bar{p}^{2}}+\tilde{G}_{n}(p^{2}).$
(15)

where $\tilde{G}_{n}(p^{2})$ is the remainder after the first $n$ terms have
been isolated. At zeroth order, we find

$\displaystyle\tilde{G}_{1}(p^{2})$ $\displaystyle=G(p^{2})-G(\bar{p}^{2})$
$\displaystyle=(-p^{2}+\bar{p}^{2})\int_{0}^{\infty}\frac{2\omega\rho(\omega)\differential{\omega}}{(\omega^{2}+\bar{p}^{2})(\omega^{2}+p^{2})}$
(16)

Through induction, we obtain

$\tilde{G}_{n}(p^{2})=(-p^{2}+\bar{p}^{2})^{n}\int_{0}^{\infty}\frac{2\omega\rho(\omega)\differential{\omega}}{(\omega^{2}+\bar{p}^{2})^{n}(\omega^{2}+p^{2})}$
(17)

This allows us to state that

$\displaystyle G(p^{2})$
$\displaystyle=\sum_{k=0}^{n-1}a_{k}(p^{2}-\bar{p}^{2})^{k}+(-p^{2}+\bar{p}^{2})^{n}\int_{0}^{\infty}\frac{2\omega\tilde{\rho}(\omega)\differential{\omega}}{\omega^{2}+p^{2}}$
$\displaystyle=\sum_{k=0}^{n-1}a_{k}(p^{2}-\bar{p}^{2})^{k}+(-p^{2}+\bar{p}^{2})^{n}\int_{-\infty}^{\infty}\frac{\tilde{\rho}(\omega)\differential{\omega}}{\omega-
ip}$ (18)

with

$\displaystyle a_{n}$
$\displaystyle=\frac{1}{n!}\partialderivative[n]{G(p^{2})}{(p^{2})}\Bigg{|}_{p^{2}=\bar{p}^{2}},$
(19) $\displaystyle\tilde{\rho}(\omega)$
$\displaystyle=\frac{\rho(\omega)}{(\omega^{2}+\bar{p}^{2})^{n}}.$ (20)

Because $\rho$ is an odd (real) function for bosonic degrees of freedom [35],
so is $\tilde{\rho}$. We will now calculate the Fourier transform of the
subtracted correlator, in order to obtain $C(\tau)$.

$\displaystyle\mathcal{F}\left\\{G(p)\right\\}$
$\displaystyle=\underbrace{\sum_{k=0}^{n-1}a_{k}\mathcal{F}\left\\{(p^{2})^{k}\right\\}}_{A}+\underbrace{\mathcal{F}\left\\{(-p^{2})^{n}\int_{-\infty}^{\infty}\frac{\tilde{\rho}(\omega)\differential{\omega}}{\omega-
ip}\right\\}}_{B}$ (21)

where we have set $\bar{p}=0$. For the first term $A$ we obtain

$\displaystyle A=\sum_{k=0}^{n-1}a_{k}\mathcal{F}\left\\{(p^{2})^{k}\right\\}$
$\displaystyle=2\pi\sum_{k=0}^{n-1}a_{k}(-1)^{k}\delta^{(2k)}(\tau)$ (22)

and for the second

$\displaystyle B$
$\displaystyle=\mathcal{F}\left\\{(-p^{2})^{n}\int_{-\infty}^{\infty}\frac{\tilde{\rho}(\omega)\differential{\omega}}{\omega-
ip}\right\\}$ (23)
$\displaystyle=\int_{0}^{\infty}\differential{\omega}\tilde{\rho}(\omega)\mathcal{F}\left\\{\frac{2\omega(-p^{2})^{n}}{\omega^{2}+p^{2}}\right\\}$
(24)
$\displaystyle=2\pi(-1)^{n}\int_{0}^{\infty}\differential{\omega}\tilde{\rho}(\omega)\left[-2\sum_{k=2}^{n}\omega^{2k-3}\delta^{(2(n-k)+2)}(\tau)+\omega^{2n}e^{-|\tau|\omega}\right]$
$\displaystyle=4\pi(-1)^{n+1}\sum_{k=2}^{n}\delta^{(2(n-k)+2)}(\tau)\int_{0}^{\infty}\differential{\omega}\omega^{2k-3}\tilde{\rho}(\omega)+\mathcal{L}\left\\{\tilde{\rho}(\omega)\omega^{2n}\right\\}(|\tau|)$

As a result the total transform is given by

$\displaystyle\mathcal{F}\left\\{G(p)\right\\}$
$\displaystyle=2\pi\sum_{k=0}^{n-1}a_{k}(-1)^{k}\delta^{(2k)}(\tau)+\mathcal{L}\left\\{\rho(\omega)\right\\}(|\tau|)$
(25)
$\displaystyle\quad+4\pi(-1)^{n+1}\sum_{k=2}^{n}\delta^{(2(n-k)+2)}(\tau)\int_{0}^{\infty}\differential{\omega}\omega^{2k-3}\tilde{\rho}(\omega)$

An important consequence of this is that all the subtractions are turned into
(derivatives of) delta functions, which end up at $\tau=0$. Therefore, taking
only $\tau>0$, we can say that

$\displaystyle C(\tau)=\mathcal{L}\left\\{\rho(\omega)\right\\}(\tau).$

##  Appendix B Numerical dependence on $N_{\omega}$

To produce the glueball spectra given in section 3.2, $\omega$ was sampled
evenly in logarithmic space between $[10^{-5},10^{5}]$
$\text{\,}\mathrm{GeV}$, in $N_{\omega}=10000$ steps, when constructing the
Laplace kernel $\boldsymbol{K}$. However, this number of steps $N_{\omega}$ is
unnecessarily large, and was chosen merely to generate aesthetically pleasing
graphs.

To demonstrate the convergence of the maxima of the spectral functions, fig. 6
shows the $\omega$ coordinate of the first and second maximum in the
$\beta=2.1$ dataset for increasing $N_{\omega}$. In both the linearly spaced
and logarithmically evenly spaced scenarios, the ground state converges
rapidly, and the second maximum soon follows suit. In addition, fig. 7 shows
the full spectrum in both scenarios, normalized such that the first maximum
has $\rho(\omega)=1$. There is excellent agreement between the two methods for
the determination of the ground state mass, but the convergence of the second
maximum is better when using logarithmically evenly spaced samples.

(a) Results for logarithmic sample of $\omega$ in $[10^{-5},10^{5}]$
$\text{\,}\mathrm{GeV}$.

(b) Results for linear sample of $\omega$ in $[0,100]$
$\text{\,}\mathrm{GeV}$.

Figure 6: The position of the first and second maxima for the $\beta=2.1$ data
set, as a function of the number of samples $N_{\omega}$ in the Laplace kernel
$\boldsymbol{K}$.

(a) $\omega$ sampled uniformly in logarithmic space between $[10^{-5},10^{5}]$
$\text{\,}\mathrm{GeV}$.

(b) $\omega$ sampled uniformly in linear space between $[0,100]$
$\text{\,}\mathrm{GeV}$.

Figure 7: $\rho(\omega)$ for the $\beta=2.1$ data set, as a function of the
number of samples $N_{\omega}$ in the Laplace kernel $\boldsymbol{K}$.
$\rho(\omega)$ was normalized such that the first maximum has
$\rho(\omega)=1$.

## References

  * [1] N. Yamanaka, H. Iida, A. Nakamura, and M. Wakayama.  Glueball scattering cross section in lattice SU(2) Yang-Mills theory.  Phys. Rev. D, 102(5):054507, 2020.
  * [2] P. A. Zyla et al.  Review of Particle Physics.  PTEP, 2020(8):083C01, 2020.
  * [3] M. Teper.  Glueball masses and other physical properties of SU(N) gauge theories in D = (3+1): A Review of lattice results for theorists.  12 1998.
  * [4] C. J. Morningstar and M. J. Peardon.  The Glueball spectrum from an anisotropic lattice study.  Phys. Rev. D, 60:034509, 1999.
  * [5] Y. Chen et al.  Glueball spectrum and matrix elements on anisotropic lattices.  Phys. Rev. D, 73:014516, 2006.
  * [6] S. Narison.  Masses, decays and mixings of gluonia in QCD.  Nucl. Phys. B, 509:312–356, 1998.
  * [7] R. C. Brower, S. D. Mathur, and C.-I. Tan.  Glueball spectrum for QCD from AdS supergravity duality.  Nucl. Phys. B, 587:249–276, 2000.
  * [8] N. Ishii, H. Suganuma, and H. Matsufuru.  Glueball properties at finite temperature in SU(3) anisotropic lattice QCD.  Phys. Rev. D, 66:094506, 2002.
  * [9] A. P. Szczepaniak and E. S. Swanson.  The Low lying glueball spectrum.  Phys. Lett. B, 577:61–66, 2003.
  * [10] D. Dudal, M. S. Guimaraes, and S. P. Sorella.  Glueball masses from an infrared moment problem and nonperturbative Landau gauge.  Phys. Rev. Lett., 106:062003, 2011.
  * [11] H. Sanchis-Alepuz, C. S. Fischer, C. Kellermann, and L. von Smekal.  Glueballs from the Bethe-Salpeter equation.  Phys. Rev. D, 92:034001, 2015.
  * [12] M. Q. Huber, C. S. Fischer, and H. Sanchis-Alepuz.  Spectrum of scalar and pseudoscalar glueballs from functional methods.  Eur. Phys. J. C, 80(11):1077, 2020.
  * [13] V. Mathieu, N. Kochelev, and V. Vento.  The Physics of Glueballs.  Int. J. Mod. Phys. E, 18:1–49, 2009.
  * [14] V. Crede and C. A. Meyer.  The Experimental Status of Glueballs.  Prog. Part. Nucl. Phys., 63:74–116, 2009.
  * [15] W. Ochs.  The Status of Glueballs.  J. Phys. G, 40:043001, 2013.
  * [16] M. P. Hertzberg, F. Rompineve, and J. Yang.  Decay of Boson Stars with Application to Glueballs and Other Real Scalars.  Phys. Rev. D, 103(2):023536, 2021.
  * [17] K. Ishikawa, M. Teper, and G. Schierholz.  The Glueball Mass Spectrum in QCD: First Results of a Lattice Monte Carlo Calculation.  Phys. Lett. B, 110:399–405, 1982.
  * [18] I. Montvay and G. Munster.  Quantum fields on a lattice.  Cambridge Monographs on Mathematical Physics. Cambridge University Press, 3 1997.
  * [19] M. E. Peskin and D. V. Schroeder.  An Introduction to quantum field theory.  Addison-Wesley, Reading, USA, 1995.
  * [20] B. Blossier, M. Della Morte, G. von Hippel, T. Mendes, and R. Sommer.  On the generalized eigenvalue method for energies and matrix elements in lattice field theory.  JHEP, 04:094, 2009.
  * [21] A. Athenodorou and M. Teper.  The glueball spectrum of SU(3) gauge theory in 3 + 1 dimensions.  JHEP, 11:172, 2020.
  * [22] J. W. Negele and H. Orland.  Quantum Many Particle Systems.  Redwood City, USA: Addison-Wesley. (Frontiers in Physics, 68), 1988.
  * [23] H. B. Meyer.  Transport Properties of the Quark-Gluon Plasma: A Lattice QCD Perspective.  Eur. Phys. J., A47:86, 2011.
  * [24] D. Dudal and M. S. Guimaraes.  On the computation of the spectral density of two-point functions: complex masses, cut rules and beyond.  Phys. Rev. D, 83:045013, 2011.
  * [25] D. Dudal, O. Oliveira, and P. J. Silva.  Källén-Lehmann spectroscopy for (un)physical degrees of freedom.  Phys. Rev. D, 89(1):014010, 2014.
  * [26] E. L. Solis, C. S. R. Costa, V. V. Luiz, and G. Krein.  Quark propagator in Minkowski space.  Few Body Syst., 60(3):49, 2019.
  * [27] D. Dudal, D. M. van Egmond, M. S. Guimarães, O. Holanda, B. W. Mintz, L. F. Palhares, G. Peruzzo, and S. P. Sorella.  Some remarks on the spectral functions of the Abelian Higgs Model.  Phys. Rev. D, 100(6):065009, 2019.
  * [28] D. Binosi and R.-A. Tripolt.  Spectral functions of confined particles.  Phys. Lett. B, 801:135171, 2020.
  * [29] F. Siringo.  Analytic structure of QCD propagators in Minkowski space.  Phys. Rev. D, 94(11):114036, 2016.
  * [30] A. K. Cyrol, J. M. Pawlowski, A. Rothkopf, and N. Wink.  Reconstructing the gluon.  SciPost Phys., 5:065, 2018.
  * [31] S. Weinberg.  The Quantum Theory of Fields, volume 1.  Cambridge University Press, 1995.
  * [32] O. Oliveira, D. Dudal, and P. J. Silva.  Glueball spectral densities from the lattice.  PoS, LATTICE2012:214, 2012.
  * [33] D. Dudal, O. Oliveira, M. Roelfs, and P. Silva.  Spectral representation of lattice gluon and ghost propagators at zero temperature.  Nucl. Phys. B, 952:114912, 2020.
  * [34] C. L. Lawson, Society for Industrial Hanson, R. J., and Applied Mathematics.  Solving least squares problems.  Philadelphia : SIAM, [rev. ed.] edition, 1995.  ”This SIAM edition is an unabridged, revised republication of the work first published by Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1974”–T.p. verso.
  * [35] M. Asakawa, T. Hatsuda, and Y. Nakahara.  Maximum entropy analysis of the spectral functions in lattice QCD.  Prog. Part. Nucl. Phys., 46:459–508, 2001.
  * [36] G. Aarts, C. Allton, M. B. Oktay, M. Peardon, and J.-I. Skullerud.  Charmonium at high temperature in two-flavor QCD.  Phys. Rev. D, 76:094513, 2007.
  * [37] A. Jakovac, P. Petreczky, K. Petrov, and A. Velytsky.  Quarkonium correlators and spectral functions at zero and finite temperature.  Phys. Rev. D, 75:014506, 2007.
  * [38] A. Rothkopf.  Bayesian inference of nonpositive spectral functions in quantum field theory.  Phys. Rev. D, 95(5):056016, 2017.
  * [39] H.-T. Ding, O. Kaczmarek, S. Mukherjee, H. Ohno, and H. T Shu.  Stochastic reconstructions of spectral functions: Application to lattice QCD.  Phys. Rev. D, 97(9):094503, 2018.
  * [40] R.-A. Tripolt, P. Gubler, M. Ulybyshev, and L. Von Smekal.  Numerical analytic continuation of Euclidean data.  Comput. Phys. Commun., 237:129–142, 2019.
  * [41] S. Schlichting, D. Smith, and L. von Smekal.  Spectral functions and critical dynamics of the $O(4)$ model from classical-statistical lattice simulations.  Nucl. Phys. B, 950:114868, 2020.
  * [42] R. Fournier, L. Wang, O. V. Yazyev, and Q. Wu.  Artificial neural network approach to the analytic continuation problem.  Phys. Rev. Lett., 124:056401, 2020.
  * [43] A. Kirsch.  An Introduction to the Mathematical Theory of Inverse Problems.  Springer-Verlag, Berlin, Heidelberg, 1996.
  * [44] N. Yamanaka, H. Iida, A. Nakamura, and M. Wakayama.  Dark matter scattering cross section and dynamics in dark Yang-Mills theory.  Phys. Lett. B, 813:136056, 2021.
  * [45] A. L. Kataev, N. V. Krasnikov, and A. A. Pivovarov.  Two Loop Calculations for the Propagators of Gluonic Currents.  Nucl. Phys. B, 198:508–518, 1982.  [Erratum: Nucl.Phys.B 490, 505–507 (1997)].


# Faces: AI Blitz XIII Solutions

Andrew Melnik Eren Akbulut Jannik Sheikh Kira Loos Michael Büttner Tobias
Lenze

(authors are listed alphabetically
Bielefeld University, Germany)

###### Abstract

AI Blitz XIII Faces challenge hosted on www.aicrowd.com platform consisted of
five problems:
Sentiment Classification, Age Prediction, Mask Prediction, Face Recognition,
and Face De-Blurring.
Our team GLaDOS took second place. Here we present our solutions and results.
Code implementation: https://github.com/ndrwmlnk/ai-blitz-xiii

##  1 Sentiment Classification: Classify Facial Expressions

###  1.1 Problem Statement and Dataset

The goal of the Sentiment Classification problem [1] of the AI Blitz XIII
Faces challenge [2] was to train a model that can classify faces into three
sentiment categories: negative, neutral and positive. The dataset doesn’t
contain real images but 512-dimensional embedding vectors generated from
images of faces by model such as ResNet18 or VGG16. The embedding vectors were
extracted from a specific layer of such a model. The original images or the
exact model from which the embedding vectors were derived were not specified.

Each face embedding consists of a 512-dimensional vector (see Table 1). The
values of the embedding vectors are between 0 and 11.24, while the mean value
of the embedding vectors is 0.73. The dataset is divided into three parts: a
training set consisting of 5000 samples, a validation set consisting of 2000
samples, and a test set consisting of 3001 samples. The predictions on the
test set were evaluated. The metrics used for evaluation are F1 score [3]
(average=weighted) as a primary score and accuracy score [4] as a secondary
score.

Table 1: Example of 5 samples from the training dataset. Each embedding has
512 dimensions. ⬇ id embedding (512-dim) label 0 [0.32, 0.98, 1.04, …]
[positive]  1 [0.05, 1.07, 0.60, …] [negative]  2 [0.41, 0.45, 1.39, …]
[negative]  3 [0.43, 0.19, 0.83, …] [positive]  4 [0.63, 0.83, 0.39, …]
[neutral]
---

###  1.2 Methods

See more details on our implementation in the Sentiment Classification
notebook [5]. We tested the following approaches: Random Forest Classifier
(RFC), K-Nearest Neighbours, Gaussian Naive Bayes, LightGBM, Neural Networks
and Support Vector Machines (SVM). For LightGBM and SVM we used Randomized
Grid Search with cross-validation. An error-based method for automatic
assignment of data subsets to the ensemble NNs using the loss profile of the
ensemble is another method for potential improvement of the prediction results
[6]. Out of all the models, Support Vector Classifier (SVC) gave the best
results. It has the following parameters: polynomial kernel with a degree of
3, a C-value of 1, and a scaled gamma with a value of
$1/(n_{features}*X.var())$. The multi-class classification is performed using
a one-vs-rest approach. Normalization with l1 norm for predictions was the
preprocessing method that improved results, albeit by a small amount (see
Table 2).

Table 2: F1 score and accuracy of the baseline Random Forest Classifier and the Support Vector Classifier without preprocessing and with normalization with l1 norm. Trained on the training set and tested on the validation set. Dataset | F1-Score | Accuracy
---|---|---
RFC | SVC | RFC | SVC
without preprocessing | 0.668 | 0.808 | 0.673 | 0.800
with l1-norm normalization | 0.679 | 0.807 | 0.687 | 0.808

###  1.3 Results

For the predictions, we used the SVC with l1 normalization. We used tenfold
cross-validation by combining the training set and the validation set,
shuffling them and normalizing them. In this way, we obtained ten estimators.
Each of the estimators can make predictions on the test set. We calculated the
mean for each sample and rounded it to the nearest label. Omitting certain
estimators has an effect on the predictions and consequently on the score. We
therefore tried different combinations of the estimators until we found a
combination that gave the best results. The final F1 and Accuracy Score on the
test set was 0.806 and 0.806 respectivly.

##  2 Age Prediction: Predict Age From A Picture

###  2.1 Problem Statement and Dataset

The goal of the Age Prediction problem [7] of the AI Blitz XIII Faces
challenge [2] was to predict the age of a human face from an input image. The
labels are categorical with 10 classes (0-10 years old, 10-20, 20-30, 30-40,
40-50, 50-60, 60-70, 70-80, 80-90, 90-100). See Figure 1 for examples. The
data is divided into training, validation, and testing sets. The train set
contains 4000, the validation set 2000 and the test set 3000 images. The
labels were given in training and validation parts of the dataset. The data
distribution is shown on Figure 2.

Images of human faces in the dataset where generated by a StyleGAN3 model [8],
which was trained using the Flickr-Faces-HQ dataset. All images are colour
images (RGB) of size 512x512x3. Since the images were generated artificially,
some images have unnatural artifacts of human faces. Figures 3(e) and 3(f)
give you an example.

Figure 1: Images from train set for each label, ascending. Upper right
starting with label 0-10. Lower right starting with label 50-60. Figure 2:
Distribution of train data, after merging

(a) age label 80-90, bounding box

(b) age label 80-90, cropped, bounding box

(c) age label 30-40, bounding box

(d) age label 30-40, cropped, bounding box

(e) age label 40-50, cartoonish quality

(f) age label 90-100, low quality

Figure 3: Example pictures from the dataset with provided age labels.

###  2.2 Evaluation Criteria

For the evaluation of models on the leader board [2], weighted average F1
score [3] was used as the primary score and the accuracy score [4] as the
secondary score. The best score for the F1 score is 1 and the worst score is
0.

$F1=2*\frac{precision*recall}{precision+recall}.$

###  2.3 Initial experiments

In this section, we present some of our initial experiments before discussing
our best model in 2.4. We started with Random Forest Classifier with
RandomizedSearchCV for hyperparamter tuning as well as Light Gradient Boosting
Classifier without proper hyperparamter tuning for a subset of the train and
test data. The results were not that promising. We achieved a F1 score of 0.33
on the test set with the RFC and even less with the LGBMClassifier.

###  2.4 Methods

See more details on our implementation in the Age Prediction notebook [5]. We
have achieved good results in age prediction using the You Only Look Once
(YOLO) architecture [9, 10, 11]. Although YOLO is an object detection
algorithm, we mainly utilized its class prediction capabilities and abandoned
the bounding box prediction functionality. The code implementation of YOLOv5
can be downloaded from GitHub [12]. Since we had no information regarding the
bounding boxes, we create fixed bounding boxes 410x410 around the center of
the images (x, y = 0.5; width, height = 0.8). Figures 3(a) and 3(c) give you a
visualization. We trained the model to predict the age class. We also merged
the training and validation sets to get more images for training.
Additionally, for every image we created a cropped fixed size version,
doubling our training set. We took central 280x280 patch of 512x512 images and
resized these 280x280 patches back to 512x512. Examples can be found on Figure
3(b) and 3(d). We trained YOLOv5 with network size L for 200 epochs and used a
batch size of 16. Network size L has 46.5M parameters. The different network
sizes and characteristics can be found on Github [12]. The training time was
about 9 minutes per epoch using Google Colab Pro with a TESLA GPU with 16GB
RAM.

###  2.5 Results

Our result for the Age Prediction problem with the mentioned method and
hyperparamters was a F1 score of 0.870 and an accuracy score of 0.871.

##  3 Mask Prediction: Predict Mask & Bounding Box From Images

###  3.1 Problem Statement and Dataset

Figure 4: Input: a picture of a person with a mask. Output: the bounding box
coordinates and the mask type. Figure 5: Samples of the dataset for each mask
type. L.t.r.: N95, cloth, surgical, KN95.

The goal of the Mask Prediction problem [13] of the AI Blitz XIII Faces
challenge [2] was to predict the type and bounding box of the mask that a
person is wearing in an image (fig.4). Trained on the provided dataset, we
achieved an accuracy of 99.3%, which was the top score in this problem.

The provided dataset was divided into three different sets. Train, validation
and test sets. Each set contains 5000, 2000, and 3000 512x512pixel images,
respectively. Each image contains one of 4 mask types: Surgical, N95, KN95 and
Cloth (fig.5). CSV files were also provided for the training and validation
images, which contained information about the mask type for each image and the
corresponding bounding box of the mask in pixel coordinates (fig.6).

Figure 6: Data structure of CSV files: ImageID, bounding box and mask type.
Also the format used for the bounding box on the right.

###  3.2 Methods

See more details on our implementation in the Mask Prediction notebook [5].
You Only Look Once (YOLO) [9] was chosen as the prediction method. YOLOv5 [12]
is the fifth iteration of YOLO, a state-of-the-art, real-time object detection
system capable of creating bounding boxes around detected objects in images.
YOLO requires a specific format and file hierarchy. For this the following
folder structure has to be implemented, containing …/images/train,
…/images/val, …/images/test images, and …/labels/train, …/labels/val labels.

Labels in the YOLO format must be created and moved to the respective labels
folders. These labels are .txt files for each image, with the same name as the
corresponding images. For example, the image …/images/train/abc.jpg must have
a counterpart …/labels/train/abc.txt. Each label contains the object class (a
numeric value starting with 0), and the bounding box in the following format:

<object-class> <x> <y> <width> <height>

It should also be noted that YOLO uses values that are image size related and
not pixel related. To convert the values from the CSV files to YOLO format,
the following formulas were used:

Training in YOLOv5 can be started by a command with parameters that specify
the image size in pixel, the batch size, number of epochs, the size of YOLOv5
network on which to train, and the path to the mask_prediction.yml [5]. The
mask_prediction.yml file contains the paths to the training and validation
images as well as the number and names of classes [5]. The training time per
epoch for the training set of 3000 images and a batch size of 16 is about 15
minutes on a Tesla K80 via Google Colab.

To get the bounding box coordinates, the –save-txt parameter is required to
generate labels in YOLO format. These bounding box coordinates have to be
converted into the format of the CSV files to be submitted for scoring. The
scoring is based on the Average Precision of mask type and bounding box ( @
IoU=0.50:0.50 )

###  3.3 Results

The prediction with the highest confidence score for each image was selected
as the predicted object. The best result, with 99.3% accuracy, was achieved
after training for 100 epochs, with a batch size of 16 images and with the
YOLOv5l model.

##  4 Face Recognition: Find The Correct Face In Crowd

###  4.1 Problem Statement and Dataset

The goal of the Face Recognition problem [14] of the AI Blitz XIII Faces
challenge [2] was to find the missing face in crowd (target image). The
missing image contains a single face that has to be found in the corresponding
target image. Each target image shows 100 faces. The unlabeled dataset [14]
consists of 1000 missing and target image pairs. The faces were all generated
using a StyleGAN3 model. The person in the missing image differ in some
attributes in the corresponding target image. See examples in Figure 7.

Figure 7: Missing image (left), target image with 100 faces (middle) and
prediction (right) for a set of images from the dataset.

###  4.2 Methods

See more details on our implementation in the Face Recognition notebook [5].
Since the provided dataset is unlabeled, we could not train the face detection
or face recognition ourselves. We used pre-trained models and methods for
face-recognition [15][16].

In our first approach, we used face-alignment [17] to compute 68 3D landmarks.
We computed the distances between each possible pair of landmarks for each
given face. Afterwards we compared these distances between the missing face
and every face seen in the target image, selecting the one with the smallest
error as the prediction. The distances are normalized by the jaw width, which
we defined as the distance between landmark 1 and 17. This got us an accuracy
of 22%, which could not be improved by using only a subset of distances to
compare, nor by weighting some distances higher.

In our second approach, we used Google’s MediaPipe [18][19] to compute 468 3D
landmarks in a given face. Initially, we only used the landmark groups [20]
which were provided by Google’s MediaPipe, but by using our own groups in
addition to that, we could increase the accuracy to 50%. Our own groups
provided landmark paths for the nose and the facial structure.

In our third approach we decided to use a face-recognition method [15]. It
allowed us to use dlib [16] library that provides an encoding function to find
an embedding for a face, and a distance function to calculate the difference
between two encodings. We used two different encoding models to find the
encodings: small uses 5 landmarks and large uses 68 landmarks. To predict the
target images, we used the mean over the distances calculated by those two
models. This increased our initial accuracy for this approach by 2%, resulting
in our final accuracy: 78.8%.

###  4.3 Results

For the third approach, we got an accuracy of 78.8%. Selection of a random
image would give us an accuracy of 1%, since there are 100 possible target
images. The prediction time for our notebook was 2.5 hours using Google Colab
Pro computation resources.

##  5 Face De-Blurring: Convert Blurred image Into Clear Image

###  5.1 Problem Statement and Dataset

The goal of the Face Deblurring problem [21] of the AI Blitz XIII Faces
challenge [2] was to generate a sharp image of a human face from a blurred
version of the same human face. The dataset used for this challenge are images
of human faces, which were generated by a StyleGAN3 model [8], which was
trained on the Flickr-Faces-HQ dataset. The data is split into train,
validation & test sets. The train set contains 5000 pairs of blurred and
original images, the validation set 2000 pairs of blurred and original images,
and the test set 3000 blurred images. All images are colour images (RGB) of
size 512x512x3. See example image pairs in Figure 8.

Figure 8: Examples of training pairs: Original (Top) and Blurred (Bottom).
Figure 9: Examples of our results: test blurred image (Bottom) and deblurred
image (Top) Figure 10: SRCNN Architecture Figure 11: MSE vs. SSIM [22]

####  5.1.1 Evaluation Criteria

The main evaluation criteria was Structural Similarity (SSIM) and for 2 images
"x" and "y" how it’s being calculated is shown below.

$\mathit{SSIM}(x,y)=\frac{(2\mu_{x}\mu_{y}+C_{1})+(2\sigma_{xy}+C_{2})}{(\mu_{x}^{2}+\mu_{y}^{2}+C_{1})(\sigma_{x}^{2}+\sigma_{y}^{2}+C_{2})}$
(1)

To solve potential conflicts in case of having an identical SSIM score a
second metric was utilized, Peak signal-to-noise ratio (PSNR).

$\mathit{PSNR}(x,y)=\frac{10\log_{10}[\max(\max(x),\max(y))]^{2}}{\lvert
x-y\rvert^{2}}$ (2)

For such image reconstruction problems SSIM is able to give more meaningful
results than Mean Squared Error (MSE) as you can see in 11.

###  5.2 Methods

####  5.2.1 SRCNN

See more details on our implementation in the Face De-Blurring notebook [5].
We mainly focused on the SRCNN [23] approach throughout this challenge with a
modified loss function. Other approaches [24] [25] [26] [27] required large
computing resources, which we did not have. SRCNN is an architecture is end-
to-end learning approach that maps low resolution images into high resolution
images. This approach is based on a CNN layers and therefore can be adopted to
all input sizes of images.

The original SRCNN paper proposes to used a MSE loss only. We combined MSE and
L1 losses with a 3rd loss function from SSIM family, Multi-scale structural
similarity (MSSSIM) [28]. The original loss function was

$\mathit{Loss}(x,y)=MSE(x,y)$ (3)

and modified into

$\mathit{Loss}(x,y)=MSE(x,y)+L1(x,y)+MSSIM(x,y)$ (4)

Too avoid MSSIM loss to return NaN from time to time, at each scale of the
MSSIM calculation we normalized the output using a Rectified Linear Unit
(ReLu) normalizer.

The network only contains 3 CNN layers (see Figure 10), the first CNN layer
with the kernel size 9 takes a lower resolution image and extracts a set of
feature maps, then the second CNN layer with the kernel size 1 maps there
features non-linearly to a higher resolution patch representations and the
final CNN layer with the kernel size 5 combines the predictions within a
spatial neighbourhood to produce the final high-resolution image. See our
results in Figure 9.

The average training time on a laptop was around 8-12 hours for 120-180 epochs
with a batch size of 1, 2 and 4. AdamW with a learning rate of 0.001 used for
the final submission. To optimize the training process we used a scheduler
with the patience of 4 epochs and factor 0.5.

###  5.3 Results

In this problem we got SSIM of 0.748 and an PSNR of 26.627. See examples of
our results in Figure 9.

## References

  * [1] AIcrowd, “AI Blitz XIII Faces: Sentiment Classification - Classify Facial Expressions.” https://www.aicrowd.com/challenges/ai-blitz-xiii/problems/sentiment-classification; accessed April 4, 2022.
  * [2] AIcrowd, “AI Blitz XIII Faces.” https://www.aicrowd.com/challenges/ai-blitz-xiii; accessed April 4, 2022\.
  * [3] “sklearn.metrics.f1score.” https://scikit-learn.org/stable/modules/generated/sklearn.metrics.f1_score.html.  accessed April 4, 2022.
  * [4] “sklearn.metrics.accuracy_score.” https://scikit-learn.org/stable/modules/generated/sklearn.metrics.accuracy_score.html.  accessed April 4, 2022.
  * [5] A. Melnik, E. Akbulut, J. Sheikh, K. Loos, M. Büttner, and T. Lenze, “Faces: AI Blitz XIII Solutions.” https://github.com/ndrwmlnk/ai-blitz-xiii; accessed April 4, 2022.
  * [6] N. Bach, A. Melnik, F. Rosetto, and H. Ritter, “An error-based addressing architecture for dynamic model learning,” in International Conference on Machine Learning, Optimization, and Data Science, pp. 617–630, Springer, 2020\.
  * [7] AIcrowd, “AI Blitz XIII Face: Age Prediction - Predict Age From A Picture.” https://www.aicrowd.com/challenges/ai-blitz-xiii/problems/age-prediction; accessed April 4, 2022.
  * [8] T. Karras, M. Aittala, S. Laine, E. Härkönen, J. Hellsten, J. Lehtinen, and T. Aila, “Alias-free generative adversarial networks,” in Proc. NeurIPS, 2021.
  * [9] C. Limberg, A. Melnik, A. Harter, and H. J. Ritter, “YOLO - You only look 10647 times,” arXiv preprint arXiv:2201.06159, 2022.
  * [10] P. Jędrzej Świeżewski, “Introduction to yolo algorithm and yolo object detection.” https://appsilon.com/object-detection-yolo-algorithm/.  accessed February 27, 2022.
  * [11] J. Redmon, S. Divvala, R. Girshick, and A. Farhadi, “You only look once: Unified, real-time object detection,” 2016.
  * [12] “YOLOv5.” https://github.com/ultralytics/yolov5/.
  * [13] AIcrowd, “AI Blitz XIII Face: Mask Prediction - Predict Mask & Bounding Box From Images.” https://www.aicrowd.com/challenges/ai-blitz-xiii/problems/mask-prediction; accessed April 4, 2022.
  * [14] AIcrowd, “AI Blitz XIII Face: Face Recognition - Find The Correct Face In Crowd.” https://www.aicrowd.com/challenges/ai-blitz-xiii/problems/face-recognition; accessed April 4, 2022.
  * [15] GitHub, “face-recognition.” https://github.com/ageitgey/face_recognition; accessed April 4, 2022.
  * [16] GitHub, “dlib.” https://github.com/davisking/dlib; accessed April 4, 2022\.
  * [17] A. Bulat, I. Toubal, J.-Y. Lee, S. Bogdanov, J. J. Lewis, S. Jadhav, S. Rezaei, P. Grunt, and A. Solanki, “face-alignment.” https://github.com/1adrianb/face-alignment, 2021.
  * [18] C. Lugaresi, J. Tang, H. Nash, C. McClanahan, E. Uboweja, M. Hays, F. Zhang, C.-L. Chang, M. G. Yong, J. Lee, W.-T. Chang, W. Hua, M. Georg, and M. Grundmann, “Mediapipe: A framework for building perception pipelines,” tech. rep., Google Research, 2019.
  * [19] GitHub, “Mediapipe.” https://github.com/google/mediapipe; accessed April 4, 2022.
  * [20] GitHub, “Mediapipe landmark groups.” https://github.com/google/mediapipe/blob/master/mediapipe/python/solutions/face_mesh_connections.py; accessed April 4, 2022.
  * [21] AIcrowd, “AI Blitz XIII Face: Face De-Blurring - Convert Blurred image Into Clear Image.” https://www.aicrowd.com/challenges/ai-blitz-xiii/problems/face-de-blurring; accessed April 4, 2022.
  * [22] H. Yeganeh, Cross Dynamic Range And Cross Resolution Objective Image Quality Assessment With Applications.  PhD thesis, University of Waterloo, 07 2014.
  * [23] C. Dong, C. C. Loy, K. He, and X. Tang, “Image super-resolution using deep convolutional networks,” IEEE transactions on pattern analysis and machine intelligence, vol. 38, no. 2, pp. 295–307, 2015.
  * [24] A. Melnik, “Deep Face Generation and Editing: A Survey.”  https://github.com/ndrwmlnk/deep-face-generation-and-editing-a-survey; accessed April 4, 2022.
  * [25] X. Mao, Y. Liu, W. Shen, Q. Li, and Y. Wang, “Deep residual fourier transformation for single image deblurring,” arXiv preprint arXiv:2111.11745, 2021.
  * [26] Z. Wang, X. Cun, J. Bao, and J. Liu, “Uformer: A general u-shaped transformer for image restoration,” arXiv preprint arXiv:2106.03106, 2021.
  * [27] S. W. Zamir, A. Arora, S. Khan, M. Hayat, F. S. Khan, and M.-H. Yang, “Restormer: Efficient transformer for high-resolution image restoration,” arXiv preprint arXiv:2111.09881, 2021.
  * [28] Z. Wang, E. Simoncelli, and A. Bovik, “Multiscale structural similarity for image quality assessment,” in The Thrity-Seventh Asilomar Conference on Signals, Systems Computers, 2003, vol. 2, pp. 1398–1402 Vol.2, 2003.


# AI-based Wildfire Prevention, Detection and Suppression System

Prisha Shroff Prisha Shroff is with Hamilton High School, Chandler, AZ.email:
<EMAIL_ADDRESS>
###### Abstract

Wildfires pose a serious threat to the world’s environment. The global
wildfire season length has increased by 19% and severe wildfires have besieged
nations around the world. Every year, forests are burned by these wildfires,
causing vast amounts of carbon dioxide to be released into the atmosphere,
contributing to climate change. Therefore, there is a need for a system which
prevents, detects, and suppresses wildfires. The AI-based Wildfire Prevention,
Detection and Suppression System (WPDSS) is a novel, fully-automated, end-to-
end, AI-based solution to effectively predict hotspots (areas where wildfires
can start) and detect wildfires, deploy drones to spray fire retardant,
preventing and suppressing wildfires. WPDSS consists of four steps. 1) Pre-
processing: WPDSS loads real-time satellite data from NASA and meteorological
data from NOAA of vegetation, temperature, precipitation, wind, soil moisture,
and land cover for prevention. For detection, it loads the real-time data of
Land Cover, Humidity, Temperature, Vegetation, Burned Area Index, Ozone, and
CO2. It uses the process of masking to eliminate not-hotspots and not-
wildfires such as water bodies, and rainfall. 2) Learning: The AI model
consists of a random forest classifier, a type of supervised machine learning,
which uses a series of decision trees to make an accurate decision. It is
trained using a labeled dataset of hotspots/wildfires and not-hotspots/not-
wildfires. 3) Identification of hotspots and wildfires: WPDSS runs the real-
time data through the model to automatically identify hotspots and wildfires.
4) Drone deployment: The drone flies to the identified hotspot/wildfire
location. WPDSS attained a 98.6% accuracy in identifying hotspots and a 98.7%
accuracy in detecting wildfires. WPDSS will reduce the impacts of climate
change, protect ecosystems and biodiversity, avert huge economic losses, and
save human lives. The power of WPDSS developed can be applied to any location
globally to prevent and suppress wildfires, reducing climate change.

###### Index Terms:

machine learning, wildfires, prevention, detection, suppression

##  I Introduction

###  I-A Problem

Wildfires pose a serious threat to the world’s environment. Forests comprise
more than 31% of the world’s land surface and contribute to the continuity of
ecological balance, and they store large amounts of carbon. When forests burn,
they release carbon dioxide into the atmosphere, contributing to climate
change.

In the last few years, severe wildfires have besieged nations around the
world, including the United States, Canada, Australia, Greece, Spain,
Portugal, Chile, Russia, China, Siberia, South Korea, Israel and Brazil. For
example, the 2019–2020 wildfire season in Australia caused the burning of over
46 million acres of forest while destroying over 10,000 structures, resulting
in over $100 billion in damages.

According to the World Wildlife Fun, the length of the global fire season has
increased on average by 19% . As stated by the National Interagency Fire
Center, there were around 58,250 wildfires and more than 10.3 million acres
burned in the US in 2020. The costs for 2020 wildfires in the Western US alone
totaled between $130-$150 billion according to AccuWeather. The lasting
effects of wildfires include ecosystem and biodiversity loss, forest and soil
degradation, air pollution, economic losses, destruction of watersheds, and
negative impacts to human health. [1].

As the world continues to get warmer and drier, experts believe that the
intensity of wildfires will increase in areas like Australia, California, and
the Mediterranean. Recently, fires above the Arctic Circle have burned vast
swaths of typically frozen, carbon-rich peat soils, releasing a record 244
megatons of carbon into the atmosphere. This is more than Spain releases from
burning fossil fuels in an entire year, according to Copernicus Atmosphere
Monitoring Service.

There is a need for a system that can effectively predict hotspots and detect
wildfires, to prevent and suppress these wildfires.

###  I-B Current Solutions

The current solutions in the market do not have wildfire prevention and
detection capabilities. A couple of the solutions are described as follows:
1\. Fire Urgency Estimator in Geosynchronous Orbit. This system detects actual
fires using satellite images, and uses drones to track the fire’s progress.
The solution doesn’t prevent the wildfire from occurring. 2\. Wildfire
Detection System InsightFD + Insight Globe. This system consists of a network
of robots that are constantly monitoring areas with thermal imaging and
cameras. Once it detects a wildfire, it alerts and notifies the fire
department. The solution doesn’t prevent wildfires from occurring and requires
robots to be installed at each site. 3\. BurnMonitor: An Early Wildfire
Detection IoT Solution. This system consists of a network of sensors that can
detect if a fire is happening nearby. It sends a signal to a gateway which
transmits the location to the fire department. The solution doesn’t prevent
wildfires from starting but instead detects the fire. 4\. San Diego Gas and
Electric iPredict. This system uses an AI based predictive model to increase
the accuracy of weather forecasts and predicts powerline caused wildfires. It
only prevents these powerline caused wildfires, and does not detect wildfires.
5\. Smart Wildfire Sensor Device. This uses sensors that are installed on the
site to predict wildfires. It analyzes the amount of dead fuel, the moisture
of the fuel and other datasets to predict the wildfires, but the biggest
problem with this solution is that it has to be manually installed in every
location, which increases the cost exponentially. It also only uses biomass as
a factor for prediction, and does not currently detect wildfires. 6\. Specific
tailor-made wildfire prevention algorithms have been created in regions such
as Spain [2], France [3], Turkey [4], the US [5], and in Europe [6]. However,
these systems can predict wildfires, but they are not end-to-end solutions
that take action to prevent the wildfires from starting (removing fuel,
spraying fire retardant, etc.). They are region-specific (use only that
region’s data) and have no wildfire detection capabilities. After looking at
the current solutions developed by the scientific community, it is clear that
there is a need for a system that can predict hotspots, detect wildfires and
send a drone to the location to prevent and suppress wildfires. Prediction and
detection of these wildfires can help prevent and mitigate the disastrous
effects.

##  II Design

Figure 1: High-level overview of WPDSS.

###  II-A Design Criteria

The AI-based Wildfire Prevention, Detection and Suppression System should be
able to use real-time satellite data and artificial intelligence to:
automatically identify hotspots, automatically identify current wildfires,
attain higher than 90% accuracy in hotspot and wildfire identification and
automatically notify a drone in the area to fly to the hotspot or current
wildfire.

###  II-B Artificial Intelligence

The Wildfire Prevention, Detection and Suppression System identifies hotspots
and wildfires using artificial intelligence. A few terms regarding artificial
intelligence are clarified. 1\. Machine Learning: Machine learning is an
intelligent system with a human-like brain that learns through data. A branch
of machine learning is supervised machine learning. Supervised learning
involves giving the algorithm training datasets with labels. Once
familiarized, it is given real-time data which it labels on its own. 2\.
Masking: Masking is the process of eliminating areas where wildfires cannot
and will not occur. These specific areas are not included in the AI Random
Forest Model. Masking is performed on areas with rainfall, ice, oceans, etc.
3\. Underfitting: Underfitting is when the model performs poorly on the
training data because the model is unable to capture the relationship between
the input examples (X) and the target values (Y). 4\. Overfitting: Overfitting
is when the model performs well on the training data but does not perform well
on the evaluation data. This is because the model is memorizing the data it
has seen and focuses on patterns relative to the dataset. When the model is
tested on a new dataset, it cannot reapply the patterns it originally
formulated. 5\. Evaluation Metrics: Precision shows out of the number of times
a model predicted 0 or 1, how often was it correct? Recall shows out of the
total actual 0 or 1, how many of them were correctly predicted as 0 or 1?

###  II-C Parameters

The Wildfire Prevention, Detection and Suppression System uses a series of
model parameters to accurately predict hotspots and identify wildfires. To
predict hotspots, WPDSS used data of land cover, wind speed, precipitation,
soil moisture, temperature, and vegetation as input to the model. The
satellite and meteorological data was loaded from the National Land Cover
Dataset (NCLD), NCEP GFS Weather (NOAA), and MODIS Aqua/Terra Satellite (NASA)
of parameters. 1) Land Cover is the type of land surface. The dataset provides
data on thematic class (for example, urban, agriculture, and forest), percent
impervious surface, and perfect tree canopy cover. WPDSS uses land cover to
mask areas with open water, perennial ice/snow, developed land (open space,
low, medium and high intensity), barren land (rock, sand, clay), moss, woody
wetlands and emergent herbaceous wetlands. 2) Wind Speed shows the speed of
the wind in meters per second. Wind plays a major role in the determination of
hotspots and wildfires as strong gusts of wind can spread wildfire flames
quickly. 3) Precipitation rate is the amount of rainfall in an area measured
in kg/(m2̂ s). Precipitation is shown to have a link with wildfires, as the
lack of precipitation has correlated with increased wildfire risk. [7] 4) Soil
Moisture is the volumetric soil moisture content around 0.0-0.1 m depth below
the land layer level. Soil moisture affects wildfires by controlling biomass
growth, fire fuel availability, and determines vegetation moisture content (O,
S, 2020). 5) Temperature measures the amount of hotness or coldness. This
temperature dataset consists of 2 datasets which are combined for accurate
depictions of temperature. 6) Vegetation is the amount of greenery of an area.
In this model, the Normalized Difference Vegetation Index (NDVI) is a dataset
that estimates the density of green on an area of land. The dataset formula is
[ndvi = (nir - red) / (nir + red)]. Whereas, to detect wildfires, WPDSS uses
data of land cover, humidity, temperature, vegetation, ozone, CO2 and burned
area as input to the model. The satellite and meteorological data was loaded
from the GOES 17 Satellite (NASA and NOAA), National Land Cover Dataset
(NCLD), NCEP GFS Weather (NOAA), and MODIS Aqua/Terra Satellite (NASA) of
parameters. The model parameters of land cover, vegetation, and temperature
are the same as the prevention model. The other 4 parameters are described as
follows: 1) Humidity is the amount of water vapor present in the air. The
dataset consists of specific humidity at the level of 2 meters above ground.
2) Burned area shows the burnt land. This consists of 2 datasets Near Infrared
and short-wave infrared combined using the formula bai = (nir - swir) / (nir +
swir). 3) Ozone is a product of wildfires and can be used to detect wildfires.
4) CO2, carbon dioxide, is emitted through wildfires. When forests burn, all
the CO2 stored in the trees gets released into the atmosphere.

###  II-D Determination of Model

To determine the type of method used for the prediction of hotspots and
detection of wildfires, three different models were tested, and the one with
the highest accuracy was used in the classifier. One of the first methods
tested used OpenCV and image processing to identify hotspots and wildfires.
When using open CV and image processing, an accuracy of 33.3% was obtained.
Logistic regression is a statistical analysis method used for binary
classification based on a dataset. A logistic regression model classifies the
category by analyzing the relationship between the model parameters. Based on
the historical data and previous hotspot and wildfire example, it can
accurately classify new data. When tested with the logistic regression model,
an accuracy of 74.8% was attained. A random forest classifier is a binary
classification model that uses a series of decision trees to make the most
accurate decision. Decision trees are a method used by an algorithm in which
it uses a simple set of decisions to classify the input. The random forest
model uses multiple decision trees and chooses the answer that has appeared
the most times to get the most accurate result. To get the most accurate
results, parameters can be tuned such as n-estimators (the number of decision
trees in the classifier) and max-depth (the number of sequential branches on
the decision tree). When tested with the random forest classifier, an accuracy
of 96.6% was attained. After these tests, it was concluded that the Random
forest classifier was the best method for identification of hotspots and
wildfires.

###  II-E Training of WPDSS

To create and train WPDSS, the labeled training/testing dataset was created by
manually identifying hotspots, not-hotspots, wildfires, and not wildfires
using the model parameters. (one-time step for supervised machine learning)
Then, the dataset is split into training (80%) and testing (20%). The training
dataset is supplied as input to the Random Forest Classifier to learn the
hotspot and wildfire occurrence patterns. The testing dataset is supplied as
input to the random forest classifier to classify areas as hotspots, not-
hotspots, wildfires and not wildfires. Then, tune the n-estimators (the number
of decision trees in the random forest classifier) and max_depth (the number
of branches in a decision tree) parameters to increase classifier accuracy.
This is completed through a series of 3 steps: 1) modify n-estimators and
max_depth from 1-15 and store the test and train accuracy, 2) determine the
highest model accuracy from the data collected and then 3) evaluate the
classifier using accuracy, confusion matrix, precision and recall.

###  II-F WPDSS Methods

WPDSS is a combined system that consists of 2 separate random forest
classifiers. One of these algorithms predicts hotspots and the other detects
wildfires. It performs 4 main steps: preprocessing, learning, identification
of hotspots and wildfires, and drone deployment. First, WPDSS loads the real-
time NASA satellite data and NOAA meteorological data of the parameters to
enhance accuracy and precision of hotspot and wildfire detection. It uses the
process of masking to eliminate not-hotspots such as water bodies and rainfall
to limit the dataset size. WPDSS then runs the real-time data (with the
masking) through the model to identify hotspots and wildfires. It removes only
the hotspots where wildfire retardant is not allowed to be sprayed. It sends
the hotspot and wildfire coordinates to the drone. Then, a prototype tello
drone is deployed to the hotspots.

Figure 2: Flow Chart

###  II-G Drones

The solution improves upon current technology for prevention and suppression
because rather than using humans and helicopters, the system will be using
drones. Drones are automated and more cost-efficient than firefighters
spraying fire retardant in a location. Also, there are helicopters which can
spray fire retardant, but require a pilot, cannot spray in the night, costs
millions of dollars and cannot fly low.

##  III ANALYSIS

###  III-A Tuning of Random Forest Classifier Parameters

Figure 3: Tuning of N_estimators

There were a few patterns noticed in the data as seen in Fig. 3. When the
number of n-estimators ¡ 7, the model is underfitting the data. When the
number of n-estimators ¿ 7, the model is overfitting the data. From graphs 1
and 2, it was analyzed that for the highest model accuracy, the number of
n-estimators should be 7.

Figure 4: Tuning of max_depth

There were some patterns noticed in these 2 graphs as seen in Fig. 4. When the
max_depth ¡ 5, the model is underfitting the data. When the max_depth ¿ 5, the
model is overfitting the data. From graphs 3 and 4, it was analyzed that for
the highest model accuracy, the max_depth should be 5.

###  III-B Statistical Methods

In total, 378 labeled coordinates were used for the training and testing of
the wildfire prevention model, 189 being hotspots and 189 being not-hotspots.
The entire dataset was split into a training (80% → 302) and a testing dataset
(20% → 76). Whereas, for wildfire detection, 158 labeled coordinates were used
for the training and testing. In this dataset, there were considerably fewer
training data points of wildfires, and more that were not wildfires since this
mimics the real-life scenario. In real life, the majority of the areas will
not be wildfires, and only a small portion will be wildfires. This dataset was
split into a training dataset (80% → 126) and a testing dataset (20% → 32).
The statistical methods used to analyze the data were accuracy (training,
testing and cross-validation), confusion matrix, precision and recall. These
are the standard methods for analyzing a random forest AI classifier.

###  III-C Accuracy

Figure 5: Accuracy

WPDSS was able to attain 98.6% testing accuracy in identifying the hotspots
for wildfire prevention, when the n-estimators is 7 and max depth is 5. WPDSS
was able to attain 98.7% testing accuracy in detecting wildfires. Cross
Validation accuracy for both wildfire prevention and detection is 96%. The
formula for accuracy used is: TP+TN/TP+TN+FP+FN.

Cross validation is a metric used to evaluate machine learning models with a
limited training dataset. It has a parameter called k, which is the number of
groups that the dataset is split into and k is dependent on the training
dataset size. The process of testing the data on unseen data to estimate how
the model will perform is called k-fold cross validation.

###  III-D Confusion Matrix

Figure 6: Confusion Matrix

Out of the 76 sample coordinates of the testing dataset for wildfire
prevention, 75 were classified correctly. WPDSS predicted the 36 true
positives (identified as a hotspot, and it is a hotspot) and 39 true negatives
(not identified as a hotspot and is not a hotspot). It predicted 1 false
positive (not a hotspot but identified as a hotspot), and 0 false negatives
(hotspots being missed). Out of the 79 sample coordinates of the testing
dataset for wildfire detection, 78 were classified correctly. WPDSS predicted
the 70 true positives (identified as a wildfire, and it is a wildfire) and 8
true negatives (not identified as a wildfire and is not a wildfire). It
predicted 1 false negative (wildfire being missed), and 0 false positives (not
a wildfire but identified as a wildfire).

###  III-E Precision and Recall

Figure 7: Precision and Recall

WPDSS was 98.5% precise while identifying a hotspot or not hotspot and was
able to recall its previous accurate decisions 98.5% of the time. WPDSS was
99.5% precise while identifying a wildfire and was able to recall its previous
accurate decisions 94.5% of the time.

###  III-F Drone Analysis

Analysis based on use of DJI Agras MG-1 drone to spray fire retardant at a
hotspot coordinates

Figure 8: Comparison of Drones

Total distance round trip to spray a hotspot coordinate
= 1.0 degree latitude (111 m) x 1.0 degree longitude (111 m)
= 12321 m² / 4047 (conversion to acres: 1 acre = 4047 m²)
= 3.04 acres
Drone’s Spraying Capacity (If drone is carrying 10 kg payload, flight time is
reduced to 10 minutes)
= 10 acres per hour (1 hour = 60 min)
= 1.66 acres per 10 minutes

To conclude, number of drones to spray 1 hotspot
= 3.04 acres / 1.66 acres
= 2 drones for spraying to be completed in 10 minutes OR 1 drone taking 2
round trips to complete spraying in 30 minutes
For Wildfire Prevention and Detection/Suppression,* the above analysis is used
for drone deployment.
* For detection/suppression rather than using the area in part a, leverage the perimeter to prevent the wildfire from spreading.

##  IV CONCLUSION

The AI-based Wildfire Prevention, Detection and Suppression System
successfully meets the design criteria. It uses real-time satellite data from
NASA and NOAA of soil moisture, precipitation, vegetation, wind speed, land
cover, and temperature to identify hotspots using artificial intelligence with
a 98.6% accuracy. It uses real-time satellite data from NASA and NOAA of
humidity, CO2, ozone, burned area, vegetation, land cover, and temperature to
detect wildfires using artificial intelligence with a 98.7% accuracy. This
system notifies drones to fly to the hotspot coordinates. It was also
theoretically proved that the usage of a DJI Agras MG-1 Drone is sufficient
enough to spray fire retardant at 1 hotspot coordinate. The AI-based Wildfire
Prevention, Detection and Suppression System can save precious lives, preserve
and protect ecosystems and biodiversity and avert huge economic losses (e.g.
in the US alone, 10.3 million acres from burning and $130-$150 billion
annually) and help rescue the impacts of climate change,. The power of the AI-
based Wildfire Prevention, Detection and Suppression System developed can be
applied to any location in the world to identify hotspots and detect
wildfires, deploy drones to the hotspots and wildfires and therefore prevent
wildfires and suppress wildfires.

##  V  FUTURE RESEARCH

After discussions with Stanford professor Eric Appel, who has built an
environmentally friendly fire retardant called fortify which can prevent fire
ignition for 6 months, Fortify will be used in WPDSS as the fire retardant.
Also, Intel has demonstrated the ability to deploy and control two thousand
drones simultaneously. Using the capability of controlling Intel’s drones
simultaneously, the wildfires can be prevented and suppressed at a faster
speed. Integrating all these efforts and creating a low cost system using
these AI based algorithms, wildfire can be prevented and suppressed.

## References

  * [1] A., M. Gill, Scott, L., Stephens, Geoffrey, J., and Cary, “The worldwide ”wildfire” problem.” _Ecological applications : a publication of the Ecological Society of America_ , vol. 23 2, pp. 438–54, 2013. [Online]. Available: https://api.semanticscholar.org/CorpusID:28398271
  * [2] A. Alonso-Betanzos, O. Fontenla-Romero, B. Guijarro-Berdiñas, E. Hernández-Pereira, M. I. Paz-Andrade, E. Jimenez, J. L. Legido, and T. Carballas, “An intelligent system for forest fire risk prediction and fire fighting management in galicia,” _Expert Syst. Appl._ , vol. 25, pp. 545–554, 2003. [Online]. Available: https://api.semanticscholar.org/CorpusID:16204189
  * [3] J. F. i Jové, P. F. i Casas, A. G. Petit, and J. Casanovas, “Firefight: A decision support system for forest fire containment,” 2014. [Online]. Available: https://api.semanticscholar.org/CorpusID:166367775
  * [4] M. U. Gumusay and K. Şahin, “Visualization of forest fires interactively on the internet,” _Scientific Research and Essays_ , vol. 4, pp. 1163–1174, 2009. [Online]. Available: https://api.semanticscholar.org/CorpusID:86614530
  * [5] E. K. Noonan-Wright, T. S. Opperman, M. Finney, G. T. Zimmerman, R. C. Seli, L. M. Elenz, D. E. Calkin, and J. R. Fiedler, “Developing the us wildland fire decision support system,” _Journal of Combustion_ , vol. 2011, pp. 1–14, 2011. [Online]. Available: https://api.semanticscholar.org/CorpusID:55117674
  * [6] L. P. Corgnati, M. Gabella, and G. Perona, “Firecast system - previsional fire danger index computation system for alpine regions,” 2008. [Online]. Available: https://api.semanticscholar.org/CorpusID:128890454
  * [7] F. Chen, S. Niu, X. Tong, J. Zhao, Y. Sun, and T. He, “The impact of precipitation regimes on forest fires in yunnan province, southwest china,” _The Scientific World Journal_ , vol. 2014, 2014. [Online]. Available: https://api.semanticscholar.org/CorpusID:18500053