Message passing on factor graphs is a powerful framework for probabilistic inference, which finds important applications in various scientific domains.
The most wide-spread message passing scheme is the sum-product algorithm (SPA) which gives exact results on trees but often fails on graphs with many small cycles. We search for an alternative message passing algorithm that works particularly well on such cyclic graphs. Therefore, we challenge the extrinsic principle of the SPA, which loses its objective on graphs with cycles. We further replace the local SPA message update rule at the factor nodes of the underlying graph with a generic mapping, which is optimized in a data-driven fashion. These modifications lead to a considerable improvement in performance while preserving the simplicity of the SPA. We evaluate our method for two classes of cyclic graphs: the $2 \times 2$ fully connected Ising grid and factor graphs for symbol detection on linear communication channels with inter-symbol interference.
To enable the method for large graphs as they occur in practical applications, we develop a novel loss function that is inspired by the Bethe approximation from statistical physics and allows for training in an unsupervised fashion.

§ INTRODUCTION

Message passing on graphical models is a powerful framework to efficiently solve inference and optimization problems. The most prominent message passing algorithm is the SPA, also known as BP [Pearl, 1988], which implements exact inference on tree-structured graphs [Kschischang et al., 2001]. Due to its simplicity, the SPA is often applied to cyclic graphs where it becomes an iterative and approximate algorithm. While this works surprisingly well for various applications, such as decoding of low-density parity-check codes [Gallager, 1963], a class of error-correcting codes, the SPA performs poorly on frustrated systems, i.e., on graphs with many cycles and strong coupling between the nodes.

The seminal work of Yedidia et al., 2000 revealed a connection between the SPA and free energy approximations of statistical physics, in particular, the fixed points of BP correspond to stationary points of the Bethe free energy.
Based on this insight, alternative message passing methods were proposed which directly minimize the Bethe free energy [Yuille, 2002, Welling and Teh, 2013]. These algorithms are guaranteed to converge to an extremum of the Bethe free energy but are computationally more demanding than plain BP.
Wainwright et al., 2003 proposed tree-reweighted BP as a message passing algorithm on the “convexified” Bethe free energy, which is guaranteed to have a global minimum. While this algorithm has stronger convergence guarantees compared to BP, it involves the selection and optimization of so-called edge appearance probabilities, a graph-specific problem that is often non-trivial for practical applications.
Yedidia et al., 2000 proposed “generalized BP” as an algorithm that passes messages between regions of nodes instead of single nodes. Larger regions will generally improve the quality of the approximation, however, they also increase the computational complexity.

Recently, model-based deep learning has shown great potential to empower various suboptimal algorithms, such as the SPA on cyclic graphs. Neural BP, proposed by Nachmani et al., 2016, unfolds the iterations of the SPA on its underlying graph and equips the resulting deep network  with tunable weights. The GAP algorithm of Schmid and Schmalen, 2022 varies the observation model by preprocessing, thereby shaping a graph with more favorable properties  with respect to BP performance.
Satorras and Welling, 2021 extend GNN to factor graphs and propose a hybrid model where BP runs conjointly to a GNN which is structurally identical to the original factor graph but has fully parametrized message updates.
All these works have in common that they are based on the SPA as a core concept which is vigorously improved using machine learning in order to compensate for its shortcomings on graphs with cycles.
In this work, we follow an alternative approach and directly search for alternative message passing algorithms that perform especially well on graphs with cycles, where the SPA tends to fail. To this end, we replace the well-known SPA message update rule with a compact NN, which is optimized to find a superior local message update rule.
Furthermore, we discuss the role of the extrinsic information principle which was originally introduced for tree-structured graphs.
Based on the close connection of BP to the Bethe approximation, we propose a novel end-to-end loss function that allows unsupervised and application-agnostic training of new message passing schemes.

§ BACKGROUND

We briefly introduce factor graphs and the SPA as a widespread framework for probabilistic inference on graphical models. We refer the reader to [Kschischang et al., 2001] for an excellent in-depth treatment of the topic.

§.§ Factor Graphs

Let $f(\mathcal{X})$ be a multivariate function of ${\mathcal{X}=\{x_1,\ldots,x_N \}}$ which factors into a product of local functions $f_j$:

\begin{equation}
f(\mathcal{X}) = \frac1Z \prod\limits_{j=1}^J f_j(\mathcal{X}_j),\quad \mathcal{X}_j \subseteq \mathcal{X}.
\label{eq:factorization}
\end{equation}

A factor graph visualizes the factorization in (<ref>) as a bipartite graph.
Every variable $x_n$ is represented by a unique vertex, a so-called variable node, which we draw as a circle in the graph. Factor nodes represent the local functions $f_j$ and are visualized by squares. The undirected edges of the graph connect a factor node $f_j(\mathcal{X}_j)$ with a variable node $x_n$ if and only if $f_j$ is a function of $x_n$, i.e., if $x_n \in \mathcal{X}_j$. From a graphical perspective, $\mathcal{X}_j$ thus corresponds to the set of adjacent variable nodes to the factor node $f_j$. Similarly, we define $\mathcal{N}(x_n)$ to be the set of adjacent factor nodes to the variable node $x_n$.

In this work, we restrict the variables ${x_n \in \{+1,-1\}}$ to be binary and the local factors $f_j$ to be either functions of a singular variable $x_n$ or functions of pairs $(x_n,x_m)$ such that the factorization becomes

\begin{equation}
f(\mathcal{X}) = \frac1Z \prod\limits_{n=1}^N \psi_n(x_n) \prod_{\mathclap{(n,m)\in \mathcal{E}}} \psi_{n,m}(x_n,x_m),
\label{eq:mrf}
\end{equation}

where $\mathcal{E}$ is the set of edges in the graph. Figure <ref> shows an exemplary factor graph.

fn = [draw, very thick, rectangle, inner sep=1em, rounded corners]
vn = [draw, very thick, circle, inner sep=.8em]
jump/.style args=(#1) to (#2) over (#3) by #4
insert path=
let 1=($(#1)-(#3)$), 1=veclen(1,1),
2=atan2(1,1), 3=abs(#4), 4=#4>0 ?180:-180  in
(#1) – ($(#1)!\n1-\n3!(#3)$)
arc (2:2+4:3) – (#2)

[auto, node distance=5.5em and 5.5em, thick]
[vn, right=of c1,label=center:$s_3$](c3);
[vn, below=of c1,label=center:$s_2$](c2);
[vn, right=of c2,label=center:$s_4$](c4);
[fn, color=KITblue, anchor=center, label=center:$\psi_{31}$](I31) at ($(c1)!0.5!(c3)$) ;
[fn, color=KITblue, anchor=center, label=center:$\psi_{21}$](I21) at ($(c1)!0.5!(c2)$) ;
[fn, color=KITblue, anchor=center, label=center:$\psi_{43}$](I43) at ($(c3)!0.5!(c4)$) ;
[fn, color=KITblue, anchor=center, label=center:$\psi_{42}$](I42) at ($(c2)!0.5!(c4)$) ;
[fn, color=KITblue, anchor=center, label=center:$\psi_{32}$, rotate=45](I32) at ($(c2)!0.3!(c3)$) ;
[fn, color=KITblue, anchor=center, label=center:$\psi_{41}$, rotate=45](I41) at ($(c1)!0.3!(c4)$) ;

[fn, color=KITred, above left=0.6em and 0.6em of c1, label=center:$\psi_{1}$](F1);
[fn, color=KITred, below left=0.6em and 0.6em of c2, label=center:$\psi_{2}$](F2);
[fn, color=KITred, above right=0.6em and 0.6em of c3, label=center:$\psi_{3}$](F3);
[fn, color=KITred, below right=0.6em and 0.6em of c4, label=center:$\psi_{4}$](F4);

[-] (c1) – (F1.center);
[-] (c2) – (F2.center);
[-] (c3) – (F3.center);
[-] (c4) – (F4.center);
[-,black] (c1) – (I31);
[-,black] (c3) – (I31);
[-,black] (c2) – (I42);
[-,black] (c4) – (I42);
[-,black] (c1) – (I21);
[-,black] (c2) – (I21);
[jump=(c3) to (I32) over ($(c2)!0.5!(c3)$) by 3pt];
[-,black] (c2) – (I32);
[-,black] (c3) – (I43);
[-,black] (c4) – (I43);
[-,black] (c1) – (I41);
[-,black] (c4) – (I41);

[fn, color=KITred, fill=white, above left=0.6em and 0.6em of c1, label=center:$\psi_{1}$](F1);
[fn, color=KITred, fill=white, below left=0.6em and 0.6em of c2, label=center:$\psi_{2}$](F2);
[fn, color=KITred, fill=white, above right=0.6em and 0.6em of c3, label=center:$\psi_{3}$](F3);
[fn, color=KITred, fill=white, below right=0.6em and 0.6em of c4, label=center:$\psi_{4}$](F4);

Factor graph representation of (<ref>) with factor nodes of degree 2 (blue) and degree 1 (red). This graph also models the $2 \times 2$ fully connected Ising graph of Sec. <ref>.

§.§ Sum-product Algorithm

The SPA is a message passing algorithm that operates in a factor graph and attempts to determine the marginals of the multivariate function $f(\mathcal{X})$.
Messages are propagated between the nodes of the factor graph along its edges and represent interim results of the marginalization.
Let $m_{f_j \rightarrow x_n}(x_n)$ denote a message sent from a factor node $f_j$ along an edge to a variable node $x_n$ and let $m_{x_n \rightarrow f_j}(x_n)$ denote a message on the same edge, but sent in the opposite direction.
If the factor graph visualizes a probabilistic model, i.e., if the variable nodes represent random variables, a message $m_{f_j \rightarrow x_n}(x_n)$ can be interpreted as a probabilistic statement from node $f_j$ about the random variable $x_n$ to be in one of its possible states [Yedidia et al., 2005].
The SPA defines the updates of the propagating messages at the nodes of the factor graph according to the simple rules [Kschischang et al., 2001]:

\begin{align}
m_{x \rightarrow f_j}(x) &= \prod_{f_i \in \mathcal{N}(x) \setminus f_j} m_{f_i \rightarrow x}(x) \label{eq:VN_update_generic} \\
m_{f_j \rightarrow x}(x) &= \sum_{\sim \{ x \}} \left( f_j(\mathcal{X}_j) \prod_{x' \in \mathcal{X}_j \setminus x} m_{x' \rightarrow f_j}(x')  \right). \label{eq:FN_update_generic}
\end{align}

The summary operator $\sum_{\sim \{x\}}$ denotes the marginalization over all variables in $\mathcal{X}_j$ except for $x$.
One key property of the SPA is the extrinsic information principle which states that the update of an outgoing message $m_{\text{A} \rightarrow \text{B}}$ at node $\text{A}$ destined to node $\text{B}$ does not depend on the incident message $m_{\text{B} \rightarrow \text{A}}$ which travels on the same edge but in opposite direction.
For the special case of degree-2 factor nodes $\psi_{n,m}(x_n,x_m)$, the SPA update rule (<ref>) thereby simplifies to

\begin{equation*}
m_{\psi_{n,m} \rightarrow x_n}(x_n) = \sum_{x_m}\psi_{n,m}(x_n,x_m) \cdot m_{x_m \rightarrow \psi_{n,m}}(x_m).
\end{equation*}

Messages at factors nodes $\psi_n(x_n)$ with degree 1 are not updated at all.

Initially, all messages are set to some unbiased state before they are iteratively updated according to a certain schedule.
For tree-structured graphs, the messages converge after they have once traveled forward and backward through the entire graph. The result of the SPA, i.e., the marginal functions $f(x_n)$, are finally obtained by a combination of all messages incident to the respective variable nodes:

\begin{equation*}
f(x_n) = \prod_{f_i \in \mathcal{N}(x_n)} m_{f_i \rightarrow x_n}(x_n).
\end{equation*}

Since the SPA makes no reference to the topology of the factor graph and the message updates are local, the SPA may also be applied to factor graphs with cycles [Yedidia et al., 2003].
On graphs with cycles, the SPA only yields an approximation of the exact marginals. While this approximation works surprisingly well in many cases, even including particular classes of graphs with many small cycles, there are also cases where the results are quite poor or where the SPA does not converge at all [Murphy et al., 1999].

Relation to the Bethe Approximation
In their seminal work, Yedidia et al., 2000 showed a revealing connection between the SPA and free energy approximations in statistical physics.
From a variational perspective, probabilistic inference can be seen as an optimization problem

\begin{equation}
q = \arg\min_{q \in \mathbb{M}} D_\text{KL}(q || p),   \label{eq:prob_inference}
\end{equation}

where we want to find the distribution $q$ from the set $\mathbb{M}$ of all globally valid probability distributions, known as the marginal polytope [Wainwright and Jordan, 2008]. Since the KL divergence $D_\text{KL}(q||p)$ is always non-negative and zero if and only if ${q=p}$, we reach the minimum exactly for ${q=p}$.
Obviously, optimizing over all possible probability distributions ${q\in \mathbb{M}}$ is generally intractable.
Based on some general assumption, free energy methods simplify the problem in (<ref>) to the minimization of a variational free energy term. We refer the reader to [Yedidia et al., 2005] for a detailed elaboration on this topic.

The Bethe approximation restricts the distribution $q(\mathcal{X})$ to be a product of univariate distributions $b_n(x_n)$ and joint distributions $b_{n,m}(x_n,x_m)$ between pairs ${(n,m) \in \mathcal{E}}$:

\begin{equation*}
q_\text{Bethe}\left( \mathcal{X} \right) := \prod_{n=1}^N b_n(x_n) \prod_{\mathclap{(n,m)\in \mathcal{E}}} b_{n,m}(x_n,x_m). \label{eq:bethe_approx}
\end{equation*}

This simplification leads to the Bethe free energy

\begin{align*}
F_\text{Bethe} = &\sum_{\mathclap{(n,m)\in \mathcal{E}}} \quad \sum_{{x_n,x_m}} b_{n,m}(x_n,x_m) \log \mleft( \frac{b_{n,m}(x_n,x_m)}{\phi_{n,m}(x_n,x_m)}\mright) \\
&-\sum_{n=1}^N (|\mathcal{X}_n| - 1) \sum_{x_n} b_n(x_n) \log \mleft( \frac{b_n(x_n)}{\psi_n(x_n)} \mright),
\end{align*}

with ${\phi_{n,m}(x_n,x_m)} := {\psi_n(x_n) \psi_{n,m}(x_n,x_m) \psi_m(x_m)}$.
Moreover, the Bethe approximation relaxes the search space in (<ref>) from the marginal polytope $\mathbb{M}$ to the local polytope

\begin{align*}
\mathbb{L} = \Big \{ \hspace{-2pt} \{&b_n(x_n), b_{n,m}(x_n,x_m): \forall n \in [1;N], (n,m) \in \mathcal{E} \} \hspace{-2pt} :  \\
& \sum_{x_n} b_{n,m}(x_n,x_m) = b_m(x_m), \, \sum_{x_n} b_n(x_n) = 1 \\
& \sum_{x_m} b_{n,m}(x_n,x_m) = b_n(x_n)\Big \}.
\end{align*}

This means that the distributions ${b_n(x_n)}$ and ${b_{n,m}(x_n,x_m)}$ only need to locally fulfill consistency in a pairwise sense.
In summary, the Bethe approximation converts (<ref>) into the optimization problem

\begin{equation}
q_\text{Bethe} = \argmin2_{{\left\{b_n, b_{n,m} \right\} \in \mathbb{L}}} \, F_\text{Bethe} \mleft( \left\{b_n, b_{n,m} \right\} \mright). \label{eq:Bethe_minimization}
\end{equation}

Yedidia et al., 2000 showed that the fixed points of BP applied to a factor graph correspond to the stationary points of the respective Bethe free energy. Seen in this light, BP is a suboptimal algorithm to minimize $F_\text{Bethe}$. The approximative nature in this sense is twofold:
First, there may exist multiple fixed points of the SPA for the same factor graph, i.e., the solution of the (converged) BP might correspond to an extremum of $F_\text{Bethe}$ other than the global minimum in $\mathbb{L}$  [Knoll et al., 2018].
Second, the beliefs only fulfill the pairwise consistency constraints at the fixed points of BP. This means that the solution only lies within the local polytope $\mathbb{L}$ after BP has converged. However, BP does not necessarily converge and failure of convergence is a major error mode [Yuille, 2002].

Various methods to directly solve (<ref>) or variants thereof were proposed (see [Yedidia et al., 2005] and references therein).
Yuille, 2002 proposed to decompose the Bethe free energy into concave and convex parts which enables the application of a CCCP. That algorithm consists of a double loop where the outer loop iteratively minimizes $F_\text{Bethe}$ and the inner loop ensures that the pairwise consistency constraints are fulfilled. Due to the CCCP, the algorithm provably converges to an extremum of the Bethe free energy.

§.§ Examples

For the remainder of this section, we introduce two important classes of factor graphs which are the basis for the numerical experiments in Sec. <ref>.

Example 1 - Ising Graphs
We consider factor graphs with ${N=M^2}$ variable nodes, arranged in a square 2D lattice, in which pairs of adjacent variable nodes $(x_n,x_m)$ are symmetrically coupled by the weights $J_{n,m}$ via factor nodes ${\psi_{n,m}(x_n,x_m)} = {\exp \mleft( J_{n,m} x_n x_m \mright)}$.
Additionally, each variable node $x_n$ has local evidence in the form of a degree-1 factor node ${\psi_n(x_n) = \exp \left( \theta_n x_n \right)}$.
The Ising model originates from statistical physics where the binary variables ${x_n \in \{+1,-1\}}$ represent the orientation of elementary magnets in a lattice [Peierls, 1936]. Each magnet is exposed to a local field $\theta_n$ and is influenced by its neighbors via an assigned pairwise coupling $J_{n,m}$.
Besides its fundamental significance in statistical physics, the Ising model is a universal mathematical model and finds applications in many other scientific domains such as image processing [Besag, 1986] and modeling of social networks [Banerjee and El Ghaoui, 2008, Wainwright and Jordan, 2008].

Following [Yedidia et al., 2005, Mooij and Kappen, 2007, Knoll et al., 2018], we study the fully connected $2 \times 2$ Ising model, i.e., ${M=2}$ and ${N=4}$, where every pair of variable nodes is connected.
A factor graph representation of this model is given in Fig. <ref>.
With more cycles than variable nodes and a girth of $3$, this graph can be parametrized to a highly frustrated system and is thus able to highlight the weaknesses of the SPA [Yedidia et al., 2005]. In particular, we consider the Ising spin glass, where the parameters $\theta_n$ and $J_{n,m}$ are iid random variables, sampled from a uniform distribution $\mathcal{U}[-S,+S]$ with ${S\in\mathbb{R}^+}$.
We are interested in the computation of the marginal functions

\begin{equation}
f(x_n) = \sum_{\sim \{x_n\}} f(x_1,x_2,x_3,x_4), \quad n=1,2,3,4, \label{eq:Ising_APP}
\end{equation}

which correspond to marginal probability distributions ${p(x_n)=f(x_n)}$ if the Ising graph represents a probabilistic model.
While the direct computation of (<ref>) is still feasible for our example with ${N=4}$, the number of summations grows exponentially with $N$, which calls for alternative methods with lower complexity. Applying the SPA on the factor graph in Fig. <ref> yields the single beliefs $b_n(x_n)$ as an approximation of $f(x_n)$ with a complexity that only grows quadratically with $N$.

Example 2 - Symbol Detection
We study the problem of symbol detection in a digital communication system [Proakis and Salehi, 2007]. A transmitter sends a sequence of $N$ independent and uniformly distributed symbols ${c_n \in \{+1,-1\}}$ over a linear channel with memory, impaired by AWGN. The receiver observes the sequence

\begin{equation}
\vek{y} =
\underbrace{
\begin{pmatrix}
h_0   &     &  &  \\
\vdots& h_0  & \scalebox{1.3}{$\vek{0}$} &  \\
h_L   &\vdots& \ddots & \\
& h_L  &        & h_0\\
&   \scalebox{1.3}{$\vek{0}$}  & \ddots  & \vdots \\
&      &  & h_L\\
\end{pmatrix}
}_{=: \vek{H}}
\underbrace{
\begin{pmatrix}
c_1 \\ c_2 \\ \vdots \\ \\ c_N
\end{pmatrix}
}_{=: \vek{c}}
\underbrace{
\begin{pmatrix}
w_1 \\ w_2 \\ \vdots \\ \\ w_{N+L}
\end{pmatrix}
}_{=: \vek{w}}
\end{equation}

where ${\vek{h} \in \mathbb{R}^{L+1}}$ describes the impulse response of the channel of length ${L+1}$ and ${w_k \sim \mathcal{CN}(0,\sigma^2)}$ are independent noise samples from a complex circular Gaussian distribution.
Applying Bayes' theorem, the posterior distribution $p(\vek{c}|\vek{y})$
can be expressed in terms of the likelihood:

\begin{equation*}
p(\vek{c}|\vek{y}) = \frac{1}{Z} p(\vek{y}|\vek{c}) = \frac{1}{Z} \exp \mleft ( -\frac{\left( \vek{y-Hc} \right) ^2}{\sigma^2} \mright ).
\end{equation*}

In the context of symbol detection, we want to infer the transmit symbols $c_n$ based on the channel observation $\vek{y}$, i.e., we are interested in the marginal distributions $p(c_n|\vek{y})$.
Based on an observation model by Ungerboeck, 1974

\begin{equation*}
p(\vek{y}|\vek{c}) \propto \exp\mleft( {2{\text{Re}}\left\{ \bm{c}^{\textrm H} \bm{H}^{\textrm H} \bm{y} \right\} -\bm{c}^{\textrm H}\bm{H}^{\textrm H}\bm{H}\bm{c}\over \sigma^2}  \mright),
\end{equation*}

we can factorize the likelihood

\begin{equation}
p(\vek{y}|\vek{c}) = \frac1Z \prod\limits_{n=1}^{N} \left [ {F_n(c_n) }
\prod\limits_{\substack{m=1 \\ m < n}}^N  {I_{n,m}(c_n,c_m) } \right ]
\label{eq:Ungerboeck_factorization}
\end{equation}

into the factors

\begin{align*}
F_n(c_n) &:=  \exp \mleft( \frac{1}{\sigma^2} {\text{Re}} \mleft\{ 2 x_n c_n^\star - G_{n,n}|c_n|^2 \mright\} \mright) \label{eq:f_fn} \\
I_{n,m}(c_n,c_m) &:= \exp \mleft(  -\frac{2}{\sigma^2} {\text{Re}} \mleft\{  G_{n,m} c_m c_n^\star  \mright\}  \mright),
\end{align*}

where ${\vek{x}:=\vek{H}^\text{H} \vek{y}}$ and ${\vek{G}:=\vek{H}^\text{H}\vek{H}}$ are the matched filtered versions of the observation and the channel matrix, respectively.
Modeling a factor graph based on (<ref>) and applying the SPA yields a low-complexity symbol detection algorithm, originally proposed by Colavolpe et al., 2011.

§ MESSAGE PASSING FOR CYCLIC GRAPHS

Despite its drawbacks on cyclic graphs, the amazing success of the SPA lies in its simplicity and generality: it is only defined by a local message update rule which can be applied to any generic factor graph based on a suitable message update schedule.
Driven by this elegant concept, we are interested in finding message passing algorithms that perform well on graphs with many cycles where the SPA fails.
More specifically, we ask the following questions:

* If the SPA fails to converge, does an alternative local message update rule exist that converges (possibly to an extremum of the Bethe free energy) and which provides better results than the SPA?

* If the SPA converges to an extremum of the Bethe free energy, is there a local message update rule which yields superior performance, either because the SPA converges to a fixed point which only corresponds to a local instead of global minimum of the Bethe free energy, or because the Bethe approximation itself is a bad approximation in this case?

§.§ On Message Update Rules

A message update rule defines a mapping from one or multiple incident messages to one outgoing message, which is applied locally at the variable or factor nodes of a factor graph. Besides the initialization of the messages and their update schedule, these mappings fully define a graph-based inference algorithm.
The SPA update rule at the variable nodes (<ref>) is simply the product of all extrinsic messages. We adopt this quite intuitive aggregation principle and focus on finding a message update rule for the factor nodes, i.e., an alternative to (<ref>).
For factor nodes of degree 2, such as in (<ref>), the update rule simplifies to a mapping from one single incident message to one outgoing message:

\begin{equation}
\text{FN}_\text{e}(\psi_{n,m}): m_{x_n \rightarrow \psi_{n,m}}(x_n) \mapsto m_{\psi_{n,m} \rightarrow x_m}(x_m). \label{eq:extrinsic_update}
\end{equation}

If the pairwise factors $\psi_{n,m}(x_n,x_m)$ are symmetric with regard to $x_n$ and $x_m$, and follow the exponential form

\begin{equation*}
\psi_{n,m}(x_n,x_m) = \exp \mleft( E_{n,m} x_n x_m \mright), \; x_n,x_m \in \{+1,-1\},
\end{equation*}

we can distill the dependency from the function $\psi_{n,m}$ to the scalar parameter ${E_{n,m}\in\mathbb{R}}$, which quantifies the repulsive (${E_{n,m}<0}$) or attractive (${E_{n,m}>0}$) coupling between the nodes $x_n$ and $x_m$. This directly coincides with the pairwise coupling weights ${J_{n,m} = E_{n,m}}$ of the Ising model in Example 1. The factor nodes $I_{n,m}$ of Example 2 can be reduced to the coupling parameters ${E_{n,m} = -2G_{n,m} / \sigma^2}$.

Challenging the Extrinsic Principle
Most of the existing message passing algorithms follow the extrinsic information principle. For instance in turbo decoding, it is known to be an important property of good message passing decoders [Richardson and Urbanke, 2001].
Ensuring that only extrinsic messages are received, it prevents backcoupling of intrinsic information in tree-structured graphs, which would otherwise lead to a self-enhancement of the messages, also known as “double counting”. Thereby, it guarantees that the SPA is exact on trees [Kuck et al., 2020].
We argue that this is in general not valid for cyclic graphs where backcoupling of messages is inevitable due to the very nature of the cycles.
Therefore, we propose a second message update rule which operates contradictory to the extrinsic principle: instead of ignoring the intrinsic message, the message update should rather actively leverage this additional information, e.g., to ensure that local consistency between neighboring nodes is fulfilled.

Without the extrinsic principle, we need to reconsider the messages from degree-1 factor nodes which are then also subject to iterative updates. To avoid an increase in complexity due to additional message updates at the degree-1 factor nodes, we apply a clustering approach similar to [Rapp et al., 2022]. We split up the single factors $\psi_n(x_n)$ into $|\mathcal{X}_n|$ parts ${\Psi_n(x_n) := \left(\psi_n(x_n)\right)^{\frac{1}{|\mathcal{X}_n|}}}$ and merge them into the adjacent pairwise factors $\psi_{n,m}(x_n,x_m)$, such that the new clustered factors are

\begin{equation*}
\Psi_{n,m}(x_n,x_m) := \Psi_n\mleft(x_n\mright) \psi_{n,m}\mleft(x_n,x_m\mright) \Psi_m\mleft(x_m\mright).
\end{equation*}

The overall factorization (<ref>) simplifies to

\begin{equation*}
f(x_1,\ldots,x_N) = \frac1Z \prod_{{(n,m)\in \mathcal{E}}} \Psi_{n,m}(x_n,x_m),
\end{equation*}

which leads to the non-extrinsic mapping

\begin{equation}
\text{FN}(\Psi_{n,m}):
\begin{pmatrix}
m_{x_n \rightarrow \Psi_{n,m}}(x_n) \\ m_{x_m \rightarrow \Psi_{n,m}}(x_m)
\end{pmatrix}
\mapsto m_{\Psi_{n,m} \rightarrow x_m}(x_m). \label{eq:nonextrinsic_update}
\end{equation}

If the single factors are in exponential form

\begin{equation*}
\Psi_n(x_n) = \exp \left( E_n x_n \right), \quad x_n \in \{ +1,-1\},
\end{equation*}

the clustered factors $\Psi_{n,m}(x_n,x_m)$ are fully characterized by the three scalars $E_n, E_{n,m}$ and $E_m$.

§.§.§ Neural Networks as Function Approximators

Finding suitable mappings (<ref>) or (<ref>) such that the overall message passing algorithm performs well is generally non-trivial.
We employ feed-forward NN, known to be efficient universal function approximators [Hornik et al., 1989], to reduce the search space of all possible mappings to a set of weights and biases $\mathcal{P}$, which fully parametrize the NN. At a factor node $f_j$, the network accepts $N_\text{in}$ inputs and produces the updated outgoing message $m_{f_j \rightarrow x_n}$.
For factor graphs with binary variables $x_n$, the messages $m_{f_j \rightarrow x_n}(x_n)$ can be expressed in scalar LLR

\begin{equation*}
L_{f_j \rightarrow x_n} := \log \mleft( \frac{m_{f_j \rightarrow x_n}(x_n=+1)}{m_{f_j \rightarrow x_n}(x_n=-1)} \mright).
\end{equation*}

A similar definition holds for the LLR $L_{x_n \rightarrow f_j}$ based on the messages $m_{x_n \rightarrow f_j}(x_n)$.
For the extrinsic update (<ref>), there are ${N_\text{in}=2}$ inputs: the LLR of the incoming extrinsic message and the coupling parameter $E_{n,m}$ of the local factor node. Without the extrinsic principle, the NN furthermore accepts the LLR of the intrinsic message as well as $E_n$ and $E_m$, i.e., in total $N_\text{in}=5$ inputs.
Since we only approximate a local mapping from a few scalar inputs to a single output, we can choose a very compact NN structure with a single hidden layer and $7$ neurons, as summarized in Table <ref>.

NN Architecture

Layer (linear)     Activation     Dimension

Input       ReLU       $(N_\text{in},7)$

Hidden      Tanh       $(7,7)$

Output      Linear     $(7,1)$

Having set up the NN structure, we are able to define a convenient message update rule by appropriately tuning the parameterization $\mathcal{P}$ of the NN.
We are interested in a local update rule such that the overall message passing performs well. To this end, we optimize $\mathcal{P}$ with respect to an objective function that evaluates the end-to-end performance of the inference task.
Therefore we apply a fixed number of message passing iterations and back-propagate the gradient of the objective function in order to iteratively optimize $\mathcal{P}$ using gradient descent based on a representative set of examples.
Note that this data-driven approach inevitably leads to a specialization of the learned message update to the data. However, we expect the result to be fairly generic and to have good generalization capabilities since we only optimize very few parameters in an otherwise model-aware system.
Moreover, despite the end-to-end optimization, we only use a single message update rule for the entire factor graph, i.e., we employ the same instance of the NN for the message updates at all factor nodes and in each iteration[As a consequence, the training procedure of the NN is not entirely local because the local copies of the NN at each factor node must be globally synchronized during optimization. However, the local nature of the message updates is still retained.].

We note that our approach can be interpreted as a special instance of a GNN as, e.g., described by Yoon et al., 2019.
In comparison, our model passes scalar messages instead of high-dimensional vectors and does not use any hidden states or embeddings at the variable nodes. For this reason, we do not require a second NN with a gated recurrent unit, as used in [Yoon et al., 2019] to update the hidden states based on the aggregated messages. Furthermore, we do not require a third NN which implements a trainable readout function to interpret the final node embeddings.

§.§ End-to-end Objective Functions

In the generic context of marginal inference, we hope to find a good approximation of the true marginals. A convenient objective function is the KL divergence which measures a type of statistical distance between the beliefs $b_n(x_n)$ and the exact marginal distributions ${p(x_n) = \sum_{\sim\{x_n\}} p(x_1,\ldots,x_N)}$:

\begin{equation}
\mathcal{L}_\text{KL} := D_\text{KL}\left(b_n(x_n) \middle\| p(x_n) \right). \label{eq:KL_to_APP}
\end{equation}

For large graphs, the computation of $p(x_n)$ might be infeasible, and $\mathcal{L}_\text{KL}$ becomes impractical. Therefore, we propose alternative loss functions in what follows.

The training of a symbol detector as in Example 2 is a typical supervised learning scenario where the labels are given by the transmitted symbols $c_n$. An appropriate performance measure for symbol detection is the BMI
which is an achievable information rate[In our case, where the symbols $c_n$ follow a Rademacher distribution, the BMI is equivalent to the mutual information.] for our scenario [Guillén i Fàbregas et al., 2008].
By a sample mean estimate over $D$ labeled examples ${(\vek{c},\vek{y})}$ from the data batch $\mathcal{D}$, the BMI can be approximated by

\begin{equation*}\label{bmi_est}
\text{BMI} \approx 1 -
\frac{1}{D N}
\sum\limits_{n=1}^{N}
\sum\limits_{(\vek{c}, \vek{y})\in \mathcal{D}}
\log_2 \mleft (  \text{e}^{-c_n L_{n}(\vek{y})} + 1  \mright ),
\end{equation*}

where ${L_n(\vek{y})}$ denotes the LLR from the belief ${b_n(c_n)}$ [Alvarado et al., 2018].

Other applications such as the Ising model in Example 1 relate to the class of unsupervised problems if the true marginals are not accessible.
For such scenarios, we consider a novel and application-agnostic objective function in the following. Inspired by the Bethe approximation, which is known to yield excellent results for many applications, even in cases where the SPA performs poorly [Yuille, 2002], we propose a regularized minimization of the Bethe free energy:

\begin{equation}
\mathcal{L}_\text{Bethe} := F_\text{Bethe} + \alpha \mathcal{L}_{\mathbb{L}}, \quad \alpha \in \mathbb{R}^+.
\label{eq:constrainedBethe}
\end{equation}

To ensure local consistency, we introduce the Bethe consistency distance

\begin{align*}
\mathcal{L}_{\mathbb{L}} := \; &D_\text{KL}\left( \sum_{x_m} b_{n,m}(x_n,x_m) \middle\| b_n(x_n) \right) \nonumber \\
&+ D_\text{KL}\left( \sum_{x_n} b_{n,m}(x_n,x_m) \middle\| b_m(x_m) \right)
\end{align*}

as a type of distance measure between the solution of the approximative inference ${\{ b_n,b_{n,m} \}}$ and the local polytope $\mathbb{L}$.
The weight $\alpha$ in (<ref>)  is a hyperparameter that controls how strictly the local consistency is enforced. With this penalty term $\mathcal{L}_\mathbb{L}$, we hope to suppress oscillations in the message passing, as they occur in the SPA for graphs with strong coupling.

§ EXPERIMENTS

We consider the examples of Sec. <ref> for numerical evaluation. To enable a deeper analysis, we fix the number of variable nodes to ${N=4}$ such that the computation of the true marginals is feasible. Despite this rather small extent, these models lead to factor graphs with a high density of short cycles and are thus expressive examples to highlight the weaknesses of the SPA.
Furthermore, we fix the global settings of the message passing to standard choices: all LLR messages are initialized with zero and we perform 10 iterations of a parallel schedule, i.e., each iteration comprises the parallel update of all messages at the factor nodes followed by message updates at all variable nodes.

A common technique to improve the performance of the SPA on graphs with cycles is the use of “momentum”, i.e., replacing a message $L^{(t)}$ of the SPA in iteration $t$ with the weighted average $(1-\mu) L^{(t)} + \mu L^{(t-1)}$ [Murphy et al., 1999]. By choosing ${0< \mu < 1}$, the idea is to improve the convergence behavior of the message passing scheme compared to the original SPA (${\mu=0}$) while retaining the same fixed points. As in [Murphy et al., 1999], we set ${\mu=0.1}$ and use this variant of the SPA as an additional baseline in the following experiments, where we refer to it as SPA$_\mu$.

Besides the SPA, we similarly apply message passing based on the newly proposed update rule (<ref>). We call the resulting inference algorithm cycBP (BP for cyclic graphs). If we use the extrinsic update rule (<ref>), we denote the algorithm with $\text{cycBP}_\text{e}$. We also consider the CCCP for the Bethe free energy as defined in [Yuille, 2002], since it gives interesting insights into the quality of the Bethe approximation. For the double loop, we apply $25$ outer iterations, each comprising $25$ inner iterations.

Ising model
We study the $2 \times 2$ fully connected spin glass model of Example 1 for ${S=2}$, i.e., all parameters $\theta_n$ and $J_{n,m}$ are independently sampled from a uniform distribution $\mathcal{U}[-2,+2]$. Table <ref> evaluates the behavior of all discussed inference schemes, averaged over $10^5$ different graphs. $\hat{\sigma}_{\mathcal{L}_\text{KL}}$ denotes the empirical standard deviation of ${\mathcal{L}}_\text{KL}$ of the individual graphs from the empirical mean.

Behavior of the Novel Message Passing Algorithm cycBP for the $2 \times 2$ Spin Glass, Averaged over $10^5$ Graphs

Algo.     Loss     ${\mathcal{L}}_\text{KL}$     $\hat{\sigma}_{\mathcal{L}_\text{KL}}$     $F_\text{Bethe}$     $\mathcal{L}_\mathbb{L}$

SPA     0.-     $0.087$     $0.265$     $-7.50$     $0.30$

SPA$_{\mu}$     0.-     $0.035$     $0.113$     $-7.49$     $0.12$

CCCP     0.-     $0.044$     $0.094$     $-7.24$     $2 \cdot 10^{-6}$

$\text{cycBP}_\text{e}$     $\mathcal{L}_\text{KL}$     $0.040$     $0.068$     $-7.37$     $0.17$

cycBP     $\mathcal{L}_\text{KL}$     $0.014$     $0.023$     $-7.37$     $0.48$

$\text{cycBP}_\text{e}$     $\mathcal{L}_\text{Bethe}$     $0.030$     $0.054$     $-7.38$     $0.11$

cycBP     $\mathcal{L}_\text{Bethe}$     $0.027$     $0.057$     $-7.47$     $0.027$

We can observe that the SPA does not leverage the full potential of the Bethe approximation, since the average loss ${\mathcal{L}_\text{KL}=0.087}$ of the SPA is twice as large compared to ${\mathcal{L}_\text{KL}=0.044}$ for the CCCP. Although the SPA reaches on average a smaller $F_\text{Bethe}$ than the CCCP, the beliefs of the SPA show local inconsistencies with ${\mathcal{L}_\mathbb{L}=0.3}$ due to non-convergent behavior. Using “momentum” in the SPA message updates can help to mitigate this behavior: the SPA$_\mu$ shows improved pairwise consistency ${\mathcal{L}_\mathbb{L}=0.12}$ and also yields in average a better approximation of the true marginals (${\mathcal{L}_\text{KL}=0.035}$).
The CCCP has a vanishing Bethe consistency distance $\mathcal{L}_\mathbb{L}$, i.e., the results of the CCCP lie within the local polytope $\mathbb{L}$.
We search for alternative message update rules, by optimizing $\mathcal{P}$ of the NN-based mappings towards minimal $\mathcal{L}_\text{KL}$. The training batches are sampled from a spin glass model with ${S=3}$ to put more emphasis on graphs with strong coupling, where the SPA is known to be susceptible to convergence errors. The results in Tab. <ref> show that there indeed exist superior message update rules to the SPA for this class of cyclic graphs. Using the extrinsic update rule (<ref>), the $\text{cycBP}_\text{e}$ algorithm reaches ${\mathcal{L}_\text{KL}=0.04}$ and thereby outperforms the original SPA as well as the CCCP.

We visualize the message update rule of the $\text{cycBP}_\text{e}$ algorithm in Fig. <ref> by plotting the optimized mapping (<ref>) from the incoming LLR message $L_{x_n \rightarrow \psi_{n,m}}$ to the outgoing LLR message $L_{\psi_{n,m} \rightarrow x_m}$. Similar to the SPA, the mapping is point-symmetric to the origin. The major difference is the behavior for incident LLR messages with high magnitudes ${|L_{x_n \rightarrow \psi_{n,m}}|>8}$, where the outgoing messages are heavily attenuated. Intuitively, this behavior reduces the potential of oscillation in graphs with strong coupling $E_{n,m}$.
We can further improve the inference performance by disabling the extrinsic principle in the message passing procedure. The resulting algorithm cycBP can be interpreted as a generalization of $\text{cycBP}_\text{e}$ and outperforms the latter with ${\mathcal{L}_\text{KL}=0.014}$, as reported in Tab. <ref>. It also yields a superior approximation of the true marginals compared to the momentum-based SPA$_\mu$, although the Bethe consistency distance ${\mathcal{L}_\mathbb{L}=0.48}$ is relatively high in this case.

[baseline, trim axis right]
grid style=gray!30,
xlabel= $L_{x_n \rightarrow \psi_{n,m}}$,
ylabel= $L_{\psi_{n,m} \rightarrow x_m}$,
y label style=yshift=-5pt,
enlarge x limits=false,
enlarge y limits=false,
xmin = -25,
xmax = +25,
ymin = -4.5,
ymax = +4.5,
point meta min=-2,
point meta max=2,
colormap=c_colorsamples of colormap=7 of viridis,
colorbar horizontal,
colorbar style=height=7pt, at=(0,1.15),anchor=north west, xtick=-2,-1.33,-0.67,0.0,0.67,1.33,2,xticklabel pos=upper,ylabel=$E_{n,m}$,
axis y line*=right,
ylabel style=rotate=-90,
width=.87 */pgfplots/parent axis width,

[index of colormap=0, line width=1.3pt, dotted] table[x=LLR_in, y=LLR_out_BP J0 ,col sep=comma] num_results/Ising/table_2x2Ising_LLR_in_out.csv;
[index of colormap=1, line width=1.3pt, dotted] table[x=LLR_in, y=LLR_out_BP J1 ,col sep=comma] num_results/Ising/table_2x2Ising_LLR_in_out.csv;
[index of colormap=2, line width=1.3pt, dotted] table[x=LLR_in, y=LLR_out_BP J2 ,col sep=comma] num_results/Ising/table_2x2Ising_LLR_in_out.csv;
[index of colormap=3, line width=1.3pt, dotted] table[x=LLR_in, y=LLR_out_BP J3 ,col sep=comma] num_results/Ising/table_2x2Ising_LLR_in_out.csv;
[index of colormap=4, line width=1.3pt, dotted] table[x=LLR_in, y=LLR_out_BP J4 ,col sep=comma] num_results/Ising/table_2x2Ising_LLR_in_out.csv;
[index of colormap=5, line width=1.3pt, dotted] table[x=LLR_in, y=LLR_out_BP J5 ,col sep=comma] num_results/Ising/table_2x2Ising_LLR_in_out.csv;
[index of colormap=6, line width=1.3pt, dotted] table[x=LLR_in, y=LLR_out_BP J6 ,col sep=comma] num_results/Ising/table_2x2Ising_LLR_in_out.csv;
[index of colormap=0, line width=1.3pt] table[x=LLR_in, y=LLR_out_NMP J0 ,col sep=comma] num_results/Ising/table_2x2Ising_LLR_in_out.csv;
[index of colormap=1, line width=1.3pt] table[x=LLR_in, y=LLR_out_NMP J1 ,col sep=comma] num_results/Ising/table_2x2Ising_LLR_in_out.csv;
[index of colormap=2, line width=1.3pt] table[x=LLR_in, y=LLR_out_NMP J2 ,col sep=comma] num_results/Ising/table_2x2Ising_LLR_in_out.csv;
[index of colormap=3, line width=1.3pt] table[x=LLR_in, y=LLR_out_NMP J3 ,col sep=comma] num_results/Ising/table_2x2Ising_LLR_in_out.csv;
[index of colormap=4, line width=1.3pt] table[x=LLR_in, y=LLR_out_NMP J4 ,col sep=comma] num_results/Ising/table_2x2Ising_LLR_in_out.csv;
[index of colormap=5, line width=1.3pt] table[x=LLR_in, y=LLR_out_NMP J5 ,col sep=comma] num_results/Ising/table_2x2Ising_LLR_in_out.csv;
[index of colormap=6, line width=1.3pt] table[x=LLR_in, y=LLR_out_NMP J6 ,col sep=comma] num_results/Ising/table_2x2Ising_LLR_in_out.csv;

The mapping (<ref>) of the $\text{cycBP}_\text{e}$ algorithm, optimized for the $2 \times 2$ spin glass with ${S=3}$ (solid), in comparison with the mapping of the SPA (dotted), plotted for different coupling $E_{n,m}$.

Moreover, we consider unsupervised training towards the proposed loss function $\mathcal{L}_\text{Bethe}$. For the $\text{cycBP}_\text{e}$ algorithm, the unsupervised training leads to a smaller loss ${\mathcal{L}_\text{KL}=0.03}$ compared to the supervised training, i.e., the loss function $\mathcal{L}_\text{Bethe}$ is better suited for the optimization via stochastic gradient descent than the loss function $\mathcal{L}_\text{KL}$ in this case.
In the unsupervised training towards $\mathcal{L}_\text{Bethe}$, we observe substantial differences between the two variants $\text{cycBP}_\text{e}$ and cycBP: while the optimization of $\text{cycBP}_\text{e}$ converges reliably, the training of cycBP is unstable and the optimization needs to run multiple times with different initializations for $\mathcal{P}$ until a reasonable result is obtained. This behavior is also reflected in the results in Tab. <ref>, where the cycBP algorithm shows a degraded performance with ${\mathcal{L}_\text{KL}=0.027}$, compared to the supervised training (${\mathcal{L}_\text{KL}=0.014}$). We conjecture that this is accounted for by the local consistency constraint $\mathcal{L}_\mathbb{L}$, which can be directly enforced at the message update at the factor nodes if the intrinsic message also takes part in the update. Optimization of the hyperparameter $\alpha$ did not lead to considerable changes in this behavior and we used ${\alpha=25}$ for all presented results.

[trim axis right]

ylabel shift = -10 pt,
ylabel style=rotate=-90,
xticklabel style=white,
xlabel shift = 10pt,
colorbar horizontal,
colorbar style=
width=1.93*/pgfplots/parent axis width,
height = 7 pt,
at=(0,1.15),anchor=north west,
ylabel style=rotate=-90,
xticklabel pos=upper,
axis y line*=right,
point meta min=0, point meta max=1.05,
only marks, mark=square*, mark options=scale=,solid,
table [col sep=comma, x=theta, y=J, z=KL] num_results/Ising/table_J_over_theta_BP_ref.csv;

at=(SPA.right of south east), anchor=left of south west,
ylabel style=rotate=-90,
xticklabel style=white,
xlabel shift = 10pt,
point meta min=0, point meta max=1.05,
only marks, mark=square*, mark options=scale=,solid]
table [col sep=comma, x=theta, y=J, z=KL] num_results/Ising/table_J_over_theta_CCCP2_ref.csv;

at=(SPA.below south west), anchor=above north west,
ylabel style=rotate=-90,
ylabel shift = -10 pt,
point meta min=0, point meta max=1.05,
only marks, mark=square*, mark options=scale=,solid]
table [col sep=comma, x=theta, y=J, z=KL] num_results/Ising/table_J_over_theta_trainedKL_NONextrinsic_feedDGR1msgs.csv;

at=(NMP_KL.right of south east), anchor=left of south west,
ylabel style=rotate=-90,
point meta min=0, point meta max=1.05,
only marks, mark=square*, mark options=scale=,solid]
table [col sep=comma, x=theta, y=J, z=KL] num_results/Ising/table_J_over_theta_trainedBethe_extrinsic.csv;

Approximation error $\mathcal{L}_\text{KL}$ on the $2 \times 2$ Ising model with constant parameters $\theta$ and $J$: SPA (top left), CCCP (top right), cycBP optimized on spin glasses w.r.t. $\mathcal{L}_\text{KL}$ (bottom left), and $\text{cycBP}_\text{e}$ trained with $\mathcal{L}_\text{Bethe}$ (bottom right).

To analyze the convergence properties on highly frustrated systems, we consider the Ising model with constant parameters $\theta$ and $J$. Note that we do not specifically optimize the models for this scenario, but rather use the previous parametrization $\mathcal{P}$ which is optimized for spin glasses with ${S=3}$. Figure <ref> plots $\mathcal{L}_\text{KL}$ over $\theta$ and $J$ for different inference algorithms. Knoll et al., 2018 showed that the Bethe free energy has a unique minimum in the complete antiferromagnetic domain (${J<0}$) and in large parts of the ferromagnetic case (${J>0}$), except for a region around ${\theta=0}$, where $F_\text{Bethe}$ has two  minima in $\mathbb{L}$. This coincides with our findings of the approximation error $\mathcal{L}_\text{KL}$ for the CCCP in Fig. <ref>. The SPA shows failure to converge in large parts of the antiferromagnetic region with ${J<-1}$, where it does not converge to the unique fixed point and produces large approximation errors. The extrinsic message passing scheme $\text{cycBP}_\text{e}$, optimized towards $\mathcal{L}_\text{Bethe}$, shows an improved behavior. However, in the antiferromagnetic case with strong repellings ($J<-1.5$), there are still considerable approximation errors. The non-extrinsic message passing algorithm cycBP shows good inference capabilities over the complete considered region if it is optimized towards $\mathcal{L}_\text{KL}$. The unsupervised training with respect to $\mathcal{L}_\text{Bethe}$ leads to a similar performance as the CCCP, however, the training procedure is again relatively unstable in this case.

Symbol Detection

[baseline, trim axis right]
grid style=gray!30,
xlabel= $E_\text{b}/N_0$ (dB),
ylabel= ${1-\text{BMI}}$,
ymode = log,
enlarge x limits=false,
enlarge y limits=false,
xmin = 2,
xmax = 14,
ymin = 0.00001,
ymax = 0.2,
legend style=at=(.03,0.03),anchor=south west, font=,
legend cell align=left,
[densely dotted, color=gray, line width=1.3pt] table[x=Eb/N0, y expr=1-4,col sep=comma] num_results/symbol_detection/BER+BMI_over_ebno_N4_BPSK_ref+nonextrinsic.csv;
[color=KITred, line width=1.3pt] table[x=Eb/N0, y expr=1-6,col sep=comma] num_results/symbol_detection/BER+BMI_over_ebno_N4_BPSK_ref+nonextrinsic.csv;
[dashed, color=KITred, line width=1.3pt] table[x=Eb/N0, y expr=1-2,col sep=comma] num_results/symbol_detection/BER+BMI_over_ebno_N4_BPSK_momentumBP.csv;
[color=KITorange, line width=1.3pt] table[x=Eb/N0, y expr=1-8 ,col sep=comma] num_results/symbol_detection/BER+BMI_over_ebno_N4_BPSK_extrinsic.csv;
[color=KITpalegreen, line width=1.3pt] table[x=Eb/N0, y expr=1-2 ,col sep=comma] num_results/symbol_detection/BER+BMI_over_ebno_N4_BPSK_extrinsic.csv;
[dotted, color=KITpalegreen, line width=1.3pt] table[x=Eb/N0, y expr=1-2 ,col sep=comma] num_results/symbol_detection/BER+BMI_over_ebno_N4_BPSK_extrinsic_constrainedBethe.csv;
[color=KITblue, line width=1.3pt] table[x=Eb/N0, y expr=1-2 ,col sep=comma] num_results/symbol_detection/BER+BMI_over_ebno_N4_BPSK_ref+nonextrinsic.csv;

[dashed, color=KITred, line width=1.3pt] table[x=Eb/N0, y expr=1-2,col sep=comma] num_results/symbol_detection/BER+BMI_over_ebno_N4_BPSK_momentumBP.csv;

True marginals
$\text{cycBP}_\text{e}$ ($\text{BMI}$)
$\text{cycBP}_\text{e}$ ($\mathcal{L}_\text{Bethe}$)
cycBP ($\text{BMI}$)

Detection performance of the proposed cycBP algorithm, averaged over $10^7$ random channels.

We further consider the factor graphs of Example 2 for approximate symbol detection on linear channels with memory ${L=2}$.
To generate random channels, we independently sample each tap $h_\ell$ of the channel impulse response for every example from a Gaussian distribution with zero mean and unit variance and subsequently normalize each channel to unit energy ${\lVert \bm{h} \rVert_2 = 1}$.
Figure <ref> evaluates the detection performance of the considered inference algorithms in terms of the BMI over the signal-to-noise ratio ${E_\text{b}/N_0 = 1/\sigma^2}$. Both, the SPA and the CCCP run into an error floor for high $E_\text{b}/N_0$, where the graphs tend to have strong coupling via the factor nodes $I_{n,m}(c_n,c_m)$.
The momentum-based message updates of the SPA$_\mu$ enhance the original SPA in the entire ${E_\text{b}/N_0}$ range under consideration and close the performance gap to the CCCP.
During the optimization of the new message update rules, the $E_\text{b}/N_0$ in dB was sampled from ${\mathcal{U}[0,16]}$ for each batch element independently.
To help the update rule adapt to different channel realizations, i.e., different $E_\text{b}/N_0$ and different channel taps $\bm{h}$, we feed $E_\text{b}/N_0$ and $\bm{h}$ as additional inputs to the NN.
Figure <ref> shows that the optimized algorithms $\text{cycBP}_\text{e}$ and cycBP clearly outperform the SPA and the CCCP, especially for high $E_\text{b}/N_0$. Consistently with our findings on the Ising graphs, the cycBP algorithm performs better than the extrinsic variant $\text{cycBP}_\text{e}$. The latter can also be trained towards $\mathcal{L}_\text{Bethe}$ without degrading the detection performance. This is particularly surprising since it thereby clearly outperforms the CCCP. The optimization of the cycBP algorithm towards $\mathcal{L}_\text{Bethe}$ does not converge and is therefore not shown in Fig. <ref>. However, since training towards the BMI is feasible for large $N$, the supervised training of the cycBP algorithm yields an attractive algorithm with low complexity and superior performance, which can be highly relevant for practical applications.

§.§ Discussion

To investigate the two central questions which we formulated in Sec. <ref> and to show the potential of our method, we investigated compact models with ${N=4}$. These models are expressive examples for analysis since they have a high density of short cycles and because the true marginals are available as ground truth data.
However, verifying the capability of the proposed cycBP algorithm for practical applications requires extensive numerical evaluation on larger graphs and varying graph structures. This is ongoing work and our preliminary results are promising.

§ CONCLUSION

This work considered message passing for approximate inference and showed the existence of message update rules which perform especially well on cyclic graphs where the SPA fails. We challenged the extrinsic information principle for cyclic graphs and proposed an alternative message update rule which also takes intrinsic information into account. The gain was demonstrated by numerical experiments on two exemplary classes of factor graphs. The learned message update rules generalize well and training is extremely fast since the update rule is defined by a very compact NN that is reused at all factor nodes. We furthermore proposed a novel unsupervised and application-agnostic loss function that follows the idea of the Bethe approximation.

This work has received funding in part from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 101001899) and in part from the German Federal Ministry of Education and Research (BMBF) within the project Open6GHub (grant agreement 16KISK010).

[Alvarado et al., 2018]
Alex Alvarado, Tobias Fehenberger, Bin Chen, and Frans M. J. Willems.
Achievable information rates for fiber optics: Applications and
J. Lightw. Technol., 360 (2):0 424–439, 2018.

[Banerjee and El Ghaoui, 2008]
Onureena Banerjee and Laurent El Ghaoui.
Model selection through sparse maximum likelihood estimation for
multivariate gaussian or binary data.
J. Mach. Learn. Re., 9:0 485–516, 2008.

[Besag, 1986]
Julian Besag.
On the statistical analysis of dirty pictures.
J. R. Stat. Soc. Ser. B. Stat. Soc., 480 (3):0
259–302, 1986.

[Colavolpe et al., 2011]
Giulio Colavolpe, Dario Fertonani, and Amina Piemontese.
SISO detection over linear channels with linear complexity in the
number of interferers.
IEEE J. Sel. Topics Signal Process., 50 (8), 2011.

[Gallager, 1963]
Robert G. Gallager.
Low-density parity check codes.
PhD thesis, MIT Press, 1963.

[Guillén i Fàbregas et al., 2008]
Albert Guillén i Fàbregas, Alfonso Martinez, and Giuseppe Caire.
Bit-interleaved coded modulation.
In Found. Trends Commun. Inf. Theory, volume 5. Now
Publishers, Delft, NL, 2008.

[Hornik et al., 1989]
Kurt Hornik, Maxwell Stinchcombe, and Halbert White.
Multilayer feedforward networks are universal approximators.
Neural Networks, 20 (5):0 359–366, January

[Knoll et al., 2018]
Christian Knoll, Dhagash Mehta, Tianran Chen, and Franz Pernkopf.
Fixed points of belief propagation – an analysis via polynomial
homotopy continuation.
IEEE Trans. Pattern Anal. Mach. Intell., 400
(9):0 2124–2136, 2018.

[Kschischang et al., 2001]
Frank R. Kschischang, Brendan J. Frey, and Hans-Andrea Loeliger.
Factor graphs and the sum-product algorithm.
IEEE Trans. Inf. Theory, 470 (2):0 22, 2001.

[Kuck et al., 2020]
Jonathan Kuck, Shuvam Chakraborty, Hao Tang, Rachel Luo, Jiaming Song, Ashish
Sabharwal, and Stefano Ermon.
Belief propagation neural networks.
Adv. Neural Inf. Process. Syst., 33:0 667–678, 2020.

[Mooij and Kappen, 2007]
Joris M. Mooij and Hilbert J. Kappen.
Sufficient conditions for convergence of the sum-product algorithm.
IEEE Trans. Inf. Theory, 530 (12):0
4422–4437, 2007.

[Murphy et al., 1999]
Kevin Murphy, Yair Weiss, and Michael I. Jordan.
Loopy belief propagation for approximate inference: An empirical
In Proc. Innov. Appl. Artif. Intell. Conf., 1999.

[Nachmani et al., 2016]
Eliya Nachmani, Yair Be'ery, and David Burshtein.
Learning to decode linear codes using deep learning.
In Proc. Annu. Allerton Conf. Commun., Control, Comput.,
Monticello, IL, 2016.

[Pearl, 1988]
Judea Pearl.
Probabilistic Reasoning in Intelligent Systems:
Networks of Plausible Inference.
Morgan Kaufmann, 1988.

[Peierls, 1936]
Rudolf Peierls.
On Ising's model of ferromagnetism.
Math. Proc. Camb. Philos. Soc., 32:0 477–481, 1936.

[Proakis and Salehi, 2007]
John Proakis and Massoud Salehi.
Digital Communications.
McGraw Hill, 5 edition, 2007.

[Rapp et al., 2022]
Lukas Rapp, Luca Schmid, Andrej Rode, and Laurent Schmalen.
Structural optimization of factor graphs for symbol detection via
continuous clustering and machine learning.
arXiv:2211.11406, 2022.

[Richardson and Urbanke, 2001]
Thomas J. Richardson and Rüdiger L. Urbanke.
The capacity of low-density parity-check codes under message-passing
IEEE Trans. Inf. Theory, 470 (2):0 599–618,

[Satorras and Welling, 2021]
Victor Garcia Satorras and Max Welling.
Neural enhanced belief propagation on factor graphs.
In Proc. Int. Conf. on Artificial Intelligence and
Statistics (AISTATS), pages 685–693, San Diego, CA, USA, 2021.

[Schmid and Schmalen, 2022]
Luca Schmid and Laurent Schmalen.
Low-complexity near-optimum symbol detection based on neural
enhancement of factor graphs.
IEEE Trans. Commun., 700 (11):0 7562–7575,
November 2022.

[Ungerboeck, 1974]
Gottfried Ungerboeck.
Adaptive maximum-likelihood receiver for carrier-modulated
data-transmission systems.
IEEE Trans. Commun., COM-220 (5):0 624–636,

[Wainwright and Jordan, 2008]
Martin J. Wainwright and Michael I. Jordan.
Graphical models, exponential families, and variational
Found. Trends Mach., 10 (1–2):0 1–305,
Number: 1–2.

[Wainwright et al., 2003]
Martin J. Wainwright, Tommi S. Jaakkola, and Alan S. Willsky.
Tree-reweighted belief propagation algorithms and approximate ML
estimation by pseudo-moment matching.
Proc. Int. Workshop Artificial Intelligence and Statistics,
R40 (PMLR R4):0 308–315, 2003.

[Welling and Teh, 2013]
Max Welling and Yee Whye Teh.
Belief optimization for binary networks: A stable alternative to
loopy belief propagation.
In Proc. Conf. on Uncertainty in Artificial
Intelligence (UAI), pages 554–561, 2013.

[Yedidia et al., 2000]
Jonathan S. Yedidia, William T. Freeman, and Yair Weiss.
Generalized belief propagation.
Advances in Neural Information Processing Systems, 13, 2000.

[Yedidia et al., 2003]
Jonathan S. Yedidia, William T. Freeman, and Yair Weiss.
Understanding belief propagation and its generalizations.
In Exploring Artificial Intelligence in the New Millennium,
volume 8, pages 236–239. Morgan Kaufmann Publishers Inc., 2003.

[Yedidia et al., 2005]
Jonathan S. Yedidia, William T. Freeman, and Yair Weiss.
Constructing free-energy approximations and generalized belief
propagation algorithms.
IEEE Trans. Inform. Theory, 510 (7):0
2282–2312, 2005.

[Yoon et al., 2019]
KiJung Yoon, Renjie Liao, Yuwen Xiong, Lisa Zhang, Ethan Fetaya, Raquel
Urtasun, Richard Zemel, and Xaq Pitkow.
Inference in probabilistic graphical models by graph neural networks.
In Proc. Asilomar Conf. Signals Syst. Comput., 2019.

[Yuille, 2002]
Alan L. Yuille.
CCCP algorithms to minimize the Bethe and Kikuchi free
energies: Convergent alternatives to belief propagation.
Neural Computation, 140 (7):0 1691–1722,


# UAV-CROWD: Violent and non-violent crowd activity simulator from the
perspective of UAV

Mahieyin Rahmun1, Tonmoay Deb1, Shahriar Ali Bijoy1, Mayamin Hamid Raha1
1Department of Electrical and Computer Engineering, North South University,
Dhaka, Bangladesh
{mahieyin.rahman, tonmoay.deb, shahriar.ali<EMAIL_ADDRESS>
###### Abstract

Unmanned Aerial Vehicle (UAV) has gained significant traction in the recent
years, particularly the context of surveillance. However, video datasets that
capture violent and non-violent human activity from aerial point-of-view is
scarce. To address this issue, we propose a novel, baseline simulator which is
capable of generating sequences of photo-realistic synthetic images of crowds
engaging in various activities that can be categorized as violent or non-
violent. The crowd groups are annotated with bounding boxes that are
automatically computed using semantic segmentation. Our simulator is capable
of generating large, randomized urban environments and is able to maintain an
average of 25 frames per second on a mid-range computer with 150 concurrent
crowd agents interacting with each other. We also show that when synthetic
data from the proposed simulator is augmented with real world data, binary
video classification accuracy is improved by 5% on average across two
different models.

###### Index Terms:

Crowd, Crowd violence, UAV, Simulation, Unreal Engine, Airsim

(a) Kicking

(b) Using weapons

(v) Punching

(d) Chasing

(e) Having conversation

(f) Walking

(g) Dispersion

(h) Talking over phone

Figure 1: Simulated activities of crowd agents. (a) - (d) and (e)-(h)
represent violent and non-violent activities, respectively

##  I Introduction

Crime rate has been rising in recent the years [1]. With technological
advancements, both law enforcement forces and criminals have become more
effective. As such, automated surveillance techniques for constantly
monitoring key locations is necessary. An important step in this automation is
to recognize the type of action from video footage, which is an active topic
of research [2], [3]. Human Activity Recognition (HAR) requires an extensive
quantity of data. Furthermore, the usage of videos has gone up significantly
in recent times. Added to it is the fact that the most common form of
surveillance is done by surveillance cameras. The aforementioned factors
combined pave the path towards revered video datasets focusing on action and
attribute recognition, namely [4, 5, 6]. Additionally, surveillance cameras
tend to provide a stationary point of view (PoV), resulting limited field of
view (FoV).

Unmanned Aerial Vehicles (UAV) have become popular recently due to their
precise and rapid movement capabilities. Although there are some datasets
captured from aerial PoV such as UCF-ARG [7] and Okutama-action [8], they
focus on general human activity. In a recent work that focuses on violent
human activity recognition [9], the researchers use dataset that was collected
by the themselves, which further proves the lack of such data. Besides,
collecting real life data on violence can be threatening to the person trying
to capture the activity. Hence, we would like to propose our solution to the
aforementioned problems in the form of a simulator which is capable of:

  * •

Generating photo-realistic RGB images captured from the perspective of an UAV,

  * •

Using procedural method for generating randomized urban environment,

  * •

Generating various forms of ground truth images,

  * •

Generating bounding box around crowd groups using semantic segmentation,

  * •

Providing control over the velocity and altitude of the UAV via input e.g.
joystick or RC.

The proposed crowd simulator is developed using Unreal Engine111Unreal Engine:
https://unrealengine.com, which is an open source game development engine. In
the upcoming sections, we describe the relevant literature to our work, our
methodology and some qualitative and quantitative results to demonstrate the
effectiveness of the proposed simulator.

##  II Related Work

The aim of our proposed simulator is two-fold. Firstly, we want to be able to
simulate violent and non-violent crowd activity in a randomized manner so that
we can obtain data related to it. And secondly, we want to be able to obtain
the data from the perspective of an UAV, which will enable us to train models
that can be applied in violence detection using UAVs as a medium of
surveillance in key locations.

###  II-A Existing crowd simulation techniques

There has been a lot of related work in the field of crowd simulation, but
most of them only involve neutral pedestrians [10, 11]. [12] deals with
modeling how crowd should move in response to various triggers and changes in
their surroundings by making use of Behavioral Finite State Machines. The
focus of the authors mainly lies in developing realistic crowd movement by
making use of goal selection, trajectory estimation and obstacle avoidance.
These are built into Unreal Engine by default as part of its core library, as
a result of which an Non-player character (NPC) can reach locations within the
world from a starting location in a realistic way when controlled by the AI.
[13] is a framework for crowd simulation interconnecting robots, yet it is
limited ground robots and is based on Menge framework for the crowd simulation
process. In [14], the authors develop a simulation environment using Unity,
which is also a video game development. They aim to generate semantic
segmentation images from real life urban scenarios captured in the form of RGB
images. They show that the synthetic dataset improves the classification
accuracy of the Deep Convolutional Neural Network (DCNN) when trained to
generate semantic segmentation images. [15] publish a dataset that is
generated using the video game Grand Theft Auto V published by Rockstar
Games222Rockstar Games: https://www.rockstargames.com/. Rockstar Games allows
third-party mods to be installed and run inside the game environment, which
allowed the researchers to instantiate environment with crowd.

Figure 2: The pipeline for designing our proposed simulator

###  II-B Existing simulators and datasets

In existing datasets related to crowd, the PoV and quantity is limited. They
mostly rely on RGB images only [16], [17]. Moreover, in our particular case,
the existence of videos from aerial viewpoint is scarce, and there aren’t many
datasets that deal with violent behavior of crowd. [18] propose a dataset
which has 123 sequences captured from an UAV perspective and simulator built
with Unreal Engine, but their focus is mostly single object tracking. [19]
propose a framework for generating synthetic data from Unreal Engine by making
use of environment modeling, activity modeling and using the generated data
for training machine learning models and data augmentation. But it is limited
to a specific type of activity recognition which is pertaining to human
interaction with cars. [20] propose a simulator designed using Unreal Engine
that is capable of generating real enough synthetic data of cars with a view
to solving the domain shift problem [21]. [22] propose a framework for
generating environments with minimal user interaction and a pipeline for
collecting image and ground truth data. [23] use Unreal Engine and character
models from Mixamo333Mixamo: https://www.mixamo.com/ to train a deep
reinforcement learning model to capture frontal faces of characters.

None of the above works focus specifically on modelling crowd behavior in
violent and non-violent contexts in a procedurally generated simulation
environment, which is the aim of our proposed simulator. We are inspired by
[23] in terms of their level design technique, hence we extend it by
introducing artificially intelligent characters and leverage procedural
randomized level generation to meet the goals of our simulator.

##  III Simulator Design

In this section, we discuss our approach in developing the simulator. The main
workflow is in Unreal Engine using Airsim [24] plugin to support UAV flight
and image capture. We use AI Behavior Toolkit444AI Behavior Toolkit:
https://www.unrealengine.com/marketplace/en-US/product/ai-behavior-toolkit
available in the Unreal Engine Marketplace555Unreal Engine Marketplace:
https://www.unrealengine.com/marketplace/en-US/store to model artificially
intelligent Non-player Characters (NPC). For creating environments, we make
use of various Unreal Engine Asset Packs available in the Unreal Engine
Marketplace. The complete pipeline of the proposed system is shown in Figure
2.

###  III-A Level Generation

When generating the simulated world, we adopt a procedural level generationa
approach which is an active research topic in the gaming community [25, 26].
Our proposed framework is quite similar to existing grid based approaches
using hexagon shaped blocks. The associated pseudocode for level generation
and character spawning can be found in the supplementary material.

###  III-B Crowd Agents, Appearance and Behavior

By crowd agents, we refer to the NPCs that are part of the simulation world
and have the ability to navigate within the simulated world. We want the crowd
agents to be fully dynamic and be able to interact intelligently with their
surroundings. The interaction can vary depending on the context e.g. violent
or non-violent. We accomplish such behavior modeling using AI Behavior
Toolkit. To make the characters realistic we use Adobe Fuse666Adobe Fuse:
https://www.adobe.com/products/fuse.html and Mixamo for character animations.

###  III-C Data and Semantic Ground Truth Generation

Once the simulator is launched, we use the Python API provided by Airsim to
spawn an UAV in the simulated world, which can be controlled using API or
joystick. We make API calls to request images from the render thread that
captures the specified image types (RGB, Segmentation, Depth-Perspective etc.)
and save them locally. The necessary calculations for generating segmentation
and depth images are handled by the Airsim plugin. We include a demo video of
the simulator tool in the supplementary materials for interested readers. We
are also able to generate ground truth boxes around crowds by assigning a
fixed color with RGB code (81, 13, 36) to the NPCs in the segmentation images,
which are then processed to obtain ground truth for the corresponding RGB
frames.

Figure 3: Ground truth is computed from the segmentation image (left column) and reapplied to the corresponding RGB frame (right column) TABLE I: A comparison of the proposed simulator with other available crowd datasets and simulators Name of Dataset |  No. of Videos |  No. of Frames | Resolution | Is simulator? |  Has Segmentation Image? |  Has Depth Image? |  Contains Violence?
---|---|---|---|---|---|---|---
Crowd Collectiveness [27] | 413 | 40,796 | 670x1000 | No | No | No | No
Data-driven crowd [28] | 212 | 121,626 | 720x480 | No | No | No | No
Violent Flows [29] | 246 | 22,074 | 320x240 | No | No | No | Yes
WWW [30] | 10,000 | 8M | 640x360 | No | No | No | Yes
UCF [31] | 46 | 18,196 | Variable | No | No | No | No
Okutama-action [8] | 43 | 77,365 | 3840x2160 | No | No | No | No
LCrowdV [32] |  $>$1M |  $>$20M | Any | Yes | No | No | No
Our Approach | Any | Any | Any | Yes | Yes | Yes | Yes

##  IV Evaluation Metrics

We evaluate the proposed simulator in a three-fold manner. First, we try to
evaluate how many concurrent active NPCs can be rendered on screen without
impacting the framerate too much. Next, we do visual comparison between some
real world images and simulator-generated images. Finally, we perform
quantitative comparison where we try to augment the simulator generated data
with real world data and see if it improves model performance in binary video
classification.

###  IV-A Framerate

Our simulator is capable of maintaining 25 frames/second when rendering 150
characters simultaneously. The benchmark is done on a mid-range PC with an
Intel Core i5-6500 @3.20 GHz, a RX 470 with 4GB VRAM and 16GB RAM. Figure 4
shows the corresponding graph that depicts the relationship between the number
of characters and framerate.

Figure 4: Average frame rate per second vs NPC count in simulation

###  IV-B Qualitative Evaluation

Figure 5 presents a qualitative comparison of the frames produced by the
simulator with real world images. It is to be noted that the goal of the
simulator is not to mimic every possible detail but to be able to simulate the
crowd behavior as closely as possible. Character models with greater detail,
more realistic and complex animations along with an increase in quantity of
both will yield better and even more realistic synthetic data.

Figure 5: A visual comparison of real world images (top row) with synthetic
data (bottom row).

###  IV-C Quantitative Evaluation

In this section, we lay down the methodology in which we evaluate the data
from the proposed simulator. We go over the datasets we use for our evaluation
strategy and discuss the effects of augmenting the simulator data with the
datasets in an endeavor to improve video classification performance.

####  IV-C1 Dataset

As we have mentioned earlier, even though there are various video datasets
that have violent and non-violent sequences, crowd activity datasets are rare
and most of those datasets are not captured from a UAV point-of-view. [9]
introduce a dataset called Aerial Violent Individual (AVI), however, that
dataset deals with images of violent individuals and not videos. [8] propose a
high resolution dataset that contains sequences captured from a UAV
perspective, but it deals with human action recognition and contains no
violent sequences. As such, we use the datasets mentioned below for
evaluation. It is to be noted that most of the sequences in the real world
video datasets are not captured from a UAV point-of-view.

##### Violent Flows [29]

The Violent Flows dataset is a classic dataset that has 246 total video
sequences of various resolution and length, 123 of which are violent while the
rest are non-violent.

##### Movie Fights [33]

This is yet another classic dataset that is composed of 200 video sequences
equally divided into violent and non-violent categories collected from action
movies.

##### Automatic Violence Detection in Videos [34]

This is a more recent dataset that has a total of 350 video clips. 120 among
them are non-violent, while the rest are violent sequences.

##### Synthetic Dataset

This is the dataset that is generated from our simulator. There are a total of
240 video sequences divided equally into violent and non-violent categories.
All videos are captured from the perspective of an UAV at an approximately 3.5
meters above ground level in real world units within the Unreal Engine
simulation. The violent sequences include punching, kicking, shooting and
chasing, while non-violent sequences include talking, walking, dispersing and
dancing. We augment this data with the real world datasets to find out whether
adding it benefits classification performance.

####  IV-C2 Preprocessing the datasets

All videos, whether real world or synthetic are normalized to have a
resolution of 640x480 at 30 frames per second, and are no longer than 10
seconds. Next, equal numbers of violent and non-violent sequences are chosen
from the total dataset to ensure class balance. They are then divided into
train, validation and test sets ensuring no data leakage among the sets. The
test set is held out and not used in training or validation. For each dataset,
20% of the data is kept for test set and another 20% from the remaining data
is kept as validation set. Table II shows the clip distribution across all
datasets after preprocessing.

TABLE II: Number of clips per dataset after preprocessing Name of Dataset | Train | Validation | Test
---|---|---|---
Violent Flows | 158 | 38 | 50
Movie Fights | 128 | 32 | 40
Automatic Violence Detection in Videos | 154 | 38 | 48
Synthetic | 154 | 38 | 48

####  IV-C3 Network Architecture

For training and inference, we use two models. All models are trained for 30
epochs on each set of data and then evaluated on the held out test data. All
the models were trained using the MMAction2 toolbox [35].

##### Temporal Segment Networks (TSN) [36]

TSN is a state of the art model for action recognition which is based on the
idea of formulating a temporal structure using temporal sampling and video-
level supervision. We train the network on the datasets using a ResNet50 [37]
backbone. A learning rate of $7.8\text{\times}{10}^{-5}\text{\,}\mathrm{i}$s
used.

##### Inflated 3D ConvNet (I3D) [38]

I3D is another state of the art model that operates by expanding the dimension
of filters and pooling layers of 2D ConvNet layers to 3D. Similar to TSN, we
use a ResNet50 backbone and a learning rate of
$7.8\text{\times}{10}^{-5}\text{\,}$

##  V Experimental Results

In this section we discuss the results obtained from our experiments. Moving
forward, we refer to the datasets using the following acronyms: Violent Flows
- VF, Movie Fights - MF, Automatic Violence Detection in Videos - AVD and
Synthetic dataset - S.

###  V-A Individual datasets

We use each of the datasets to establish a baseline and see how well the
models perform when it is just trained on a single dataset and tested on test
set from the same dataset. Results are shown in Table III. It is seen that
almost all models perform quite well on the datasets.

TABLE III: Baseline performances of the models on the 4 datasets Name of Dataset | Model | Test Set |  Test Accuracy
---|---|---|---
VF |  TSN I3D | VF |  88.0% 83.8%
MF |  TSN I3D | MF |  92.5% 88.0%
AVD |  TSN I3D | AVD |  91.7% 87.5%
S |  TSN I3D | S |  98.0% 100.0%

###  V-B Augmenting synthetic data with real world data

Next, we start augmenting the synthetic data with real world data to see if it
helps in improving classification performance. In Table IV, we see that adding
synthetic data nets us a few more predicitons correct across most datasets. We
see a maximum gain of 8.2% and a minimum gain of 2.1%.

TABLE IV: Results after augmenting synthetic data with real world data Name of Dataset | Model | Test Set |  Test Accuracy
---|---|---|---
VF + S |  TSN I3D | VF |  94.0% 92.0%
MF + S |  TSN I3D | MF |  95.0% 92.5%
AVD + S |  TSN I3D | AVD |  93.8% 93.8%

###  V-C Augmenting real world data with real world data

After that, we combine data from the datasets together to see how much
improvement occurs if real world data is used instead of synthetic data. This
is to evaluate if the models are performing well simply because of having more
training data from the synthetic dataset. Table V shows us that the
performance is worse in some instances and see no improvements in others. This
is because, the real world videos have random camera movements and jerky
motion compared to the data from the simulator, which has more stable camera
position. This contributes to the classifier failing to make predictions
correctly, since the quality of the data is not as good.

TABLE V: Results after increasing real world available data and using no synthetic data Name of Dataset | Model | Test Set |  Test Accuracy
---|---|---|---
VF + MF |  TSN I3D | VF |  88.0% 96.0%
VF + MF |  TSN I3D | MF |  92.5% 87.0%
AVD + VF |  TSN I3D | VF |  88.0% 85.0%
AVD + MF |  TSN I3D | AVD |  93.8% 95.8%
AVD + VF |  TSN I3D | MF |  88.0% 90.0%

###  V-D Testing domain difference

Finally, we switch up the domains of the data. That is, we use models trained
on one dataset on test data from another dataset. Table VI illustrates the
results. This is where the domain difference of the simulated data and real
world data is evident, as all models perform really poorly. We intend to
address this issue by experimenting with domain adaptation techniques to see
if the results can be improved in a future work.

TABLE VI: Results of domain shift Name of Dataset | Model | Test Set |  Test Accuracy
---|---|---|---
VF |  TSN I3D | S |  43.0% 37.8%
MF |  TSN I3D | S |  65.0% 53.8%
AVD |  TSN I3D | S |  48.0% 45.5%
S |  TSN I3D | VF |  50.0% 56.8%
S |  TSN I3D | MF |  48.0% 52.0%
S |  TSN I3D | AVD |  50.0% 50.0%

Overall, we are able to confirm that the data from the simulator is benefiting
the classification accuracy in most instances. The associated Receiver
Operating Characteristics curves for TSN and I3D models are shown in Figure 6.

Figure 6: ROC curves of the models when trained with various combinations of
the mentioned datasets

##  VI Limitations, Conclusion and Future Work

We propose a simulator that is capable of producing photo-realistic image
sequences very near to real world environment in the context of violent and
non-violent crowd activity captured from the perspective of an UAV. Leveraging
a procedural level generation approach, we are able to create randomized urban
environment, deploy AI-controlled NPCs in it and start generating data. Using
segmentation images, we are able to provide bounding box annotations for the
crowd in a particular frame automatically without the need of human
intervention. Our work is presented as a baseline which is open for a wide
range of extensions. Currently, the simulator lacks in the number of character
models and diversified animations, which could be a point of argument that the
data from the simulator is simplistic. If the number of character models and
animations could be increased, it can greatly add to the realism of the
simulated environment. Adding vehicles, destructible and flammable objects
will open up room for more complex crowd behavior modelling. Additionally, we
would like to extend our work in the future by introducing instance
segmentation to the simulator, which will allow us to obtain bounding box data
for individual NPCs, as well as generate videos from the perspective of a
moving UAV. We also would like to experiment with domain adaptation techniques
and see if the classification results can be further improved.

## References

  * [1] R. D. Harp and K. B. Karnauskas, “Global warming to increase violent crime in the united states,” _Environmental Research Letters_ , vol. 15, no. 3, p. 034039, 2020.
  * [2] T. Deb, A. Arman, and A. Firoze, “Machine cognition of violence in videos using novel outlier-resistant vlad,” in _2018 17th IEEE International Conference on Machine Learning and Applications (ICMLA)_. IEEE, 2018, pp. 989–994.
  * [3] P. Zhou, Q. Ding, H. Luo, and X. Hou, “Violence detection in surveillance video using low-level features,” _PLoS one_ , vol. 13, no. 10, 2018.
  * [4] S. Abu-El-Haija, N. Kothari, J. Lee, P. Natsev, G. Toderici, B. Varadarajan, and S. Vijayanarasimhan, “Youtube-8m: A large-scale video classification benchmark,” _arXiv preprint arXiv:1609.08675_ , 2016.
  * [5] G. Castanon, N. Shnidman, T. Anderson, and J. Byrne, “Out the window: A crowd-sourced dataset for activity classification in surveillance video,” _arXiv preprint arXiv:1908.10899_ , 2019.
  * [6] M. Monfort, A. Andonian, B. Zhou, K. Ramakrishnan, S. A. Bargal, T. Yan, L. Brown, Q. Fan, D. Gutfreund, C. Vondrick _et al._ , “Moments in time dataset: one million videos for event understanding,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 42, no. 2, pp. 502–508, 2019.
  * [7] A. Nagendran, D. Harper, and M. Shah, “New system performs persistent wide-area aerial surveillance,” _SPIE Newsroom_ , vol. 5, pp. 20–28, 2010\.
  * [8] M. Barekatain, M. Martí, H.-F. Shih, S. Murray, K. Nakayama, Y. Matsuo, and H. Prendinger, “Okutama-action: An aerial view video dataset for concurrent human action detection,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_ , 2017, pp. 28–35.
  * [9] A. Singh, D. Patil, and S. Omkar, “Eye in the sky: Real-time drone surveillance system (dss) for violent individuals identification using scatternet hybrid deep learning network,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_ , 2018, pp. 1629–1637.
  * [10] D. Kaup, T. L. Clarke, L. Malone, and R. Oleson, “Crowd dynamics simulation research,” _Simulation Series_ , vol. 38, no. 4, p. 365, 2006.
  * [11] C. Loscos, D. Marchal, and A. Meyer, “Intuitive crowd behavior in dense urban environments using local laws,” in _Proceedings of Theory and Practice of Computer Graphics, 2003._ IEEE, 2003, pp. 122–129.
  * [12] S. Curtis, A. Best, and D. Manocha, “Menge: A modular framework for simulating crowd movement,” _Collective Dynamics_ , vol. 1, pp. 1–40, 2016.
  * [13] A. Aroor, S. L. Esptein, and R. Korpan, “Mengeros: A crowd simulation tool for autonomous robot navigation,” in _2017 AAAI Fall Symposium Series_ , 2017\.
  * [14] G. Ros, L. Sellart, J. Materzynska, D. Vazquez, and A. M. Lopez, “The synthia dataset: A large collection of synthetic images for semantic segmentation of urban scenes,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 3234–3243.
  * [15] Q. Wang, J. Gao, W. Lin, and Y. Yuan, “Learning from synthetic data for crowd counting in the wild,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2019, pp. 8198–8207.
  * [16] J. C. van Gemert, C. R. Verschoor, P. Mettes, K. Epema, L. P. Koh, and S. Wich, “Nature conservation drones for automatic localization and counting of animals,” in _European Conference on Computer Vision_. Springer, 2014, pp. 255–270.
  * [17] M. A. Olivares-Mendez, C. Fu, P. Ludivig, T. F. Bissyandé, S. Kannan, M. Zurad, A. Annaiyan, H. Voos, and P. Campoy, “Towards an autonomous vision-based unmanned aerial system against wildlife poachers,” _Sensors_ , vol. 15, no. 12, pp. 31 362–31 391, 2015.
  * [18] M. Mueller, N. Smith, and B. Ghanem, “A benchmark and simulator for uav tracking,” in _European conference on computer vision_. Springer, 2016, pp. 445–461.
  * [19] T. S. Kim, M. Peven, W. Qiu, A. Yuille, and G. D. Hager, “Synthesizing attributes with unreal engine for fine-grained activity analysis,” in _2019 IEEE Winter Applications of Computer Vision Workshops (WACVW)_. IEEE, 2019, pp. 35–37.
  * [20] R. Khirodkar, D. Yoo, and K. Kitani, “Domain randomization for scene-specific car detection and pose estimation,” in _2019 IEEE Winter Conference on Applications of Computer Vision (WACV)_. IEEE, 2019, pp. 1932–1940.
  * [21] S. J. Pan and Q. Yang, “A survey on transfer learning,” _IEEE Transactions on knowledge and data engineering_ , vol. 22, no. 10, pp. 1345–1359, 2009.
  * [22] M. Müller, V. Casser, J. Lahoud, N. Smith, and B. Ghanem, “Sim4cv: A photo-realistic simulator for computer vision applications,” _International Journal of Computer Vision_ , vol. 126, no. 9, pp. 902–919, 2018.
  * [23] A. Tzimas, N. Passalis, and A. Tefas, “Leveraging deep reinforcement learning for active shooting under open-world setting,” in _2020 IEEE International Conference on Multimedia and Expo (ICME)_. IEEE, 2020, pp. 1–6.
  * [24] S. Shah, D. Dey, C. Lovett, and A. Kapoor, “Airsim: High-fidelity visual and physical simulation for autonomous vehicles,” in _Field and service robotics_. Springer, 2018, pp. 621–635.
  * [25] O. Roglà Pujalt, N. Pelechano Gómez, and G. A. Patow, “Procedural semantic cities,” in _CEIG 2017: XXVII Spanish Computer Graphics Conference: Sevilla, Spain, June 28-30, 2017_. European Association for Computer Graphics (Eurographics), 2017, pp. 113–120.
  * [26] B. Maxim, D. Ciugurean, and D. Gorgan, “Enhancing player experience using procedural level generation.” in _RoCHI_ , 2018, pp. 141–146.
  * [27] B. Zhou, X. Tang, and X. Wang, “Measuring crowd collectiveness,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2013, pp. 3049–3056.
  * [28] M. Rodriguez, J. Sivic, I. Laptev, and J.-Y. Audibert, “Data-driven crowd analysis in videos,” in _2011 International Conference on Computer Vision_. IEEE, 2011, pp. 1235–1242.
  * [29] T. Hassner, Y. Itcher, and O. Kliper-Gross, “Violent flows: Real-time detection of violent crowd behavior,” in _2012 IEEE Computer Society Conference on Computer Vision and Pattern Recognition Workshops_. IEEE, 2012, pp. 1–6.
  * [30] J. Shao, K. Kang, C. Change Loy, and X. Wang, “Deeply learned attributes for crowded scene understanding,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2015, pp. 4657–4666.
  * [31] S. Ali and M. Shah, “A lagrangian particle dynamics approach for crowd flow segmentation and stability analysis,” in _2007 IEEE Conference on Computer Vision and Pattern Recognition_. IEEE, 2007, pp. 1–6.
  * [32] E. Cheung, T. K. Wong, A. Bera, X. Wang, and D. Manocha, “Lcrowdv: Generating labeled videos for simulation-based crowd behavior learning,” in _European Conference on Computer Vision_. Springer, 2016, pp. 709–727.
  * [33] E. B. Nievas, O. D. Suarez, G. B. García, and R. Sukthankar, “Violence detection in video using computer vision techniques,” in _International conference on Computer analysis of images and patterns_. Springer, 2011, pp. 332–339.
  * [34] M. Bianculli, N. Falcionelli, P. Sernani, S. Tomassini, P. Contardo, M. Lombardi, and A. F. Dragoni, “A dataset for automatic violence detection in videos,” _Data in Brief_ , vol. 33, p. 106587, 2020.
  * [35] D. L. Yue Zhao, Yuanjun Xiong, “Mmaction,” https://github.com/open-mmlab/mmaction, 2019.
  * [36] L. Wang, Y. Xiong, Z. Wang, Y. Qiao, D. Lin, X. Tang, and L. Van Gool, “Temporal segment networks: Towards good practices for deep action recognition,” in _European conference on computer vision_. Springer, 2016, pp. 20–36.
  * [37] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” _arXiv preprint arXiv:1512.03385_ , 2015.
  * [38] J. Carreira and A. Zisserman, “Quo vadis, action recognition? a new model and the kinetics dataset,” in _proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2017, pp. 6299–6308.