# Improved Modeling of the Discrete Component of the Galactic Interstellar
$\gamma$-ray Emission and Implications for the Fermi–LAT Galactic Center
Excess

Christopher M. Karwin<EMAIL_ADDRESS>Department of Physics and
Astronomy, Clemson University, Clemson, South Carolina 29634, USA  Department
of Physics and Astronomy, University of California, Irvine, California 92697,
USA  Alex Broughton<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of California, Irvine, California 92697, USA  Simona Murgia
<EMAIL_ADDRESS>Department of Physics and Astronomy, University of California,
Irvine, California 92697, USA  Alexander Shmakov<EMAIL_ADDRESS>Mohammadamin Tavakoli<EMAIL_ADDRESS>Pierre Baldi<EMAIL_ADDRESS>Department of Computer Science, University of California, Irvine, California
92697, USA

###### Abstract

The aim of this work is to improve models for the $\gamma$-ray discrete or
small-scale structure related to $\ce{H_{2}}$ interstellar gas. Reliably
identifying this contribution is important to disentangle $\gamma$-ray point
sources from interstellar gas, and to better characterize extended
$\gamma$-ray signals. Notably, the Fermi–LAT Galactic center (GC) excess,
whose origin remains unclear, might be smooth or point-like. If the data
contain a point-like contribution that is not adequately modeled, a smooth GC
excess might be erroneously deemed to be point-like. We improve models for the
$\ce{H_{2}}$-related $\gamma$-ray discrete emission for a 50${}^{\circ}\times
1^{\circ}$ region along the Galactic plane via $\ce{H_{2}}$ proxies better
suited to trace these features. We find that these are likely to contribute
significantly to the $\gamma$-ray Fermi–LAT data in this region, and the
brightest ones are likely associated with detected Fermi–LAT sources, a
compelling validation of this methodology. We discuss prospects to extend this
methodology to other regions of the sky and implications for the
characterization of the GC excess.

††preprint: APS/PRL

The Galactic $\gamma$-ray interstellar emission (IE) traces interactions of
cosmic rays (CRs) with the interstellar medium, and constitutes most of the
$\gamma$-ray emission observed by Fermi–LAT. Uncertainties in modeling the IE
are large and difficult to constrain, and they impact the study of other
$\gamma$-ray sources in the Fermi–LAT data, point-like as well as extended. We
focus on modeling the small-scale structure of the IE which, if not robustly
captured by the model, confuses the determination of point sources, especially
along the Galactic plane [1]. It was shown in [2] that a large fraction of the
point sources detected by Fermi–LAT in the Galactic Center (GC) region could
be misidentified gas structure, and strongly dependent in number and spatial
distribution on the IE model employed to analyze the data. This result
indicates the presence of point-like emission in the data arising from
unmodeled structure in the interstellar gas, and it underlines the importance
of accurately modeling this component to reliably identify point sources in
the $\gamma$-ray data.

More reliable modeling of the small scale gas structure could also impact the
characterization of extended sources. A prominent example is the Fermi–LAT GC
excess [3, 4, 5, 6, 7, 8, 2]. Striking features of this excess are its spatial
morphology and spectrum which are compatible with annihilating dark matter
(DM). Alternative explanations have been proposed, with the leading one
attributing the signal to a collection of discrete emitters such as an
unresolved point source population of millisecond pulsars. The origin of the
Fermi-LAT GC excess remains a debated topic [9]. This debate could be settled
by determining whether the spatial morphology of the excess is consistent with
a smooth distribution, as predicted for DM, or with the cumulative emission of
a collection of point-like emitters. Statistical techniques have been employed
to accomplish this [10, 11, 12, 13], including the non-Poissonian template fit
(NPTF) which can detect upward fluctuations in the photon statistics above
Poisson noise which are associated with a collection of point sources.
However, the results of the NPTF technique strongly depend on the modeling of
the IE, and an uncontroversial resolution to the origin of the GC excess has
not yet been reached. We posit that if the interstellar gas is more structured
and point-like than current models predict, the unresolved point source
contribution in the data could be erroneously inflated by the fainter
component (also below detection threshold) of the small scale gas features. In
particular, statistical methods such as the NPTF might attribute this
component to a smooth GC excess and conclude it is point-like. These
uncertainties might therefore hinder the ability to distinguish the smooth
versus point-like duplicity of the excess. We note that wavelet decomposition
is another statistical technique that has been employed to resolve the GC
excess. While these studies are not as impacted by the IE model, the related
results on the nature of the GC excess remain uncertain nonetheless [10, 14,
15].

In this work, we present a novel approach to improve modeling of the small-
scale structure in the interstellar gas. The IE is due primarily to CRs
interacting with the interstellar hydrogen gas (and radiation field), in
molecular ($\ce{H_{2}}$), atomic ($\ce{H\textsc{i}}$), and ionized
($\ce{H\textsc{ii}}$) forms. $\ce{H_{2}}$ and $\ce{H\textsc{i}}$ are highly
structured compared to $\ce{H\textsc{ii}}$, and this structure is traced by
the related $\gamma$-ray emission. In this work, we focus on $\ce{H_{2}}$ gas
because of its high degree of structure and our methodology hinges on the
availability of additional proxies better suited to model it. The impact of
the other components of the $\gamma$-ray IE (specifically $\ce{H\textsc{i}}$
and dark gas [16]), and related uncertainties, is not considered in this study
and will be the focus of later work. Since $\ce{H_{2}}$ does not emit at a
characteristic radio frequency, other molecules are used to trace its
distribution. In particular, carbon monoxide ($\mathrm{{}^{12}CO}$, or
$\ce{CO}$ hereafter) is typically employed as a proxy. Radio surveys trace the
distribution of $\ce{CO}$ across the sky and the $\ce{H_{2}}$ column density
is inferred by scaling the $\ce{CO}$ content by a conversion factor (${\rm
X_{CO}}$) which gives the ratio of the integrated $\ce{CO}$ line emission to
the $\ce{H_{2}}$ column density. The bulk of the $\ce{H_{2}}$ is traced
following this methodology, and the survey of the $\ce{CO}$ $J$=1–0 transition
line from [17] has been widely employed. However, $\ce{CO}$ is typically
optically thick in the denser cores of molecular clouds and it underestimates
the total $\ce{H_{2}}$ column density there. $\gamma$-ray IE models that
employ $\ce{CO}$ to trace $\ce{H_{2}}$ gas may therefore underestimate its
finer structure. This limitation can be addressed by exploiting $\ce{CO}$
isotopologues, $\ce{{}^{13}CO}$ and $\ce{C^{18}O}$, also found in $\ce{H_{2}}$
clouds. Although rarer, these isotopologues are not as optically thick and
therefore more reliable to probe dense cores.

We employ the data from the Mopra Southern Galactic Plane CO Survey (data
release 3) [18] to trace the denser $\ce{H_{2}}$ regions. The survey covers a
50 square degree region, spanning Galactic longitudes $l$ = 300∘-350∘ and
latitudes $|b|\leq$ 0.5∘, and it targets the $J$ = 1–0 transitions of
$\ce{CO}$, $\ce{{}^{13}CO}$, $\ce{C^{18}O}$, and $\ce{C^{17}O}$. The data
reduction process, which involves six main stages of processing to perform the
band-pass correction and background subtraction, yields data cubes of the
brightness temperature for each spectral line, as a function of Galactic
coordinates and local standard of rest velocity. Mopra’s observations are
carried out in $1^{\circ}\times 1^{\circ}$ segments, which we combine into a
mosaic for the full region. To ensure the highest quality data in the rare
isotopologues, only pixels for which the brightness temperature exceeds the
$1\sigma$ noise provided by Mopra are used. This step is necessary since the
noise could otherwise predict spurious $\gamma$-ray emission. In this analysis
we utilize the $\ce{CO}$ and $\ce{{}^{13}CO}$ data, as the $\ce{C^{18}O}$
emission is extremely sparse, and the $\ce{C^{17}O}$ emission was too faint to
be detected in the survey. The spatial resolution is $0.6^{\prime}$ and the
spectral resolution is 0.1 km $\mathrm{sec^{-1}}$. For the region it surveys,
the Mopra data provide a sharper view of the $\ce{CO}$ emission compared to
[17], in addition to probing the rarer $\ce{H_{2}}$ tracers.

We calculate the $\ce{H_{2}}$ column density corresponding to Mopra’s
$\ce{CO}$ and $\ce{{}^{13}CO}$, referred to as
$N(\mathrm{H_{2}})_{\mathrm{CO}}$ and $N(\mathrm{H_{2}})_{\mathrm{CO13}}$,
respectively. Following the method of [19], the gas is separated into 17
Galactocentric radial bins based on its velocity and therefore corresponding
to different distances from the GC, assuming uniform circular motion about the
GC. The $N(\mathrm{H_{2}})_{\mathrm{CO}}$ is determined as:

$N(\mathrm{H_{2}})_{\mathrm{CO}}=W(\mathrm{CO})\times X_{\mathrm{CO}},$ (1)

where $W(\mathrm{CO})$ is the line strength of the CO gas, and we adopt the
radially-dependent $X_{\mathrm{CO}}$ from [19] for the 17 radial bins.

$N(\mathrm{H_{2}})_{\mathrm{CO13}}$ is evaluated as:

$N(\mathrm{H_{2}})_{\mathrm{CO13}}=N(\mathrm{CO13})\times\Big{[}\frac{\mathrm{H_{2}}}{\mathrm{CO13}}\Big{]}.$
(2)

Following [20, 18], we derive the $\ce{{}^{13}CO}$ column density
$N(\mathrm{CO13})$ assuming that the gas is in local thermodynamic equilibrium
at a fixed excitation temperature of 10 K. For the abundance ratio
$[\frac{\mathrm{H_{2}}}{\mathrm{CO13}}]$, we choose the upper bound of the
range (2.7 - 7.5) $\times 10^{5}$ provided in [21], in an attempt to assess
the maximal impact on the $\gamma$-ray data, although it is not excluded that
larger values are possible.

The CR propagation code GALPROP (v56) [22, 23, 24, 25, 26, 27, 28, 29, 30, 31,
32] is used to evaluate the $\gamma$-ray intensity maps for the
$\mathrm{H_{2}}$-related emission traced by $\ce{CO}$ and $\ce{{}^{13}CO}$.
GALPROP self-consistently calculates spectra and abundances of Galactic CR
species and associated diffuse emissions ($\gamma$-rays, but also radio,
X-rays) in 2D and 3D. We adopt the same GALPROP input parameters as described
in [33]. GALPROP returns radially dependent $\gamma$-ray intensity all-sky
maps (in 17 Galactocentric annuli) which allow us to determine, for the region
observed by Mopra, the additional contribution in $\gamma$-rays from the
$\mathrm{H_{2}}$ traced by $\ce{{}^{13}CO}$, and assess its significance in
simulated Fermi–LAT data.

Figure 1: The Modified Excess Template, calculated as the difference between
the Modified Map and the baseline $\ce{CO}$ map. The color scale shows the
$\gamma$-ray intensity. The solid cyan circles show unassociated $\gamma$-ray
point sources in version 4FGL-DR2 of the fourth Fermi–LAT catalog [34], and
the dashed green circles show new sources, not in the 4FGL-DR2, that we find
in this work (see text). The radius of each source corresponds to the 95%
localization uncertainty.

The difference between $N(\mathrm{H_{2}})_{\mathrm{CO13}}$ and
$N(\mathrm{H_{2}})_{\mathrm{CO}}$ gives an estimate of the $\ce{H_{2}}$ that
is missed in dense regions when only CO is used as a tracer. To calculate the
corresponding intensity in $\gamma$-rays, we define the “Modified Map”. This
is calculated in the following way: for pixels with
$N(\mathrm{H_{2}})_{\mathrm{CO13}}$ $>$ $N(\mathrm{H_{2}})_{\mathrm{CO12}}$,
the value of $N(\mathrm{H_{2}})_{\mathrm{CO12}}$ is replaced with the value of
$N(\mathrm{H_{2}})_{\mathrm{CO13}}$. We define the “Modified Excess Template”
as the difference between the Modified Map and the baseline $\ce{CO}$ map,
which accounts for the additional $\mathrm{H_{2}}$-related $\gamma$-ray
emission not included in current IE models, shown in Fig. 1.

We evaluate the significance of the Modified Excess Template in the Fermi–LAT
data by simulating the data collected between 2008 August 04 to 2020 November
11 ($\sim$12 years). The simulated events have energies in the range $1-100$
GeV and are binned in 8 energy bins per decade, for event class P8R3_CLEAN
(FRONT+BACK). The analysis is performed using Fermipy (v0.19.0), which
utilizes the Fermitools (v1.2.23). In these simulations, we only focus on the
$\mathrm{H_{2}}$-related $\gamma$-ray emission, and exclude all other
components of the $\gamma$-ray sky. The goal is to assess the significance of
the Modified Excess Template, i.e. the contribution of the newly modeled
$\mathrm{H_{2}}$ fine structure, in the optimistic scenario where all other
components are satisfactorily modeled. The simulated events trace the
$\mathrm{H_{2}}$-related $\gamma$-ray emission modeled with the Modified Map.
The simulated data are then fit based on a binned maximum likelihood method to
a model that includes two components, the $\gamma$-rays traced by the baseline
$\ce{CO}$ map and the Modified Excess Template. The latter is assigned the
spectrum determined by GALPROP, and its normalization is free to vary in the
fit. The normalization of the $\ce{CO}$ map is also free to vary and its
spectrum constrained to that calculated by GALPROP. As mentioned above, the
$\gamma$-ray flux is calculate in 17 radial bins, since the predicted
$\mathrm{H_{2}}$-related $\gamma$-ray emission depends on the CR density,
which is a function of Galactocentric radius. In the simulations, the total
emission is integrated along the line-of-sight. Moreover, the individual maps
have a high level of degeneracy. We therefore combine the maps into 4 radial
bins. Specifically, we combine bins 1-6, 7-10, 11-13, and 14-17, which we
refer to as A1, A2, A3, and A4.

We simulate 1000 realizations of Fermi–LAT data, and calculate the Test
Statistics (TS) for the nested models ($\mathrm{-2log(L_{0}/L)}$, where
$\mathrm{L_{0}}$ corresponds to value of the likelihood function for the null
hypothesis (CO baseline), and $\mathrm{L}$ to the alternative hypothesis (CO
baseline and Modified Excess Template.) The statistical significance is
approximated by $\sigma\approx\sqrt{\mathrm{TS}}$. The distribution of the
$\sqrt{TS}$ for the 1000 simulations is shown in Fig. 2. The mean of the
distribution is 48.30$\pm$1.02, for the 50 squared degree Mopra region, and
therefore very significant in a scenario where other components contributing
to the Fermi–LAT data are perfectly modeled, and if the $\gamma$-ray emission
traced by $\ce{{}^{13}CO}$ is at the high end of the range we have considered
(with the caveats discussed above.)

Figure 2: Distribution of the statistical significance
($\sigma\approx\sqrt{\mathrm{TS}}$) of the Modified Excess Template for the
50${}^{\circ}\times$1∘ region covered by Mopra for 1000 realizations of
$\sim$12 years of Fermi–LAT data. A fit with a Gaussian function is overlaid.
Figure 3: Fractional count residuals as a function of energy. The green and
yellow bands show the 68% and 95% confidence regions from 1000 simulations,
respectively. As an example, we also plot the results for a single simulation,
shown with black data points.

The fractional residuals as a function of energy for the 1000 simulations are
shown in Fig. 3. They are consistent with zero, as expected. Fig. 4 shows the
distributions of the flux of each model component, including the Modified
Excess Template. In the Mopra region, the mean integrated flux is $(7.2\pm
0.2)\times 10^{-8}\ \mathrm{ph\ cm^{-2}\ s^{-1}}$. Overall, the Modified
Excess Template accounts for a fair fraction (15.4%) of the total
$\ce{H_{2}}$-related $\gamma$-ray emission in the Mopra region. For
comparison, the integrated GC excess flux in a 15${}^{\circ}\times$15∘ region
around the GC, which is 4$\times$ larger than the region observed by Mopra, is
in the range $18.3-25.0\times 10^{-8}$ $\mathrm{ph\ cm^{-2}\ s^{-1}}$ (from
[2]). In intensities, the GC excess corresponds to $2.67-3.65\times
10^{-6}~{}\mathrm{ph\ cm^{-2}\ s^{-1}\ sr^{-1}}$ compared to $4.72\times
10^{-6}~{}\mathrm{ph\ cm^{-2}\ s^{-1}\ sr^{-1}}$ for the Modified Excess
Template. Because of the different spatial morphology of the GC excess and
$\ce{H_{2}}$-related $\gamma$-ray emission traced by $\ce{{}^{13}CO}$, a
direct comparison is unwarranted and we do not expect the latter to account
for the majority of the GC excess. However, the $\mathrm{H_{2}}$ emission
extends beyond the latitudes considered in this study and, albeit dimmer at
higher latitudes, the estimates provided here do not indicate its small scale
component to be negligible. This contribution could therefore confuse the GC
excess morphology, depending on the exact $\ce{H_{2}}$ emission outside of
$\pm$0.5∘ in latitude, because the techniques that track the point-like
fluctuations could erroneously ascribe unaccounted point-like emission
(originating from $\ce{H_{2}}$ gas in this case) to the GC excess [12].
Considerations based on the spectrum of the unresolved source population could
be powerful to settle this degeneracy, however a robust determination is
required.

Figure 4: The distributions of fluxes for each component of the best-fit model
using the Modified Excess Template over the full region for 1000 realizations
of $\sim$ 12 years of Fermi–LAT data. Note that one of the annuli (A4) has
negligible contribution and the normalization was fixed in the fitting process
(corresponding to a flux of $4.98\times 10^{-9}$
$\mathrm{p}\mathrm{h}\,\mathrm{c}\mathrm{m}^{-2}\,\mathrm{s}^{-1}$).

We identify the structure traced by $\ce{{}^{13}CO}$ that is bright enough to
be detected as a Fermi–LAT point source and compare it to the unassociated
sources in version 4FGL-DR2 of the fourth Fermi–LAT catalog of $\gamma$-ray
sources [34]. A significant overlap would indicate that its brightest
contribution has already been detected by Fermi–LAT and validate the
methodology to trace this component. To this end, we perform a likelihood fit
where, instead of including the Modified Excess Template in the model, we only
consider the baseline $\ce{CO}$ map and find additional point sources using
the fermipy find_srcs function. The method calculates TS maps, and identifies
point sources based on peaks in the TS. We use a test point source modeled
with a power law spectrum with spectral index=2, and a minimum TS threshold of
TS $\geq$ 16\. An index of 2.7 is compatible with CRs interacting with gas,
and if used, it yields consistent results. In total we find 23 new point
sources with TS$\geq$16, which are shown with green circles in Fig. 1. For
comparison, we also overlay 4FGL-DR2 unassociated sources. We find that 8/23
($35\%$) of the new sources are spatially consistent with an unassociated
4FGL-DR2 source (based on an overlap of the 95% localization errors), which
accounts for 8/46 ($17\%$) of the total unassociated sources in the region. To
quantify the probability that the associations happen by chance, we perform
1000 simulations, where for each realization we randomly distribute 23 sources
in the Mopra region. For each source the Galactic coordinates are drawn from a
uniform distribution and we assign an error radius drawn from a Gaussian
distribution, with the mean and standard deviation determined from the 23
detected sources. From the 1000 simulations we find the average number of
random associations to be 3.16 $\pm$ 1.72, which corresponds to a p-value
(i.e. N($\geq 8$)/1000) of 0.012. It is therefore likely that the brightest of
the point-like emission traced by Mopra $\ce{{}^{13}CO}$ have already been
detected in the Fermi–LAT data.

The results presented here demonstrate that there likely is significant
structure in the $\ce{H_{2}}$-related $\gamma$-ray emission detected by
Fermi–LAT that is not currently included in the IE models. Since its spatial
morphology has point-like features, it directly affects the detection of
resolved and, potentially, unresolved point source populations in the
$\gamma$-ray data. More specifically, we have discussed the impact on the
4FGL-DR2 catalog and on the interpretation of the GC excess. Although in this
analysis we have focused on the region covered by Mopra, similar conclusions
would hold elsewhere because of the same limitations in tracing
$\mathrm{H_{2}}$. The methodology described here however can only be readily
applied to very small regions of the sky because of the paucity of
observations of the rare $\ce{H_{2}}$ tracers (e.g. [35]). Also, we cannot
simply extrapolate the Modified Excess Template with the less ambitious goal
of providing only an estimate of this emission since the Mopra data on which
it is based is tightly confined, especially in latitude. To address the
limited available observations, we have resorted to machine learning to
develop a methodology that predicts the distribution of the small scale
$\ce{H_{2}}$-related $\gamma$-ray emission for other regions of the sky based
on the existing Mopra data. Because of its complexity, we describe the machine
related work in a companion paper [36].

Conclusions.— We assess whether there is point-like emission in the
$\gamma$-ray data associated with the interstellar gas and not currently
included in models of the Galactic interstellar $\gamma$-ray emission. We
employ the data from the Mopra Southern Galactic Plane CO Survey, which
includes tracers of the small scale structure of the $\mathrm{H_{2}}$-related
$\gamma$-rays, to improve models by more accurately describing the point-like
features in the gas. We find that significant point-like emission originates
from $\mathrm{H_{2}}$ gas and the significance could be as high as
$\sqrt{TS}\sim$ 48 (depending on assumptions) in a 50${}^{\circ}\times$1∘
region covered by Mopra, corresponding to $\sim$ 15% of the modeled
$\ce{H_{2}}$-related $\gamma$-ray emission in the region. We also show that
this newly found point-like component may account for some fraction
($\lesssim$17%) of $\gamma$-ray point sources detected by Fermi–LAT in the
Galactic plane whose origin is so far unknown. That a significant number of
unidentified 4FGL-DR2 sources along the Galactic plane originate from
unmodelled structure in the gas is not unexpected, but here we develop a
robust and reliable methodology to identify the contribution of both bright
and dimmer components of the $\ce{H_{2}}$ gas discrete emission. Given the
significance of this emission, its contribution in the GC region could
introduce a significant systematic uncertainty in determining whether the GC
excess is smooth or point-like, and therefore more consistent with a dark
matter or pulsar interpretation, specifically when relying on statistical
techniques such as the NPTF. Identifying this component in the GC region could
be a crucial step to settle the origin of the GC excess and to reliably
determine whether it is smooth and therefore a signal from DM.

##  I Acknowledgements

We thank Troy Porter for many helpful discussions and insights. The work of
AS, MT, and PB is in part supported by grants NSF NRT 1633631 and ARO 76649-CS
to PB.

## References

  * Abdollahi et al [2020] S. Abdollahi et al, Fermi large area telescope fourth source catalog, The Astrophysical Journal Supplement Series 247, 33 (2020).
  * Ajello  _et al._ [2016] M. Ajello, A. Albert, W. Atwood, G. Barbiellini, D. Bastieri, K. Bechtol, R. Bellazzini, E. Bissaldi, R. Blandford, E. Bloom, _et al._ , Fermi-lat observations of high-energy $\gamma$-ray emission toward the galactic center, The Astrophysical Journal 819, 44 (2016).
  * Goodenough and Hooper [2009] L. Goodenough and D. Hooper, Possible evidence for dark matter annihilation in the inner milky way from the fermi gamma ray space telescope (2009), arXiv:0910.2998 [hep-ph] .
  * Hooper and Goodenough [2011] D. Hooper and L. Goodenough, Dark Matter Annihilation in The Galactic Center As Seen by the Fermi Gamma Ray Space Telescope, Phys. Lett. B 697, 412 (2011), arXiv:1010.2752 [hep-ph] .
  * Abazajian [2011] K. N. Abazajian, The consistency of fermi-LAT observations of the galactic center with a millisecond pulsar population in the central stellar cluster, Journal of Cosmology and Astroparticle Physics 2011 (03), 010.
  * Abazajian  _et al._ [2014] K. Abazajian, N. Canac, S. Horiuchi, and M. Kaplinghat, Astrophysical and dark matter interpretations of extended gamma ray emission from the galactic center, Physical Review D 90 (2014).
  * Calore  _et al._ [2015] F. Calore, I. Cholis, and C. Weniger, Background model systematics for the fermi GeV excess, Journal of Cosmology and Astroparticle Physics 2015 (03), 038.
  * Daylan  _et al._ [2016] T. Daylan, D. P. Finkbeiner, D. Hooper, T. Linden, S. K. Portillo, N. L. Rodd, and T. R. Slatyer, The characterization of the gamma-ray signal from the central milky way: A case for annihilating dark matter, Physics of the Dark Universe 12, 1 (2016).
  * Murgia [2020] S. Murgia, The fermi-lat galactic center excess: Evidence of annihilating dark matter?, Annual Review of Nuclear and Particle Science 70, 455 (2020), https://doi.org/10.1146/annurev-nucl-101916-123029 .
  * Bartels  _et al._ [2016] R. Bartels, S. Krishnamurthy, and C. Weniger, Strong support for the millisecond pulsar origin of the Galactic center GeV excess, Phys. Rev. Lett. 116, 051102 (2016), arXiv:1506.05104 [astro-ph.HE] .
  * Lee  _et al._ [2016] S. K. Lee, M. Lisanti, B. R. Safdi, T. R. Slatyer, and W. Xue, Evidence for unresolved $\gamma$-ray point sources in the inner galaxy, Phys. Rev. Lett. 116, 051103 (2016).
  * Leane and Slatyer [2019] R. K. Leane and T. R. Slatyer, Revival of the Dark Matter Hypothesis for the Galactic Center Gamma-Ray Excess, Phys. Rev. Lett. 123, 241101 (2019), arXiv:1904.08430 [astro-ph.HE] .
  * Zhong  _et al._ [2020] Y.-M. Zhong, S. D. McDermott, I. Cholis, and P. J. Fox, Testing the Sensitivity of the Galactic Center Excess to the Point Source Mask, Phys. Rev. Lett. 124, 231103 (2020), arXiv:1911.12369 [astro-ph.HE] .
  * Balaji  _et al._ [2018] B. Balaji, I. Cholis, P. J. Fox, and S. D. McDermott, Analyzing the gamma-ray sky with wavelets, Phys. Rev. D 98, 043009 (2018), arXiv:1803.01952 [astro-ph.HE] .
  * Zhong  _et al._ [2020] Y.-M. Zhong, S. D. McDermott, I. Cholis, and P. J. Fox, Testing the Sensitivity of the Galactic Center Excess to the Point Source Mask, Phys. Rev. Lett. 124, 231103 (2020), arXiv:1911.12369 [astro-ph.HE] .
  * Grenier  _et al._ [2005] I. A. Grenier, J.-M. Casandjian, and R. Terrier, Unveiling Extensive Clouds of Dark Gas in the Solar Neighborhood, Science 307, 1292 (2005).
  * Dame  _et al._ [2001] T. M. Dame, D. Hartmann, and P. Thaddeus, The Milky Way in Molecular Clouds: A New Complete CO Survey, Astrophys. J. 547, 792 (2001), arXiv:astro-ph/0009217 [astro-ph] .
  * Braiding  _et al._ [2018] C. Braiding, G. F. Wong, N. I. Maxted, D. Romano, M. G. Burton, R. Blackwell, M. D. Filipović, M. S. R. Freeman, B. Indermuehle, and J. Lau, et al, The mopra southern galactic plane co survey—data release 3, Publications of the Astronomical Society of Australia 35, 10.1017/pasa.2018.18 (2018).
  * Ackermann  _et al._ [2012] M. Ackermann, M. Ajello, W. Atwood, L. Baldini, J. Ballet, G. Barbiellini, D. Bastieri, K. Bechtol, R. Bellazzini, B. Berenji, _et al._ , Fermi-lat observations of the diffuse $\gamma$-ray emission: implications for cosmic rays and the interstellar medium, The Astrophysical Journal 750, 3 (2012).
  * Wilson  _et al._ [2009] T. L. Wilson, K. Rohlfs, and S. Hüttemeister,  _Tools of Radio Astronomy_ (2009).
  * Pineda  _et al._ [2008] J. E. Pineda, P. Caselli, and A. A. Goodman, CO Isotopologues in the Perseus Molecular Cloud Complex: the X-factor and Regional Variations, Astrophys. J. 679, 481 (2008), arXiv:0802.0708 [astro-ph] .
  * Moskalenko and Strong [1998] I. V. Moskalenko and A. W. Strong, Production and propagation of cosmic ray positrons and electrons, Astrophys. J. 493, 694 (1998), arXiv:astro-ph/9710124 [astro-ph] .
  * Moskalenko and Strong [2000] I. V. Moskalenko and A. W. Strong, Anisotropic inverse compton scattering in the galaxy, Astrophys. J. 528, 357 (2000), arXiv:astro-ph/9811284 [astro-ph] .
  * Strong and Moskalenko [1998] A. W. Strong and I. V. Moskalenko, Propagation of cosmic-ray nucleons in the galaxy, Astrophys. J. 509, 212 (1998), arXiv:astro-ph/9807150 [astro-ph] .
  * Strong  _et al._ [2000] A. W. Strong, I. V. Moskalenko, and O. Reimer, Diffuse continuum gamma-rays from the galaxy, Astrophys. J. 537, 763 (2000), [Erratum: ApJ 541,1109(2000)], arXiv:astro-ph/9811296 [astro-ph] .
  * Ptuskin  _et al._ [2006] V. S. Ptuskin, I. V. Moskalenko, F. C. Jones, A. W. Strong, and V. N. Zirakashvili, Dissipation of Magnetohydrodynamic Waves on Energetic Particles: Impact on Interstellar Turbulence and Cosmic-Ray Transport, Astrophys. J. 642, 902 (2006), astro-ph/0510335 .
  * Strong  _et al._ [2007] A. W. Strong, I. V. Moskalenko, and V. S. Ptuskin, Cosmic-ray propagation and interactions in the Galaxy, ARNPS 57, 285 (2007), arXiv:astro-ph/0701517 [astro-ph] .
  * Vladimirov  _et al._ [2011] A. E. Vladimirov, S. W. Digel, G. Johannesson, P. F. Michelson, I. V. Moskalenko, P. L. Nolan, E. Orlando, T. A. Porter, and A. W. Strong, GALPROP WebRun: an internet-based service for calculating galactic cosmic ray propagation and associated photon emissions, Comput. Phys. Commun. 182, 1156 (2011), arXiv:1008.3642 [astro-ph.HE] .
  * Jóhannesson  _et al._ [2016] G. Jóhannesson _et al._ , Bayesian analysis of cosmic-ray propagation: evidence against homogeneous diffusion, Astrophys. J. 824, 16 (2016), arXiv:1602.02243 [astro-ph.HE] .
  * Porter  _et al._ [2017] T. A. Porter, G. Johannesson, and I. V. Moskalenko, High-energy gamma rays from the milky way: Three-dimensional spatial models for the cosmic-ray and radiation field densities in the interstellar medium, Astrophys. J. 846, 23pp (2017).
  * Johannesson  _et al._ [2018] G. Johannesson, T. A. Porter, and I. V. Moskalenko, The Three-Dimensional Spatial Distribution of Interstellar Gas in the Milky Way: Implications for Cosmic Rays and High-Energy Gamma-Ray Emissions, Astrophys. J. 856, 45 (2018), arXiv:1802.08646 [astro-ph.HE] .
  * Génolini  _et al._ [2018] Y. Génolini, D. Maurin, I. V. Moskalenko, and M. Unger, Current status and desired precision of the isotopic production cross sections relevant to astrophysics of cosmic rays: Li, be, b, c, and n, PhRvC 98, 034611 (2018).
  * Karwin  _et al._ [2019] C. M. Karwin, S. Murgia, S. Campbell, and I. V. Moskalenko, Fermi-lat observations of $\gamma$-ray emission toward the outer halo of m31, The Astrophysical Journal 880, 95 (2019).
  * Ballet  _et al._ [2020] J. Ballet, T. H. Burnett, S. W. Digel, and B. Lott, Fermi Large Area Telescope Fourth Source Catalog Data Release 2, arXiv e-prints , arXiv:2005.11208 (2020), arXiv:2005.11208 .
  * Blackwell  _et al._ [2016] R. Blackwell, M. Burton, and G. Rowell, Mopra central molecular zone carbon monoxide survey status, Proceedings of the International Astronomical Union 11, 164–165 (2016).
  * Shmakov  _et al._ [2022] A. Shmakov, M. Tavakoli, P. Baldi, C. Karwin, A. Broughton, and S. Murgia, Deep learning models of the discrete component of the galactic interstellar $\gamma$-ray emission (2022).


# Lyapunov-driven deep reinforcement learning for edge inference empowered by
Reconfigurable Intelligent Surfaces

###### Abstract

In this paper, we propose a novel algorithm for energy-efficient, low-latency,
accurate inference at the wireless edge, in the context of 6G networks endowed
with reconfigurable intelligent surfaces (RISs). We consider a scenario where
new data are continuously generated/collected by a set of devices and are
handled through a dynamic queueing system. Building on the marriage between
Lyapunov stochastic optimization and deep reinforcement learning (DRL), we
devise a dynamic learning algorithm that jointly optimizes the data
compression scheme, the allocation of radio resources (i.e., power,
transmission precoding), the computation resources (i.e., CPU cycles), and the
RIS reflectivity parameters (i.e., phase shifts), with the aim of performing
energy-efficient edge classification with end-to-end (E2E) delay and inference
accuracy constraints. The proposed strategy enables dynamic control of the
system and of the wireless propagation environment, performing a low-
complexity optimization on a per-slot basis while dealing with time-varying
radio channels and task arrivals, whose statistics are unknown. Numerical
results assess the performance of the proposed RIS-empowered edge inference
strategy in terms of trade-off between energy, delay, and accuracy of a
classification task.

Index Terms—  Edge intelligence, inference, reconfigurable intelligent
surfaces, Lyapunov optimization, reinforcement learning.

##  1 Introduction

The development of the next generation of wireless communication systems,
known as 6G, is still at its infancy. The main challenge of 6G is to provide
an Artificial Intelligence (AI) and Machine Learning (ML) native communication
infrastructure. This concept is known as Edge ML/Edge AI [1]. Edge AI comes
with a twofold perspective, including both the benefits of AI/ML algorithms
exploited for network optimization and orchestration, and the benefit of a
powerful and efficient communication platform to process and distill large
volumes of data collected by heterogeneous devices. The final aim is to
accurately perform learning tasks within low end-to-end (E2E) delays, in the
most efficient way from different perspectives entailing energy, communication
overhead, etc. Edge AI will strongly benefit from Multi-access Edge Computing
(MEC) [2], thanks to the deployment of distributed computing resources close
to end users, namely in Mobile Edge Hosts (MEHs) that are, e.g. co-located
with radio Access Points (APs). This will allows mobile end devices to access
computing resources, albeit limited with respect to central clouds, in a fast,
secure, reliable, and possibly sustainable manner. In this work, we focus on
both roles of edge AI in future systems, devising a resource allocation
framework to enable dynamic energy efficient edge classification on data
collected by end devices, with target E2E delays and inference reliability
constraints. Besides the paradigm of communications for AI, in this work, the
twofold role of edge AI is explored through the marriage between model-based
Lyapunov stochastic network optimization, and Deep Reinforcement Learning
(DRL), with the latter playing the role of an AI-based optimization algorithm,
able to compensate the limitations of model-based optimization, namely
complexity and/or lack of system modeling. Previous works successfully
attempted to merge such tools [3, 4], however never focusing on the edge AI
use cases, i.e. not taking into account application performance. Instead, the
goal of this paper is to focus on the performance (in terms of accuracy) of an
edge classification task, along with typical MEC performance including End-to-
End (E2E) delay and energy consumption.

Related works. As illustrated in the overview papers [5, 1], an efficient
design of edge-inference hinges on several aspects, such as memory footprint
optimization [6], adaptive model selection [7], or goal-oriented optimization,
for instance adapting frame rate and resolution of video streaming for
efficient inference [8], [9]. The authors of [10] consider a video analytics
task, maximizing the average accuracy under a frame rate and delay constraint.
Also, the work in [11] proposed a joint management of radio and computation
resources for edge ML, hinging on data quantization to control the accuracy of
several learning tasks. Finally, the closest reference [12] proposes a dyanmic
joint optimization or radio, computing resources, and JPEG data compression,
fully based on Lyapunov optimization, assuming full knowledge of the system
model.

All the aforementioned works enabled edge computing and/or inference
considering the presence of a suitable wireless propagation environment.
Moving toward millimeter wave communications (and beyond), the performance of
MEC-oriented systems can be severely reduced due to poor channel conditions
and blocking events. In this context, a strong performance improvement can be
obtained exploiting Reconfigurable Intelligent Surfaces (RISs) [13, 14, 15,
16, 17], which are programmable surfaces made of hundreds of nearly passive
reflective elements controlled to realize dynamic transformations of the
wireless propagation environment, both in indoor and outdoor scenarios. The
inclusion of RISs in MEH systems offer a two-fold benefit: i) it alleviates
the effect of blockages, which leads to low offloading rates and poor
performance; ii) it enables better exploitation of the computing resources of
the edge server thanks to the improved offloading capabilities. Several works
in the literature have already exploited RISs to empower wireless
communications [18, 19, 20, 21]) and, very recently, also for computation
offloading [22, 23]. Finally, preliminary results on edge learning empowered
by RISs appear in [24], which considered static resource allocation for edge
inference, and in the work [25] that instead focused on adaptive federated
learning.

Contributions. In this paper, we focus on an edge inference task, aimed at
performing classification of data collected, distilled, and processed at the
edge of a wireless network endowed with MEC capabilities and RISs. Our design
aims at striking the minimum average power necessary by the system to perform
low-latency inference with accuracy constraints. To this aim, we propose an
online method based on the marriage of Lyapunov stochastic optimization DRL,
to dynamically optimize: i) data compression scheme; ii) users’ transmission
precoding and power; iii) RISs’ reflectivity parameters; iv) MEC resources.
Our method leads to the convergence of model-based (i.e., Lyapunov) and data-
driven (i.e., DRL) stochastic optimization into a single holistic framework,
thus paving the way towards fully reconfigurable networks endowed with
effective and efficient edge AI capabilities. The main novelties with respect
to the closest work [12] are the following: i) the function linking accuracy
and compression scheme is supposed to be unknown, and optimized through a DRL-
based approach; ii) the RIS is not present in [12], and it is optimized here
through the DRL-based approach, differently from [26], in which it is
optimized though a projected gradient descent method. The computation part is
handled as in [12, 26], thus it does not represent a novelty in this paper.
Finally, numerical results assess the performance of the proposed method.
Notation. Bold upper case letters denote matrices, $\textrm{Tr}(\cdot)$ denote
the trace operator, $|\cdot|$ denotes the determinant, and the superscript
$(\cdot)^{H}$ denotes the hermitian operator; $\mathbf{I}_{N}$ and
$\mathbf{0}_{N}$ denote the $N\times N$ identity matrix and the all-zeros
matrix, respectively. Finally, given a variable $X(t)$,
$\overline{X}=\lim_{T\to\infty}\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\\{X(t)\\}$
if the sequence converges, or $\overline{X}=\lim\ {\rm
sup}_{T\to\infty}\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}\\{X(t)\\}$ otherwise.

##  2 System Model

Let us consider a wireless setup with a set $\mathcal{K}$ of Mobile Devices
(MDs), each one equipped with $N_{u}$ antennas, an $N_{a}$ antennas Access
Point (AP), with a co-located MEH able to process tasks through a pre-trained
and pre-uploaded ML model (e.g. a deep neural network). Finally, an RIS with
$M$ reflecting unit elements is assumed to be available, and it can be
dynamically reconfigured to enhance MEC service performance. As in [12], we
model the inference process as a flow of request arrivals, which are first
buffered locally at each MD before transmission and, after transmission,
buffered remotely before computation. Time is organized in slots $t$ of equal
duration $\tau$. Each device has its own buffer and all MDs compete for radio
and computing resources. Then, uplink communication, along with computation,
is the focus of this work.

###  2.1 Channel model

An RIS-aided wireless channel can be modeled by two components: i) a direct
path from MD $k$ to the AP, whose time-varying coefficient are the elements of
a matrix $\mathbf{H}_{k,d}(t)\in\mathbb{C}^{N_{a}\times N_{u}}$, and ii) an
indirect path, created by the reflection of the RIS. The latter includes a
channel matrix $\mathbf{H}_{k,r}(t)\in\mathbb{C}^{M\times N_{u}}$ from MD $k$
to the RIS, and a channel matrix
$\mathbf{H}_{r,a}(t)\in\mathbb{C}^{N_{a}\times M}$ from the RIS to the AP.
Finally, the RIS is composed of passive elements, whose phases can be
dynamically and opportunistically reconfigured, and its response can be
written as a diagonal matrix $\mathbf{\Phi}(t)$, with diagonal elements
$\\{\mathbf{\Phi}_{i,i}(t)=e^{j\phi_{i}(t)}\\}_{i=1}^{M}$, where $\phi_{i}(t)$
denotes the (reconfigurable) phase shift of elements $i$. Then, the overall
channel is [20]

$\mathbf{H}_{k}(t)=\mathbf{H}_{k,d}(t)+\mathbf{H}_{r,a}(t)\mathbf{\Phi}(t)\mathbf{H}_{k,a}(t),\quad\forall
k\in\mathcal{K}.$ (1)

We assume that, at each slot, every MD selects a transmit precoding strategy,
based on current connect-compute system conditions. Then, denoting by
$\mathbf{F}_{k}(t)$ the input covariance matrix of MD $k$, the instantaneous
data rate can be written as follows:

$R_{k}(t)=W\log_{2}\left|\mathbf{I}_{N_{u}}+\frac{1}{\sigma^{2}}\mathbf{H}_{k}(t)\mathbf{F}_{k}(t)\mathbf{H}_{k}^{H}(t)\right|,\;\;\forall
k\in\mathcal{K}$ (2)

where $\sigma^{2}=N_{0}W$, with $N_{0}$ the noise power spectral density, and
$W$ bandwidth. All users are coupled by the RIS reflection.

###  2.2 Communication and computation queuing models

Similarly to [12], buffers are intended as units of patterns (or data unit).
At each slot, a generic MD $k$ accepts $A_{k}(t)$ new patterns into its
communication buffer, while transmitting previously buffered patterns through
the (RIS-aided) wireless connection with the AP, at a rate $R_{k}(t)$ (cf.
(2)). Then, denoting by $Q_{k}^{l}(t)$ the buffer size at time $t$, the
communication queue evolves as follows:

$Q_{k}^{l}(t+1)=\max\left(0,Q_{k}^{l}(t)-\left\lfloor\frac{\tau
R_{k}(t)}{n_{k,b}(c_{k}(t))}\right\rfloor\right)+A_{k}(t),$ (3)

where $n_{k,b}(c_{k}(t))$ is the number of bits encoding all patterns
transmitted at time $t$, which is a function of the data compression level
$c_{k}(t)\in\mathcal{C}_{k}$ at time slot $t$. The transmitted patterns join a
remote computation queue, which is drained by the MEH’s processing, i.e. the
issuing of inference results. Then, denoting by $f_{k}(t)$ the MEH’s CPU cycle
frequency assigned to user $k$, the queue evolves as:

$\displaystyle Q_{k}^{r}(t+1)=$
$\displaystyle\,\max\left(0,Q_{k}^{r}(t)-\left\lfloor\frac{\tau
f_{k}}{w_{k}}\right\rfloor\right)$
$\displaystyle\quad+\min\left(Q_{k}^{l}(t),\left\lfloor\frac{\tau
R_{k}(t)}{n_{k,b}(c_{k}(t))}\right\rfloor\right),$ (4)

where $w_{k}$ is the computation load, i.e. the number of CPU cycles needed to
output one inference result.

###  2.3 Inference performance indicators

Effective and efficient edge inference entails timing (i.e., the delay from
request issuing at the device until its treatment at the MEH), accuracy (e.g.,
correctly classified patterns), and energy. The average delay entails
communication and computation phases. Also, if the queues are strongly stable,
i.e. $\overline{Q}_{k}^{l(r)}<\infty$, the average E2E delay is finite and can
be written in closed form thanks to Little’s law:
$\overline{D}_{k}=\tau(\overline{Q}_{k}^{l}+\overline{Q}_{k}^{r})/\overline{A}_{k}$,
where $\overline{A}_{k}$ denotes the average number of arrivals per slot. In
this paper, we consider a generic metric $\mathcal{G}_{k}(c_{k})$ of inference
reliability, which we mildly assume to be a function of the employed
compression scheme $c_{k}$. Such generic function could be the accuracy of a
classification/regression/estimation task, or any other measure reflecting
effective operation of an ML model running in the MEH (e.g., classification
confidence). The aim will be to keep the long-term average of the accuracy
metric $\overline{\mathcal{G}}_{k}(c_{k})$ above a predefined threshold, set a
priori as an application requirement. Finally, as source of MDs’ power
consumption we consider the long-term average transmit power, which can be
expressed as $\overline{\sum_{k\in\mathcal{K}}\textrm{Tr}(\mathbf{F}_{k})}$.

##  3 Problem formulation

In this work, our aim is to guarantee energy-efficient edge inference with a
minimum level of accuracy and a given E2E delay. The problem is formulated as
follows:

$\displaystyle\underset{\\{\mathbf{F}(t)\\}_{k,t},\\{f_{k}(t)\\}_{k,t},\\{\phi(t)\\}_{i,t},\\{c_{k}(t)\\}_{k,t}}{\min}\;\;\overline{\sum_{k\in\mathcal{K}}{\textrm{Tr}(\mathbf{F}_{k})}}$
(5) subject to $\displaystyle(a)\;\overline{Q}_{k}^{l}<\infty,\;\forall
k\in\mathcal{K}\qquad\qquad(b)\;\overline{Q}_{k}^{r}<\infty,\;\forall
k\in\mathcal{K}$
$\displaystyle(c)\;\overline{\mathcal{G}}_{k}(c_{k})\geq\mathcal{G}_{k,\textrm{th}},\;\forall
k\in\mathcal{K}\quad\;\;(d)\;\mathbf{F}_{k}(t)\succcurlyeq 0,\;\forall
k\in\mathcal{K}$ $\displaystyle(e)\;\textrm{Tr}(\mathbf{F}_{k}(t))\leq
P_{k},\;\forall k\in\mathcal{K}\quad(f)\;c_{k}\in\mathcal{C}_{k},\;\forall
k\in\mathcal{K}$ $\displaystyle(g)\;\;\phi_{i}(t)\in[0,2\pi],\;i=1,\ldots M$
$\displaystyle(h)\;f_{k}(t)\geq
0\quad(i)\;\sum\nolimits_{k\in\mathcal{K}}f_{k}(t)\leq f_{\max}$

Besides queues’ stability in $(a)$ and $(b)$, the constraints of (5) have the
following meaning: $(c)$ the average accuracy is higher than a predefined
threshold; $(d)$ the covariance matrix of each user is positive semidefinite;
$(e)$ the instantaneous transmit power of each user is lower than a threshold;
$(f)$ the compression scheme belongs to a discrete set $\mathcal{C}_{k}$;
$(g)$ the RIS elements’ phases are chosen between $0$ and $2\pi$; $(h)$ the
CPU cycle frequency assigned to each MD is non negative; $(i)$ the sum of the
CPU cycle frequencies assigned to all MDs does not exceed the maximum MEH’s
frequency.

##  4 Lyapunov-driven DRL: The convergence of model-based and data-driven
optimization

Due to the lack of knowledge of channels and arrivals statistics, we now apply
Lyapunov stochastic optimization tools to decouple the long-term problem in
(5) into a sequence of simpler per-slot problems, whose solution only requires
the instantaneous observation of such context parameters. As it will be
explained in the sequel, part of the per-slot problem can be solved with
closed-form expressions and fast algorithms, while part of it is solved
through DRL, which can efficiently handle the inherent complexity and the lack
of knowledge of the involved functions. First of all, to handle constraint
$(c)$, we introduce a virtual queue $Z_{k}(t),\forall k\in\mathcal{K}$, whose
scope is to detect instantaneous violations of the constraint, and take
consequent actions to guarantee the long-term desired performance (i.e. the
accuracy above a threshold). In particular, the virtual queue evolves as
follows:

$Z_{k}(t+1)=\max(0,Z_{k}(t)-\epsilon(\mathcal{G}_{k}(c_{k}(t))-\mathcal{G}_{k,\textrm{th}})),\;\forall
k\in\mathcal{K}.$ (6)

By definition, virtual queue $Z_{k}$ grows whenever the constraint is not met
instantaneously, and it is drained otherwise. Also, if the mean rate stability
of $Z_{k}$ is guaranteed, i.e. $\lim_{T\to\infty}\mathbb{E}\\{Z(T)\\}/T=0$,
constraint $(c)$ is guaranteed [27]. To guarantee physical queue and virtual
queue stability, we first introduce the Lyapunov Function, a scalar measure of
the system’s congestion state [27]:

$L((t))=\frac{1}{2}\sum_{k\in\mathcal{K}}\left[Q_{k}^{l}(t)^{2}+Q_{k}^{r}(t)^{2}+Z_{k}^{2}(t)\right],$
(7)

with $\mathbf{\Lambda}(t)=[\\{Q_{k}^{l}(t)\\}_{\forall
k\in\mathcal{K}},\\{Q_{k}^{r}(t)\\}_{\forall
k\in\mathcal{K}},\\{Z_{k}(t)\\}_{\forall k\in\mathcal{K}}]$. From (7), we can
define the drift-plus-penalty function, i.e. the conditional expected
variation of the Lyapunov function over two successive time slots, penalized
by the objective function of problem (5):

$\Delta_{p}(\mathbf{\Lambda}(t))=\mathbb{E}\\{L(\mathbf{(}t+1))-L((t))+V\sum_{k\in\mathcal{K}}\textrm{Tr}(\mathbf{F}_{k}(t))|\mathbf{\Lambda}(t)\\}.$
(8)

Now, based on the theoretical findings in [27], our method proceeds by
instantaneously minimizing a suitable upper bound of (8) at each slot, based
on instantaneous observations of time-varying context parameters, and state
variables (physical and virtual queues). The proposed upper bound (whose
derivations are omitted due to the lack of space), built using [27, Eq.
($4.46$),(4.47)] reads as follows:

$\displaystyle\Delta_{p}(\mathbf{\Lambda}(t))\leq
B+\mathbb{E}\bigg{\\{}\sum_{k\in\mathcal{K}}\bigg{(}(Q_{k}^{r}(t)-Q_{k}^{l}(t))\frac{\tau
R_{k}(t)}{n_{k,b}(c_{k}(t))}$ $\displaystyle-Q_{k}^{r}(t)\tau
f_{k}(t)/w_{k}-Z_{k}(t)(\mathcal{G}_{k}(c_{k}(t))-\mathcal{G}_{k}^{\textrm{th}})\bigg{)}\bigg{|}\mathbf{\Lambda}(t)\bigg{\\}},$
(9)

where $B>0$ is a finite constant, whose expression is omitted due to the lack
of space. Minimizing (4) in a per-slot basis leads to two sub-problems
involving communication and computation variables, respectively, whose
solution is illustrated in the next paragraphs.

###  4.1 Computation sub-problem

As already mentioned in the introduction, the computation sub-problem is the
same as [26, section III.B, problem $(10)$], i.e. it does not represent a
novelty in this work. The formulation involves the variables
$\\{f_{k}\\}_{k\in\mathcal{K}}$, and its solution consists on iteratively
assigning computing resources to the users with the highest computation buffer
load. The novelty consists in the formulation and solution of the
communication sub-problem that follows.

###  4.2 Communication sub-problem

The communication sub-problem involves: i) users’ covariance matrices
$\mathbf{F}_{k},\forall k\in\mathcal{K}$, ii) RIS reflectivity parameters
$\phi_{i},i=1,\ldots,M$ and iii) the compression scheme $c_{k}$. The problem
is formulated as follows (the time index $t$ is omitted to ease the notation):

$\displaystyle\underset{\\{\mathbf{F}_{k}\\}_{k},\\{\phi\\}_{i},\\{c_{k}\\}_{k}}{\min}\;\;\mathcal{J}$
(10) subject to $\displaystyle(a)\;\mathbf{F}_{k}\succcurlyeq 0,\;\forall
k\in\mathcal{K}\quad(b)\;\textrm{Tr}(\mathbf{F}_{k}(t))\leq P_{k},\;\forall
k\in\mathcal{K}$ $\displaystyle(c)\;c_{k}\in\mathcal{C}_{k},\;\forall
k\in\mathcal{K}\quad(d)\;\;\phi_{i}\in[0,2\pi],\;i=1,\ldots M$

where
$\mathcal{J}=\sum_{k\in\mathcal{K}}\bigg{[}(Q_{k}^{r}-Q_{k}^{l})\frac{\tau
R_{k}}{n_{k,b}(c_{k})}+V\textrm{Tr}(\mathbf{F}_{k})-Z_{k}\mathcal{G}_{k}(c_{k})\bigg{]}.$
Despite the dramatic complexity reduction with respect to the original problem
(5), (10) is challenging for two main reasons: i) it is a mixed-integer non-
convex program, and ii) the function $\mathcal{G}_{k}(c_{k})$ is generally
unknown. However, once $\\{\phi_{i}\\}_{\forall i}$ and $\\{c_{k}\\}_{\forall
k}$ are fixed, the problem boils down to a convex problem that can be
efficiently solved through a water-filling procedure that is presented in
[26]. Therefore, we propose to first select RIS parameters and compression
schemes through a DRL-based algorithm, to get rid of the complexity introduced
by discrete non-convex function and lack of model. Then, once RIS parameters
and compression schemes have been set, we solve the remaining part, which only
involves the transmit covariance matrices, using a model-based approach and
classical tools from convex optimization [28].

####  4.2.1 DRL-based selection of RIS reflectivity and compression

Reinforcement Learning aims to learn a parameterized policy function (i.e. a
neural network) that maps from environmental observations to available
actions, so that the expected cumulative reward signal $r(t)$ is maximized.
For the communication sub-problem, the agent being trained at time $t$
observes a vector containing the system’s past state variables and current
channels in vectorized forms:

$\displaystyle\mathbf{s}(t)=$
$\displaystyle\Large[\\{Q_{k}^{l}(t-1)\\},\\{Q_{k}^{r}(t-1)\\},\\{Z_{k}(t-1)\\},\\{\mathcal{G}_{k}(c_{k}(t-1))\\},$
$\displaystyle\\{f_{k}(t-1)\\},\\{\mathbf{H}_{k,d}(t)\\},\\{\mathbf{H}_{r,a}(t)\\},\\{\mathbf{H}_{k,a}(t)\\}\Large]$
(11)

of dimensionality $5K+2KN_{a}N_{u}+2KMN_{u}2N_{a}M$ for which the real and
imaginary parts of the channel vectors are received separately. Based on
$\mathbf{s}(t)$, the agent selects the action
$\mathbf{a}(t)=[\\{c_{k}(t)\\}_{\forall k},\\{\phi_{m}(t)\\}_{\forall m}]$,
which fixes the current compression and RIS profiles, allowing for the optimal
values for $\\{\mathbf{F}_{k}(t)\\}_{\forall k}$ to be computed (see below)
and consequently, for the objective value of (10) to be calculated. Since the
RL problem is posed as maximization, we simply set $r(t)=-\mathcal{J}(t)$ to
attain the equivalent optimization problem of (10). The system’s state then
proceeds to $t+1$.

It is important to note that all component variables of $\mathbf{s}(t)$
appearing in (4.2.1), as well as the objective value of (10) at time $t$
depend exclusively on the values of $\mathbf{s}(t-1)$ and $\mathbf{a}(t)$.
This remark ensures that the system evolves as a Markov Decision Process (MDP)
[20], allowing for theoretically optimal policies to be derived without the
need for knowledge of its evolution from $t-2$ and backward. Besides, the DRL
formulation is algorithmic-agnostic, which allows for a wide-variety of DRL
agents to be applied with different component mechanisms. Our methodology
involves the use of the well-established Proximal Policy Evaluation (PPO)
algorithm [29], although we omit its detailed description due to space
restrictions. To make PPO applicable to the system’s mechanics, we configure
its policy network to output $\mathbf{a}(t)$ as vectors of continuous values
in $[0,1]^{K+M}$ and we then construct $c_{k}(t)$ by discretizing each of the
first $K$ elements to the allowed compression levels and $\phi_{m}(t)$ by
multiplying the remaining elements by $2\pi$, as proposed in [30].

####  4.2.2 Water-filling for covariance matrix optimization

Once RIS reflectivity and compression schemes are fixed, the problem to find
the optimal covariance matrix can be solved optimally. Indeed, first of all,
once the RIS parameters are fixed, the problem can be decoupled across
different users. Then, for all users
$k\in\widetilde{\mathcal{K}}=\\{\widetilde{k}\in\mathcal{K}:Q_{\widetilde{k}}^{l}\leq
Q_{\widetilde{k}}^{r}\\}$, the optimal solution is
$\mathbf{F}=\mathbf{0}_{N_{u}}$, a the first two terms in (10) are non-
decreasing functions of the user transmit power, as in [26, 12]. Also, for all
users $k\in\mathcal{K}\setminus\widetilde{\mathcal{K}}$, the problem is convex
and can be solved through the water-filling procedure described in [26,
Algorithm 2]. The overall solution seeks to balance between data rate
(weighted by physical queues, i.e. delay), accuracy (weighted by virtual
queues), and objective function (weighted by the trade-off parameter $V$).

##  5 Numerical Evaluation

System parameters. To illustrate the effectiveness our proposed methodology we
conceive an edge inference scenario of image recognition, in the presence of
an RIS, and an AP. A ResNet32 [31] model is deployed on a MEH machine and
trained to classify CIFAR-10 images, after being trained to approximately
$92\%$ accuracy (over the uncompressed images), offering a reasonable
compromise between performance and E2E delay. For the shake of clarity of the
results, we consider a single mobile user with no direct wireless link to the
AP. The UE performs image compression on device using the JPEG protocol, so
that $c_{k}$ values correspond to the compression quality. In terms of
inference reliability metric, we consider the average accuracy in order to
provide interpretable thresholds and performance evaluations. While indeed,
the accuracy of the network in unlabelled images cannot be known by
definition, we have numerically quantified its average metric values for all
possible JPEG compression levels of dataset’s test set, to be used for
performance evaluation. The rest of the parameter values are presented next:
Classifier parameters: $4.7\cdot 10^{5}$, $w_{k}=5.6$GHz, $\bar{A_{k}}=4$
arrivals per time step, $c_{k}\in\\{1,\dots,100\\}$, $f_{\text{max}}=3.6$ GHz,
$P_{k}=100$ mW, $\tau=0.01$ s, $W=100$ MHz,
$\mathcal{G}_{k,\textrm{th}}=0.85$, Rice factor: $25$ dB, $\sigma^{2}=-120$
dB, operating frequency: $5$ GHz, UE-RIS attenuation: $62.60$ dB, AP-RIS
attenuation: $66.34$ dB, Maximum user movement displacement: $5$ m. A
frequency-division multiplexing transmission scheme is employed, where the
power allocation among the narrowband frequency bins is kept equal and is not
treated as part of the optimization problem. This assumption is necessary to
integrate a realistic transmission paradigm without increasing the complexity
of the problem.

Fig. 1: Evaluation of the achieved system objectives for different trade-off
values of $V$.

Evaluation: To solve the problem (5), we jointly apply the partial
methodologies presented in Section 3. The water-filling and CPU scheduling
optimization routines solve their respective sub-problems optimally at every
time step. Therefore, to assess the overall formulation we numerically
evaluate performance of the DRL’s selections of RIS reflections and
compression schemes. As baselines, we introduce the policies of always
performing (i) maximum compression (which offers the minimum delay), (ii) no
compression (offering maximum accuracy) and (iii) compression at a uniformly
random level, while employing the rest of the optimization components and
controlling the RIS reflections at random. Upon fixing a trade-off value for
$V$ from $\\{1\cdot 10^{5},2\cdot 10^{5},3\cdot 10^{5},4\cdot 10^{5},5\cdot
10^{5},3\cdot 10^{6}\\}$, we train a PPO instance with $5$ layers of $32$
neurons for its policy and value networks for $10^{6}$ time steps, resetting
the system at randomized episodes of length $1500$, using the final episode as
evaluation.

The training phases resulted in objective values orders of magnitudes lower
than the baselines, although the scales of each objective instance is heavily
influenced by the choice of $V$. The resulting accuracy offered by PPO ranges
from $0.85$ to $0.91$ denoting that the DRL component always satisfies the
desired constraint. The maximum and random compression schemes offer accuracy
scores of $0.20$ and $0.69$, respectively. Fig. 1 presents the achieved delay
and mean transmit power offered by each method. Compared to the no-compression
approach, DRL can achieve up to $47\%$ reduction in the average UE power
consumption and up to $59\%$ reduction in the maximum E2E delay. At the same
time, it results from $45\%$ to up to $2.5$ times (in the most delay-sensitive
cases) higher power consumption and $5-8$ms added delay in comparison to the
full compression policy, but the latter falls extremely short of any practical
accuracy constraints. Clearly, the system endowed with they Lyapunov-based
techniques and DRL is able to automatically offer a balance between
performance (accuracy and E2E delay) and power consumption.

##  6 Conclusion

In this paper, the problem of accurate, fast, and low-power edge inference has
been investigated. A RIS-empowered MEC system has been proposed and a power
minimization problem under E2E delay and accuracy constraints was formulated.
The problem was solved through a combination of Lyapunov-driven optimization
and DRL tools, showing how edge AI will play the twofold role of an efficient
optimization tool and a service enabled by edge computing resources in 6G. A
numerical evaluation in a RIS-empowered wireless scenario illustrated the
capabilities of the methodology of achieving desired thresholds between
accuracy, delay, and power consumption, striking typical the trade-off of edge
AI native wireless networks.

## References

  * [1] K. B. Letaief et al.,  “Edge artificial intelligence for 6G: Vision, enabling technologies, and applications,”  IEEE J. Sel. Areas Commun., vol. 40, no. 1, pp. 5–36, 2022.
  * [2] S. Kekki et al.,  “MEC in 5G networks,” ETSI White Paper, no. 28, 2018.
  * [3] Suzhi Bi, Liang Huang, Hui Wang, and Ying-Jun Angela Zhang,  “Lyapunov-guided deep reinforcement learning for stable online computation offloading in mobile-edge computing networks,”  IEEE Trans. Wireless Commun., vol. 20, no. 11, pp. 7519–7537, 2021\.
  * [4] M. Sana et al.,  “Energy efficient edge computing: When lyapunov meets distributed reinforcement learning,”  in Proc. IEEE ICC, 2021, pp. 1–6.
  * [5] Z. Zhou et al.,  “Edge intelligence: Paving the last mile of artificial intelligence with edge computing,”  Proc. IEEE, vol. 107, no. 8, pp. 1738–1762, 2019.
  * [6] S. Han et al.,  “Learning both weights and connections for efficient neural networks,”  in Proc. Neural Information Processing Systems (NIPS) - Volume 1, Cambridge, MA, USA, 2015, NIPS’15, p. 1135–1143, MIT Press.
  * [7] B. Taylor et al.,  “Adaptive deep learning model selection on embedded systems,”  SIGPLAN Not., vol. 53, no. 6, pp. 31–43, Jun. 2018.
  * [8] J. Jiang et al.,  “Chameleon: Scalable adaptation of video analytics,”  in Proc. ACM SIGCOMM, Budapest, Hungary, 2018, p. 253–266.
  * [9] X. Ran et al.,  “DeepDecision: A mobile deep learning framework for edge video analytics,”  in Proc. IEEE INFOCOM, 2018, pp. 1421–1429.
  * [10] A. Galanopoulos et al.,  “AutoML for video analytics with edge computing,”  in Proc. IEEE INFOCOM, 2021, pp. 1–10.
  * [11] M. Merluzzi et al.,  “Wireless edge machine learning: Resource allocation and trade-offs,”  IEEE Access, vol. 9, pp. 45377–45398, 2021.
  * [12] M. Merluzzi et al.,  “Energy-efficient classification at the wireless edge with reliability guarantees,”  in Proc. IEEE ICC, 2022, pp. 109–114.
  * [13] M. Di Renzo et al.,  “Smart radio environments empowered by reconfigurable AI meta-surfaces: An idea whose time has come,”  EURASIP J. Wireless Commun. Netw., vol. 2019, no. 1, pp. 1–20, 2019\.
  * [14] G. C. Alexandropoulos et al.,  “Reconfigurable intelligent surfaces and metamaterials: The potential of wave propagation control for 6G wireless communications,”  IEEE ComSoc TCCN Newslett., vol. 6, no. 1, pp. 25–37, Jun. 2020\.
  * [15] E. Calvanese Strinati et al.,  “Wireless environment as a service enabled by reconfigurable intelligent surfaces: The rise-6G perspective,”  in Proc. EuCNC & 6G Summit, Porto, Portugal, 2021, pp. 562–567.
  * [16] E. Calvanese Strinati et al.,  “Reconfigurable, intelligent, and sustainable wireless environments for 6g smart connectivity,”  IEEE Commun. Mag., vol. 59, no. 10, pp. 99–105, 2021.
  * [17] M. Jian et al.,  “Reconfigurable intelligent surfaces for wireless communications: Overview of hardware designs, channel models, and estimation techniques,”  Intell. Converged Netw., vol. 3, no. 1, pp. 1–32, 2022.
  * [18] C. Huang et al.,  “Reconfigurable intelligent surfaces for energy efficiency in wireless communication,”  IEEE Trans. Wireless Commun., vol. 18, no. 8, pp. 4157–4170, 2019\.
  * [19] Q. Wu and R. Zhang,  “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,”  IEEE Trans. Wireless Commun., vol. 18, no. 11, pp. 5394–5409, 2019\.
  * [20] G. C. Alexandropoulos et al.,  “Pervasive machine learning for smart radio environments enabled by reconfigurable intelligent surfaces,”  Proc. IEEE, vol. 110, no. 9, pp. 1494–1525, 2022.
  * [21] K. Stylianopoulos et al.,  “Deep contextual bandits for orchestrating multi-user MISO systems with multiple RISs,”  in Proc. IEEE ICC, 2022, pp. 1556–1561.
  * [22] T. Bai et al.,  “Latency minimization for intelligent reflecting surface aided mobile edge computing,”  IEEE J. Sel. Areas Commun., vol. 38, no. 11, pp. 2666–2682, 2020\.
  * [23] P. Di Lorenzo et al.,  “Dynamic edge computing empowered by reconfigurable intelligent surfaces,”  EURASIP J. Wireless Commun. Netw., vol. 2022, no. 122, Dec 2022\.
  * [24] S. Huang et al.,  “Reconfigurable intelligent surface assisted mobile edge computing with heterogeneous learning tasks,”  IEEE Trans. Cognitive Commun. Netw., vol. 7, no. 2, pp. 369–382, Jun. 2021.
  * [25] C. Battiloro et al.,  “Lyapunov-based optimization of edge resources for energy-efficient adaptive federated learning,”  IEEE Trans. Green Commun. Netw., 2022.
  * [26] F. E. Airod et al.,  “Reconfigurable intelligent surface aided mobile edge computing over intermittent mmwave links,”  in Proc. IEEE SPAWC, 2022, pp. 1–5.
  * [27] M. J. Neely,  Stochastic Network Optimization with Application to Communication and Queuing Systems,  Morgan and Claypool, 2010.
  * [28] S. Boyd and L. Vandenberghe,  Convex Optimization,  Cambridge university press, 2004.
  * [29] J. Schulman et al.,  “Proximal policy optimization algorithms,”  arXiv [cs.LG/1707.06347], 2017.
  * [30] K. Stylianopoulos and G. C. Alexandropoulos,  “Online RIS configuration learning for arbitrary large numbers of 1-bit phase resolution elements,”  in Proc. IEEE SPAWC, 2022, pp. 1–5.
  * [31] K. He et al.,  “Deep residual learning for image recognition,”  in Conference on Computer Vision and Pattern Recognition (CVPR), Jun. 2016.