# The mysterious magnetic ground state of Ba14MnBi11 is likely altermagnetic

Po-Hao Chang<EMAIL_ADDRESS>Igor I. Mazin<EMAIL_ADDRESS>Department of
Physics and Astronomy, George Mason University, Fairfax, VA 22030, USA
Quantum Science and Engineering Center, George Mason University, Fairfax, VA
22030, USA

###### Abstract

Mn-based transition metal Zintl compounds in the “14-1-11” phase are known to
host complex atomic and magnetic structures owing to their intricate crystal
structure. Among this family of compounds, Ba14MnBi11 stands out as one of the
least understood compounds, with experimental measurements and theoretical
findings largely inconsistent. Following up on the earlier attempt [D.
Sánchez-Portal et al., Phys. Rev. B 65, 144414 (2002)] at establishing a
connection between metallicity and magnetism through a DFT-based analysis, our
work aims to provide additional insights to resolve the existing
contradictions. Our key findings is that the magnetic ground state is very
susceptible to charge doping (band filling). Density functional calculations
for stoichiometric Ba14MnBi11 give a rather stable ferromagnetic metallic
ground state. However, by adding exactly one additional electron per Mn, the
system becomes semiconducting as expected, and consequently the magnetic
ground state becomes weakly antiferromagnetic (AF). On the other hand, upon
small hole doping the system transitions to a special type of AF state known
as altermagnetic ordering. The observed trends suggest that hole and electron
doping-induced phase transitions likely result from different underlying
mechanisms, influencing various exchange pathways. Additionally, our projected
density-of-states along with bandstructure analyses indicate that, besides the
largest hole contribution coming from the tetrahedral unit of Bi, the isolated
Bi sites also play a substantial role and the dispersive bands near VBM
suggest a rather complex hybridization network involving both Bi band
characters. Through a comprehensive comparison of available data and our
analysis, we propose that the inconsistency in magnetic states between
experimental findings and DFT calculations is not a failure of DFT, but rather
indicates nonstoichiometric effects, likely impurities or defects,
contributing to the observed discrepancies. A possibility of stabilizing,
through doping, an altermagnetic state, is exciting.

##  I Introduction

Zintl compounds with the general formula A14MnPn11, where A is an alkaline
earth, an alkaline metal, or a trivalent rare earth, and Pn is a pnictogen,
have been attracting attention as potential thermoelectrics (especially
Yb14MnPn11) for a few decades. These were intensively studied in 1992-2002,
but larger forgotten after that. Yet, their physics, and especially their
magnetic properties remain very poorly understood, which is especially true
about Ba14MnBi11, as summarized in Ref. [1]. Below we list the most enigmatic
experimental and computational findings [2, 3, 1, 4]

1\. Magnetic susceptibility is clearly antiferromagnetic (AF), but DFT
calculations, with large confidence, give ferromagnetism (FM) as the ground
state.

2\. Yet, the Curie-Weiss temperature is positive (ferromagnetic). This is a
common phenomenon in A-type antiferromagnets consisting of strongly coupled FM
planes, with a weak AF interaction between them. However, the crystal
structure of Ba14MnBi11, in its space group I41acd, is rather three-
dimensional and, at first glance, highly symmetric: each Mn has 10 nearest
neighbors in all directions, with four Mn-Mn bonds forming a tetrahedron,
another four exactly the same tetrahedron, but flipped, and two more forming a
linear structure. Geometrically, it seems that the structure cannot be easily
partitioned into weakly coupled FM sublattices.

3\. The Curie-Weiss effective magnetic moment is 4.8 $\mu_{B}$, consistent
with the spin $S\approx 2$ ($M=4$ $\mu_{B}$), but the saturation magnetization
was claimed from the experiment in the field up to 5 T is $\approx 3$
$\mu_{B}$. At the same time, DFT calculations, again with a great degree of
confidence show Mn to be in a divalent high-spin state, i.e., $M=5$ $\mu_{B}$.

4\. Calculations give a metallic ground state, with a few Bi-derived bands at
the Fermi level, and no Mn bands near $E_{F}$, consistent with Mn2+. It was
claimed to be metallic in the experiment as well, but in the only transport
measurement [3] the resistivity is barely dependent on temperature, with
Residual Resistivity Ratio (RRR) about 1.3, and the resistivity being weakly
metallic below the Neél transition.

5\. Last but not least, the linear specific heat coefficient corresponds to
the density of states of 44 states/eV$\cdot$formula, while the calculations
give, for a sizeable range of Fermi energies, less than 7.5
states/eV$\cdot$formula.

So, in and by itself it is an intriguing compound with many unresolved
mysteries. In addition, there is no experimental data on the magnetic pattern,
except that it is clearly antiferromagnetic. As pointed out in Ref. [1], there
are three distinct-symmetry magnetic patterns even for $q=0$, which they
called AFM(a-c). In retrospect, one of these patterns, namely AFM(b) is
altermagnetic (AM), a novel type of ordered collinear magnetism that has been
intensively discussed in the last few years[5, 6, 7, 8]

In this paper, we will address some of the unresolved issues either exact or
speculative, assuming that the experimentally available samples have a small
hole doping of the order of 1 hole per 4 formula units (F.U.), for instance,
due to 1% vacancies on Ba site. Interestingly, at the same time, doping
stabilized the AM order. The paper is organized as follows. Section II
provides an overview of the crystal structure. The computational details are
described in Section III. Section IV presents our data for the doping effect
on the ground states and how two different AF orderings can be induced.
Section V discusses the doping effect on the exchange coupling using two
different methods. Section VI analyzes the band structure and the projected
density-of-states near the Fermi energy to establish a connection between band
characters and the doping effect on the magnetic structure. Finally, the
conclusions are summarized in Section VIII

##  II Crystal structure

Figure1:  The crystal structure views with different emphases: (a) The
schematic view of the first three nearest exchange paths defined in the
unitcell of four formula unit. (b) Depiction of the model magnetic cell
showing four magnetic sites and how two inequivalent exchange coupling
constants $J_{1}$ and $J_{1^{\prime}}$ are spatially oriented. (c) Extended
structural view revealing four inequivalent Bi sites organized into
tetrahedral units, trimers, and isolated Bi sites. (d) Additional symmetry
feature: The view illustrates that $J_{1}$ bonds have an inversion center at
the midpoints and $J_{1^{\prime}}$ do not.

The compound has a body-center-tetragonal space group I41acd unit cell with
four formula units as depicted in Fig. 1 (a) with additional Mn sites (blue)
included to show $J_{2}$. There are four inequivalent Bi sites shown in
different colors as defined in the figure. The lattice parameters
$a=18.665(3)$ and $c=24.429$ Å are taken from the original experimental work
[3].

While all Mn are symmetry equivalent, there are two inequivalent 1st NN
exchange coupling interactions: $J_{1}$ (red) that passes through two vertices
from two tetrahedra and $J_{1^{\prime}}$ (blue) that do not. More importantly,
as Fig. 1 (d) shows, the $J_{1}$ bonds have an inversion center at the
midpoints and $J_{1^{\prime}}$ do not.

A more complete view of how $J_{1}$ and $J_{1^{\prime}}$ are spatially
oriented for each Mn site is presented in Fig. 1 (b). All four Mn sites are
labeled with numbers to indicate their order for later discussion.

Fig. 1 (c) shows the extended view of the structure that exhibits a
complicated tetrahedra network. Due to the large number of atoms in the unit
cell, for better visibility, Ba sites are not displayed. It becomes clear that
the crystal structure per F.U. consists of one linear trimer
(Bi${}_{3}^{7-}$), four isolated Bi3- and a tetrahedral unit
(MnBi${}_{4}^{9-}$). The tetrahedra, as explained in Ref. [1], are
interconnected through either $J_{1}$ and $J_{1}^{\prime}$, two inequivalent
first NN exchange coupling interactions, which form two different networks.

##  III Methods

The small energy scale of the magnetic interactions and the large crystal
structure require considerable computational flexibility; since different DFT
codes are best at addressing different aspects of the calculations, we
employed three DFT codes of different types of basis sets to ensure accuracy
as well as provide more insight. As a reference, the electronic structure
calculations were first performed to obtain the total energies using Vienna ab
initio Simulation Package (VASP) [9] within projector augmented wave (PAW)
method.[10, 11] The Perdew-Burke-Enzerhof (PBE) [12] generalized gradient
approximation was employed to describe exchange-correlation effects. To
improve the description for localized d-electrons in Mn2+ ion to be strongly
correlated, we added a Hubbard $U$ correction with the fully localized limit
double-counting recipe[13, 14]. The effective parameter $U-J=5$ eV was used.

Additionally, as an alternative to SIESTA [15] used in Ref [1], a similar
numerical-orbital-based [16] DFT code OpenMX [17], designed for large-scale
simulation was used for the calculations of magnetic properties. In these
calculations, core electrons are replaced with norm-conserving pseudopotential
[18, 19] and a sufficiently large numerical atomic orbital set is chosen as
Mn6.0-$s3p2d2f1$, Ba8.0-$s3p2d2$ and Bi8.0-$s3p2d2f1$ to accurately determine
the band structure. As Fig. 7 shows, the band structures calculated using both
codes are in good agreement.

The exchange coupling parameters are calculated using two different methods.
We first use the standard energy-mapping approach following Ref. [1] to
calculate exchange coupling constants. Due to the computation cost, and large
distances between Mn sites, only the interactions up to 2nd nearest neighbors
(NN) are considered. In this method, the DFT energies for the four
configurations (i.e. FM and three AF) are fitted onto the following spin
Hamiltonian

$H=-\sum_{\left\langle ij\right\rangle_{1}}J_{1}s_{i}s_{j}-\sum_{\left\langle
ij\right\rangle_{1^{\prime}}}J^{\prime}_{1}s_{i}s_{j}-\sum_{\left\langle
ij\right\rangle_{2}}J_{2}s_{i}s_{j},$ (1)

where $s$ are the normalized moments, $(\mathbf{s}=\mathbf{S}/|S|$,
$|\mathbf{s}|=1$), $J_{1}$, $J_{1^{\prime}}$ and $J_{2}$ are the exchange
parameters of two inequivalent 1st NN, and 2nd NN defined in Fig. 1 (a), and
summation is over all different bonds. Note that in our definition, $J>0$
indicates a FM interaction. This method will be applied exclusively to the
data obtained using VASP and will be referred to as the “energy difference
method” in the subsequent discussion.

To gain more insight, we also employed the Green’s function method [20, 21]
implemented in OpenMX 3.9[22]. In this approach, the exchange interaction
between any pair of given magnetic sites can be directly calculated from a
single magnetic state for any given interionic distance. Naturally, this
approach can only be directly applied to OpenMX as it requires the Hamiltonian
to be in local orbital representation.

Finally, for historical reasons, as well as in order to compare with an all-
electron method, we used the Linear Augmented Wave method Wien2K [23] for the
Fermi surface properties, such as the Fermi surface cuts and the plasma
frequencies as a function of the Fermi energy.

##  IV Doping effect on total energy

Figure2: Relative energy of AF configurations compared to FM as a function of
the charge doping per unit cell (4 F.U.). Negative charge represents hole
doping.

Fig. 2 shows the relative energies compared to the FM state $uuuu$ as a
function of doped charge ranging from light hole doping up to exactly one
electron per Mn (four in total) for $udud$, $uudd$ and $uddu$ three AF
orderings, where $u$ $(d)$ represents the up (down) spin moment on the Mn site
and the structure and the order of Mn are defined in Fig. 1. The positive
(negative) charge in the $x$-axis indicates electron (hole) doping. This
convention is also applied to the rest of the discussion.

Generally, the relative energies for $uudd$ and $uddu$, in contrast to that of
$udud$, share a very similar pattern and both are very sensitive to doping. In
the light doping region near stoichiometry ($-1$ to 2), the trends of all
three curves are rather consistent. Both $uudd$ and $uddu$ vary monotonically
and almost linearly with very steep slopes as the energies change from around
$-10$ to 30 meV within the region while $udud$ does not show any obvious
dependence and the energy stays at roughly $9$ meV above FM.

Without doping, all three states are very close in energy and FM is the ground
state. The energy differences then begin to widen rapidly as electrons or
holes are introduced into the system. With light electron doping, the ground
state stays FM. However, when a small fraction of holes is introduced into the
system, the GS quickly transitions into $uudd$, a very interesting type of AF,
known as AM. [5, 6, 7, 8]. Indeed, as Fig. 1(d) shows, if the $J_{1}$ bond is
antiferromagnetic, there is an inversion operation connecting the two
sublattices, which is a signature of a trivial AF state [24] , while if the AF
bonds are $J_{1}^{\prime}$, the structure becomes AM.

As more electrons are added to the system, a qualitative change of behaviors
can be seen in all three cases around $e=2$ (0.5 $e$ per Mn). Above this
point, the energies for all three states begin to decrease monotonically until
$e=4$, where the Fermi level reaches the gap. The system hence becomes
semiconducting and transitions into a very weakly coupled AF state in $udud$
phase.

In this gapped system all four magnetic states (including FM) are nearly
degenerate with only about 1–2 meV difference in total energy despite the
large supercell. This is due to the fact that the original exchange pathways
mediated by long-range RKKY interaction in the metallic state are no longer
available. As a result the magnetic interactions between Mn sites separated by
a large distance are now instead governed mainly by the much shorter-range
superexchange mechanism, which tends to favor AF. This $udud$ phase induced by
electron doping, as previous suggested, corresponds to the configuration that
is ferromagnetic within the same network formed by the tetrahedra units are
connected through $J_{1}$ and antiferromagnetic with adjacent network
connected by $J_{1}^{\prime}$[1].

##  V Doping effect on exchange coupling

Figure3:  The first few nearest exchange coupling parameters as functions of
charge doping calculated using (a) Greens function method and (b) total energy
difference.

Fig. 3 shows the exchange coupling parameters as functions of doping for the
first few nearest NN using both (a) GF function method and (b) total energy
difference method. For the latter, due to the computation limitation by the
system size, only up to $J_{2}$ are calculated.

Generally, some differences between the two methods are expected, as in the
case of the energy difference method, the further NN interactions which could
still contribute in the RKKY-driven regime are combined into the nearest three
exchange parameters, while GF is a perturbation method on a single determinant
which allows the calculation between any given pair of any distance. Also,
these two methods differ in that the perturbative Green function method
corresponds to small deviations of collinearity and the total energy method
assumes full spin flips; the two may physically differ if the magnetic
Hamiltonian deviates from the Heisenberg form

Despite the difference, several key features are rather consistently predicted
by both methods and offer useful insights into understanding the effect of
doping. Most importantly, both suggest that $J_{1}$ is the most susceptible
component to doping and its sign and strength vary greatly depending on the
type of doping while $J_{1}^{\prime}$ and $J_{2}$ change more gradually. Other
qualitative behaviors such as the crossover between the $J_{1}^{\prime}$ and
$J_{2}$ as well as the sign change of $J_{1}$ in the hole doping region are
also captured.

Both $uudd$ and $uddu$ respond to doping in a very similar way suggest that
they share the same dominant mechanism that is very different from that of
$udud$. Indeed, from both Figs 3 (a) and (b), one can see that the qualitative
behaviors of $uudd$ and $uddu$ closely resemble that of $J_{1}$ which connect
opposite spins in both states indicating that $J_{1}$ is the main driving
mechanism of the dependency.

Although other major components such as $J_{1}^{\prime}$ and $J_{2}$ are also
significant in terms of the strength, the effects from $J_{1}^{\prime}$ and
farther NN interaction $J_{2}$ roughly cancel out. This cancellation is mainly
due to the fact that, while the strength of $J_{2}$ is about twice bigger than
$J_{1}^{\prime}$ but, for any given Mn site, the number of its 2nd NN is only
of half its first NN connected through $J_{1}^{\prime}$. This explains why
$uudd$ and $uddu$ nearly degenerate in electron doping cases.

However, in the hole doping case, this cancellation no longer holds, as
$J_{1}^{\prime}$ continues to increase with more hole carriers introduced into
the system, while $J_{2}$ reaches the peak and begin to decline. As a result,
in Fig. 2, one can see the energy gap between $uudd$ and $uddu$ begins to
widen and $uudd$ becomes the ground state. It is interesting to note that,
despite some discrepancy in the qualitative behavior, both Figs 3 (a) and (b)
exhibit the same cancellation in charge doping region and the cancellation is
removed due to hole doping.

On the other hand, in $udud$, the energy term that involves $J_{1}$, the most
sensitive component to doping, has the same sign as that of the FM ordering,
therefore the doping dependence of the energy difference between the two, is
mainly influenced by other exchange paths, such as less sensitive
$J_{1}^{\prime}$ and $J_{2}$ .

Figure4: Exchange coupling constants as functions of distance for different
amounts of doping calculated using Green’s function method.

Fig. 4 shows the exchange coupling parameters as functions of distance for
Bi14MnBi11 with different amounts of doping, calculated using the Green’s
function method. The exchange coupling parameters in all three metallic states
($0$, $0.5$ and $1$) possess long-range (there are still significant
contributions beyond 20 Å) and sign alternating behavior which are both
signatures of the RKKY interaction where the localized $3d$ Mn moments are
coupled through conduction electrons. This is likely the reason for the
unusually high $T_{N}$ given such large spacing between the magnetic Mn sites.

For $e=+4$, the system in a semiconducting state as mentioned earlier, likely
governed by super-exchange, has only much weaker and shorter range
interactions that favor AF, as shown in Fig. 4(d) all $J_{1-3}$ are much
smaller, and essentially negligible beyond $J_{3}$.

However, it is important to note that, while there are still significant
interactions beyond 3NN ($\geq 14$ Å) in the metallic states, most of the
farther neighbors are rather insensitive to doping as shown in Figs. 4(a)-(c).
It is therefore reasonable to assume that the qualitative behaviors due to
doping are mainly dictated by $J_{1}$, $J_{1}^{\prime}$ and $J_{2}$.

Figure5:  Mean-field $T_{CW}$ and $T_{N}$ as functions of doping using (a)
total energy (VASP) and (b) Green’s function method (OpenMX)

###  $T_{CW}$ and $T_{N}$

|  $J_{1}$ |  $J_{1^{\prime}}$ |  $J_{2}$
---|---|---|---
CW |  $+$ |  $+$ |  $+$
FM |  $+$ |  $+$ |  $+$
$udud$ |  $+$ |  $-$ |  $-$
$uudd$ |  $-$ |  $+$ |  $-$
$uddu$ |  $-$ |  $-$ |  $+$
Table 1: The sign ($\sigma=s_{i}s_{j}$) indicates the alignment between the
magnetic moments connected by a given $J$ for each magnetic state.

All Mn sites are equivalent, and each Mn is connected to eight first NN: four
through $J_{1}$ and four through $J_{1^{\prime}}$, as well as two second NN
through $J_{2}$. The critical temperatures $T_{c}$ (c=CW and N for Curie-Weiss
and Nèel temperature respectively) in the mean-field approximation are
directly related to the strength of exchange interactions $J_{i}$ and can be
expressed in the following simple form

$T_{c}=\frac{2}{3k_{B}}(4J_{1}\sigma_{1}+4J_{1^{\prime}}\sigma_{1^{\prime}}+2J_{2}\sigma_{2}),$
(2)

where $k_{B}$ is the Boltzmann constant and $\sigma=\bm{s_{i}s_{j}}$ are the
signs, as summarized in TABLE 1, indicate the alignment between the magnetic
moments connected by the given $J_{i}$ in each magnetic state. Figure 5 shows
the mean-field $T_{N}$ for the three AF states and $T_{CW}$ estimated from the
exchange parameters obtained using (a) Green’s function method and (b) total
energy difference method as presented in Fig. 3. It is worth noting that the
Curie temperature $T_{C}$ for $uuuu$ has an identical expression as $T_{CW}$
in the mean-field approximation. Despite some small deviations in $J_{2}$,
both methods agree reasonably well suggesting that the Mn moments are fairly
localized. Generally, the state with the highest critical temperature ($T_{N}$
or $T_{C}$) corresponds to the most stable and energetically favorable state.
This becomes clear when compared with Figs. 2 and 3.

Light electron doping stabilizes the FM state leading to a positive and
increasing $T_{CW}$. At $e=4$ all temperatures become very small as a result
of weak couplings in semiconducting states, with $udud$ being the ground state
and $T_{N}^{udud}$ having the largest value. On the other hand, light hole
doping has the opposite effect promoting $uudd$ while destabilizing the FM
state. Given the close resemblance between $uudd$ and $udud$ shown in Fig. 2,
it is not surprising that their $T_{N}$ values also have the same dependence
on doping, despite $uudd$ being more stable.

Most importantly, near the hole doping concentration where $uudd$ (AM) is
predicted to be the ground state (i.e. $e=1.0$ in Fig. (a) and $e=0.4$ in Fig.
(b)), two curves $T_{CW}$ and $T_{N}^{uudd}$ intersect at around $T=15$ K.
This naturally implies that both temperatures have the same sign and are very
close in magnitude ($\sim$15 K) in the proximity of the intersection, as
indicated by the dashed circles in both Figs. 5 (a) and (b). This is indeed
supported by the experimental values showing $T_{N}=15$ K and $T_{CW}=17$ K
[3]. This finding is rather intriguing because a positive $T_{CW}$ suggests a
FM state which contradicts a positive $T_{N}$ that indicates an AF state. As
introduced earlier, this phenomenon is common in the antiferromagnets
consisting of strongly coupled FM planes with a weak AF interaction between
them.

Although, at first glance, the highly symmetric crystal structure of
Ba14MnBi11 appears three-dimensional and does not allow such type of planer
coupling, the magnetic moments, instead of planar, can in fact be visualized
as two interpenetrating 3D networks [1].

The discrepancy can then be easily explained if the Mn moments within the same
network, coupled through $J_{1^{\prime}}$, strongly favor the FM ordering,
meanwhile a weak AF interaction, coupled through $J_{1}$, exists between
networks. This explanation is consistent with our data shown in Fig. 3 further
supporting the argument made in previous work [1]. However, it is important to
note that, the definition of Mn network here is different from that in Ref.
[1], where they consider the Mn tetrahedra complexes connected by $J_{1}$ to
belong to the same network.

##  VI Exchange pathway and band characters

Due to the large distances that separate the Mn atoms, the behaviors of the
exchange couplings are expected to be heavily influenced by the detailed
electronic structure of Bi-$p$ bands as they are the major components near the
Fermi level [1]. To gain more insight into the underlying mechanism, it is
helpful to first look into the exchange pathway and identify the band
characters associated with it.

Fig. 1(a) shows the definitions of the three dominant exchange parameters
$J_{1}$, $J_{1^{\prime}}$ and $J_{2}$ along with the relavant surrouding
elements displayed. One can see a clear distinction between the three. For
$J_{1}$, the path directly goes through two Bi3 sites while $J_{1^{\prime}}$
does not overlap with any sites directly. However, its exchange pathway likely
involves two Bi4 atoms near the path that form a plane with the two Mn sites
connected by $J_{1^{\prime}}$. The Next NN $J_{2}$ has a direct overlap with
the Bi1 at the midpoint which is likely to involve the entire trimer (Bi1 and
Bi2). As the Bi band consists of four different inequivalent Bi-sites as
introduced in Section II, one can immediately notice that each pathway can be
associated with at least one distinct type of Bi character that is directly
involved.

Figure6: Projected density of state of four inequivalent Bi sites near the
Fermi level.

Fig.6 shows the projected density of states (pDOS) for all four inequivalent
Bi sites near the Fermi level and valence band maximum (VBM). Generally, the
Bi sites from the tetrahedral units have the largest contribution around the
Fermi level and host the biggest portion of the hole concentration (near VBM).
This explains why $J_{1}$ is the most susceptible parameter to band filling.
In addition, there is also a significant amount of isolated Bi (Bi2). In
comparison, the trimer characters are less significant near the Fermi level.
This configuration of orbital characters is not uncommon. A similar
arrangement is also seen in other Zintel 14-11-1 compounds, as the
bandstructures closely resemble one another [25].

It is worth mentioning that, our estimation of DOS for FM state at $E_{F}$ is
approximately 7.0 states/eV$\cdot$formula which is very close to the earlier
theoretical value (7.8 states/eV$\cdot$formula) [1]. The large difference
compared to the experimental value ($\sim 44$ states/eV$\cdot$formula) [26]
has been attributed to the effect of hyperfine splitting not being considered
However, the mechanism behind the unexpectedly large observed DOS remains
unclear.

Figure7:  The bandstructure of Ba14MnBi11 in FM states for (a) majority and
(b) minority spins using VASP (red) and OpenMX (blue).

Fig. 7 shows the bandstructure near the VBM for Ba14MnBi11 in FM state for
both (a) majority and (b) minority spins [calculated using VASP (red) and
OpenMX (black) as an additional check. The two methods yield remarkable
agreement].

Similar to the previous findings [25, 1], near the Fermi level a few rather
dispersive light bands emerge from the flatter valence at $\Gamma$ point both
seen in (a) and (b), indicates a p-type conductor. As shown in the pDOS
analysis Fig. 6, these bands near VBM are mainly made up of the Bi characters
from the tetrahedra unit with significant contribution from the isolated Bi
sites. The relative high discursiveness of the hole pockets around $\Gamma$
point suggests a more complicated hybridization network most likely between
the two dominant characters.

##  VII Electronic and transport properties

In the transport measurement, the temperature dependence of RRR indicates
metallic behavior, as resistivity decreases with the reducing temperature,
despite the high resistivity observed. In the original work, this high
resistivity is attributed to the polycrystal sample [3]. However, our
transport calculations suggest that the compound inherently has poor
conductivity regardless of sample quality. Figure 8 shows the plasma frequency
as a function of Fermi energy. The values calculated along both $x$ (red) and
$z$ (green) directions are generally smaller than 0.5 eV2, which is at least
an order of magnitude smaller than in typical metals. In the Drude
approximation, the static conductivity $\sigma_{0}$ is determined by

$\sigma_{0}=\omega^{2}_{p}{\tau}/{4\pi}$ (3)

where $\tau$ and $\omega_{p}$ are the relaxation time and the plasma
frequency, respectively. Such a small $\omega_{p}^{2}$ inherently leads to a
small conductivity, although a short relaxation time can also contribute.

Figure8:  The plasma frequency as a function of the Fermi level calculated
along $x$ and $z$ directions.

To better understand magnetism in Ba14MnBi11, it is also important to
highlight the intricate connection between the magnetic and electronic
structures. As explained in detail in Ref. [1], a comparison with the
electronic structure of a similar compound Ca14GaAs11 (where Ga is identified
to be trivalent) [27], reveals that if Mn were in a trivalent state, the
system would have been semiconducting [27], which contradicts the transport
measurements [3]. To be consistent with a metallic band structure, Mn should
instead be in the divalent state (Mn2+) $d^{5}$ configuration, which is
consistent with the DFT calculations. The deviation from the experimentally
observed net magnetization of $\sim 4$ $\mu_{B}$ is attributed to the opposite
polarization of the Bi atoms that form tetrahedral units with Mn. Therefore,
the entire MnBi${}_{4}^{-}$ complex should be treated as one magnetic unit in
order to explain the discrepancy. This conclusion is further supported by our
Voronoi charge analysis, which shows a clear $d^{5}$ Mn state ($\sim
4.8\mu_{B}$) and the opposite polarization on Bi sites ($\sim 0.03-0.05$
$\mu_{B}$/Bi), consistent with the conclusions of Ref.[1].

Intuitively, the FM order obtained from DFT is related to the metallicity of
the stoichiometric compound; indeed, adding one additional electron per F.U to
the system opens a semiconducting gap and the ground states transitions into a
weakly coupled AF state that corresponds to the state proposed in Ref. [1].
However, this does not account for the observed metallicity (albeit this
observation may be questioned) and the weak interactions hardly account for
the $T_{N}$ observed.

On the other hand, we found that the magnetic states are susceptible to band
filling particularly near stoichiometry. So, through a small amount of hole
doping, we were able to obtain another type of AF state known as
altermagmetic, while preserving the metallic conductivity. We therefore
believe that the discrepancy in the magnetic state between DFT calculations
and experimental findings stems from nonstoichiometry and intrinsic small hole
doping and the AF state observed in the experiment is in fact altermagnetic.

Figure 9 shows the cross-section of Fermi surfaces at $k_{z}=0$ of the AM
state. The red (blue) color represents the up (down) spin. The band splits in
a compensated magnetic order is the signature of altermagnetism.

Figure9:  The cross-section of Fermi surfaces at $k_{z}=0$ in the
altermagnetic state. The red and blue contours represent up and down spins,
respectively.

##  VIII Conclusions and considerations regarding altermagetism

In this study, we conducted an extensive and comprehensive examination of the
magnetic properties of Ba14MnBi11 using several different Density Functional
Theory (DFT) codes. We have successfully established a connection between the
electronic and magnetic structures, and were able to explain some of the
unresolved issues. Our investigation revealed a significant sensitivity of the
first few nearest exchange coupling interactions to band filling. With either
electron or hole doping, we observed the emergence of distinct
antiferromagnetic (AF) states. By adding exactly one electron per F.U., which
fills the valence hole pocket, the system transitions into an AF state where
the moments are ferromagnetically coupled through $J_{1}$ and
antiferromagnetically coupled through $J_{1^{\prime}}$. Conversely, the
introduction of light hole doping results in a unique type of AF state known
as altermagnetism, which simultaneously also preserves the metallic state
observed in the experiment. In this case, the Neél and Curie-Weiss
temperatures obtained from mean-field approximation are both in excellent
agreement with the experimental data. Based on our analysis, the AF observed
state is expected to be altermagnetic and the discrepancy in the magnetic
ground state between the DFT calculations and experimental observation is
likely due to deviations from the stoichiometry such as Ba vacancies.

Importantly, the predicted altermagnetic state has a number of interesting
properties suggesting that if better quality samples can be manufactured they
can provide considerable advantage over existing altermagnetic candidates.
First, it belongs to the $d-$wave class, per classification of Ref. [8], which
reduces the number of vertical nodal lines (compared to the two other
convincingly established semiconducting/metallic altermagnets, MnTe and CrSb).
Most importantly, despite its relatively high tetragonal symmetry, it has non-
zero anomalous transport for all possible orientations of the Néel vector.
Specifically, in the case of the easy-axis anisotropy, $\mathbf{M}||c$,
anomalous conductivity $\sigma_{xz}=\sigma_{xz}\neq 0$ is allowed, and for the
easy plane anisotropy, $\mathbf{M}\perp c$, all three components are non-zero,
$\sigma_{xz}=\sigma_{xz}\neq\sigma_{xy}\neq 0$. This property is quite unique
among altermagnetic candidates. The low Néel temperature of 15 K precludes
practical application, but this material remains an interesting candidate to
study altermagnetic physics.

###### Acknowledgements.

This work was supported by the Army Research Office under Cooperative
Agreement Number W911NF-22-2-0173.

## References

  * Sánchez-Portal  _et al._ [2002] D. Sánchez-Portal, R. M. Martin, S. M. Kauzlarich, and W. E. Pickett, Bonding, moment formation, and magnetic interactions in Ca14MnBi11 and Ba14MnBi11, Phys. Rev. B 65, 144414 (2002).
  * Webb  _et al._ [1991] D. J. Webb, T. Y. Kuromoto, and S. M. Kauzlarich, New ternary magnets (Ca, Sr, Ba)14MnBi11, Journal of Magnetism and Magnetic Materials 98, 71 (1991).
  * Kuromoto  _et al._ [1992] T. Y. Kuromoto, S. M. Kauzlarich, and D. J. Webb, Structure and properties of the transition-metal zintl compounds: A14MnBi11 (A = Ca, Sr, Ba), Chem. Mater. 4, 435 (1992).
  * Kauzlarich  _et al._ [2004] S. M. Kauzlarich, A. C. Payne, and D. J. Webb, Magnetism and Magnetotransport Properties of Transition Metal Zintl Isotypes, in  _Magnetism: Molecules to Materials_, edited by J. S. Miller and M. Drillon (Wiley, 2004) 1st ed., pp. 37–62.
  * Mazin [2023] I. I. Mazin, Altermagnetism in MnTe: Origin, predicted manifestations, and routes to detwinning, Phys. Rev. B 107, L100418 (2023).
  * Mazin and The PRX Editors [2022] I. Mazin and The PRX Editors, Editorial: Altermagnetism—A New Punch Line of Fundamental Magnetism, Phys. Rev. X 12, 040002 (2022).
  * Šmejkal  _et al._ [2022a] L. Šmejkal, J. Sinova, and T. Jungwirth, Emerging Research Landscape of Altermagnetism, Phys. Rev. X 12, 040501 (2022a).
  * Šmejkal  _et al._ [2022b] L. Šmejkal, J. Sinova, and T. Jungwirth, Beyond Conventional Ferromagnetism and Antiferromagnetism: A Phase with Nonrelativistic Spin and Crystal Rotation Symmetry, Phys. Rev. X 12, 031042 (2022b).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996), publisher: American Physical Society.
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994), publisher: American Physical Society.
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999), publisher: American Physical Society.
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. 77, 3865 (1996), publisher: American Physical Society.
  * Liechtenstein  _et al._ [1995] A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Density-functional theory and strong interactions: Orbital ordering in Mott-Hubbard insulators, Phys. Rev. B 52, R5467 (1995), publisher: American Physical Society.
  * Dudarev  _et al._ [1998] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study, Phys. Rev. B 57, 1505 (1998), publisher: American Physical Society.
  * Artacho  _et al._ [1999] E. Artacho, D. S?nchez-Portal, P. Ordej?n, A. Garc?a, and J. Soler, Linear-Scaling ab-initio Calculations for Large and Complex Systems, phys. stat. sol. (b) 215, 809 (1999).
  * Ozaki [2003] T. Ozaki, Variationally optimized atomic orbitals for large-scale electronic structures, Phys. Rev. B 67, 155108 (2003), publisher: American Physical Society.
  * [17] OpenMX, The code OPENMX pseudoatomic basis functions, and pseudopotentials are available under terms of the GNU-GPL on a web site, http://www.openmx-square.org/ .
  * Vanderbilt [1990] D. Vanderbilt, Soft self-consistent pseudopotentials in a generalized eigenvalue formalism, Phys. Rev. B 41, 7892 (1990), publisher: American Physical Society.
  * Morrison  _et al._ [1993] I. Morrison, D. M. Bylander, and L. Kleinman, Nonlocal Hermitian norm-conserving Vanderbilt pseudopotential, Phys. Rev. B 47, 6728 (1993), publisher: American Physical Society.
  * Katsnelson and Lichtenstein [2000] M. I. Katsnelson and A. I. Lichtenstein, First-principles calculations of magnetic interactions in correlated systems, Phys. Rev. B 61, 8906 (2000).
  * Antropov  _et al._ [1997] V. Antropov, M. Katsnelson, and A. Liechtenstein, Exchange interactions in magnets, Physica B: Condensed Matter 237-238, 336 (1997).
  * Terasawa  _et al._ [2019] A. Terasawa, M. Matsumoto, T. Ozaki, and Y. Gohda, Efficient Algorithm Based on Liechtenstein Method for Computing Exchange Coupling Constants Using Localized Basis Set, J. Phys. Soc. Jpn. 88, 114706 (2019).
  * P. Blaha  _et al._ [2002] P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka, and J. Luitz, WIEN2K (2002), pages: 2201 Publication Title: J. Phys. Chem. Sol. Volume: 63.
  * Smolyanyuk  _et al._ [2024] A. Smolyanyuk, L. Šmejkal, and I. I. Mazin, A tool to check whether a symmetry-compensated collinear magnetic material is antiferro- or altermagnetic (2024), arXiv:2401.08784 [cond-mat].
  * Liu  _et al._ [2021] Y. Liu, M. Y. Toriyama, Z. Cai, M. Zhao, F. Liu, and G. Jeffrey Snyder, Finding the order in complexity: The electronic structure of 14-1-11 zintl compounds, Applied Physics Letters 119, 213902 (2021).
  * Siemens  _et al._ [1992] D. P. Siemens, J. Del Castillo, W. Potter, D. J. Webb, T. Y. Kuromoto, and S. M. Kauzlarich, Specific heat of the ternary Zintl compounds (Sr14MnBi11 and Ba14MnBi11), Solid State Communications 84, 1029 (1992).
  * Gallup  _et al._ [1992] R. F. Gallup, C. Y. Fong, and S. M. Kauzlarich, Bonding properties of calcium gallium arsenide, Ca14GaAs11: a compound containing discrete GaAs4 tetrahedra and a hypervalent As3 polyatomic unit, Inorg. Chem. 31, 115 (1992).


# Variational Inference of Parameters in Opinion Dynamics Models

Jacopo Lenti1,2 Fabrizio Silvestri1&Gianmarco De Francisci Morales2 1Sapienza
University
2CENTAI
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Despite the frequent use of agent-based models (ABMs) for studying social
phenomena, parameter estimation remains a challenge, often relying on costly
simulation-based heuristics. This work uses variational inference to estimate
the parameters of an opinion dynamics ABM, by transforming the estimation
problem into an optimization task that can be solved directly.

Our proposal relies on probabilistic generative ABMs (PGABMs): we start by
synthesizing a probabilistic generative model from the ABM rules. Then, we
transform the inference process into an optimization problem suitable for
automatic differentiation. In particular, we use the Gumbel-Softmax
reparameterization for categorical agent attributes and stochastic variational
inference for parameter estimation. Furthermore, we explore the trade-offs of
using variational distributions with different complexity: normal
distributions and normalizing flows.

We validate our method on a bounded confidence model with agent roles (leaders
and followers). Our approach estimates both macroscopic (bounded confidence
intervals and backfire thresholds) and microscopic ($200$ categorical, agent-
level roles) more accurately than simulation-based and MCMC methods.
Consequently, our technique enables experts to tune and validate their ABMs
against real-world observations, thus providing insights into human behavior
in social systems via data-driven analysis.

##  1 Introduction

Agent-based models (ABMs) are computational frameworks that simulate the
actions and interactions of autonomous agents (individuals or collective
entities) within a system. Their goal is to understand how complex phenomena
and properties of a system emerge from the individual behavior of agents.
These models integrate elements from various disciplines, including game
theory, complex systems, and computational sociology; and they are popular
across diverse scientific domains, such as biology, ecology, economics, and
the social sciences.

In ABMs, the agents possess individual attributes and act following predefined
rules of behavior and interaction, and system outcomes can be validated
against real-world data. For example, in the field of opinion dynamics,
Fortunato and Castellano (2007) show that a simple voting model can reproduce
the distribution of votes received by political candidates across several
countries. Slightly more complex models of the same nature have been used to
analyze the negotiations around the Paris Agreement on climate change
(Bernardo et al., 2021) and to study the emergence of ideological spaces in US
(Baumann et al., 2021). However, tuning the parameters of these models is
still a considerable challenge (Flache et al., 2017).

Especially in the context of opinion dynamics, both the global (macroscopic)
and the individual (microscopic) level ABM parameters are of great interest,
and this is the focus of the current work. While ABMs offer great flexibility
and a clear representation of the causal mechanisms that govern the system,
they pose significant challenges in the inference phase. The obstacles come
from several factors, such as the considerable computational effort required
to generate data (a complete simulation of the model), and the high
dimensionality of microscopic parameters, which scale with the number of
agents. Additionally, the heterogeneity of the agents and the diversity of the
rules in the ABM literature are an additional hindrance in the tractability of
these models (Platt, 2020).

Figure 1: The proposed parameter estimation pipeline. First, we translate the
agent-based model into a probabilistic generative agent-based model. Then, we
apply variational inference to get an approximate posterior of the target
parameters within a given dataset.

The current practice for parameter estimation in ABMs mostly focuses on
simulation-based methods. These methods compute parameter estimates by running
a large amount of computationally expensive simulations of the model with
different values of the parameters, and comparing the observed and simulated
traces of the system state (a process called _calibration_) (Fagiolo et al.,
2019; Lux and Zwinkels, 2018). Simulation-based approaches rely on summary
statistics of the system state, primarily focusing on the replication of some
macroscopic property, rather than comparing the microscopic behaviors of the
agents. These summary statistics cause an inevitable loss of information at a
granular scale, and potentially cause identification issues since different
parameterizations of the ABM can lead to indistinguishable summary statistics
(Gallo et al., 2022). Moreover, when the parameters are high-dimensional, an
exhaustive exploration of the parameter space is not feasible. For this
reason, simulation-based approaches have been mostly used to estimate low-
dimensional parameters, while ignoring the microscopic parameters. Indeed,
most of the research on ABM calibration focuses on ($i$) a proper estimation
choice, following the comparison of the simulated and observed summary
statistics (Grazzini and Richiardi, 2015), ($ii$) an effective representation
of the data trajectories (Lamperti et al., 2018), and ($iii$) an efficient
exploration of the parameter space (Quera-Bofarull et al., 2023).

A recent line of work has suggested tackling the parameter estimation problem
in ABMs in a radically different way: via a likelihood-based approach (Monti
et al., 2020, 2022; Lenti et al., 2024). They propose a paradigm shift, by
translating the ABMs into probabilistic generative models, called
Probabilistic Generative ABMs (PGABMs). PGABMs explicitly define the latent
and observed variables of the system according to the available data, and
describe the conditional probabilities connecting them. This representation
allows deriving the likelihood function of the latent variables of the system
from the observed data, especially by leveraging the rigorous theory of
probabilistic graphical models to extrapolate the conditional independencies
among the variables in play. As a result, this approach is more principled,
more accurate, and faster than simulation-based alternatives (Lenti et al.,
2024).

Nevertheless, the likelihood-based approach requires a non-trivial analytical
process to derive the likelihood function. In addition, handling categorical
variables is still an open challenge: for instance, a maximum-likelihood
approach in our setting would require exploring an exponentially large
combination of agents’ parameter values (up $2^{200}$). The main contribution
of this work is a technique that circumvents the need to write an explicit
derivation of the likelihood, by using variational inference (VI) instead. VI
assumes a tractable, parametric, functional form for the approximation of the
posterior distribution of the parameters to be estimated. Then, it directly
optimizes the parameters of the approximation to minimize the distance between
the approximated and real posterior. In so doing, our VI approach can address
models with intractable likelihood functions (e.g., with a large number of
categorical variables). Figure 1 depicts the proposed pipeline.

VI approximates the posterior distribution of the parameters and thus
represents a Bayesian approach. Several Bayesian, simulation-based
methodologies exist to calibrate ABMs: some employ Monte Carlo schemes
(Grazzini et al., 2017; Ternes et al., 2022), others use neural approximation
techniques (Dyer et al., 2022). However, the computational burden of the model
simulations often limits the scalability of the inference procedure. In
particular, Grazzini et al. compare different Bayesian approaches, including
Markov Chain Monte Carlo (MCMC) and Approximate Bayesian Computation (ABC),
and lament the major obstacle of simulation time. In stark contrast to the
existing literature, our approach is _not_ simulation-based, as it relies on
the recovery of the probabilistic generative model, and the parameter
estimates are the result of a direct optimization procedure.

To demonstrate the viability of our approach, we apply it to a bounded
confidence model where each agent is parameterized by a categorical role
(Weisbuch et al., 2002). When estimating these roles, the VI-based approach is
$3\times$ more accurate than an MCMC baseline, and almost $8\times$ more
accurate than a simulation-based ABC. In addition, VI is almost an order of
magnitude faster than the alternatives.

##  2 Preliminaries

###  2.1 Variational Inference

In the typical Bayesian framework, the objective is to determine the posterior
distribution $p(\theta\mid y)$ of the model parameters $\theta$ given the
observed data $y$. Bayes’ theorem provides the relationship $p(\theta\mid
y)=p(y\mid\theta)p(\theta)/p(y)$, where $p(\theta)$ is the prior distribution
of the model parameters which encodes our initial beliefs about $\theta$, and
$p(y\mid\theta)$ is the likelihood of observing the data according to the
underlying probabilistic model. However, as $p(y)=\int p(\theta,y)\,d\theta$
is generally intractable, and closed-form solutions for $p(\theta\mid y)$ are
limited, Bayesian inference often relies on approximations.

MCMC is one of the main techniques: it uses a sampling routine for $\theta$
based on the prior distribution and is guaranteed to converge to the correct
posterior, asymptotically. However, MCMC algorithms do not scale well to large
datasets and often struggle to approximate multimodal posteriors effectively
(Blei et al., 2017). Variational Inference (VI) is an alternative to MCMC
algorithms for approximating a target density via an optimization step (Blei
et al., 2017). VI approximates the target posterior $p(\theta\mid y)$ by using
a more tractable family of densities $q_{\lambda}(\theta)$, called the
_variational distribution_ , indexed by the variational parameter $\lambda$.

The optimal variational parameter $\lambda^{*}$ is the one that minimizes the
KL-divergence of $q_{\lambda}(\theta)$ from $p(\theta\mid y)$. In practice,
the search for this parameter is done by maximizing the evidence lower bound
(ELBO), denoted as $\mathcal{L}(\lambda)$, as the KL-divergence itself is
intractable

$\displaystyle\mathcal{L}(\lambda)=\mathbbm{E}_{q}[\log
p(y,\theta)-q_{\lambda}(\theta)],$ (1)

which is equivalent to

$\displaystyle\mathcal{L}(\lambda)=\log
p(y)-KL(q_{\lambda}(\theta)\;\|\;p(\theta\mid y)).$ (2)

So, maximizing the ELBO is equivalent to minimizing the KL-divergence of
$q_{\lambda}(\theta)$ from $p(\theta\mid y)$.

Stochastic Variational Inference (SVI) is a common routine for maximizing the
ELBO (Hoffman et al., 2013). SVI uses gradient ascent to optimize the ELBO,
and estimates the gradient on subsamples of the data instead of using the
entire dataset for improved speed.

Variational Inference with Normalizing Flows. The choice of the variational
distribution represents a critical step in the VI process, as it encodes the
functional form of the approximation of $p(\theta\mid y)$. The family of
normal distributions is a common choice, thanks to its simplicity and
prevalence in Bayesian statistics. However, the normal family imposes strong
constraints on the variational distribution, such as symmetry and unimodality.
A recent and more flexible solution is to adopt normalizing flows (NFs) as the
variational distribution (Rezende and Mohamed, 2015). NFs are transformations
of simple probability distributions (e.g., normal) into richer and more
complex distributions through a sequence of invertible and differentiable
mappings. The expressiveness of these transformations enables NFs to represent
arbitrarily complex distributions, which we can efficiently sample from and
compute the ELBO of (Papamakarios et al., 2021).

Categorical variables. As the ELBO is maximized via gradient ascent, the
computation is tractable only for continuous, differentiable variables.
Categorical variables pose a great challenge, both because of their discrete
nature and because there is no inherent ordering within their support. The
Gumbel-Softmax reparameterization offers a viable solution for these cases
(Maddison et al., 2016; Jang et al., 2016). Rather than sampling a variable
$z$ from a categorical distribution with probabilities
$(\pi_{1},\ldots,\pi_{k})$, we sample from a distribution that smoothly
approximates the categorical distribution. Mathematically, this is expressed
as

$\displaystyle
z=\operatorname*{soft\,max}\limits_{i}\left[g_{i}+\log\pi_{i}\right],$ (3)

where $g_{i}$ is a Gumbel noise. A temperature $\tau$ parameterizes the
softmax and controls how closely it approximates the argmax (closer for
$\tau\to 0$). If $k=2$, the Gumbel-Softmax reparameterizes a Bernoulli
distribution. This technique provides a way to relax a probabilistic
generative model with categorical random variables which allows for optimizing
the ELBO and, consequently, estimating the target parameters.

In the realm of ABMs, differentiable ABMs have been developed for keeping the
gradient of categorical variables in simulation-based inference (Andelfinger,
2021). In this case, the categorical variables are modeled as deterministic
approximations of continuous variables, equipped with the gradient of such
continuous variables. We prefer the Gumbel-Softmax for two reasons. First, it
finds broader applicability in the probabilistic machine learning community as
a versatile solution for handling categorical variables. Second, it is
supported by a more robust theoretical framework for incorporating categorical
variables into probabilistic generative models.

###  2.2 Opinion Dynamics Model

Bounded confidence model with backfire effect. We analyze the bounded
confidence model with backfire effect, a popular opinion dynamics model that
depicts the dynamics in a social network with agreement and disagreement
behaviors (Jager and Amblard, 2005). Each agent has an opinion in $[0,1]$,
where 0 and 1 represent the two polar opposites on a given topic (e.g., on the
left–right political spectrum). At each time step, an agent interacts with one
of its neighbors, thus emulating a social network. If the opinions of the two
interacting agents are closer than a threshold—the bounded confidence interval
($\varepsilon^{+}$)—they have a positive interaction ($s=1$) and their
opinions converge by a convergence rate ($\mu^{+}$). Conversely, if their
opinions are further than a backfire threshold ($\varepsilon^{-}$), then they
have a negative interaction ($s=-1$) and their opinions diverge by a
divergence rate ($\mu^{-}$). In the other cases, their opinions remain
unchanged ($s=0$).

Let $x^{t}_{u}$ be the opinion of agent $u$ at time $t$, the dynamics of the
opinion after the interaction $e_{t}=(u,v)$ are as follows

$\displaystyle\left\\{\begin{array}[]{lcl}x^{t+1}_{v}=x^{t}_{v}+\mu^{+}\Delta
x^{t}_{uv}&\text{if}&\lvert\Delta x^{t}_{uv}\rvert<\varepsilon^{+},\\\
x^{t+1}_{v}=x^{t}_{v}&\text{if}&\varepsilon^{+}\leq\lvert\Delta
x^{t}_{uv}\rvert<\varepsilon^{-},\\\ x^{t+1}_{v}=x^{t}_{v}-\mu^{-}\Delta
x^{t}_{uv}&\text{if}&\lvert\Delta x^{t}_{uv}\rvert\geq\varepsilon^{-}.\\\
\end{array}\right.$ (7)

After each update, the opinions are clamped to $[0,1]$,

$\displaystyle x^{t+1}_{u}=\max(0,\min(x^{t+1}_{u},1)).$ (8)

Different choices of $\varepsilon^{+}$ and $\varepsilon^{-}$ lead to different
final configurations of the opinions at the end of the simulation (time step
$T$), which can capture consensus, polarization, or fragmentation (Jager and
Amblard, 2005).

Bounded confidence model with backfire effect and roles. We extend the model
to represent a social context where the agents are divided into non-
overlapping roles, similar to the model by Weisbuch et al. (2002). Each agent
has a role ($r_{u}$) which can be either _leader_ ($r_{u}=L$) or _follower_
($r_{u}=F$). Leaders and followers adopt different behaviors: followers are
more prone to change their opinions compared to leaders.

An interaction $(u,v)$ in this model is asymmetric since only the opinion of
$v$ gets updated and the dynamics depend on the role of $v$. It follows a
bounded confidence model with backfire effect with parameters
$\varepsilon^{+}_{F}$ and $\varepsilon^{-}_{F}$ if $v$ is a follower, and with
parameters $\varepsilon^{+}_{L}$ and $\varepsilon^{-}_{L}$ if $v$ is a leader.
We assume that followers are more inclined to change their opinions, and thus
$\varepsilon^{+}_{F}\geq\varepsilon^{+}_{L}$,
$\varepsilon^{-}_{F}\leq\varepsilon^{-}_{L}$, $\mu^{+}_{F}\geq\mu^{+}_{L}$,
and $\mu^{-}_{F}\geq\mu^{-}_{L}$. For conciseness, we define
$\varepsilon=(\varepsilon^{+}_{F},\varepsilon^{+}_{L},\varepsilon^{-}_{F},\varepsilon^{-}_{L})$.
In summary, at each time step, the model evolves according to these rules:

1: Sample two agents $(u,v)$ uniformly at random

2: Let $r_{v}\in\\{F,L\\}$ be the role of $v$

3: if $\lvert\Delta x^{t}_{uv}\rvert<\varepsilon^{+}_{K}$ then

4: $s\coloneqq 1$

5: $x^{t+1}_{v}\coloneqq x^{t}_{v}+\mu^{+}_{r_{v}}\Delta x^{t}_{uv}$

6: else if $\lvert\Delta x^{t}_{uv}\rvert\geq\varepsilon^{-}_{r_{v}}$ then

7: $s\coloneqq-1$

8: $x^{t+1}_{v}\coloneqq x^{t}_{v}-\mu^{-}_{r_{v}}\Delta x^{t}_{uv}$

9: else

10: $s\coloneqq 0$

11: $x^{t+1}_{v}\coloneqq x^{t}_{v}$

12: end if

13: $x^{t+1}_{v}\coloneqq\max(0,\min(x^{t+1}_{v},1))$

##  3 Probabilistic Generative ABMs

Translating the ABM into its probabilistic generative counterpart is a crucial
step in the inference process, as it defines the generative model we wish to
approximate. We split the variable $s$ into two variables, $s^{+}=1$ and
$s^{-}=0$ if $s=1$, else $s^{+}=0$ and $s^{-}=1$ if $s=-1$, otherwise
$s^{+}=s^{-}=0$. Then, we define the conditional probabilities of observing
the interaction outcomes $s^{+}$ and $s^{-}$

$\displaystyle P(s^{+}_{j}=1\mid
x_{t},e_{j},\varepsilon,r_{v})=\sigma(\rho\cdot(\varepsilon^{+}_{r_{v}}-\lvert\Delta
x^{t}_{uv}\rvert))$ (9) $\displaystyle P(s^{-}_{j}=1\mid
x_{t},e_{j},\varepsilon,r_{v})=\sigma(-\rho\cdot(\varepsilon^{-}_{r_{v}}-\lvert\Delta
x^{t}_{uv}\rvert)),$ (10)

where $\sigma(\cdot)$ is the sigmoid function and $\rho$ is its steepness
(Monti et al., 2020). This way, we have a probabilistic model where we observe
$s^{+}=1$ with a high probability if $\lvert\Delta
x^{t}_{uv}\rvert<\varepsilon^{+}_{r_{v}}$, and $s^{-}=1$ if $\lvert\Delta
x^{t}_{uv}\rvert\geq\varepsilon^{-}_{r_{v}}$. As $\rho\rightarrow+\infty$ such
probabilities go to 1, and the model converges to the one in Equation 7.111In
our experiments we set $\rho$ = 32, for which the probabilities of
interactions outsides the bounds are already negligible ($\leq 10^{-5})$. Note
that the samples of $s^{+}$ and $s^{-}$ are independent, so at each
interaction, we can have both a positive and a negative interaction, resulting
in a simultaneous convergence and divergence of the opinions. This process
ensures that any sample from the parameter space can be associated with a non-
zero probability of having generated the data.

Then, we use the probabilistic graphical model (represented in Figure 2) to
capture the conditional independencies between the random variables and the
repeated patterns. To this end, we identify the observed variables, the latent
variables, and the parameters of the PGM, together with the relationships that
connect them. In our setting, we assume that the pairs of interacting agents
($e$), the outcomes of the interactions ($s^{+}$, $s^{-}$), the initial
opinions ($X_{0}$), and the convergence and divergence rates
($\mu=(\mu^{+}_{L},\mu^{+}_{F},\mu^{-}_{L},\mu^{-}_{F})$) are observed.
$X_{t}$ is a deterministic variable, and it is computable from the previous
observations. Notice that $\mu$ and $e$ only affect observed variables, so we
consider them as parameters of the PGM (we omit $\mu$ from Figure 2). Instead,
$\varepsilon$ and the agent roles $r$ are latent. In particular, we consider
$r$ as a micro-parameter, as it is a vector of size equal to the number of
agents $N$, while $\varepsilon$ has dimension 4. So we need to estimate $N+4$
latent variables.

$\mathbf{X}_{t}$$\mathbf{X}_{0}$$\mathbf{X}_{t+1}$$s$$e$$\varepsilon$$\mathbf{r}$

$T$

Figure 2: Probabilistic Graphical Model associated with the Bounded Confidence
model with backfire effect and latent roles. Circles represent stochastic
variables, diamonds deterministic variables, and letters without enclosures
are given PGM parameters. Shaded variables are observed and white ones are
latent. $\boldsymbol{X_{t}}$ is the opinion vector at time $\boldsymbol{t}$,
$\boldsymbol{e}$ is the vector of the interacting agents, $\boldsymbol{s}$
encodes the interactions outcomes, $\boldsymbol{r}$ are the roles of the
agents, and $\boldsymbol{\varepsilon}$ is the latent vector
($\boldsymbol{\varepsilon^{+}_{F}}$, $\boldsymbol{\varepsilon^{+}_{L}}$,
$\boldsymbol{\varepsilon^{-}_{F}}$, $\boldsymbol{\varepsilon^{-}_{L}}$) of ABM
parameters. In our framework, we estimate the latent variables,
$\boldsymbol{r}$ and $\boldsymbol{\varepsilon}$, given the observations on
$\boldsymbol{s}$ and $\boldsymbol{X}$.

A note on nomenclature. In the field of ABMs, and dynamical systems in
general, a _variable_ depends on time while a _parameter_ is independent of
time. Conversely, the nomenclature for PGMs is more concerned about
randomness: a variable can assume different values, according to either a
stochastic or deterministic process. An observed variable without parents in
the graphical model is called a _parameter_. It should be clear from the
context when we are referring to a latent ABM parameter that we are estimating
vs. an observed PGM parameter that is given. Similarly, the same objects are
treated as ABM _parameters_ and PGM _latent variables_ , and this distinction
should be clear from the context.

##  4 Inference

After having derived the PGM, we need to infer the target latent variables. We
compare 4 different methods to do so: two methods based on variational
inference that represent our proposal, an MCMC method that is a more classic
take on the concept of PGABM, and a Bayesian simulation-based method.

Variational Inference. The first two methods are based on VI and differ only
in the choice of the family for the variational distribution. We cast the
model into a probabilistic programming language (NumPyro in our case). Since
we infer a $4+N$-dimensional variable, the variational distribution
$q_{\lambda}$ has the same dimensions. The two distinct VI methods use
different $q_{\lambda}$, respectively a multivariate normal distribution and
normalizing flows (see Section 2.1).

Upon sampling from $q_{\lambda}$, we need to restrict the parameters in the
appropriate domains, since $q_{\lambda}$ has support in $\mathbbm{R}^{4+N}$.
We compute $\hat{\varepsilon}$ and $\hat{r}$ from $\theta\sim q_{\lambda}$ by
using a sigmoid transform.

  * $\bullet$

$\hat{\varepsilon}^{+}_{F}=\sigma(\theta_{1})/2$

  * $\bullet$

$\hat{\varepsilon}^{+}_{L}=\sigma(\theta_{2})/2$

  * $\bullet$

$\hat{\varepsilon}^{-}_{F}=\sigma(\theta_{3})/2+1/2$

  * $\bullet$

$\hat{\varepsilon}^{-}_{L}=\sigma(\theta_{4})/2+1/2$

This step ensures that $\hat{\varepsilon}^{+}_{F}$ and
$\hat{\varepsilon}^{+}_{L}$ are bounded within $[0,0.5]$ and
$\hat{\varepsilon}^{-}_{F}$ and $\hat{\varepsilon}^{-}_{L}$ within $[0.5,1]$.

The procedure for the role vector $r$ is slightly more complex. We first
define $\hat{\phi}_{u}=\sigma(\theta_{u+4})$, which represents the posterior
probability that $u=L$. Instead of sampling the role with a Bernoulli
extractions of $\hat{\phi}$, we use the Gumbel-Softmax relaxation, with
temperature 0.1. Consequently, we define $\tilde{r}_{u}\sim\text{Gumbel-
Softmax}(\phi_{u},1-\phi_{u})$. Finally, we can express the bounded confidence
interval for the node as the expectation over its role
$\tilde{\varepsilon}_{u}=\tilde{r}_{u}\cdot\hat{\varepsilon}_{L}+(1-\tilde{r}_{u})\cdot\hat{\varepsilon}_{F}$,
in a way that is similar to a linear relaxation of an integer program. Note
that this expression simply represents a superposition of the beliefs of the
model given the current evidence for the role of an agent, but does not
require changing the model to define an intermediate chase ‘in-between’ roles.
This formulation allows determining $\tilde{\varepsilon}_{u}$ based on the
sampled Gumbel-Softmax probabilities and builds an interpolation that is
easier to handle.

Then, we define the process used to generate the data. Since our goal is to
estimate the ELBO, we need to sample from $q_{\lambda}$ and from
$p(s,\theta)$, as outlined in Equation 1. From the PGABM, we have closed form
solutions for $P(s^{+}_{j}=1\mid\varepsilon,X,r)$ and
$P(s^{-}_{j}=1\mid\varepsilon,X,r)$, representing $p(s\mid\theta)$. Hence, we
can directly sample $p(s,\theta)$, via $p(\theta)$. Thus, this formulation
allows us to automatically maximize the ELBO by seeking the optimal
approximation $q_{\lambda^{*}}$ of the posterior distribution of $\varepsilon$
and $r$.

Markov Chain Monte Carlo. An alternative approach to estimate the latent
variables of the PGM is sampling-based inference via an MCMC algorithm. We use
No-U Turn Sampler (NUTS), an adaptive Hamiltonian Monte Carlo method that is
known for its high efficiency and flexibility (Hoffman et al., 2014). The PGM
is the same as in the variational case.

Approximate Bayesian Computation. While the previous methods are based on
PGMs, Approximate Bayesian Computation (ABC) is a likelihood-free, simulation-
based method (Csilléry et al., 2010). ABC samples a set of parameters from a
prior distribution, and for each sample it runs the entire simulation of the
opinion dynamics model. It then considers summary statistics of the generated
data trace to see how ‘close’ it is to the original data. In this case, we use
the number of positive interactions and the number of negative interactions
per time step as summary statistics. If the distance between the summary
statistics of the observed and simulated trajectories is lower than an
acceptance threshold, then the simulation is accepted; otherwise it is
rejected. Finally, it derives the posterior distribution of the target
parameters by counting the number of acceptances and rejections from the
different regions of the parameter space via Bayes’ rule.

##  5 Results

Figure 3: Comparison between actual values $\varepsilon$ (x-axis) and
estimates $\boldsymbol{\hat{\varepsilon}}$ (y-axis) for each macroscopic
parameter (rows) and method (columns). Each experiment samples
$\varepsilon^{+}_{F}$ and $\varepsilon^{+}_{L}$ in
$\boldsymbol{\\{0.05,0.15,0.25,0.35,0.45\\}}$, such that $\varepsilon^{+}_{F}$
$\boldsymbol{\geq}$ $\varepsilon^{+}_{L}$, and $\varepsilon^{-}_{F}$ and
$\varepsilon^{-}_{L}$ in $\boldsymbol{\\{0.55,0.65,0.75,0.85,0.95\\}}$, such
that $\varepsilon^{+}_{F}$ $\boldsymbol{\leq}$ $\varepsilon^{+}_{L}$. Points
on the diagonals represent exact estimates.

Experimental Setting. We run a broad set of experiments to compare the
estimates of $\varepsilon^{+}_{F}$, $\varepsilon^{+}_{L}$,
$\varepsilon^{-}_{F}$, $\varepsilon^{-}_{L}$, and $r$ when using the four
Bayesian methods mentioned in Section 4.

  * $\bullet$

SVI with normal distribution as variational distribution. The ELBO is
maximized via Adam, with learning rate $0.01$, running for $20\,000$ epochs.

  * $\bullet$

SVI with NFs, transformed via Block Neural Autoregressive flow (De Cao et al.,
2020), with 2 flows. The ELBO is maximized via Adam, with learning rate
$0.01$, running for $10\,000$ epochs.

  * $\bullet$

Markov Chain Monte Carlo (MCMC) sampling using No U-Turn Sampler (NUTS). We
use $5000$ simulations as burn-in and $5000$ simulations for the estimate.

  * $\bullet$

Approximate Bayesian Computation (ABC), a well-known method in the literature
of ABM calibration, similar to the work by Grazzini et al. (2017). The
distance between the simulated time series and the data summary is the $L^{2}$
distance. We set an acceptance threshold of $5\sqrt{T}$, and we launch $1000$
simulations.

We run a grid of $912$ experiments where we sample $\varepsilon$ with the
constraint that $\varepsilon^{+}_{F}\geq\varepsilon^{+}_{L}$, both in
$[0,0.5]$, and $\varepsilon^{-}_{F}\leq\varepsilon^{-}_{L}$, both in
$[0.5,1]$. We set $\mu^{+}_{F}$ = $\mu^{+}_{L}$ = $\mu^{-}_{F}$ =
$\mu^{-}_{L}$ = $0.02$ and use 10 interactions per time step. We vary the
number of agents $N$, the length of the simulation $T$, and the proportion of
leaders to have all the possible combinations with $N\in\\{50,100,200,400\\}$,
$T\in\\{128,512,2048,8192\\}$, and leaders proportion
$\in\\{0.01,0.02,0.04,0.1,0.2\\}$. For each combination of parameters, we
simulate the ABM, then estimate the target parameters with the four methods.
For practical reasons, we abort an experiment after $3$ hours.

Since the methods under scrutiny are all Bayesian, the parameter estimates are
distributions. To ease understanding of the results, we compare the average
values of $200$ samples from each posterior distribution.222The code to
reproduce these results is open-source and available at
https://anonymous.4open.science/r/learning_micromacro-1B9C

Parameter Estimates. Figure 3 compares the estimates of the _macroscopic_
parameters with their actual values for all the experiments. Summarizing the
results, we obtain an average root-mean-square error (RMSE) of $0.044$ for SVI
with normal, of $0.036$ for SVI with NFs, $0.051$ for MCMC and $0.125$ for
ABC. The improved performance of the SVI-based methods is evident even by
inspecting the figure visually.

The main hurdle in the inference is discriminating between leaders and
followers, both because of the high dimensionality of the parameter and its
categorical nature. Since most of the users are followers, the leader’s
parameters are the hardest to estimate: when the methods are not able to
identify the leaders they tend to set
$\hat{\varepsilon}^{+}_{L}=\hat{\varepsilon}^{+}_{F}$ and
$\hat{\varepsilon}^{-}_{L}=\hat{\varepsilon}^{-}_{F}$. Since
$\varepsilon^{+}_{F}$ $\geq$ $\varepsilon^{+}_{L}$, several observations in
Figure 3 (row 3) are above the diagonal. Analogously, since
$\varepsilon^{-}_{F}\leq\varepsilon^{-}_{L}$, most of the estimates in Figure
3 (row 4) are below the diagonal.

Figure 4: Specific errors of $\varepsilon^{+}_{F}$, $\varepsilon^{+}_{L}$,
$\varepsilon^{-}_{F}$, $\varepsilon^{-}_{L}$, and $\mathbf{r}$ as functions of
$\boldsymbol{N}$ (left), $\boldsymbol{T}$ (center), and proportion of leaders
(right). The error bars represent the standard errors.  Figure 5: Estimation
time as a function of $\boldsymbol{N}$ (left) and $\boldsymbol{T}$ (right).
The error bars represent standard errors. Figure 6: Root-mean-square error
(RMSE) of $\varepsilon$ against some features of the ABMs trajectories, (i)
number of positive interactions of followers, (ii) number of positive
interactions of leaders, (iii) absolute difference between
$\varepsilon^{+}_{F}$ and $\varepsilon^{+}_{L}$, (iv) absolute difference
between $\varepsilon^{-}_{F}$ and $\varepsilon^{-}_{L}$, (v) variance of the
opinion at the end of the simulation. The regression lines are in logarithmic
scale.

Estimate Robustness. Figure 4 shows the average error for the four macroscopic
parameters of the ABM, and the error rate for the microscopic ones (the
roles), i.e., the proportion of roles that are not correctly estimated, as
functions of the three hyperparameters that we vary in the experiments—the
number of agents $N$, the length of the data trace $T$, and the proportion of
leaders. The main insights from Figure 4 are:

  * $\bullet$

SVI with NFs, followed by SVI with normal, outperforms the other competitors
in each scenario.

  * $\bullet$

The three methods based on PGMs (the two SVI methods and MCMC) show better
performances compared to the simulation-based method, ABC. In particular, ABC
is not able to guess the roles of the agents, with error rates on the role
microparameters above $40\%$.

  * $\bullet$

As $T$ grows, we have more data, thus the estimates of SVI and MCMC improve.

  * $\bullet$

Coherently to Figure 3, the estimates of leaders’ parameters,
$\varepsilon^{+}_{L}$ and $\varepsilon^{-}_{L}$, are worse than followers’
ones, $\varepsilon^{+}_{F}$ and $\varepsilon^{-}_{F}$. As mentioned, when the
roles are not correctly estimated the methods favor estimating
$\varepsilon^{+}_{F}$ and $\varepsilon^{-}_{F}$ correctly, which affects the
majority of the agents’ behaviors.

  * $\bullet$

The error rate of the roles, and consequently all the other errors, grows as
$N$ increases. Indeed, as $N$ grows we have both an increased number of target
parameters and a reduced number of observations per agent.

  * $\bullet$

A higher proportion of leaders leads to more accurate estimations. Indeed,
more leaders imply a larger number of interactions involving
$\varepsilon^{+}_{L}$ and $\varepsilon^{-}_{L}$, which helps in estimating the
correct roles by differentiating the agents.

  * $\bullet$

The VI methods are the only ones able to estimate the role of microparameters
with high accuracy. The average error rates for SVI with NFs, SVI with normal,
MCMC, and ABC, are $0.06$, $0.09$, $0.22$, and $0.47$, respectively. The
proportion of experiments with 100% accuracy is, respectively, $0.54$, $0.38$,
$0.29$, and $0$.

Running time. The estimation times vary considerably among the tested methods
(Figure 5). First, both the VI methods scale well with the length of the ABMs
trajectories ($T$). The scalability of VI arises from its optimization-based
approach, which avoids the need to explicitly evaluate the entire dataset
during inference. This computational advantage becomes more evident as the
dataset size grows (Blei et al., 2017). Second, the estimation time of SVI
with NFs grows linearly with $N$. This reflects the increased number of
parameters underlying the neural networks implementing the NFs. Conversely,
the estimation time of MCMC increases linearly with $T$, reflecting the number
of samples required to replicate larger datasets. Overall, SVI with normal is
the most efficient method: it reduces the estimation time by a factor of
$56.0$ compared to SVI with NFs, $36.9$ compared to MCMC, and $36.5$ compared
to ABC. On average, the estimation time of SVI with normal in the largest
experiments ($N=400$, $T=8192$) is $61$ seconds. The experiments having SVI
with NFs with $N=400$, and the ones with ABC and MCMC with $T=8192$ ran over
the time limit.

Error analysis. Figure 6 shows the error rate of the roles as a function of
some summary statistics of the generated data trace. This analysis provides
new insights into how to estimate the parameters based on the available
observations.

The first two panels show that an increase in the number of positive outcomes
of the interactions of the followers and of the leaders leads to a lower
error. Higher values of $\lvert\varepsilon^{+}_{F}-\varepsilon^{+}_{L}\rvert$
and $\lvert\varepsilon^{-}_{L}-\varepsilon^{-}_{F}\rvert$ are associated with
lower errors. Such differences determine how distinguishable are the behaviors
of the leaders and the followers. In the extreme case where
$\varepsilon^{+}_{F}$ = $\varepsilon^{+}_{L}$ and $\varepsilon^{-}_{F}$ =
$\varepsilon^{-}_{L}$ the leaders and followers are indistinguishable. The
last plot of Figure 6 compares the variance of the opinion at the end of the
simulation and the error rate. A low variance of the opinions at the final
time $T$ corresponds to a convergence of the opinions: a scenario of
consensus. In such scenarios, the probability of observing negative
interactions is extremely low and the probability of having positive
interactions is extremely high, which hinders discriminating between different
values of $\varepsilon$. Moreover, all the agents are likely to have similar
behavior despite the different roles.

These observations help to identify which settings provide learnable
experiments and help the modeler determine the properties of the opinion
dynamics model that affect the outcomes of the estimates. For instance, if we
focus on the experiments with
$\lvert\varepsilon^{+}_{L}-\varepsilon^{+}_{F}\rvert>0$,
$\lvert\varepsilon^{-}_{L}-\varepsilon^{-}_{F}\rvert>0$, and $\sum s_{L}>10$
the error rate of the roles with SVI with NFs drops from $0.06$ 0 to $0.029$.
Moreover, if we restrict to experiments with $var(X_{T})>0.01$, this error
rate decreases to $0.018$.

##  6 Discussion

This work explores the use of AI technology in synergy with agent-based models
(ABMs). On the one hand, ABMs are behavioral models that put the human at the
center of the system, by focusing on the macroscopic phenomena arising from
individual, microscopic choices. This perspective provides a mechanistic
explanation of social phenomena. On the other hand, the AI and ML literature
offers a rich toolbox to incorporate data into models. Moreover, these
technologies can assist model design by providing new insights into the
features that allow learning a model. They enable a model refinement phase
that iteratively adjusts the ABM rules to better fit the specific data
context. Furthermore, model design can benefit from model selection, where
different models are compared in the analysis of the same dataset to compare
different hypotheses.

Until now, the applications of AI in ABMs have primarily centered on
developing intelligent agents (Zhang et al., 2021), or representing ABMs with
neural meta-models (van der Hoog, 2019). By properly estimating the parameters
of a popular opinion dynamics model with variational inference, we pave the
way to the application of the broad theory of probabilistic AI in the realm of
inference on ABMs. Parameter estimation in PGABMs has been tackled via
expectation maximization (Monti et al., 2020, 2022) or by direct maximization
of the likelihood (Lenti et al., 2024). While effective, these approaches may
be limited in their applicability as they require an analytical derivation of
the likelihood. In this work, we offer a proof-of-concept of a general and
flexible approach able to estimate high-dimensional, categorical, and
continuous parameters in ABMs.

The model we study focuses on polarizing human behavior in social media. This
area is of great importance, as polarization has been identified as a major
societal risk, together with
misinformation.333https://www.weforum.org/publications/global-risks-
report-2024 These phenomena are related, as misinformation spreads faster in
polarized environments (Nikolov et al., 2021). A wide range of opinion
dynamics models for these phenomena has been proposed (Del Vicario et al.,
2016; Proskurnikov et al., 2015; Törnberg, 2018; Azzimonti and Fernandes,
2023). However, they have found limited applicability in real-world scenarios
due to the lack of a methodological framework to estimate and validate them.
Our methodology offers a data-driven perspective on opinion dynamics, thus
facilitating the exploration of different scenarios, and the evaluation of
potential policy interventions.

More generally, we emphasize the versatility of PGABMs for applications to
ABMs in different fields, such as human mobility, economics, and epidemiology.
The connection between PGABMs and real-world data can foster interdisciplinary
collaborations between AI specialists and domain experts, thus advancing
knowledge of human behavior and its consequences.

Limitations and Future Works. Although our study focuses only on one specific
model, we underline that the latter is not only popular but also challenging
because of the high dimensionality and categorical nature of the latent
parameters. Indeed, it would be unfeasible to visit all the $2^{200}$
combinations of roles within a simulation-based or likelihood-based approach.
Additionally, the presented pipeline required only minimal customization in
the translation from the ABM to the PGABM, its implementation in NumPyro, and
the optimization of the ELBO. Applying it to other models should be easy given
its flexibility and generality.

While the model we study aims to capture wide social behavior, the scale of
this study is far smaller. Real-world social media studies often involve tens
of thousands of users. However, our experiments were constrained to $3$ hours
of running time to accommodate the presentation of hundreds of experiments. In
practice, computational efforts could be allocated differently on a single
dataset. Additionally, probabilistic programming languages support GPU
accelerators and can be optimized for greater efficiency. Moreover, if
computational time is an issue, SVI with normal represents a great option as
it reached good estimates in the most demanding experiments in only $61$
seconds and scales well.

Finally, our work focuses on methodological development, so a closer
inspection of the results was beyond the scope of this paper. We acknowledge
that a deeper analysis of the posterior distributions returned by the models
would be an interesting research direction. Leveraging the Bayesian nature of
VI, future analyses could explore multimodal posterior distributions to
address identifiability problems.

Conclusions. We presented a novel methodology for estimating both the
macroscopic and the microscopic parameters of a bounded confidence model with
backfire effect with leaders and followers. Our approach can estimate high-
dimensional, numerical and categorical, global and individual parameters of
ABMs with high accuracy. To achieve this result, we rely on stochastic
variational inference with normalizing flows as variational distribution,
which outperformed other Bayesian baselines. Using a normal variational
distribution achieved good results at a fraction of the time.

## References

  * Andelfinger (2021) Philipp Andelfinger.  Differentiable agent-based simulation for gradient-guided simulation-based optimization.  In _Proceedings of the 2021 ACM SIGSIM Conference on Principles of Advanced Discrete Simulation_ , pages 27–38, 2021.
  * Azzimonti and Fernandes (2023) Marina Azzimonti and Marcos Fernandes.  Social media networks, fake news, and polarization.  _European Journal of Political Economy_ , 76:102256, January 2023.  ISSN 0176-2680.  doi: 10.1016/j.ejpoleco.2022.102256.  URL https://www.sciencedirect.com/science/article/pii/S0176268022000623.
  * Baumann et al. (2021) Fabian Baumann, Philipp Lorenz-Spreen, Igor M Sokolov, and Michele Starnini.  Emergence of polarized ideological opinions in multidimensional topic spaces.  _Physical Review X_ , 11(1):011012, 2021.
  * Bernardo et al. (2021) Carmela Bernardo, Lingfei Wang, Francesco Vasca, Yiguang Hong, Guodong Shi, and Claudio Altafini.  Achieving consensus in multilateral international negotiations: The case study of the 2015 paris agreement on climate change.  _Science Advances_ , 7(51):eabg8068, 2021.
  * Blei et al. (2017) David M Blei, Alp Kucukelbir, and Jon D McAuliffe.  Variational inference: A review for statisticians.  _Journal of the American statistical Association_ , 112(518):859–877, 2017.
  * Csilléry et al. (2010) Katalin Csilléry, Michael GB Blum, Oscar E Gaggiotti, and Olivier François.  Approximate bayesian computation (abc) in practice.  _Trends in ecology & evolution_, 25(7):410–418, 2010.
  * De Cao et al. (2020) Nicola De Cao, Wilker Aziz, and Ivan Titov.  Block neural autoregressive flow.  In _Uncertainty in artificial intelligence_ , pages 1263–1273. PMLR, 2020.
  * Del Vicario et al. (2016) Michela Del Vicario, Alessandro Bessi, Fabiana Zollo, Fabio Petroni, Antonio Scala, Guido Caldarelli, H. Eugene Stanley, and Walter Quattrociocchi.  The spreading of misinformation online.  _Proceedings of the National Academy of Sciences_ , 113(3):554–559, January 2016.  ISSN 0027-8424, 1091-6490.  doi: 10.1073/pnas.1517441113.  URL https://pnas.org/doi/full/10.1073/pnas.1517441113.  tex.ids= del2016spreading tex.date-added: 2019-10-28 15:35:56 +0000 tex.date-modified: 2019-10-28 15:35:56 +0000 publisher: National Acad Sciences.
  * Dyer et al. (2022) Joel Dyer, Patrick Cannon, J Doyne Farmer, and Sebastian Schmon.  Black-box bayesian inference for economic agent-based models.  _arXiv preprint arXiv:2202.00625_ , 2022.
  * Fagiolo et al. (2019) Giorgio Fagiolo, Mattia Guerini, Francesco Lamperti, Alessio Moneta, and Andrea Roventini.  Validation of agent-based models in economics and finance.  In _Computer simulation validation_ , pages 763–787. Springer, 2019.
  * Flache et al. (2017) Andreas Flache, Michael Mäs, Thomas Feliciani, Edmund Chattoe-Brown, Guillaume Deffuant, Sylvie Huet, and Jan Lorenz.  Models of Social Influence: Towards the Next Frontiers.  _Journal of Artificial Societies and Social Simulation_ , 20(4):2, 2017.
  * Fortunato and Castellano (2007) Santo Fortunato and Claudio Castellano.  Scaling and universality in proportional elections.  _Physical Review Letters_ , 99(13):138701, 2007.
  * Gallo et al. (2022) Luca Gallo, Mattia Frasca, Vito Latora, and Giovanni Russo.  Lack of practical identifiability may hamper reliable predictions in covid-19 epidemic models.  _Science advances_ , 8(3):eabg5234, 2022.
  * Grazzini and Richiardi (2015) Jakob Grazzini and Matteo Richiardi.  Estimation of ergodic agent-based models by simulated minimum distance.  _Journal of Economic Dynamics and Control_ , 51:148–165, 2015.
  * Grazzini et al. (2017) Jakob Grazzini, Matteo G Richiardi, and Mike Tsionas.  Bayesian estimation of agent-based models.  _Journal of Economic Dynamics and Control_ , 77:26–47, 2017.
  * Hoffman et al. (2013) Matthew D Hoffman, David M Blei, Chong Wang, and John Paisley.  Stochastic variational inference.  _Journal of Machine Learning Research_ , 2013.
  * Hoffman et al. (2014) Matthew D Hoffman, Andrew Gelman, et al.  The no-u-turn sampler: adaptively setting path lengths in hamiltonian monte carlo.  _J. Mach. Learn. Res._ , 15(1):1593–1623, 2014.
  * Jager and Amblard (2005) Wander Jager and Frédéric Amblard.  Uniformity, bipolarization and pluriformity captured as generic stylized behavior with an agent-based simulation model of attitude change.  _Computational & Mathematical Organization Theory_, 10(4):295–303, 2005.
  * Jang et al. (2016) Eric Jang, Shixiang Gu, and Ben Poole.  Categorical reparameterization with gumbel-softmax.  In _International Conference on Learning Representations_ , 2016.
  * Lamperti et al. (2018) Francesco Lamperti, Andrea Roventini, and Amir Sani.  Agent-based model calibration using machine learning surrogates.  _Journal of Economic Dynamics and Control_ , 90:366–389, 2018.
  * Lenti et al. (2024) Jacopo Lenti, Corrado Monti, and Gianmarco De Francisci Morales.  Likelihood-based methods improve parameter estimation in opinion dynamics models.  _Proceedings of the 17th ACM International Conference on Web Search and Data Mining_ , 2024.
  * Lux and Zwinkels (2018) Thomas Lux and Remco CJ Zwinkels.  Empirical validation of agent-based models.  In _Handbook of computational economics_ , volume 4, pages 437–488. Elsevier, 2018.
  * Maddison et al. (2016) Chris J Maddison, Andriy Mnih, and Yee Whye Teh.  The concrete distribution: A continuous relaxation of discrete random variables.  In _International Conference on Learning Representations_ , 2016.
  * Monti et al. (2020) Corrado Monti, Gianmarco De Francisci Morales, and Francesco Bonchi.  Learning opinion dynamics from social traces.  In _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, pages 764–773, 2020.
  * Monti et al. (2022) Corrado Monti, Marco Pangallo, Gianmarco De Francisci Morales, and Francesco Bonchi.  On learning agent-based models from data.  _arXiv preprint arXiv:2205.05052_ , 2022.
  * Nikolov et al. (2021) Dimitar Nikolov, Alessandro Flammini, and Filippo Menczer.  Right and left, partisanship predicts (asymmetric) vulnerability to misinformation.  _Harvard Kennedy School Misinformation Review_ , 2021.  doi: 10.37016/mr-2020-55.  URL https://misinforeview.hks.harvard.edu/?p=5702.
  * Papamakarios et al. (2021) George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan.  Normalizing flows for probabilistic modeling and inference.  _The Journal of Machine Learning Research_ , 22(1):2617–2680, 2021.
  * Platt (2020) Donovan Platt.  A comparison of economic agent-based model calibration methods.  _Journal of Economic Dynamics and Control_ , 113:103859, 2020.
  * Proskurnikov et al. (2015) Anton V Proskurnikov, Alexey S Matveev, and Ming Cao.  Opinion dynamics in social networks with hostile camps: Consensus vs. polarization.  _IEEE Transactions on Automatic Control_ , 61(6):1524–1536, 2015.
  * Quera-Bofarull et al. (2023) Arnau Quera-Bofarull, Ayush Chopra, Anisoara Calinescu, Michael Wooldridge, and Joel Dyer.  Bayesian calibration of differentiable agent-based models.  _arXiv preprint arXiv:2305.15340_ , 2023.
  * Rezende and Mohamed (2015) Danilo Rezende and Shakir Mohamed.  Variational inference with normalizing flows.  In _International conference on machine learning_ , pages 1530–1538. PMLR, 2015.
  * Ternes et al. (2022) Patricia Ternes, Jonathan A Ward, Alison Heppenstall, Vijay Kumar, Le-Minh Kieu, and Nick Malleson.  Data assimilation and agent-based modelling: towards the incorporation of categorical agent parameters.  _Open Research Europe_ , 1:131, 2022.
  * Törnberg (2018) Petter Törnberg.  Echo chambers and viral misinformation: Modeling fake news as complex contagion.  _PLOS ONE_ , 13(9):e0203958, September 2018.  ISSN 1932-6203.  doi: 10.1371/journal.pone.0203958.  URL https://dx.plos.org/10.1371/journal.pone.0203958.
  * van der Hoog (2019) Sander van der Hoog.  Surrogate modelling in (and of) agent-based models: A prospectus.  _Computational Economics_ , 53(3):1245–1263, 2019.
  * Weisbuch et al. (2002) Gérard Weisbuch, Guillaume Deffuant, Frédéric Amblard, and Jean-Pierre Nadal.  Meet, discuss, and segregate!  _Complexity_ , 7(3):55–63, 2002.
  * Zhang et al. (2021) Wei Zhang, Andrea Valencia, and Ni-Bin Chang.  Synergistic integration between machine learning and agent-based modeling: A multidisciplinary review.  _IEEE Transactions on Neural Networks and Learning Systems_ , 2021.