# Cosmic Reionization On Computers: Statistics, Physical Properties and
Environment of Lyman Limit Systems at $z\sim 6$

Jiawen Fan（樊稼问  Department of Physics; University of Michigan, Ann Arbor, MI
48109, USA  Hanjue Zhu (朱涵珏) Department of Astronomy & Astrophysics; The
University of Chicago; Chicago, IL 60637, USA  Camille Avestruz Department of
Physics; University of Michigan, Ann Arbor, MI 48109, USA  Leinweber Center
for Theoretical Physics; Department of Physics; University of Michigan, Ann
Arbor, MI 48109, USA  Nickolay Y. Gnedin  Theoretical Physics Division; Fermi
National Accelerator Laboratory; Batavia, IL 60510, USA  Kavli Institute for
Cosmological Physics; The University of Chicago; Chicago, IL 60637, USA
Department of Astronomy & Astrophysics; The University of Chicago; Chicago, IL
60637, USA

###### Abstract

Lyman Limit Systems (LLSs) are dense hydrogen clouds with high enough HI
column densities ($N_{\rm HI}$) to absorb Lyman continuum photons emitted from
distant quasars. Their high column densities imply an origin in dense
environments; however, the statistics and distribution of LLSs at high
redshifts still remain uncertain. In this paper, we use self-consistent
radiative transfer cosmological simulations from the “Cosmic Reionization On
Computers” (CROC) project to study the physical properties of LLSs at the tail
end of cosmic reionization at $z\sim 6$. We generate 3000 synthetic quasar
sightlines to obtain a large number of LLS samples in the simulations. In
addition, with the high physical fidelity and resolution of CROC, we are able
to quantify the association between these LLS samples and nearby galaxies. Our
results show that a higher fraction of higher column density LLSs are
spatially associated with nearby galaxies. Moreover, we find that LLSs that
are not near any galaxy typically reside in filamentary structures connecting
neighboring galaxies in the intergalactic medium (IGM). This quantification of
the distribution and associations of LLSs to large scale structures informs
our understanding of the IGM-galaxy connection during the Epoch of
Reionization, and provides a theoretical basis for interpreting future
observations.

galaxies — methods, numerical — cosmology

††journal: ApJ

##  1 Introduction

The recent launch of JWST opens up a new window in studying galaxy formation
during the Epoch of Reionization. Recent and forthcoming data from JWST’s
near-infrared and mid-infrared bands enable observers to study the earliest
galaxies and to probe galaxy formation to fainter magnitudes. As the
observations push beyond $z\sim 6$ (Crighton et al., 2019; Adams et al.,
2022), however, they lose an important window into the galaxy formation - the
Intergalactic Medium (IGM) as probed by the Ly-$\alpha$ forest absorption in
the spectra of distant quasars. The class of Ly-$\alpha$ absorbers that is the
hardest to study at the highest redshifts is Lyman limit systems (LLSs). The
nature of these systems is still not fully understood, and the observational
challenges in observing them at the highest redshifts places an undue burden
on theoretical modeling.

As the first stars and galaxies form, photons emitted from these luminous
sources travel through the IGM. Emitted photons encounter hydrogen atoms
residing in the diffuse IGM and in LLSs. These photons rapidly ionize the
neutral hydrogen in the diffuse IGM, but LLSs, which presumably trace high-
density regions, shield a significant fraction of enclosed neutral hydrogen.
As the ionized regions grow and merge during the process of reionization, the
enclosed LLSs alter the propagation of ionizing photons inside ionized
regions, thereby changing the Mean Free Path (MFP) of the Lyman continuum
(LyC) photons. LLSs therefore affect reionization processes and better
understanding of LLSs provide further insights in the physics of reionization.

LLSs are usually defined by the value of the hydrogen column density along the
line-of-sight. By a formal definition LLSs have column densities greater
$N_{\rm HI}>1.6\times 10^{17}\rm cm^{-2}$ (this corresponds to an optical
depth of 1) and less than $N_{\rm HI}<2\times 10^{20}\rm cm^{-2}$ (with higher
column density systems called damped Ly$\alpha$ systems), although of course
this definition is rather arbitrary and does not correspond to any physical
characterization of absorbers. LLSs contribute more to the ionizing photon
mean free path than Ly$\alpha$ forest (systems with $N_{\rm HI}\lesssim
1.6\times 10^{17}\rm cm^{-2}$), and they are much more frequently seen than
the denser but more rare damped Ly$\alpha$ systems. While many low-z
observations of the LLSs exist in the literature Prochaska et al. (2010);
Ribaudo et al. (2011); O’Meara et al. (2013), far fewer observations exist at
high-z. Observational measurements of the statistics and properties of LLSs at
high-z require more quasar spectra, particularly with complete spectral
coverage of both the Ly$\alpha$ and the Lyman break (Shukla et al., 2016).

With numerical simulations of high physical fidelity, especially with the
advancement of radiative transfer implemented in these simulations, we can now
study the LLSs at high-z to an unprecedented level. A similar study from
Kohler & Gnedin (2007) investigated LLS at z = 4 with two sets of simulations
with box length $L_{box}$ = 4 $h^{-1}$cMpc, and $L_{box}$ = 8 $h^{-1}$cMpc,
but these simulations did not include the self-consistent implementation of
galaxy formation physics accounted for in the work presented here.

Our goal in this paper is to study the LLSs at $z\sim 6$ in the Cosmic
Reionization on Computers (CROC) simulations. The structure of this paper is
as follows. We describe the CROC simulations and our methods to identify and
quantify LLSs in Section 2, present the statistics and properties of LLSs in
Section 3, and summarize our results in Section 4.

##  2 Methodology

###  2.1 “Cosmic Reionization on Computers”: CROC Simulation

CROC simulations use the Adaptive Refinement Tree code (Kravtsov et al., 1997,
2002; Rudd et al., 2008). To account for physics necessary to self-
consistently model cosmic reionization, CROC simulations include gravity, gas
dynamics, fully coupled (to gas dynamics) radiative transfer, atomic cooling
and heating processes, molecular hydrogen formation, star formation, and
stellar feedback. Full details of the simulation are described in the CROC
method paper (Gnedin, 2014).

For our analysis, we use three simulation boxes with box length $L_{box}$ = 40
$h^{-1}$cMpc. These simulation boxes have a maximum spatial resolution of 100
pc in proper units. We extract 1000 line-of-sight samples in each of the boxes
that correspond to independent random realizations. The line-of-sights are
sampled non-uniformly, fully tracing the underlying high resolution of the
adaptively refined grid. Each line-of-sight sample is 100 $h^{-1}$cMpc in
length, which allows us to find even rare absorbers.

###  2.2 Identifying Lyman Limit Systems

The observational signature of an LLS is a break at the Lyman limit (a drop in
transmitted flux at 912Å rest frame). In numerical simulations, we can simply
look for regions with high enough neutral fraction ($X_{\rm HI}$). Figure 1
shows the neutral fraction along one line-of-sight in our simulation. Along
each line of sight we search for regions with $X_{\rm HI}$ values larger than
a given threshold, e.g. the region indicated by the red rectangle in the lower
panel of Figure 1. We define the edges of these highly neutral regions to be
at the intersection points between the threshold and the line-of-sight, e.g.
the left and right sides of the red rectangle, which also corresponds to the
“piece” of the LLS identified.

We use $X_{\rm HI}=0.001$ as our fiducial threshold for highly neutral
regions. This choice is somewhat arbitrary; we tested other $X_{\rm HI}$
values ($X_{\rm HI}=\\{0.0001,0.001,0.003\\}$), which all lead to similar
distributions of column densities in the LLS range, though such a method for
determining the HI column density would not work reliably in the Ly$\alpha$
forest. We illustrate this comparison in Figure 2 and further describe the
distribution of column densities in comparisons across different CROC boxes
and also with observations in Section 3.1. The primary takeaway from Figure 2
is that the number of LLSs with $N_{\rm HI}>3\times 10^{17}\rm cm^{-2}$ does
not significantly vary with threshold choice.

Figure 1: Neutral fraction along an example line of sight. The top panel shows
two identified “highly neutral” regions in a single line-of-sight (red peaks),
and the bottom panel zooms into the second peak with the red rectangle
indicating the region of the line-of-sight above a threshold neutral fraction
of $X_{\rm HI}=0.001$ used to identify such regions. Figure 2: Distribution of
the number of LLSs at a given column densities, where we identify the LLSs at
different threshold neutral fractions, $X_{\rm HI}=\\{0.0001,0.001,0.003\\}$.
The distribution does not significantly vary with threshold choice for $N_{\rm
HI}>3\times 10^{17}{\rm cm}^{-2}$.

###  2.3 Galaxy Associations and Visualization of Lyman Limit Systems

A primary goal of this paper is to have a clearer understanding of where LLSs
are located with respect to surrounding galaxies. We use the Rockstar halo
finder (Behroozi et al., 2013) to identify dark matter halos in our
simulations. We use the identified center of the particles’ mass distribution
inside the dark matter halos as a proxy for the location of each hosted
galaxy, enabling spatial associations with LLSs.

We define the “association” between LLSs and nearby galaxies using both the
spatial locations and virial radii ($R_{vir}$) of the dark matter halos
hosting these galaxies. We quantify the strength of associations using the
quantity $f_{n}(R_{vir})$, where $n$ is an integer. With $n=1$,
$f_{1}(R_{vir})$ is the fraction of LLSs whose neutral fraction peaks are
inside the $R_{vir}$ of any galaxy. With $n=2$, $f_{2}(R_{vir})$ is the
fraction of LLSs with neutral fraction peaks within 2 $\times R_{vir}$ of some
galaxy. Higher $n$ values indicate associations at larger distances with
respect to $R_{vir}$.

We also generate gas density projection plots around every LLS in the
simulation box snapshots using the “yt” simulation analysis package (Turk et
al., 2011). In these projection plots, one can visually identify structures
surrounding the LLS including galaxies and connective filaments.

##  3 Results

###  3.1 Statistics of LLSs

An important quantity that observers have studied extensively at low-z is the
column density distribution of LLSs. Figure 3 shows the column density
distribution of LLSs. We compute the y-axis quantity by taking the LLSs number
count in each column density bin ($\Delta N_{\rm HI}$) normalized by the total
path length sampled in the simulation ($\Delta l$) and the bin size ($\Delta
N_{HI}$). We then multiply this quantity by the column density value in each
bin to reduce the dynamic range for comparison, which is equivalent to
measuring the column density distribution per log of neutral hydrogen column
density.

The left panel shows the individual distributions from each of the three
different realizations of the CROC simulations to account for the cosmic
variance. Here, we see that the column density distributions of LLSs in CROC
exhibit an overall trend where the LLS number density decreases with
increasing column density. Overall, the three simulation boxes produce similar
column density distributions.

The right hand panel shows the same column density distribution of LLSs from
the combined three CROC boxes, alongside fits to observational data from
Crighton et al. (2019) at $z=2.4$ and $z=4.4$. We see that our simulation
result and the observational fits have a relatively similar shape, with CROC
having more LLSs with $16<\rm log_{10}(N_{\rm HI})<19$. These differences are
acceptable, as simulations are not expected to be at the level of fidelity to
match all the data perfectly. It is sufficient for our purpose that they
capture the dominant trends in the data, such as a depression at $N_{\rm
HI}\sim 10^{18}{\rm cm}^{-2}$.

Historically, the observational data has been reported per unit ”absorption
distance” $\Delta X$; it is mostly a historical artifact and is defined as
follows:

$dX=\frac{H_{0}}{H(z)}(1+z)^{2}dz,$ (1)

whereas we measure the path length $\Delta l$ in the simulation in the
comoving space. The conversion between $\Delta l$ and $\Delta X$ is:

$dl/dX=\frac{(1+z)^{2}H_{0}}{c}.$ (2)

We apply this conversion to the Crighton et al. (2019) observational fits to
enable comparison between observations and simulations.

Figure 3: Left: Column density distribution of LLSs for three different 40
$h^{-1}$cMpc boxes. Right: Comparison between the observational fit at
$z=2.4$, $z=4.4$ (Crighton et al., 2019) and our simulations.

###  3.2 Association of LLSs with Galaxies

The association between LLSs and dense gas structures surrounding galaxies is
well established neither at low nor at high redshift. A primary goal of this
paper is to provide theoretical predictions of how LLSs at $z\sim 6$, near the
end of cosmic reionization, are associated with galaxies and nearby
structures.

Figure 4: The fraction of Lyman Limit Systems within $n\times R_{vir}$ of any
galaxies, with $n\in[1,2,3,4]$. Blue, orange, green, and red indicate
respective fraction of Lyman Limit Systems within that distance of any galaxy.
We see that higher column density systems have progressively stronger
associations with galaxies.

In order to quantify the association of Lyman Limit Systems and galaxies in
our simulation, we show in Figure 4 the fraction of LLSs whose neutral
fraction peaks are within $n\times R_{vir}$ of any galaxy, i.e. the
$f_{n}(R_{vir})$ quantity we discussed in Section 2.3 with
$n\in\\{1,2,3,4\\}$. The fraction is for all LLSs from each of our 40
$h^{-1}$cMpc simulation boxes. The overall trend with the neutral hydrogen
column density is very clear: the damped Ly$\alpha$ systems of $\rm
log_{10}\it N_{HI}>$20 are predominantly located inside galactic halos. At
lower column densities the association of LLSs with galaxies weakens, with a
larger fraction of these LLSs located further than $2\times R_{vir}$ from any
galaxy. As we go down to the Ly$\alpha$ forest range, $N_{\rm HI}<10^{17}{\rm
cm}^{-2}$, only 10% of these systems are within $1\times R_{vir}$ of a galaxy
center and only 40% are within $2\times R_{vir}$.

Of the population of highest column density LLSs (i.e. those with $10^{19}{\rm
cm}^{-2}<N_{\rm HI}<10^{20}{\rm cm}^{-2}$), 80% are within $2\times R_{vir}$
of some galaxy. Figure 5 shows the projected neutral hydrogen density around a
LLS with the column density of $N_{HI}=10^{19.4}{\rm cm}^{-2}$, located within
$1\times R_{vir}$ of the virial radius of a galaxy. The left panel shows the
projection in the y-z plane, and the right in the x-z plane. The depth of
projection corresponds to the estimated size of the LLS ($\sim$10 kpc).

We indicate the LLS location with a green star, and surrounding dark matter
halo centers with their virial radii respectively in black crosses and black
circles. The galaxy associated with the LLS has a dark matter halo mass and
virial radius of $M=2.8\times 10^{9}{\rm\,M}_{\odot}$ and $R_{vir}=6.8\rm\
kpc$ respectively; we show it in red.

The remaining 20% of the highest column density LLSs are outside of $2\times
R_{vir}$ of any galaxy. While these are not closely associated with any
galaxy, we find that they are associated with overdense filamentary structures
in the large-scale surrounding environments. Figure 6 shows the projected
neutral hydrogen map surrounding such a system, with similar projections and
color mapping as Figure 5. Here, we see that the LLS sits in a filamentary
bridge connecting galaxies, but lies well outside twice the virial radius of
nearby galaxies.

Figure 5: Projected neutral hydrogen maps centered around a Lyman Limit system
with $N_{HI}=10^{19.4}{\rm cm}^{-2}$. The projection depth of the graph is
$\sim$ 10 kpc, corresponding to the characteristic size of the Lyman Limit
system. Left: zy-projection plane. Right: xz-projection plane. The green star
indicates the position of the Lyman Limit system; the black circles center on
galaxies and indicate their corresponding dark matter halo virial radii; the
red circle corresponds to the galaxy that contains the Lyman Limit system
within its virial radius. Most of the high column density Lyman limit systems
reside within one virial radius of a galaxy.

Figure 6: Projected neutral hydrogen maps centered around an example Lyman
Limit system weakly associated with any galaxy. The projection depth is $\sim$
10 kpc. Markers indicate the same properties as in Figure 5. This Lyman Limit
system is well outside of 2 times the virial radius of any galaxy. For
illustration purposes, we annotate galaxies up to 5 virial radii away from the
Lyman Limit system.

###  3.3 Physical Properties of LLSs

Figure 7: Lyman Limit system optical depth vs various physical properties of
the Lyman Limit systems in CROC. From left to right: characteristic size,
neutral fraction, and gas overdensity. The top panels are hexbin plots to more
clearly visualize the distribution of data points with the color bar
indicating the number of data points in each hexbin. The bottom panels show
individual data points colored by their distance to the nearby galaxies at
fixed multiples of the galaxy virial radius. The color coding is the same as
in Figure 4: blue, orange, green, and red correspond to systems within 1, 2,
3, and 4 times the virial radius of any galaxy.

In this section, we investigate the relationship between the optical depth of
LLSs and their various physical properties in CROC simulations. The LLS
physical properties of interest here are: optical depth, characteristic size,
neutral fraction, and density.

Optical depth describes opacity of the LLSs and is directly proportional to
the column density. Optical depth is defined as:

$\tau=\sigma_{\rm LL}N_{\rm HI},$ (3)

where $\sigma_{\mathrm{LL}}=6.3\times 10^{-18}\ \rm cm^{2}$ is the hydrogen
photoionization cross section at the Lyman limit.

We define the characteristic size of a LLS as:

$\L_{\rm eff}=N_{\rm HI}/n_{\rm HI},$ (4)

where $n_{\rm HI}$ is the peak HI number density of the LLS. If the peak had a
Gaussian shape, that definition would closely correspond to the FWHM.

The top panel of Figure 7 illustrates the distribution of simulated Lyman
Limit Systems in these four quantities. We use hexbin plots to clearly
represent over-dense data points without overlapping them. The color bar on
top of each plot represents the number of data points inside each hexagon. The
horizontal dashed lines show the $\tau$ range corresponding to the LLS column
density range, although we note that the numbers are arbitrarily defined.

We see from the middle and right panels that $\tau$ increases with the peak
neutral fraction and the peak gas density of the LLSs. This means that higher
column density LLSs are more neutral and denser at the center. However, there
is not much of a trend between $\tau$ and $\L_{\rm eff}$ \- i.e., at a fixed
column density (or, equivalently, $\tau$), some LLSs are compact dense gas
clumps while others are more spatially extended but less dense. Most LLSs,
however, have sizes between 1 and $\sim$5 physical kpc. Note, that for LLSs
with $\tau\geq$ 100, their neutral fractions are close to unity. We also know
that these high column density LLSs are strongly associated with nearby
galaxies from Figure 4. This means that despite their proximity to ionizing
sources, high column density LLSs still have highly neutral centers that are
not affected by the ionizing photons.

It is also instructive to separate LLSs based on the strength of their
association to galaxies, and plot the physical quantities in the bottom panel
of Figure 7. The colors we use here are the same as in Figure 4. It is clear
that the strength of association with galaxies does not significantly affect
the physical properties of LLSs at $z\sim 6$.

##  4 Summary and Discussion

We use the state-of-the-art reionization simulation, CROC, to investigate the
statistics, environment and physical properties of lyman limit systems at
$z\sim 6$. We cast a total of 3000 lines-of-sight to simulate synthetic quasar
sight lines. Along these sight lines, we are able to identify Lyman Limit
systems in the simulated IGM.

Our main results are:

  1. 1.

The highest column density self-shielded absorbers, such as damped Ly$\alpha$
systems, have the strongest association with galaxies. These are typically
within the virial radius of a galaxy.

  2. 2.

Galaxy distribution largely dictates where we will find most of the LLSs; most
LLSs are within 1-2 virial radii of a galaxy.

  3. 3.

LLSs that are not located close to any galaxy are preferentially located along
the filament connecting more distant galaxies.

  4. 4.

The column density of LLSs strongly correlates with gas density and the
hydrogen neutral fraction. Their characteristic sizes, however, show little
correlation with other physical properties.

In future work, we plan to examine the main free path of LyC photons and
relate it to the distributions of LLSs. With such future studies, we can
develop a better understanding of what role the LLSs play in shaping the
ionization history during the epoch of reionization.

JF was supported by the Summer Undergraduate Research Experience fund from the
University of Michigan’s Physics Department. CA acknowledges support from the
Leinweber Center for Theoretical Physics. Fermilab is operated by Fermi
Research Alliance, LLC, under Contract No. DE-AC02-07CH11359 with the United
States Department of Energy. This work used resources of the Argonne
Leadership Computing Facility, which is a DOE Office of Science User Facility
supported under Contract DE-AC02-06CH11357. An award of computer time was
provided by the Innovative and Novel Computational Impact on Theory and
Experiment (INCITE) program. This research is also part of the Blue Waters
sustained-petascale computing project, which is supported by the National
Science Foundation (awards OCI-0725070 and ACI-1238993) and the state of
Illinois. Blue Waters is a joint effort of the University of Illinois at
Urbana-Champaign and its National Center for Supercomputing Applications.

## References

  * Adams et al. (2022) Adams, N. J., Conselice, C. J., Ferreira, L., et al. 2022, arXiv e-prints, arXiv:2207.11217
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, ApJ, 762, 109
  * Crighton et al. (2019) Crighton, N. H. M., Prochaska, J. X., Murphy, M. T., et al. 2019, MNRAS, 482, 1456
  * Gnedin (2014) Gnedin, N. Y. 2014, ApJ, 793, 29
  * Kohler & Gnedin (2007) Kohler, K., & Gnedin, N. Y. 2007, ApJ, 655, 685
  * Kravtsov et al. (2002) Kravtsov, A. V., Klypin, A., & Hoffman, Y. 2002, ApJ, 571, 563
  * Kravtsov et al. (1997) Kravtsov, A. V., Klypin, A. A., & Khokhlov, A. M. 1997, ApJS, 111, 73
  * O’Meara et al. (2013) O’Meara, J. M., Prochaska, J. X., Worseck, G., Chen, H.-W., & Madau, P. 2013, ApJ, 765, 137
  * Prochaska et al. (2010) Prochaska, J. X., O’Meara, J. M., & Worseck, G. 2010, ApJ, 718, 392
  * Ribaudo et al. (2011) Ribaudo, J., Lehner, N., Howk, J. C., et al. 2011, ApJ, 743, 207
  * Rudd et al. (2008) Rudd, D. H., Zentner, A. R., & Kravtsov, A. V. 2008, ApJ, 672, 19
  * Shukla et al. (2016) Shukla, H., Mellema, G., Iliev, I. T., & Shapiro, P. R. 2016, MNRAS, 458, 135
  * Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9


ifaamas [AAMAS ’23]The 14th Workshop on Optimization and Learning in
Multiagent Systems (OptLearnMAS-23) AT the 22nd International Conference on
Autonomous Agents and Multiagent Systems (AAMAS 2023)May 29 – June 2, 2023
London, United KingdomA. Ricci, W. Yeoh, N. Agmon, B. An (eds.) 2023 2023 ???
Eindhoven University of Technology Eindhoven The Netherlands Eindhoven
University of Technology Eindhoven The Netherlands Eindhoven University of
Technology Eindhoven The Netherlands Eindhoven University of Technology
Eindhoven The Netherlands

# Coordinating Fully-Cooperative Agents Using Hierarchical Learning
Anticipation

Ariyan Bighashdel<EMAIL_ADDRESS>,  Daan de Geus<EMAIL_ADDRESS>,
Pavol Jancura<EMAIL_ADDRESS>and  Gijs Dubbelman<EMAIL_ADDRESS>
###### Abstract.

Learning anticipation is a reasoning paradigm in multi-agent reinforcement
learning, where agents, during learning, consider the anticipated learning of
other agents. There has been substantial research into the role of learning
anticipation in improving cooperation among self-interested agents in general-
sum games. Two primary examples are Learning with Opponent-Learning Awareness
(LOLA), which anticipates and shapes the opponent’s learning process to ensure
cooperation among self-interested agents in various games such as iterated
prisoner’s dilemma, and Look-Ahead (LA), which uses learning anticipation to
guarantee convergence in games with cyclic behaviors. So far, the
effectiveness of applying learning anticipation to fully-cooperative games has
not been explored. In this study, we aim to research the influence of learning
anticipation on coordination among common-interested agents. We first
illustrate that both LOLA and LA, when applied to fully-cooperative games,
degrade coordination among agents, causing worst-case outcomes. Subsequently,
to overcome this miscoordination behavior, we propose Hierarchical Learning
Anticipation (HLA), where agents anticipate the learning of other agents in a
hierarchical fashion. Specifically, HLA assigns agents to several hierarchy
levels to properly regulate their reasonings. Our theoretical and empirical
findings confirm that HLA can significantly improve coordination among common-
interested agents in fully-cooperative normal-form games. With HLA, to the
best of our knowledge, we are the first to unlock the benefits of learning
anticipation for fully-cooperative games.

###### Key words and phrases:

Multi-agent reinforcement learning, Learning anticipation, Hierarchical
reasoning, Fully-cooperative games

##  1\. Introduction

One of the key characteristics of learning in multi-agent systems is the non-
stationary environment. As a result, agents should continuously interact with
each other and adapt their strategies accordingly. However, in various game
settings, these interactions commonly lead to worst-case outcomes for all
agents (Foerster et al., 2018a). A prominent example is the case of general-
sum games, particularly the Iterated Prisoner’s Dilemma (IPD), where self-
interested agents typically converge to defect-defect, which is the worst
result globally. Recently, the learning anticipation paradigm, where agents
take into account the anticipated learning of other agents, has been broadly
employed to avoid such catastrophic outcomes Foerster et al. (2018a); Letcher
et al. (2019); Zhang and Lesser (2010).

For instance, the Learning with Opponent-Learning Awareness (LOLA) method
(Foerster et al., 2018a) has proven to be successful in the IPD game.
Specifically, in LOLA, agents anticipate and shape the learning step of others
to enforce cooperation among themselves. This so-called opponent shaping in
LOLA leads to a tit-for-tat strategy in IPD (starting out cooperating and
otherwise mirroring the opponent’s last move), which achieves mutual
cooperation when chosen by both players (Axelrod and Hamilton, 1981).
Subsequent works include: a) Higher-order LOLA (HOLA) Foerster et al. (2018a)
that assumes non-naïve opponents, b) Look-Ahead (LA)Zhang and Lesser (2010);
Letcher et al. (2019) that discards opponent shaping to guarantee convergence
in games with cyclic behaviors (e.g., matching pennies), c) Stable Opponent
Shaping (SOS) Letcher et al. (2019) that guarantees convergence to the fixed
points (i.e., points with zero gradients) in differentiable games, and d)
Consistent LOLA (COLA) Willi et al. (2022) which provides a closed-form
solution of infinite-order LOLA.

Nevertheless, all of these methods are proposed to improve cooperation in
games with self-interested agents. So far, no works have explored how these
methods perform when agents are fully-cooperative, i.e., common-interested.
Current investigations on fully-cooperative games are limited to the study of
Letcher et al. (Letcher et al., 2019), proving the convergence and non-
convergence of LOLA to stable and unstable fixed points, respectively.
However, we demonstrate that in fully-cooperative normal-form games with
unstable fixed points, although the agents in both LOLA and LA do not converge
to the unstable fixed points, they are subject to miscoordination and,
consequently, worse overall rewards.

Considering the above, the key research goal of this work is to research the
influence of learning anticipation on coordination among fully-cooperative
agents. We believe that as multiple agents should interact and cooperate to
achieve a common goal, learning anticipation has the potential to improve
coordination among agents, resulting in better overall rewards.

To accomplish our goal, we first theoretically prove that in a two-agent two-
action coordination game (Claus and Boutilier, 1998), both LOLA and LA have
the tendency to lead to miscoordination among common-interested agents,
causing a worse outcome. To solve this miscoordination problem and improve the
applicability of learning anticipation to fully-cooperative games, we then
propose Hierarchical Learning Anticipation (HLA), a new learning method
explicitly developed for improving coordination in games with common-
interested agents. Specifically, HLA assigns all agents to hierarchy levels,
which define the reasoning orders of agents. We theoretically prove that HLA
can avoid miscoordination in the aforementioned coordination game.
Furthermore, we empirically show that in a two-agent three-action coordination
game (Claus and Boutilier, 1998), HLA, as opposed to LOLA and LA,
significantly improves coordination among agents, leading to better overall
rewards. Finally, we discuss the shortcomings of our study and provide future
research directions.

##  2\. Background

Our work assumes a multi-agent task that is commonly described as a Markov
Game (MG) (Littman, 1994). An MG can be defined as a tuple
$(\mathcal{N},\mathcal{S},\\{\mathcal{A}_{i}\\}_{i\in\mathcal{N}},\\{\mathcal{R}_{i}\\}_{i\in\mathcal{N}},\mathcal{T},\rho,\gamma)$,
where $\mathcal{N}$ is the set of agents ($|\mathcal{N}|=n$), $\mathcal{S}$ is
the set of states, and $\mathcal{A}_{i}$ is the set of possible actions for
agent $i\in\mathcal{N}$. Agent $i$ chooses its action
$a_{i}\in\mathcal{A}_{i}$ through the policy network
$\pi_{\theta_{i}}:\mathcal{S}\times\mathcal{A}_{i}\rightarrow[0,1]$
parameterized by $\theta_{i}$ conditioning on the given state
$s\in\mathcal{S}$. Given the actions of all agents, each agent $i$ obtains a
reward $r_{i}$ according to its reward function
$\mathcal{R}_{i}:\mathcal{S}\times\mathcal{A}_{1}\times...\times\mathcal{A}_{n}\rightarrow\mathbb{R}$.
Given an initial state, the next state is produced according to the state
transition function
$\mathcal{T}:\mathcal{S}\times\mathcal{A}_{1}\times...\times\mathcal{A}_{n}\times\mathcal{S}\rightarrow[0,1]$.
Given an episode $\tau$ of horizon $T$, the discounted return for each agent
$i$ at time step $t\leq T$ is defined by
$G^{t}_{i}(\tau)=\sum_{l=t}^{T}\gamma^{l-t}r_{i}$ where $\gamma$ is a
predefined discount factor. The expected return given the agents’ policy
parameters approximates the state value function for each agent
$V_{i}(s,\theta_{1},...,\theta_{n})=\mathbb{E}[G_{i}^{t}(\tau|s^{t}=s)]$. Each
agent $i$ aims to maximize the expected return given the distribution of the
initial state $\rho(s)$, denoted by the performance objective
$J_{i}=\mathbb{E}_{\rho(s)}V_{i}(s,\theta_{1},...,\theta_{n})$. A naïve agent
updates its policy parameters in the direction of the objective’s gradient

$\begin{split}\nabla_{\theta_{i}}J_{i}=\mathbb{E}_{\rho(s)}\nabla_{\theta_{i}}V_{i}(s,\theta_{1},...,\theta_{n}).\end{split}$
(1)

Learning With Opponent-Learning Awareness (LOLA). Unlike naïve agents, LOLA
agents modify their learning objectives by differentiating through the
anticipated learning steps of the opponents (Foerster et al., 2018a). Given
$n=2$ for simplicity, a first-order LOLA agent assumes a naïve opponent and
uses policy parameter anticipation to optimize
$V_{1}^{\text{LOLA}}(s,\theta_{1},\theta_{2}+\Delta\theta_{2})$ where
$\Delta\theta_{2}=\mathbb{E}_{\rho(s)}\eta\nabla_{\theta_{2}}V_{2}(s,\theta_{1},\theta_{2})$
and $\eta\in\mathbb{R}^{+}$ is the prediction length. Using first-order Taylor
expansion and by differentiating with respect to $\theta_{1}$, the gradient
adjustment for the first LOLA agent (Foerster et al., 2018a) is given by

$\begin{split}\nabla_{\theta_{1}}V_{1}^{\text{LOLA}}(s,\theta_{1},\theta_{2}+\Delta\theta_{2})\approx&\nabla_{\theta_{1}}V_{1}+(\nabla_{\theta_{2}\theta_{1}}V_{1})^{\intercal}\Delta\theta_{2}\\\
&+\underbrace{(\nabla_{\theta_{1}}\Delta\theta_{2})^{\intercal}\nabla_{\theta_{2}}V_{1}}_{\text{shaping}},\end{split}$
(2)

where $V_{1}=V_{1}(s,\theta_{1},\theta_{2})$. The rightmost term in the LOLA
update allows for active shaping of the opponent’s learning. This term has
been proven effective in enforcing cooperation in various games with self-
interested agents, including IPD (Foerster et al., 2018a, b). The LOLA update
can be further extended to non-naïve opponents, resulting in HOLA agents
(Foerster et al., 2018a; Willi et al., 2022).

Look Ahead (LA). LA agents assume that the opponents’ learning steps cannot be
influenced, i.e., cannot be shaped (Zhang and Lesser, 2010; Letcher et al.,
2019). In other words, agent one assumes that the prediction step,
$\Delta{\theta}_{2}$, is independent of the current optimization, i.e.,
$\nabla_{\theta_{1}}\Delta{\theta}_{2}=0$. Therefore, the shaping term
disappears, and the gradient adjustment for the first-order LA agent will be

$\begin{split}\nabla_{\theta_{1}}V_{1}^{\text{LA}}(s,\theta_{1},{\theta}_{2}+\perp\Delta{\theta}_{2})\approx\nabla_{\theta_{1}}V_{1}+(\nabla_{{\theta}_{2}\theta_{1}}V_{1})^{\intercal}\Delta{\theta}_{2},\end{split}$
(3)

where $\perp$ prevents gradient flowing from $\Delta{\theta}_{2}$ upon
differentiation.

The benefits of LOLA and LA have been frequently shown throughout the
literature in games with self-interested agents Zhang and Lesser (2010);
Foerster et al. (2016); Letcher et al. (2019). In the next section, we analyze
the effectiveness of LOLA and LA in fully-cooperative games with common-
interested agents, i.e., $\mathcal{R}_{i}=\mathcal{R}_{j}\;\forall
i,j\in\mathcal{N}$ and, consequently, $V_{i}=V_{j}\;\forall
i,j\in\mathcal{N}$.

##  3\. Miscoordination analysis in fully-cooperative games

To investigate the influence of LOLA and LA on coordination among common-
interested agents, we consider a two-agent two-action coordination game (Claus
and Boutilier, 1998) with an added miscoordination penalty. The game is
defined by a common reward matrix

$\begin{split}\mathcal{R}_{1}=\mathcal{R}_{2}=\begin{bmatrix}\alpha&k\\\
k&\alpha\end{bmatrix},\end{split}$ (4)

where $\alpha>0$ is the coordination reward and $k\leq 0$ is the
miscoordination penalty. We further define $g=\alpha-k>0$ as the
miscoordination regret. The agents are parameterized by $\theta_{1}\in[0,1]$
and $\theta_{2}\in[0,1]$, denoting the probability of choosing the first
action by agents one and two, respectively. With a joint strategy
$(\theta_{1},\theta_{2})$, the common value function of the game is

$\begin{split}V_{1}(\theta_{1},\theta_{2})=V_{2}(\theta_{1},\theta_{2})=2g\theta_{1}\theta_{2}-g(\theta_{1}+\theta
2)+\alpha,\end{split}$ (5)

The game has two equilibrium points, i.e., $(\theta_{1}=0,\theta_{2}=0)$ and
$(\theta_{1}=1,\theta_{2}=1)$, where each agent receives a reward of $\alpha$.
Furthermore, two miscoordination points of the games are
$(\theta_{1}=0,\theta_{2}=1)$ and $(\theta_{1}=1,\theta_{2}=0)$, where agents
receive a penalty of $k$. To solve the coordination game, agents iteratively
adjust their parameters in the direction of the value function gradients,
i.e., $\nicefrac{{V}}{{\theta_{1}}}$ and $\nicefrac{{V}}{{\theta_{2}}}$ in
naïve update rule, $\nicefrac{{V^{\text{LOLA}}}}{{\theta_{1}}}$ and
$\nicefrac{{V^{\text{LOLA}}}}{{\theta_{2}}}$ in LOLA, and
$\nicefrac{{V^{\text{LA}}}}{{\theta_{1}}}$ and
$\nicefrac{{V^{\text{LA}}}}{{\theta_{2}}}$ in LA. Therefore, we can analyze
the dynamics of $\theta_{1}$ and $\theta_{2}$ in LA, LOLA, and naïve agents to
investigate their behaviors.

Figure 1. Phase planes of the LA and LOLA dynamics. Left: unstable saddle
fixed point. Right: unstable fixed point.

If, in the previously defined two-agent two-action coordination game with a
miscoordination regret $g$, the agents are updated following the LA method and
a fixed prediction length $\eta$, then they can be subject to miscoordination
for $g>\tfrac{1}{2\eta}$. Proof. Given Eq. 3, the unconstrained dynamics of LA
agents can be defined by the following differential equations:

$\begin{split}\begin{bmatrix}\nicefrac{{d\theta_{1}}}{{dt}}\\\
\nicefrac{{d\theta_{2}}}{{dt}}\end{bmatrix}=\begin{bmatrix}4\eta g^{2}&2g\\\
2g&4\eta g^{2}\end{bmatrix}\begin{bmatrix}\theta_{1}\\\
\theta_{2}\end{bmatrix}-\begin{bmatrix}2\eta g^{2}+g\\\ 2\eta
g^{2}+g\end{bmatrix}.\end{split}$ (6)

This system of equations has a unique fixed point (zero gradients) at
$\theta_{1}=\theta_{2}=0.5$ (see Figure 1). The eigenvalue analysis of the
coefficient matrix yields two real eigenvalues, $\lambda_{1}=4\eta g^{2}+2g$
and $\lambda_{2}=4\eta g^{2}-2g$, and two respective diagonal and off-diagonal
eigenvectors. While $\lambda_{1}$ is always positive, the sign of
$\lambda_{2}$ depends on the values of both $\eta$ and $g$. For a fixed
prediction length, non-positive values of $\lambda_{2}$ are reached by
$g\leq\tfrac{1}{2\eta}$. In this case, the fixed point is an unstable saddle
point (or unstable line in case of $\lambda_{2}=0$), and the agents, with any
initial values of $\theta_{1}$ and $\theta_{2}$ (except on the fixed point
itself), converge to the equilibrium points (see Figure 1-Left). However, when
the miscoordination regret increases, $g>\tfrac{1}{2\eta}$, the fixed point
becomes an unstable (source) point (see Figure 1-Right). Therefore, some
initial values of $\theta_{1}$ and $\theta_{2}$ naturally lead to the
miscoordination points $(\theta_{1}=0,\theta_{2}=1)$, and Theorem 3 is proved.

If, in the previously defined two-agent two-action coordination game with a
miscoordination regret $g$, the agents are updated following the LOLA method
and a fixed prediction length $\eta$, then they can be subject to
miscoordination for $g>\tfrac{1}{4\eta}$. Proof. Given the LOLA update rule in
Eq. 2, the unconstrained dynamics can be defined as

$\begin{split}\begin{bmatrix}\nicefrac{{d\theta_{1}}}{{dt}}\\\
\nicefrac{{d\theta_{2}}}{{dt}}\end{bmatrix}=\begin{bmatrix}8\eta g^{2}&2g\\\
2g&8\eta g^{2}\end{bmatrix}\begin{bmatrix}\theta_{1}\\\
\theta_{2}\end{bmatrix}-\begin{bmatrix}4\eta g^{2}+g\\\ 4\eta
g^{2}+g\end{bmatrix}.\end{split}$ (7)

This system of equations has a unique fixed point, again at
$\theta_{1}=\theta_{2}=0.5$ (see Figure 1). The eigenvalue analysis of the
coefficient matrix yields two real eigenvalues, $\lambda_{1}=8\eta g^{2}+2g$
and $\lambda_{2}=8\eta g^{2}-2g$, and two respective diagonal and off-diagonal
eigenvectors. Similar to the case of LA agents, $\lambda_{1}$ is always
positive, and the sign of $\lambda_{2}$ depends on the values of both $\eta$
and $g$. For a fixed prediction length, non-positive values of $\lambda_{2}$
are reached by $g\leq\tfrac{1}{4\eta}$. In this case, the fixed point is an
unstable saddle point (or unstable line in case of $\lambda_{2}=0$), and the
agents, with any initial values of $\theta_{1}$ and $\theta_{2}$ (except on
the fixed point itself), converge to the equilibrium points. However, when the
miscoordination regret increases, $g>\tfrac{1}{4\eta}$, the fixed point
becomes an unstable (source) point. Therefore, some initial values of
$\theta_{1}$ and $\theta_{2}$ naturally lead to the miscoordination points,
and Theorem 3 is proved.

If, in the previously defined two-agent two-action coordination game with a
miscoordination regret $g$, the agents follow the naïve updates, then they are
never subject to miscoordination for any value of $g$. Proof. In the case of
the naïve agents, we have

$\begin{split}\begin{bmatrix}\nicefrac{{d\theta_{1}}}{{dt}}\\\
\nicefrac{{d\theta_{2}}}{{dt}}\end{bmatrix}=\begin{bmatrix}0&2g\\\
2g&0\end{bmatrix}\begin{bmatrix}\theta_{1}\\\
\theta_{2}\end{bmatrix}-\begin{bmatrix}g\\\ g\end{bmatrix}.\end{split}$ (8)

Similar to the case of LOLA and LA, this system of equations has a unique
fixed point (zero gradients) at $\theta_{1}=\theta_{2}=0.5$. The eigenvalue
analysis of the coefficient matrix yields two real eigenvalues,
$\lambda_{1}=2g$ and $\lambda_{2}=-2g$, and two respective diagonal and off-
diagonal eigenvectors. This time, however, the eigenvalues are of opposite
signs for any values of $g$, and the fixed point is always an unstable saddle
point. Therefore, any initial values of $\theta_{1}$ and $\theta_{2}$ (except
on the fixed point itself) naturally lead to the equilibrium points.

A closer inspection of LOLA and LA methods reveals two important aspects of
their fundamental ideas.

  1. (1)

Anticipating other agents’ learning is only effective when it is close to
their true future learning. Both LOLA and LA assume a reasoning order for
other agents. If this assumption is wrong, it can negatively affect the
coordination among the cooperative agents. For self-interested agents, it is
natural for them to not unveil their true reasoning orders to each other as
they have different goals. However, common-interested agents can benefit more
from this reasoning information to achieve their common goal.

  2. (2)

The idea of shaping other agents’ learning can be misleading if the other
agents do not follow, making agents more likely to suffer from
miscoordination. LOLA agents constantly underestimate each other, and each
agent intends to shape the learning of others. Letcher et al. Letcher et al.
(2019) indicated that these arrogant behaviors lead to outcomes that are
strictly worse for all agents. It is also clear from Theorems 3 and 3 that the
range of $g$ that can lead to miscoordination in LOLA ($g>\tfrac{1}{4\eta}$)
is larger than the range of $g$ in LA ($g>\tfrac{1}{2\eta}$).

Given the above discussion, we hypothesize that if agents be informed of the
reasoning orders and properly follow the shaping plans, they can improve
coordination among themselves.

##  4\. Hierarchical Learning Anticipation

In this section, we propose Hierarchical Learning Anticipation (HLA), a
methodology designed to improve coordination among fully-cooperative agents.
In contrast to LOLA and LA, HLA determines a hierarchy among the agents to
specify their reasoning orders. Specifically, we first assign $n$ agents to
$n$ hierarchy levels, where levels one and $n$ are the lowest and highest
hierarchy levels, respectively. In each hierarchy level, the assigned agent is
a leader of the lower hierarchy levels and a follower of the higher ones, with
two reasoning rules:

  1. (1)

A leader knows the reasoning orders of the followers and is one level higher

  2. (2)

A follower cannot shape the leaders and only follows their shaping plans

With these reasoning rules, we can address the two previously mentioned
shortcomings of LOLA and LA. Specifically, the first reasoning rule makes sure
that a leader has correct assumptions about the followers’ reasoning orders
and, consequently, can accurately anticipate their future behaviors. With the
second reasoning rule, we can control the shaping plans of the agents and make
sure that the shaping plans are followed. Below, we describe the update rules
in our proposed HLA.

Initialize $\theta_{i}$ $\forall i\in\mathcal{N}$, the prediction length
$\eta$, and the learning rate $\lambda$

Randomly assign the agents into $n$ hierarchy levels.

Rename the agent assigned to level $i$ as the agent $i$, $\forall
i\in\mathcal{N}$

for $\textrm{iter}=1\textrm{ to MaxIter}$ do

for level $i=n\textrm{ to }1$ do

for level $j=1\textrm{ to }i$ do

if $\;\;\,j=1\;\&\;i\neq n$ then
$\Delta\theta_{j}=\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathbb{E}_{\rho(s)}\eta\frac{\partial}{\partial\theta_{j}}V(s,\theta_{1},...,\theta_{i},\bar{\theta}_{i+1},...,\bar{\theta}_{n})$

elif $j\neq 1\;\&\;i=n$ then
$\Delta\theta_{j}=\;\;\mathbb{E}_{\rho(s)}\eta\frac{\partial}{\partial\theta_{j}}V(s,\theta_{1}+\Delta\theta_{1},...,\theta_{j-1}+\Delta\theta_{j-1},\theta_{j},...,\theta_{n})$

elif $j=1\;\&\;i=n$ then
$\Delta\theta_{j}=\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathbb{E}_{\rho(s)}\eta\frac{\partial}{\partial\theta_{j}}V(s,\theta_{1},...,\theta_{n})$

else
$\Delta\theta_{j}=\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathbb{E}_{\rho(s)}\eta\frac{\partial}{\partial\theta_{j}}V(s,\theta_{1}+\Delta\theta_{1},...,\theta_{j-1}+\;\;\;\;\;\;\;\;\Delta\theta_{j-1},\theta_{j},\bar{\theta}_{j+1},...,\bar{\theta}_{n})$

end for

Set $\bar{\theta}_{i}=\theta_{i}+\Delta\theta_{i}$

end for

Update $\theta_{i}=\theta_{i}+\lambda\frac{\Delta\theta_{i}}{\eta}$ $\forall
i\in\mathcal{N}$

end for

Algorithm 1 HLA for a set of $n$ common-interested agents ($\mathcal{N}$).

For simplicity, we set $n=2$, and we assume that agents one ($\theta_{1}$) and
two ($\theta_{2}$) are assigned to the hierarchy levels one and two,
respectively. In other words, agent one is a naïve follower, and agent two is
a first-order leader. Based on our first reasoning rule, the leader performs
first-order reasoning, and its gradient adjustment is similar to a first-order
LOLA agent:

$\begin{split}\nabla_{\theta_{2}}V^{\text{HLA-
Leader}}(s,\theta_{1}+\Delta\theta_{1},\theta_{2})\approx&\nabla_{\theta_{2}}V+(\nabla_{\theta_{1}\theta_{2}}V)^{\intercal}\Delta\theta_{1}\\\
&+(\nabla_{\theta_{2}}\Delta\theta_{1})^{\intercal}\nabla_{\theta_{1}}V,\end{split}$
(9)

where $V=V(s,\theta_{1},\theta_{2})$ is the common value function, and
$\Delta\theta_{1}=\mathbb{E}_{\rho(s)}\eta\nabla_{\theta_{1}}V$. However,
unlike first-order LOLA agents, the first-order leader knows the reasoning
level of the follower, which is a naïve agent. The shaping plan of the first-
order leader is to change its parameters as

$\begin{split}\bar{\theta}_{2}=\theta_{2}+\mathbb{E}_{\rho(s)}\eta\nabla_{\theta_{2}}V^{\text{HLA-
Leader}}(s,\theta_{1}+\Delta\theta_{1},\theta_{2}),\end{split}$ (10)

so that an optimal increase in the common value is achieved after its new
parameters are taken into account by the naïve follower. Therefore, based on
our second reasoning rule, the naïve follower must follow the plan and adjust
its parameters through

$\begin{split}\nabla_{\theta_{1}}&V^{\text{HLA-
Follower}}(s,\theta_{1},\bar{\theta}_{2})\approx\nabla_{\theta_{1}}V\\\
&+(\nabla_{\theta_{2}\theta_{1}}V)^{\intercal}\mathbb{E}_{\rho(s)}\eta\nabla_{\theta_{2}}V^{\text{HLA-
Leader}}(s,\theta_{1}+\Delta\theta_{2},\theta_{2}),\end{split}$ (11)

Since the follower can no longer shape the leader, the shaping term in Eq.
(11) is zero. Therefore, the gradient adjustment for the naïve follower in HLA
is similar to a first-order LA agent, which predicts the leader parameters as
$\theta_{2}+\mathbb{E}_{\rho(s)}\eta\nabla_{\theta_{2}}V^{\text{HLA-
Leader}}(s,\theta_{1}+\Delta\theta_{1},\theta_{2})$.

Algorithm 1 illustrates the HLA update rules for the case of $n$ agents. In
the theorem below, we show how HLA agents can effectively avoid
miscoordination in the two-agent two-action coordination game.

If, in the previously defined two-agent two-action coordination game with a
miscoordination regret $g$, the agents are updated following HLA, then they
are not subject to miscoordination for any value of $g$. Proof. Given Eqs. 9
and 11, the unconstrained dynamics of the HLA agents can be defined by the
following differential equations

$\begin{split}\begin{bmatrix}\nicefrac{{\partial\theta_{1}}}{{\partial t}}\\\
\nicefrac{{\partial\theta_{2}}}{{\partial
t}}\end{bmatrix}=\begin{bmatrix}4\eta g^{2}&2g+16\eta^{2}g^{3}\\\ 2g&8\eta
g^{2}\end{bmatrix}\begin{bmatrix}\theta_{1}\\\
\theta_{2}\end{bmatrix}-\begin{bmatrix}8\eta^{2}g^{3}+4\eta g^{2}+g\\\ 4\eta
g^{2}+g\end{bmatrix},\end{split}$ (12)

resulting in a unique fixed point at $\theta_{1}=\theta_{2}=0.5$ and two real
eigenvalues, $\lambda=6\eta g^{2}\pm 2p\sqrt{9\eta^{2}g^{2}+1}$. Unlike the
case of LOLA and LA, the eigenvalues are now of opposite signs for any values
of $g$, and the fixed point is always an unstable saddle point. Therefore, any
initial values of $\theta_{1}$ and $\theta_{2}$ (except on the fixed point
itself) naturally lead to the equilibrium points, and Theorem 4 is proved.

With Theorem 4, we have shown that HLA naturally avoids miscoordination, and
therefore, our hypothesis is correct. However, the main goal is to demonstrate
that HLA improves coordination among common-interested agents with respect to
naïve learning, which does not take into account learning anticipation at all.
If HLA does not improve the coordination, there is no clear benefit over the
naïve learners, as they also avoid miscoordination.

Figure 2. Converged results for various values of miscoordination regret in
the three-action coordination game.

To further show the benefits of HLA, we employ a standard two-agent three-
action coordination game (Claus and Boutilier, 1998). The game has a common
reward matrix

$\begin{split}\mathcal{R}_{1}=\mathcal{R}_{2}=\begin{bmatrix}10&0&k\\\
0&2&0\\\ k&0&10\end{bmatrix},\end{split}$ (13)

and we define $g=10-k$ as the miscoordination regret. Each agent
$i\in\\{1,2\\}$ is parameterized with three parameters:
$\theta_{i}=\\{\theta_{i}^{1},\theta_{i}^{2},\theta_{i}^{3}\\}$
($\theta_{i}^{j}>0\;\forall j\in\\{1,2,3\\}$ and $\sum_{j}\theta_{i}^{j}=1$),
representing the probability of taking the actions one, two, and three,
respectively. Consequently, the common value function of the game can be
defined as

$\begin{split}V_{1}=V_{2}=10(\theta_{1}^{1}\theta_{2}^{1}+\theta_{1}^{3}\theta_{2}^{3})+2\theta_{1}^{2}\theta_{2}^{2}+k(\theta_{1}^{1}\theta_{2}^{3}+\theta_{1}^{3}\theta_{2}^{1}).\end{split}$
(14)

The game has two global equilibrium points
($\theta_{1}^{1}=1,\theta_{2}^{1}=1$ and $\theta_{1}^{3}=1,\theta_{2}^{3}=1$),
one local equilibrium point ($\theta_{1}^{2}=1,\theta_{2}^{2}=1$), and two
miscoordination points ($\theta_{1}^{1}=1,\theta_{2}^{3}=1$ and
$\theta_{1}^{3}=1,\theta_{2}^{1}=1$).

In Figure 2, we depict the converged results for this game for naïve, LA,
LOLA, and HLA agents, for various values of miscoordination regret $g$. The
experiments are run 500 times until convergence, with random initializations.
For HLA, we randomly assigned the agents to the hierarchy levels (leader and
follower) in each experiment. From Figure 2, we find that both LA and LOLA
agents are subject to miscoordination for high values of $g$, which is
consistent with our findings for the two-action coordination game. However,
the most interesting aspect of this experiment is that by increasing the value
of $g$, the coordination among the naïve agents reduces, leading them to the
local equilibrium point, whereas our HLA consistently achieves the highest
reward, independently of the miscoordination regret. These results clearly
show the benefits of HLA over other methods.

##  5\. Discussion and Future Work

In this study, we extended the applicability of learning anticipation to games
with fully-cooperative agents. We demonstrated that methods such as LOLA and
LA, which heavily benefit from learning anticipation in general-sum games, can
significantly reduce coordination when the agents are fully-cooperative. We
first hypothesized that when agents know the reasoning orders of others and
properly follow the shaping plans, they could improve coordination among
themselves. To verify our hypothesis, we proposed the novel HLA method, which
incorporates a hierarchy to regulate the reasoning orders and shaping plans of
agents. Having HLA, we can now use the benefits of learning anticipation in
fully-cooperative games.

Nevertheless, there are still some unanswered questions about applying HLA to
more complex fully-cooperative games. For instance, we limit ourselves to
fully-cooperative normal-form games with differentiable objective functions
where agents can access gradients and Hessians. In many multi-agent problems,
the objective functions are non-differentiable. In these cases, the agents
must estimate the higher-order gradients with various approximation methods.
Furthermore, the agents may not have access to other agents’ exact parameters,
and they may need to infer the parameters from state-action trajectories.
Consequently, future studies on the current topic are required.

## References

  * (1)
  * Axelrod and Hamilton (1981) Robert Axelrod and William D Hamilton. 1981.  The evolution of cooperation.  _science_ 211, 4489 (1981), 1390–1396.
  * Claus and Boutilier (1998) Caroline Claus and Craig Boutilier. 1998.  The dynamics of reinforcement learning in cooperative multiagent systems.  _AAAI/IAAI_ 1998, 746-752 (1998), 2.
  * Foerster et al. (2018a) Jakob Foerster, Richard Y Chen, Maruan Al-Shedivat, Shimon Whiteson, Pieter Abbeel, and Igor Mordatch. 2018a.  Learning with Opponent-Learning Awareness. In _Proceedings of the 17th International Conference on Autonomous Agents and MultiAgent Systems_. 122–130.
  * Foerster et al. (2018b) Jakob Foerster, Gregory Farquhar, Maruan Al-Shedivat, Tim Rocktäschel, Eric Xing, and Shimon Whiteson. 2018b.  Dice: The infinitely differentiable monte carlo estimator. In _International Conference on Machine Learning_. PMLR, 1529–1538.
  * Foerster et al. (2016) Jakob N Foerster, Yannis M Assael, Nando De Freitas, and Shimon Whiteson. 2016.  Learning to communicate with deep multi-agent reinforcement learning.  _arXiv preprint arXiv:1605.06676_ (2016).
  * Letcher et al. (2019) Alistair Letcher, Jakob Foerster, David Balduzzi, Tim Rocktäschel, and Shimon Whiteson. 2019\.  Stable Opponent Shaping in Differentiable Games. In _International Conference on Learning Representations_.
  * Littman (1994) Michael L Littman. 1994\.  Markov games as a framework for multi-agent reinforcement learning.  In _Machine learning proceedings 1994_. Elsevier, 157–163.
  * Willi et al. (2022) Timon Willi, Alistair Hp Letcher, Johannes Treutlein, and Jakob Foerster. 2022.  COLA: consistent learning with opponent-learning awareness. In _International Conference on Machine Learning_. PMLR, 23804–23831.
  * Zhang and Lesser (2010) Chongjie Zhang and Victor Lesser. 2010.  Multi-agent learning with policy prediction. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 24.


# A Logistic Regression Approach to Field Estimation Using Binary Measurements

Alex S. Leong, Mohammad Zamani, and Iman Shames A. Leong and M. Zamani are
with Land Division, Defence Science and Technology Group, Melbourne,
Australia, email {alex.leong<EMAIL_ADDRESS>. I. Shames is
with the CIICADA Lab, College of Engineering and Computer Science, Australian
National University, Canberra, Australia, e-mail<EMAIL_ADDRESS>
###### Abstract

In this letter, we consider the problem of field estimation using binary
measurements. Previous work has formulated the problem as a parameter
estimation problem, with the parameter estimation carried out in an online
manner using sequential Monte Carlo techniques. In the current work, we
consider an alternative approach to the parameter estimation based on online
logistic regression. The developed algorithm is less computationally intensive
than the sequential Monte Carlo approach, while having more reliable
estimation performance.

##  I Introduction

The mapping of a contaminated area due to a chemical, biological,
radiological, or nuclear (CBRN) leakage is an important task in allowing
operations to be carried out in the affected area. The use of mobile
autonomous vehicles in carrying out this mapping is likely to be crucial, as
existing infrastructure for detecting/sensing such contaminants in the
affected area is often lacking.

Previous work in field estimation using mobile autonomous vehicles include [1,
2, 3]. The field was modelled as a sum of radial basis functions, with the
field estimation problem then reduced to estimating the parameters of the
basis functions.

Motivated by the limitations of current technologies in chemical sensors [4,
5], in [3] it was assumed that the noisy sensor measurements were coarsely
quantized, and in particular binary [6, 7]. The field parameters were
estimated using sequential Monte Carlo techniques [8], and an active sensing
mechanism [9, 10] for choosing the next measurement location given the
measurements currently collected, was proposed.

In this letter, we also consider the use of binary measurements. Instead of
the sequential Monte Carlo approach, we instead propose an approach that views
the parameter estimation problem as an online logistic regression problem
[11]. The parameters are updated recursively using a Newton-type update [12].
With further approximations, the computational load can be significantly
reduced from that of the sequential Monte Carlo approach. The active sensing
mechanism of [3] based on maximization of an information measure does not
extend in a straightforward way to the current approach, and thus an
alternative active sensing mechanism is proposed based on numerical properties
of the Hessian of the cost function.

The letter is organized as follows. Section II gives the system model. The
logistic regression approach is presented in Section III, first in an “exact”
form using online Newton’s method (ONM). Section III-B then introduces an
approximate ONM that is computationally less intensive than the exact
counterpart. An active sensing mechanism, inspired by the results of the
preceding sections, is presented in Section IV. Simulation results that
illustrates the performance of the proposed methods are presented in Section
V.

_Notation_ : The operators $\sigma_{\textnormal{min}}(\cdot)$ and
$\lambda_{\textnormal{min}}(\cdot)$ denote the minimum singular value and
minimum eigenvalue of their arguments, respectively.

##  II System Model

We assume that the field can be sufficiently well approximated [13] by

$\phi(\mathbf{x})=\sum_{i=1}^{p}\beta_{i}K_{i}(\mathbf{x})$ (1)

where $\mathbf{x}\in\mathbb{R}^{2}$ is the position, $K_{i}(\cdot)$ are basis
functions, and $\beta_{i}$ are the coefficients multiplying the basis
functions. Similar models have been used in e.g. [1, 14, 2, 3]. For the basis
functions we will use the Gaussian kernel

$K_{i}(\mathbf{x})=\exp\left(-\frac{\|\mathbf{c}_{i}-\mathbf{x}\|^{2}}{\sigma_{i}^{2}}\right)$

where $\mathbf{c}_{i}\in\mathbb{R}^{2}$ represents the center of the $i$-th
basis function. Note that this model captures a pollutant field at steady-
state in a sufficiently rich, yet tractable analytic model. In reality the
physical processes that generate such a field may be a lot more complicated.
However, this model serves as a surrogate that approximates the true field
provided that a sufficient number of basis functions are considered [13].

For the measurement model, at position $\mathbf{x}$, the sensor is first
assumed to have a noisy measurement

$y(\mathbf{x})=\phi(\mathbf{x})+v(\mathbf{x}),$ (2)

where $v(\mathbf{x})$ is noise. This noisy measurement is then thresholded to
give a binary measurement

$z(\mathbf{x})=\mathds{1}(y(\mathbf{x})>\tau),$

where $\tau$ is a positive scalar threshold, and $\mathds{1}(\cdot)$ is the
indicator function that returns 1 if its argument is true and 0 otherwise.

In the absence of knowledge of the “true” parameters, we assume that
$\mathbf{c}_{i}$ and $\sigma_{i}^{2},i=1,\dots,p$ are to be chosen by
us.111This corresponds to case (i) of [3]. Case (ii), which estimates all
$\beta_{i}$, $\mathbf{c}_{i}$, $\sigma_{i}^{2},i=1,\dots,p$, was also
considered in [3], but was found to give less reliable estimates. The
objective is to estimate the parameters
$\bm{\beta}\triangleq(\beta_{1},\dots,\beta_{p})$ from binary measurements
collected by an autonomous vehicle, where the vehicle can collect measurements
from any position in the area of interest $\mathcal{S}$.

##  III Logistic regression approach

We consider a cost function based on logistic regression. Recall the logistic
function $\ell(x)=1/(1+\exp(-\eta x))$, where $\eta>0$ is a parameter. Larger
values of $\eta$ will more closely approximate the function $\mathds{1}(x>0)$.

Let $\tilde{z}\triangleq 2z-1.$ Since $z\in\\{0,1\\}$, we have
$\tilde{z}\in\\{-1,+1\\}$. Consider the following cost function,
$J_{k}(\bm{\beta})=\sum_{t=0}^{k}g_{t}(\bm{\beta})$ with per stage cost

$g_{t}(\bm{\beta})=\log\left(1+\exp(-\eta\tilde{z}_{t}(\bm{\beta}^{T}\mathbf{K}(\mathbf{x}_{t})-\tau))\right),$

where
$\mathbf{K}(\mathbf{x})\triangleq\begin{bmatrix}K_{1}(\mathbf{x})&K_{2}(\mathbf{x})&\dots&K_{p}(\mathbf{x})\end{bmatrix}^{T}$
and $t$ is a time/iteration index. The measurement $\tilde{z}_{t}$ and
position $\mathbf{x}_{t}$ are assumed to be provided to us via a sensor and
e.g. GPS, respectively. The case where $\eta=1$ is similar to the costs used
in logistic regression problems, see e.g. [15, p.516].

We are interested in online estimation of $\bm{\beta}$ recursively rather than
solving an optimization problem in batch mode. In this setup a new measurement
may be collected after each time step. We use the online Newton’s method (ONM)
of [12]. To this aim, the gradient of $J_{k}(\bm{\beta})$ is computed as the
following:

$\begin{split}\nabla J_{k}(\bm{\beta})&=\sum_{t=0}^{k}\nabla
g_{t}(\bm{\beta})\\\
&=\sum_{t=0}^{k}\frac{-\eta\tilde{z}_{t}}{1+\exp(\eta\tilde{z}_{t}(\bm{\beta}^{T}\mathbf{K}(\mathbf{x}_{t})-\tau))}\mathbf{K}(\mathbf{x}_{t}).\end{split}$
(3)

Similarly, the Hessian of $J_{k}(\bm{\beta})$ is

$\begin{split}&\nabla^{2}J_{k}(\bm{\beta})=\sum_{t=0}^{k}\nabla^{2}g_{t}(\bm{\beta})\\\
&=\sum_{t=0}^{k}\frac{\eta^{2}\tilde{z}_{t}^{2}\exp(\eta\tilde{z}_{t}(\bm{\beta}^{T}\mathbf{K}(\mathbf{x}_{t})-\tau))}{\big{(}1+\exp(\eta\tilde{z}_{t}(\bm{\beta}^{T}\mathbf{K}(\mathbf{x}_{t})-\tau))\big{)}^{2}}\mathbf{K}(\mathbf{x}_{t})\mathbf{K}^{T}(\mathbf{x}_{t}).\end{split}$
(4)

The ONM for recursively estimating $\bm{\beta}$ is then given by

$\bm{\hat{\beta}}_{k+1}=\bm{\hat{\beta}}_{k}-\left(\nabla^{2}J_{k}(\bm{\hat{\beta}}_{k})\right)^{-1}\nabla
J_{k}(\bm{\hat{\beta}}_{k}),$ (5)

where $\nabla J_{k}(\bm{\hat{\beta}}_{k})$ and
$\nabla^{2}J_{k}(\bm{\hat{\beta}}_{k})$ can be obtained by substituting
$\bm{\hat{\beta}}_{k}$ into (3) and (4), respectively.

###  III-A Performance of online Newton’s method

Define $\bm{\beta}^{*}_{k}=\arg\min_{\bm{\beta}}J_{k}(\bm{\beta})$. We make
the following:

######  Assumption 1.

  1. (i)

There exists a positive real scalar $\bar{v}$ such that
$\|\bm{\beta}^{*}_{k+1}-\bm{\beta}^{*}_{k}\|\leq\bar{v}$ for all $k>0$.

  2. (ii)

There exist positive real scalars $h_{k}$ such that
$\|(\nabla^{2}J_{k}(\bm{\beta}^{*}_{k}))^{-1}\|\leq 1/h_{k}$.

  3. (iii)

In a ball of radius $\rho$ centered at $\bm{\beta}^{*}_{k}$, the Hessian is
Lipschitz continuous with Lipschitz constant $L_{k}$.

  4. (iv)

There exist positive integer and real scalars $\overline{k}$ and
$\underline{\lambda}$ such that
$\lambda_{\textnormal{min}}\left(\sum_{t=\tau}^{\overline{k}+\tau-1}\nabla^{2}g_{t}(\bm{\beta})\big{)}\right)\geq\underline{\lambda},\,\forall\tau.$

Assumptions (i)-(iii) are assumptions from [12], while assumption (iv) is a
“persistence of excitation” type assumption that there is sufficient richness
in our collected measurements.

We aim to bound the regret of the ONM defined by

$\displaystyle\mathrm{Reg}(T)=\sum_{k=0}^{T}\left(J_{k}(\hat{\bm{\beta}}_{k})-J_{k}(\bm{\beta}^{*}_{k})\right).$
(6)

For a bound on $\|(\nabla^{2}J_{k}(\bm{\beta}^{*}_{k}))^{-1}\|$, suppose we
use the induced 2-norm. Then

$\displaystyle\|(\nabla^{2}J_{k}(\bm{\beta}^{*}_{k}))^{-1}\|_{2}$
$\displaystyle=\frac{1}{\sigma_{\textnormal{min}}\big{(}\sum_{t=1}^{k}\nabla^{2}g_{t}(\bm{\beta}^{*}_{k})\big{)}}$
$\displaystyle=\frac{1}{\lambda_{\textnormal{min}}\big{(}\sum_{t=1}^{k}\nabla^{2}g_{t}(\bm{\beta}^{*}_{k})\big{)}}\leq\dfrac{1}{\lceil
k/\overline{k}\rceil\underline{\lambda}},$

where the last inequality holds by Assumption 1 (iv) and Corollary 4.3.15 of
[16]. Thus we can set $h_{k}\triangleq\lceil
k/\overline{k}\rceil\underline{\lambda}$.

For a Lipschitz constant $L_{k}$ of the Hessian of $J_{k}(\cdot)$, we can
obtain $L_{k}\triangleq kL_{g}$ that increases linearly with $k$, since one
has a Lipschitz constant $L_{g}$ for the Hessian of each $g_{t}(\cdot)$, and
there are $k$ such terms in the summation in $J_{k}(\cdot)$.

For any iterations such that

$\|\bm{\hat{\beta}}_{k}-\bm{\beta}^{*}_{k}\|<\min\\{\rho,\frac{2\underline{\lambda}}{3\overline{k}L_{g}}\\}\triangleq\gamma,$

then
$\|\bm{\hat{\beta}}_{k}-\bm{\beta}^{*}_{k}\|<\min\\{\rho,\frac{2h_{k}}{3L_{k}}\\}$,
since

$\displaystyle\|\bm{\hat{\beta}}_{k}-\bm{\beta}^{*}_{k}\|$
$\displaystyle<\min\\{\rho,\frac{2\underline{\lambda}}{3\overline{k}L_{g}}\\}=\min\\{\rho,\frac{2k/\overline{k}\underline{\lambda}}{3kL_{g}}\\}$
$\displaystyle\leq\min\\{\rho,\frac{2\lceil
k/\overline{k}\rceil\underline{\lambda}}{3kL_{g}}\\}=\min\\{\rho,\frac{2h_{k}}{3L_{k}}\\}.$

Assume also that
$\bar{v}\leq\gamma-(3\overline{k}L_{g}/2\underline{\lambda})\gamma^{2}.$ Then
the regret of ONM satisfies [12]:

$\mathrm{Reg}(T)\leq\dfrac{\kappa}{1-\dfrac{3\overline{k}L_{g}}{2\underline{\lambda}}\gamma}(V_{T}+\delta)$
(7)

where $\kappa$ is a positive scalar such that if
$\|\bm{\beta}-\bm{\beta}^{*}_{k}\|\leq\gamma$ then
$|J_{k}(\bm{\beta})-J_{k}(\bm{\beta}^{*}_{k})|\leq\kappa\|\bm{\beta}-\bm{\beta}^{*}_{k}\|$
for $k=1,\dots,T$. Moreover,
$V_{T}=\sum_{k=1}^{T}\|\bm{\beta}^{*}_{k}-\bm{\beta}^{*}_{k-1}\|$ and
$\delta=(3\overline{k}L_{g}/2\underline{\lambda})\left(\|\bm{\hat{\beta}}_{0}-\bm{\beta}^{*}_{0}\|^{2}-\|\bm{\hat{\beta}}_{T}-\bm{\beta}^{*}_{T}\|^{2}\right)$.

Additionally, [12, Lemma 2] yields

$\displaystyle\|\bm{\hat{\beta}}_{k}-\bm{\beta}^{*}_{k}\|$
$\displaystyle=\|\bm{\hat{\beta}}_{k}-\bm{\beta}^{*}_{k}-\bm{\beta}^{*}_{k-1}+\bm{\beta}^{*}_{k-1}\|$
$\displaystyle\leq\|\bm{\hat{\beta}}_{k}-\bm{\beta}^{*}_{k-1}\|+\|\bm{\beta}^{*}_{k-1}-\bm{\beta}^{*}_{k}\|$
$\displaystyle\leq\dfrac{3\overline{k}L_{g}}{2\underline{\lambda}}\|\bm{\hat{\beta}}_{k-1}-\bm{\beta}^{*}_{k-1}\|^{2}+\|\bm{\beta}^{*}_{k-1}-\bm{\beta}^{*}_{k}\|$
$\displaystyle\leq\dfrac{3\overline{k}L_{g}}{2\underline{\lambda}}\|\bm{\hat{\beta}}_{k-1}-\bm{\beta}^{*}_{k-1}\|^{2}+\bar{v}.$
(8)

If $\bar{v}\leq\underline{\lambda}/(6\overline{k}L_{g})$, then (8) gives

$\lim\sup_{k\rightarrow\infty}\|\bm{\hat{\beta}}_{k}-\bm{\beta}^{*}_{k}\|\leq\dfrac{\underline{\lambda}}{3\overline{k}L_{g}}.$
(9)

###  III-B Computationally efficient approximate ONM

The online Newton’s method (ONM) (5) requires computation of the gradient
$\nabla J_{k}(\bm{\hat{\beta}}_{k})$ and Hessian
$\nabla^{2}J_{k}(\bm{\hat{\beta}}_{k})$, both of which are summations where
the number of terms will increase with $k$. This can be computationally
intensive for large $k$. For more efficient numerical implementation, we use a
recursive formulation to approximately capture the accumulated cost, and
derive the approximate online Newton step for minimizing it. Under this
approximation, the computational cost per iteration does not increase over
time. Moreover, the computation required is significantly less than the
sequential Monte Carlo approach of [3], see Section V.

Write the cost function as

$J_{k}(\bm{\beta})=J_{k-1}(\bm{\beta})+g_{k}(\bm{\beta}).$

First, motivated by Assumption 1 (i), we assume that $\bm{\hat{\beta}}_{k}$
will approximately minimize $J_{k-1}(\cdot)$, so that by the first order
condition $\nabla J_{k-1}(\bm{\hat{\beta}}_{k})\approx 0$, and hence $\nabla
J_{k}(\bm{\hat{\beta}}_{k})\approx\nabla g_{k}(\bm{\hat{\beta}}_{k})$.

Next, the Hessian of $J_{k}(\cdot)$ at $\hat{\bm{\beta}}_{k}$ is approximated
as

$\begin{split}\nabla^{2}J_{k}(\hat{\bm{\beta}}_{k})&=\nabla^{2}J_{k-1}(\hat{\bm{\beta}}_{k})+\nabla^{2}g_{k}(\hat{\bm{\beta}}_{k})\\\
&\approx\nabla^{2}J_{k-1}(\hat{\bm{\beta}}_{k-1})+\nabla^{2}g_{k}(\hat{\bm{\beta}}_{k}),\end{split}$
(10)

where in the second line of (10), motivated by Assumption 1 (iii), we make the
assumption that
$\nabla^{2}J_{k-1}(\hat{\bm{\beta}}_{k})\approx\nabla^{2}J_{k-1}(\hat{\bm{\beta}}_{k-1})$.
This allows us to use the Hessian $\nabla^{2}J_{k-1}(\hat{\bm{\beta}}_{k-1})$
from the previous time step, without having to recalculate the terms
$g_{t}(\hat{\bm{\beta}}_{k}),t=0,\dots,k-1$ based on the new estimate
$\hat{\bm{\beta}}_{k}$.

We can then consider the iterations for an approximate online Newton’s method
of the form

$\bm{\hat{\beta}}_{k+1}=\bm{\hat{\beta}}_{k}-\left(H_{k}(\bm{\hat{\beta}}_{k})\right)^{-1}G_{k}(\bm{\hat{\beta}}_{k}),$
(11)

where

$\begin{split}G_{k}(\bm{\hat{\beta}}_{k})&=\nabla
g_{k}(\bm{\hat{\beta}}_{k})\\\
H_{k}(\hat{\bm{\beta}}_{k})&=H_{k-1}(\hat{\bm{\beta}}_{k-1})+\nabla^{2}g_{k}(\hat{\bm{\beta}}_{k}).\end{split}$
(12)

We next provide a numerically robust formulation of the approximate Newton
iteration (11) and (12). The first idea is to use a full rank initialization
$(H_{0}(\bm{\hat{\beta}}_{0}))^{-1}=\epsilon I$ for some $\epsilon>0$. Recall
the Hessian of the per stage cost

$\nabla^{2}g_{k}(\bm{\hat{\beta}}_{k})=\frac{\eta^{2}\tilde{z}_{k}^{2}\exp(\eta\tilde{z}_{k}(\bm{\hat{\beta}}_{k}^{T}\mathbf{K}(\mathbf{x}_{k})-\tau))}{\big{(}1+\exp(\eta\tilde{z}_{k}(\bm{\hat{\beta}}_{k}^{T}\mathbf{K}(\mathbf{x}_{k})-\tau))\big{)}^{2}}\mathbf{K}(\mathbf{x}_{k})\mathbf{K}^{T}(\mathbf{x}_{k}).$

In short this can be written as

$\nabla^{2}g_{k}(\bm{\hat{\beta}}_{k})=h(\mathbf{x}_{k},\bm{\hat{\beta}}_{k})\mathbf{K}(\mathbf{x}_{k})\mathbf{K}^{T}(\mathbf{x}_{k}),$

where $h(\mathbf{x}_{k},\bm{\hat{\beta}}_{k})\in\mathbb{R}$. Then the rank one
update (see e.g. [15, 16]) can be used to update the next and subsequent
iterations:

$\displaystyle(H_{k}(\bm{\hat{\beta}}_{k}))^{-1}=(H_{k-1}(\bm{\hat{\beta}}_{k-1}))^{-1}$
$\displaystyle\quad\times\Big{(}I-\dfrac{h(\mathbf{x}_{k},\bm{\hat{\beta}}_{k})\mathbf{K}(\mathbf{x}_{k})\mathbf{K}^{T}(\mathbf{x}_{k})(H_{k-1}(\bm{\hat{\beta}}_{k-1}))^{-1}}{1+h(\mathbf{x}_{k},\bm{\hat{\beta}}_{k})\mathbf{K}^{T}(\mathbf{x}_{k})(H_{k-1}(\bm{\hat{\beta}}_{k-1}))^{-1}\mathbf{K}(\mathbf{x}_{k})}\Big{)}.$

######  Remark 1.

The formal treatment of the question of robustness of (12) is beyond the scope
of this letter. However, in light of [17] and the numerical example presented
later in this letter, we conjecture that the proposed approximate ONM is
robust to the introduced approximation errors.

##  IV Active Sensing

In this section we consider the problem of adaptively selecting the next
measurement location $\mathbf{x}_{k+1}$, based on the measurements
$z_{0},\dots,z_{k}$ collected so far, with the aim of achieving better
estimation performance. We refer to this as _active sensing_ [9]. The active
sensing mechanism of [3] was based optimizing an information measure [10]
which does not extend in a straightforward manner to the current approach.
Based on Assumption 1 and informed by the results of Section III-A, having a
larger $\underline{\lambda}$ may be desirable in terms of convergence
behaviour. Motivated by this, in this letter we propose to choose position
$\mathbf{x}_{k+1}$ to be the one that maximizes the minimum eigenvalue of the
expected Hessian (13).

For a given future position $\mathbf{x}^{\prime}$, consider an “expected
Hessian” term:

$\small\begin{split}&\textnormal{Hess}^{+}(\mathbf{x}^{\prime})\triangleq\nabla^{2}J_{k}(\bm{\hat{\beta}}_{k+1})+\mathbb{P}(\tilde{z}_{k+1}=1)\\\
&\quad\quad\times\frac{\eta^{2}\exp(\eta(\bm{\hat{\beta}}_{k+1}^{T}\mathbf{K}(\mathbf{x}^{\prime})-\tau))}{\big{(}1+\exp(\eta(\bm{\hat{\beta}}_{k+1}^{T}\mathbf{K}(\mathbf{x}^{\prime})-\tau))\big{)}^{2}}\mathbf{K}(\mathbf{x}^{\prime})\mathbf{K}^{T}(\mathbf{x}^{\prime})\\\
&\quad+\mathbb{P}(\tilde{z}_{k+1}=-1)\\\
&\quad\quad\times\frac{\eta^{2}\exp(-\eta(\bm{\hat{\beta}}_{k+1}^{T}\mathbf{K}(\mathbf{x}^{\prime})-\tau))}{\big{(}1+\exp(-\eta(\bm{\hat{\beta}}_{k+1}^{T}\mathbf{K}(\mathbf{x}^{\prime})-\tau))\big{)}^{2}}\mathbf{K}(\mathbf{x}^{\prime})\mathbf{K}^{T}(\mathbf{x}^{\prime})\\\
&=\nabla^{2}J_{k}(\bm{\hat{\beta}}_{k+1})\\\
&\quad+\frac{\eta^{2}\exp(\eta(\bm{\hat{\beta}}_{k+1}^{T}\mathbf{K}(\mathbf{x}^{\prime})-\tau))}{\big{(}1+\exp(\eta(\bm{\hat{\beta}}_{k+1}^{T}\mathbf{K}(\mathbf{x}^{\prime})-\tau))\big{)}^{2}}\mathbf{K}(\mathbf{x}^{\prime})\mathbf{K}^{T}(\mathbf{x}^{\prime})\end{split}$
(13)

where the last equality holds since
$\exp(-y)/(1+\exp(-y))^{2}=\exp(y)/(1+\exp(y))^{2}$ and
$\mathbb{P}(\tilde{z}_{k+1}=1)+\mathbb{P}(\tilde{z}_{k+1}=-1)=1$.

The optimization problem of interest, to maximize the minimum eigenvalue of
the expected Hessian (13), is

$\mathbf{x}_{k+1}=\textnormal{arg}\max\limits_{\mathbf{x}^{\prime}\in\mathcal{X}_{k+1}}\lambda_{\min}(\textnormal{Hess}^{+}(\mathbf{x^{\prime}})),$

where $\mathcal{X}_{k+1}$ is the set of possible future positions222The set
$\mathcal{X}_{k+1}$ may, e.g., capture the set of reachable positions from the
current state of the mobile sensor platform..

In the case where some of the positions in $\mathcal{X}_{k+1}$ may be far from
the vehicle’s current position, we may instead decide to only travel a limited
distance $\rho$ towards the intended position before taking a measurement.
That is, we collect a new measurement at position
$\mathbf{x}_{k+1}=\mathbf{x}_{k}+\rho\mathbf{d}_{k}$, where
$\mathbf{d}_{k}=(\tilde{\mathbf{x}}_{k+1}-\mathbf{x}_{k})/\|\tilde{\mathbf{x}}_{k+1}-\mathbf{x}_{k}\|$
and
$\tilde{\mathbf{x}}_{k+1}=\textnormal{arg}\max\limits_{\mathbf{x}^{\prime}\in\mathcal{X}_{k+1}}\lambda_{\min}(\textnormal{Hess}^{+}(\mathbf{x^{\prime}})).$
One possibility for $\mathcal{X}_{k+1}$ with this approach is then the set of
centers of the basis functions, i.e.
$\mathcal{X}_{k+1}=\\{\mathbf{c}_{1},\dots,\mathbf{c}_{p}\\}$. We can also
smooth out the motion by taking a convex combination of the new
($\mathbf{d}_{k}$) and previous ($\tilde{\mathbf{d}}_{k-1}$) directions,
namely $\mathbf{x}_{k+1}=\mathbf{x}_{k}+\rho\tilde{\mathbf{d}}_{k}$, where
$\tilde{\mathbf{d}}_{k}=(\alpha\mathbf{d}_{k}+(1-\alpha)\tilde{\mathbf{d}}_{k-1})/\|\alpha\mathbf{d}_{k}+(1-\alpha)\tilde{\mathbf{d}}_{k-1}\|$,
for some $\alpha\in[0,1]$. We have found this smoothing for values of $\alpha$
of around 0.3-0.5 to significantly improve performance.

For the approximate ONM of Section III-B, we can similarly choose positions to
maximize the minimum eigenvalue of an expected Hessian term, but instead of
(13) we use instead

$\displaystyle\widetilde{\textnormal{Hess}}^{+}(\mathbf{x}^{\prime})\triangleq
H_{k}(\bm{\hat{\beta}}_{k})$
$\displaystyle\quad+\frac{\eta^{2}\exp(\eta(\bm{\hat{\beta}}_{k+1}^{T}\mathbf{K}(\mathbf{x}^{\prime})-\tau))}{\big{(}1+\exp(\eta(\bm{\hat{\beta}}_{k+1}^{T}\mathbf{K}(\mathbf{x}^{\prime})-\tau))\big{)}^{2}}\mathbf{K}(\mathbf{x}^{\prime})\mathbf{K}^{T}(\mathbf{x}^{\prime}).$

##  V Numerical Results

We consider a $[0,100\textnormal{m}]\times[0,100\textnormal{m}]$ area of
interest $\mathcal{S}$. We consider randomly generated true fields using the
model (1) with $p^{\textnormal{true}}=4$, with parameters
$\beta_{i}^{\textnormal{true}}\sim U(0.7,1.4)$, each component of
$\textbf{c}_{i}^{\textnormal{true}}$ $\sim U(5,95)$, and
$\sigma_{i}^{\textnormal{true}}\sim U(25,45)$, for
$i=1,\dots,p^{\textnormal{true}}$. In the measurement model (2), we assume
$v(\mathbf{x})\sim\mathcal{N}(0,\sigma_{v}^{2})$ is i.i.d. Gaussian noise,
with $\sigma_{v}^{2}=0.1$, and the sensor threshold $\tau=1$.

Figure 1: True field (left) and estimated field after 1000 measurements
(right)

First, we use the approximate online Newton’s method of Section III-B together
with the active sensing method of Section IV. The logistic function parameter
is chosen as $\eta=5$. We use $p=16$ basis functions, with centers
$\mathbf{c}_{i}\in\\{12.5,37.5,62.5,87.5\\}^{2}$ and
$\sigma_{i}=25,i=1,\dots,16$. The initial estimate $\bm{\hat{\beta}}_{0}$ is
randomly chosen with each component $\sim U(0,1)$. The full rank
initialization $(H_{0}(\bm{\hat{\beta}}_{0}))^{-1}=\epsilon I$ with
$\epsilon=0.1$. In active sensing we use
$\mathcal{X}_{k+1}=\\{\mathbf{c}_{1},\dots,\mathbf{c}_{16}$},
$\rho=5\textnormal{m}$, and $\alpha=0.4$.

A plot of a sample randomly generated field that we wish to estimate is given
on the left of Fig. 1. On the right of Fig. 1 we show the estimated field
after a single simulation run consisting of the collection of 1000
measurements. We see that the approximate ONM can capture the qualitative
shape of the true field.

We next compare the performance of approximate ONM, “exact” ONM, and the
sequential Monte Carlo (SMC) approach of [3]. We note that for this particular
example, during the initial steps, the exact ONM is not numerically stable due
to the Hessian being close to singular (which violates Assumption 1 (ii)). We
used a hybrid Newton implementation that used a damping multiplier of $0.1$
for the derivative until the Hessian was large enough, at which point we
switched to the original exact ONM. Moreover a small regularization factor of
$0.1I$ was added to the Hessian when it got close to being singular. For
performance evaluation, we use the following Mean Squared Error (MSE) type
measure:

$\frac{1}{|\tilde{\mathcal{S}}|}\sum_{\mathbf{x}\in\tilde{\mathcal{S}}}\left(\mathbb{P}(z=1|\bm{\beta};\mathbf{x})-\mathbb{P}(z=1|\bm{\hat{\beta}}_{k};\mathbf{x})\right)^{2},$
(14)

which was shown in [3] to be a useful measure in comparing the closeness of
two fields. In (14) $\tilde{\mathcal{S}}$ is a discretization of $\mathcal{S}$
where the $x$ and $y$ axes are each discretized into 32 values, giving
cardinality $|\tilde{\mathcal{S}}|=32^{2}=1024$. From the model of Section II,
we can derive the expression
$\mathbb{P}(z=1|\bm{\beta};\mathbf{x})=1-\Phi\Big{(}\big{(}\tau-\sum_{i=1}^{p}\beta_{i}K_{i}(\mathbf{x})\big{)}/\sigma_{v}\Big{)}$,
where $\Phi(\cdot)$ is the cdf of the standard normal distribution.

In Fig. 2 we show boxplots for the MSE’s at $k=1000$, obtained from simulation
runs for 100 different fields. The computations are done using a Core i7-9700
desktop with 16 Gb of RAM. The results are summarised in Table I. The
approximate ONM and SMC methods generally perform similar to each other.
However, we see that the approximate ONM performs more consistently than the
exact ONM and SMC approaches, with fewer outliers that have large MSEs.
Moreover, the running time for the approximate ONM is significantly faster
than the other two methods. In regards to the poorer performance of the exact
ONM compared to the approximate ONM, Fig. 3 plots the estimates of
$\bm{\beta}$ over time for a single run of both the approximate and exact
ONMs. We see that the parameter estimates for the exact ONM have larger
variations than the approximate ONM, especially in the earlier iterations.
This might be due the fact that the exact ONM is sensitive to initial
condition and as a result had to be modified in our implementation. A rigorous
study of this phenomenon and the dependency of the methods to the initial
condition is a potential future research direction.

TABLE I: Performance of different methods Method | median MSE | min MSE | max MSE | Time/Run (s)
---|---|---|---|---
Approx. ONM | 0.00374 | 0.00108 | 0.0121 | 1.3
Exact ONM | 0.00470 | 0.00187 | 0.0319 | 13.0
SMC [3] | 0.00391 | 0.00092 | 0.0301 | 9.5
Figure 2: Box plots corresponding to the performance of different methods
after 1000 measurements over 100 scenarios.

Figure 3: Parameter estimates for approximate ONM (left) and exact ONM (right)

##  VI Conclusion

An approach based on online logistic regression for the problem of field
estimation using binary measurements has been proposed. The developed
algorithm is less computationally intensive than an existing sequential Monte
Carlo approach, while having more reliable estimation performance.

## References

  * [1] H. M. La, W. Sheng, and J. Chen, “Cooperative and active sensing in mobile sensor networks for scalar field mapping,” _IEEE Trans. Syst., Man, Cybern., Syst._ , vol. 45, no. 1, pp. 1–12, Jan. 2015.
  * [2] R. A. Razak, S. Sukumar, and H. Chung, “Scalar field estimation with mobile sensor networks,” _Int. J. Robust Nonlinear Control_ , vol. 31, pp. 4287–4305, 2021.
  * [3] A. S. Leong and M. Zamani, “Field estimation using binary measurements,” _Signal Processing_ , vol. 194, no. 108430, May 2022.
  * [4] K. Rosser, K. Pavey, N. FitzGerald, A. Fatiaki, D. Neumann, D. Carr, B. Hanlon, and J. Chahl, “Autonomous chemical vapour detection by micro UAV,” _Remote Sensing_ , vol. 7, pp. 16 865–16 882, 2015.
  * [5] B. Ristic, A. Gunatilaka, and R. Gailis, “Achievable accuracy in Gaussian plume parameter estimation using a network of binary sensors,” _Information Fusion_ , vol. 25, pp. 42–48, 2015.
  * [6] ——, “Localisation of a source of hazardous substance dispersion using binary measurements,” _Atmospheric Environment_ , vol. 142, pp. 114–119, 2016.
  * [7] D. D. Selvaratnam, I. Shames, D. V. Dimarogonas, J. H. Manton, and B. Ristic, “Co-operative estimation for source localisation using binary sensors,” in _Proc. IEEE Conf. Decision and Control_ , Melbourne, Australia, Dec. 2017, pp. 1572–1577.
  * [8] “Sequential Monte Carlo methods in practice,” A. Doucet, N. de Freitas, and N. Gordon, Eds. New York: Springer, 2001.
  * [9] C. M. Kreucher, A. O. Hero, K. D. Kastella, and M. R. Morelande, “An information-based approach to sensor management in large dynamic networks,” _Proc. IEEE_ , vol. 95, no. 5, pp. 978–999, May 2007.
  * [10] B. Ristic, M. Morelande, and A. Gunatilaka, “Information driven search for point sources of gamma radiation,” _Signal Processing_ , vol. 90, pp. 1225–1239, 2010.
  * [11] R. Jézéquel, P. Gaillard, and A. Duri, “Efficient improper learning for online logistic regression,” in _Conf. Learning Theory_ , Graz, Austria, Jul. 2020, pp. 2085–2108.
  * [12] A. Lesage-Landry, J. A. Taylor, and I. Shames, “Second-order online nonconvex optimization,” _IEEE Trans. Autom. Control_ , vol. 66, no. 10, pp. 4866–4872, Oct. 2021.
  * [13] J. Park and I. W. Sandberg, “Approximation and radial-basis-function networks,” _Neural Computation_ , vol. 5, no. 2, pp. 305–316, 1993.
  * [14] M. R. Morelande and A. Skvortsov, “Radiation field estimation using a Gaussian mixture,” in _Proc. Intl. Conf. Inf. Fusion_ , Seattle, USA, Jul. 2009, pp. 2247–2254.
  * [15] G. Calafiore and L. El Ghaoui, _Optimization Models_. UK: Cambridge University Press, 2014.
  * [16] R. A. Horn and C. R. Johnson, _Matrix Analysis_ , 2nd ed. UK: Cambridge University Press, 2013.
  * [17] G. G. Colabufo, P. M. Dower, and I. Shames, “Newton’s method: sufficient conditions for practical and input-to-state stability,” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 6334–6339, 2020.