# HRFNet: High-Resolution Forgery Network for Localizing Satellite Image
Manipulation

###### Abstract

Existing high-resolution satellite image forgery localization methods rely on
patch-based or downsampling-based training. Both of these training methods
have major drawbacks, such as inaccurate boundaries between pristine and
forged regions, the generation of unwanted artifacts, etc. To tackle the
aforementioned challenges, inspired by the high-resolution image segmentation
literature, we propose a novel model called HRFNet to enable satellite image
forgery localization effectively. Specifically, equipped with shallow and deep
branches, our model can successfully integrate RGB and resampling features in
both global and local manners to localize forgery more accurately. We perform
various experiments to demonstrate that our method achieves the best
performance, while the memory requirement and processing speed are not
compromised compared to existing methods.

Index Terms—  Forgery Detection, High-Resolution Image, Image Manipulation,
Satellite Image

##  1 Introduction

The advent of satellites equipped with advanced imaging technology has enabled
a wide range of applications that utilize satellite imagery, such as
agricultural crop classification [1], wildlife monitoring, etc. Concurrently,
the proliferation of image and video editing software, along with advanced AI
techniques such as deepfakes [2, 3], has led to the creation of fake images.
Consequently, image forgery has emerged as a significant socio-technical
concern. Satellite images are also susceptible to manipulation using these
tools, highlighting the need for automated methods to detect tampered areas.
However, accurately identifying tampered areas with different types of forgery
(including splicing, copy-move, removal, etc.) remains highly challenging.

Typically, satellite images are high-resolution images with large spatial
sizes. On the other hand, existing general-purpose image forgery localization
methods are designed for lower-resolution images. Therefore, using these
methods directly on high-resolution satellite images is not feasible due to
resource constraints. To address this issue, several approaches [4, 5, 6, 7]
specifically targeting satellite image forgery localization have been
proposed. These approaches often involve patch extraction from the satellite
image or downsampling the high-resolution image before inputting it to the
forgery localization network for segmentation.

However, such approaches have their limitations. Chen et al. [8] have shown
that patch-based or downsampling-based training for the segmentation of high-
resolution images has major drawbacks. Specifically, patch-based training
lacks spatial contexts and neighborhood dependency information, making it
difficult to distinguish between different classes. On the other hand,
downsampling-based training suffers from “jiggling” artifacts and inaccurate
boundaries due to the missing details from downsampling. To further tackle
these challenges and facilitate high-resolution image segmentation on
satellite images, a few approaches have been proposed. For example, GLNet [8]
uses a downsampled version of the full-resolution input image in the global
branch and extracts patches from the same input image into the local branch,
with the aim that the two branches will collaborate with each other. MBNet [9]
improves upon GLNet by taking co-centered patches of various sizes, instead of
only one resolution. In addition, ISDNet [10] takes the full-resolution image
as input but uses a lightweight network to extract features. To improve the
low-level fine-grained features extracted by the lightweight network, a deep
network is also used to extract features from the downsampled full-resolution
image in a parallel branch.

Fig. 1: Overview of proposed HRFNet: The shallow branch consisting of a
lightweight model takes the full resolution satellite image as input. In
contrast, the deep branch consisting of a deep network takes the downsampled
satellite image as input. Both RGB and resampling features are extracted and
fused to output the final prediction mask.

However, methods like ISDNet and MBNet are specifically tailored for the
segmentation of high-resolution images and do not readily translate for usage
in the case of satellite image forgery localization. This is because forgery
localization, although a segmentation task, requires additional features to
accurately pinpoint the forged region within the image. Moreover, no attempt
has been made yet to significantly deviate from the traditional patch-based or
downsampling-based training approaches in this field. As a result of all the
above factors and challenges, existing satellite forgery localization methods
still suffer from drawbacks such as inaccurate boundaries and the generation
of artifacts.

Our work aims to overcome such drawbacks by leveraging training strategy from
high-resolution image segmentation literature and tailoring and translating it
to fit the satellite image manipulation localization task. Specifically,
inspired by the shallow and deep integration of ISDNet [10], we propose a
novel model called “HRFNet”, which successfully integrates resampling and RGB
features from both the shallow branch and deep branch to effectively perform
satellite image manipulation localization. The RGB branches assist in
distinguishing tampered regions from authentic regions by capturing visual
inconsistencies at tampered boundaries, such as the contrast effect between
tampered regions and authentic regions, etc. On the other hand, SRM branches
analyze the local noise features in an image, such as resampling artifacts,
etc.

Both the RGB and SRM branches are further divided into shallow and deep
branches. The shallow branch enables the extraction of features in a global
manner with a wide field of view, whereas the deep branch extracts more
detailed and hierarchical features. By fusing features from both branches, our
model successfully captures image manipulation traces and becomes effective in
localizing high-resolution satellite image manipulation in an efficient
manner. It should also be noted that we only adopt the training strategy of
ISDNet, which is the local-global collaboration strategy. However, the
underlying architecture of HRFNet is fundamentally different compared to
ISDNet. In summary, our main contributions are:

  * •

We propose a novel satellite image manipulation localization model named
HRFNet that is equipped with RGB and SRM branches of shallow and deep nature
for successfully integrating both global contextual information and fine-
grained local information to capture forgery traces from high-resolution
satellite images effectively.

  * •

Our approach deviates from the conventional patch-based or downsampling-based
methods, specifically in the case of satellite image manipulation localization
task, aiming to tackle the drawbacks inherent to the previous approaches.

  * •

We perform several experiments on satellite manipulation dataset to show that
our method achieves the highest performance compared to existing methods,
while the memory footprint and processing speed are not compromised.

Related Works:  Extensive research has been conducted on generic image forgery
localization, which is not specific to satellite images. ManTra-Net [11] is a
prominent forgery localization model that uses a VGG-based feature extractor
and an LSTM-based detection module. It is trained to detect various types of
image manipulation traces. Moreover, CFLNet [12] combines contrastive learning
with cross-entropy loss for more generalizable forgery localization. However,
these methods are designed for lower-resolution images. Therefore,
downsampling or patch extraction is necessary when applied to manipulated
satellite images, which reduces performance.

On the other hand, several research works have been proposed which are
specifically focused on satellite image manipulation localization. One
approach by [4] extracts patches from full-resolution satellite images and
uses a deep belief network (DBN) to reconstruct them. The model detects and
localizes forgery based on the reconstruction error, which becomes significant
for manipulated images. Authors in [6] employ a conditional GAN, where the
generator aims to produce the ground truth mask and the discriminator
classifies forgery using real image-mask pairs and synthetic pairs. In [7],
authors use an auto-encoder to reconstruct patches from satellite images and
apply an SVDD-based one-class classifier to detect and localize forgery.

##  2 Method

In Fig. 1, we present the overall diagram of our proposed HRFNet, where the
inputs are the full resolution RGB image $I$ and its downsampled version
$I^{\prime}$. As $I$ represents the full resolution image, using deep networks
is impractical due to memory constraints. Hence, we utilize a lightweight
network in the shallow RGB branch to extract features from $I$. However, a
shallow network performs poorly in capturing long-range and high-level
semantic cues. To complement the shallow branch features, we introduce a deep
network in the deep RGB branch, which takes $I^{\prime}$ as input. The deep
network effectively extracts high-level features from downsampled images.
However, the deep network loses spatial information due to the downsampled
input image, which is better captured by the shallow network using the full
resolution image input.

Next, we incorporate a fusion module to merge features from the shallow and
deep RGB branches, aiming to enhance segmentation performance through their
complementary nature. While sophisticated fusion methods have been proposed in
the literature [10, 13], it is important to note that our paper primarily
focuses on introducing a novel training strategy for satellite image
manipulation localization. Hence, we adopt a simpler fusion mechanism that
involves concatenating the features.

Additionally, we utilize SRM filters [14] on $I$ to generate the output
$I_{s}$. SRM filters are high pass filters that capture resampling features
and enhance the high-frequency information of the input image, which is
beneficial for forgery localization. $I_{s}$ is then downsampled to obtain
$I_{s}^{\prime}$. Similarly, $I_{s}$ and $I_{s}^{\prime}$ are passed through
the shallow and deep SRM branches respectively. Subsequently, feature
concatenation is again employed to fuse the features from both SRM branches.

Following the fusion module, we apply a conv-relu-conv layer to the final
features from the RGB and SRM branches, followed by another fusion module that
concatenates features, combining the information from both branches. An ASPP
module [15] is then utilized on the resulting fused feature maps to extract
multi-scale information. It has been reported in [16] that global context aids
in gathering more clues for manipulation detection. The ASPP module
contributes to this by extracting information at different scales, thereby
making global context and fine-grained pixel-level context available. The
output of the ASPP module is passed to the decoder, where we employ a
DeepLabv3+ style segmentation head to generate the final segmentation mask.

##  3 Experiments

Table 1: AUC Scores (in %). Methods | AUC (%)
---|---
MantraNet (CVPR’19) | 77.31
DBN (CVPRW’20) | 63.80
Vision Transformer (CVPRW’21) | 73.26
CFLNet (WACV’23) | 83.52
HRFNet (ours) | 87.36

Dataset: There is a scarcity of benchmark high-resolution satellite image
manipulation datasets. Hence, we take the benchmark satellite dataset
Deepglobe [17] and follow the process from [4] to manipulate the satellite
images. Each final manipulated image in the dataset has a resolution of
$1,000\times 1,000$, as described in [4].

Baseline Models: We compare our method with methods that are specifically
designed for satellite image manipulation localization - Vision transformer
with post processing [5] and DBN [4]. We also compare with the prominent image
forgery localization method - MantraNet [11] and the current SOTA image
forgery localization model CFLNet [12]. For MantraNet and CFLNet, we resize
input images to the size specified by their respective papers for comparisons.

Implementation Details: We use ResNet-18 as a network backbone for both deep
RGB and deep SRM branches. MobileNet-v3 is used for the shallow RGB and
shallow SRM branches. For the deep networks, we downsample the input image to
$224\times 224$ pixel-size. We train HRFNet with Adam optimizer with a
learning rate of $1e-3$. We reduce the learning rate by 20% after every 20
epochs and train for 100 epochs. Cross-entropy loss is weighted to provide the
tampered class ten times more weight [12]. We set the batch size to 4 and
train the model on NVIDIA Tesla K80 GPU.

Results: We report the AUC scores (in %) of our method and the baseline models
in Table 1. As shown, our method achieves the highest performance compared to
the other methods, yielding a $3.84\%$ AUC improvement over the current state-
of-the-art methods for generic image manipulation localization. One reason
behind the poor performance of Vision Transformer and DBN is that both the
methods generate heatmap and then threshold on the heatmap to generate the
output mask instead of using any learnable decoder.

To compare the resource requirement for HRFNet, we compute the memory and
processing time (in frames per second or FPS) of our model and the best
performing baseline models of Table 1. The objective of this experiment is to
demonstrate that HRFNet achieves the best result without requiring higher
resource requirements compared to the best baseline models. To measure the GPU
memory usage of a model, we use the command line tool “gpustat”, with the
minibatch size of 1, and avoid calculating any gradients as instructed in [8].

Table 2: Memory requirement in MB and processing speeds in FPS. Our method does not require extra resources compared to baseline models while achieving the best performance. Methods | Memory (MB) | FPS
---|---|---
MantraNet | 682 | 5.18
CFLNet | 1,740 | 43.12
HRFNet (ours) | 1,402 | 33.68

Table 2 shows the memory performance results. Our model achieves comparable
results with the best baseline models that do not take the full resolution
image as an input. Here, CFLNet requires the most amount of memory. Although
we use a two parallel architecture similar to CFLNet, the network in our two
deep branches is ResNet-18, whereas CFLNet uses two ResNet-50 based encoders.
ResNet-50 has more than twice the number of parameters than ResNet-18. Hence,
the overall memory requirement of HRFNet is less than CFLNet. Although
MantraNet requires the least amount of memory, the AUC score of MantraNet is
also the least among the three models. HRFNet achieves a balanced spot
regarding the memory requirement, while the AUC score is the highest. For low-
powered mobile devices, processing speed is also a major factor. We have also
measured the processing time in FPS for the three models in Table 2. In fact,
CFLNet achieves the best processing speed due to its simpler encoder-decoder
type architecture, while HRFNet performs better than MantraNet. Here again, it
is demonstrated that our model achieves a balanced spot and does not need
extra resources compared to best performing baseline models. Based on table 2,
we can easily conclude that the resource requirement for HRFNet is very much
comparable to the baseline models, while it achieves the best localization
performance.

Visualization: In Fig. 2, we present the sample visualizations of the
predicted masks generated by HRFNet and compare them to those of the best
performing baseline model in our experiment - CFLNet. It is evident that the
predicted masks of CFLNet have inaccurate boundaries compared to HRFNet due to
downscaling the input image. This is in accordance with [8], where it has been
claimed that downsampling high-resolution images results in inaccurate
boundaries and unwanted artifacts in the predicted masks. Additionally, when
the forged object’s shape is very small, CFLNet struggles to localize it,
which is expected because downsampling discards important spatial information.
On the other hand, the shallow branch of HRFNet extracts features that have
accurate spatial information, and the deep branch extracts high-level semantic
cues. Complementing both features yields more accurate localization results
for the forged objects using our method.

|  |  |
---|---|---|---
|  |  |
|  |  |
|  |  |
Input | GT Mask | CFLNet | Ours
Fig. 2: Comparison of the predicted mask with the best performing baseline
model - CFLNet. It is evident that the prediction of HRFNet is closer to the
ground truth mask compared to CFLNet.

##  4 Conclusion

In this paper, we propose a novel satellite image manipulation localization
model, HRFNet. Deviating from conventional training strategy that involves
patch extraction or downsampling, our method introduces a new training
strategy specifically targeted for the task of high-resolution satellite image
manipulation localization. Our method exploits shallow and deep integration to
successfully capture both global context and high level semantic clues, as
well as noise features in order to trace and localize forgery. We have
compared our method with existing satellite image localization methods, as
well as SOTA generic image manipulation localization methods and demonstrated
that our method achieves the best performance while the memory requirement and
processing time remain similar to existing methods. We hope our work will
inspire more research in this area for effectively localizing satellite image
manipulation.

## Acknowledgments

This work was partly supported by Institute for Information & communication
Technology Planning & evaluation (IITP) grants funded by the Korean government
MSIT: (No. 2022-0-01199, Graduate School of Convergence Security at
Sungkyunkwan University), (No. 2022-0-01045, Self-directed Multi-Modal
Intelligence for solving unknown, open domain problems), (No. 2022-0-00688, AI
Platform to Fully Adapt and Reflect Privacy-Policy Changes), (No.
2021-0-02068, Artificial Intelligence Innovation Hub), (No. 2019-0-00421, AI
Graduate School Support Program at Sungkyunkwan University), and (No.
RS-2023-00230337, Advanced and Proactive AI Platform Research and Development
Against Malicious deepfakes).

## References

  * [1] Han Gao, Changcheng Wang, Guanya Wang, Qian Li, and Jianjun Zhu,  “A new crop classification method based on the time-varying feature curves of time series dual-polarization sentinel-1 data sets,”  IEEE Geoscience and Remote Sensing Letters, vol. 17, no. 7, pp. 1183–1187, 2019.
  * [2] Binh M Le and Simon S Woo,  “Add: Frequency attention and multi-view based knowledge distillation to detect low-quality compressed deepfake images,”  arXiv preprint arXiv:2112.03553, 2021.
  * [3] Shahroz Tariq, Sowon Jeon, and Simon S Woo,  “Am i a real or fake celebrity? evaluating face recognition and verification apis under deepfake impersonation attack,”  in Proceedings of the ACM Web Conference 2022, 2022, pp. 512–523.
  * [4] János Horváth, Daniel Mas Montserrat, Hanxiang Hao, and Edward J Delp,  “Manipulation detection in satellite images using deep belief networks,”  in Proceedings of the IEEE/CVF conference on computer vision and pattern recognition workshops, 2020, pp. 664–665.
  * [5] János Horváth, Sriram Baireddy, Hanxiang Hao, Daniel Mas Montserrat, and Edward J Delp,  “Manipulation detection in satellite images using vision transformer,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2021, pp. 1032–1041.
  * [6] Emily R Bartusiak, Sri Kalyan Yarlagadda, David Güera, Paolo Bestagini, Stefano Tubaro, Fengqing M Zhu, and Edward J Delp,  “Splicing detection and localization in satellite imagery using conditional gans,”  in 2019 IEEE Conference on Multimedia Information Processing and Retrieval (MIPR). IEEE, 2019, pp. 91–96.
  * [7] János Horváth, David Güera, Sri Kalyan Yarlagadda, Paolo Bestagini, Fengqing Maggie Zhu, Stefano Tubaro, and Edward J Delp,  “Anomaly-based manipulation detection in satellite images,”  networks, vol. 29, pp. 21, 2019.
  * [8] Wuyang Chen, Ziyu Jiang, Zhangyang Wang, Kexin Cui, and Xiaoning Qian,  “Collaborative global-local networks for memory-efficient segmentation of ultra-high resolution images,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2019, pp. 8924–8933.
  * [9] Lianlei Shan and Weiqiang Wang,  “Mbnet: A multi-resolution branch network for semantic segmentation of ultra-high resolution images,”  in ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2022, pp. 2589–2593.
  * [10] Shaohua Guo, Liang Liu, Zhenye Gan, Yabiao Wang, Wuhao Zhang, Chengjie Wang, Guannan Jiang, Wei Zhang, Ran Yi, Lizhuang Ma, et al.,  “Isdnet: Integrating shallow and deep networks for efficient ultra-high resolution segmentation,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2022, pp. 4361–4370.
  * [11] Yue Wu, Wael AbdAlmageed, and Premkumar Natarajan,  “Mantra-net: Manipulation tracing network for detection and localization of image forgeries with anomalous features,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2019, pp. 9543–9552.
  * [12] Fahim Faisal Niloy, Kishor Kumar Bhaumik, and Simon S Woo,  “Cfl-net: Image forgery localization using contrastive learning,”  in Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision, 2023, pp. 4642–4651.
  * [13] Jun Fu, Jing Liu, Haijie Tian, Yong Li, Yongjun Bao, Zhiwei Fang, and Hanqing Lu,  “Dual attention network for scene segmentation,”  in Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, 2019, pp. 3146–3154.
  * [14] Jessica Fridrich and Jan Kodovsky,  “Rich models for steganalysis of digital images,”  IEEE Transactions on information Forensics and Security, vol. 7, no. 3, pp. 868–882, 2012.
  * [15] Liang-Chieh Chen, George Papandreou, Florian Schroff, and Hartwig Adam,  “Rethinking atrous convolution for semantic image segmentation,”  arXiv preprint arXiv:1706.05587, 2017.
  * [16] Peng Zhou, Xintong Han, Vlad I Morariu, and Larry S Davis,  “Learning rich features for image manipulation detection,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2018, pp. 1053–1061.
  * [17] Ilke Demir, Krzysztof Koperski, David Lindenbaum, Guan Pang, Jing Huang, Saikat Basu, Forest Hughes, Devis Tuia, and Ramesh Raskar,  “Deepglobe 2018: A challenge to parse the earth through satellite images,”  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops, 2018, pp. 172–181.


See pages 1-last of Wan_CVPR18.pdf


#  A visual introduction to Gaussian Belief Propagation
A Framework for Distributed Inference with Emerging Hardware.
https://gaussianbp.github.io

Joseph Ortiz1, Talfan Evans1,2, Andrew J. Davison1
1Imperial College London, 2DeepMind
<EMAIL_ADDRESS>

(June 21, 2021)

###### Abstract

In this article, we present a visual introduction to Gaussian Belief
Propagation (GBP), an approximate probabilistic inference algorithm that
operates by passing messages between the nodes of arbitrarily structured
factor graphs. A special case of loopy belief propagation, GBP updates rely
only on local information and will converge independently of the message
schedule. Our key argument is that, given recent trends in computing hardware,
GBP has the right computational properties to act as a scalable distributed
probabilistic inference framework for future machine learning systems.

This is a citable PDF version of the online interactive article hosted at:
https://gaussianbp.github.io. This version is purely for citation purposes and
the article should be read in the online format. In the online version, all
figures are interactive figures and in this PDF version we include screenshots
of the interactive figures.

Figure 1:  INTERACTIVE FIGURE – view at https://gaussianbp.github.io. Gaussian
Belief Propagation performs probabilistic inference iteratively and is
convergent even when messages are passed randomly through the graph. Here GBP
is applied to a geometric grid alignment problem. See an interactive version
of this figure later in the article.

##  1 Introduction

###  1.1 Bayesian Probabilistic Inference

Bayesian probability theory is the fundamental framework for dealing with
uncertain data, and is at the core of practical systems in machine learning
and robotics [23, 11]. A probabilistic model relates unknown variables of
interest to observable, known or assumed quantities and most generally takes
the form of a graph whose connections encode those relationships. Inference is
the process of forming the posterior distribution to determine properties of
the unknown variables, given the observations, such as their most probable
values or their full marginal distributions.

There are various possible algorithms for probabilistic inference, many of
which take advantage of specific problem structure for fast performance.
Efficient inference on models represented by large, dynamic and highly inter-
connected graphs however remains computationally challenging and is already a
limiting factor in real embodied systems. Inference on such graphs will only
become more important as general purpose perception systems strive towards
more heterogeneous and dynamic representations containing abstractions (such
as objects or concepts) with continually changing relationships [5, 11] 111 In
AI, there has been a recent trend towards sparse representations and graphs
which are well-suited for representing sparse high dimensional data. Part of
this trend has been driven by the fact that graphs are a natural
representation in certain domains, and evidence of this is the rise of graph
neural networks and graph-based probabilistic inference frameworks [15, 11].
Another trend is driven by the idea that massive sparsity is required to scale
AI compute and evidence of this is the scaling of neural networks with gated
computation [47, 20] and the development of novel AI processors specifically
for high performance on sparse representations [32]. We believe that inference
on large scale and dynamic probabilistic graphs will be an important part of
this trend towards sparse computing to scale AI compute..

###  1.2 The Bitter Lesson

As we search for the best algorithms for probabilistic inference, we should
recall the “Bitter Lesson” of machine learning (ML) research: “general methods
that leverage computation are ultimately the most effective, and by a large
margin” [51]. A reminder of this is the success of CNN-driven deep learning,
which is well suited to GPUs, the most powerful widely-available processors in
recent years.

Looking to the future, we anticipate a long-term trend towards a “hardware
jungle” [49] of parallel, heterogeneous, distributed and asynchronous
computing systems which communicate in a peer-to-peer manner. This will be at
several levels: across networks of multiple smart devices operating in the
same environment; across the many sensors, actuators and processors within
individual embodied devices; and even within single processor chips themselves
222For single processor chips, a trend towards multicore designs with
distributed on-core memory and ad-hoc communication is rapidly emerging [32,
24]. This design has the key aim of increasing performance while reducing
power usage by minimizing the ”bits x millimetres” of data movement during
computation [52]. .

To scale arbitrarily and leverage all available computation in this ”hardware
jungle”, our key hypothesis is that we need inference methods which can
operate with distributed local processing / storage and message passing
communication, without the need for a global view or coordination of the whole
model 333The opposite proposition would be a global centralized inference
method, for example via a monolithic cloud-based processor. This centralized
approach would not leverage most of the available compute resources and even
if feasible may not be desirable, for communication bandwidth or privacy
reasons..

###  1.3 Gaussian Belief Propagation

Fortunately, a well-known inference algorithm exists which has these
properties: Loopy Belief Propagation [42, 29, 39]. Belief Progagation (BP) was
invented in the 1980s [42], and is an algorithm for calculating the marginals
of a joint distribution via local message passing between nodes in a factor
graph. Factor graphs are a type of bipartite graphical model that connects
variables via factors which represent independent relationships [15, 14].

BP guarantees exact marginal computation with one sequence of forward-backward
message passing in tree-structured graphs [7], but empirically produces good
results when applied to “loopy” graphs with cycles [39]. Unlike Graph neural
networks [45, 9, 4] which learn edge and node updates that are applied over a
fixed number of message passing steps, loopy BP applies probabilistic message
passing updates with iterative convergent behaviour. Although Loopy BP has
been successful in a few applications such as error-correcting codes [34], it
has not as of yet been applied to broader machine learning problems.

One issue that has precluded the use of general Loopy BP is that it lacks
convergence guarantees. However, a second and perhaps more relevant issue
given modern hardware trends, is that its computational properties have not
fitted the dominant processing paradigms of recent decades, CPUs and GPUs.
Consequently, other factor graph inference algorithms have been preferred
which take advantage of global problem structure to operate much more rapidly
and robustly than loopy BP on a CPU.

We believe that now is the right time to re-evaluate the properties of loopy
BP given the rise of graph-like computing hardware and sparse models. Indeed
there has been a recent resurgence in interest in the algorithm [22, 44, 30,
41, 40]; in particular, one notable work [41] demonstrated that BP on novel
graph processor hardware can achieve a 24x speed improvement over standard
methods for bundle adjustment. We focus on the special case of Gaussian Belief
Propagation (GBP) [12], in which all factors are Gaussian functions of their
dependent variables, and we also infer Gaussian-distributed marginals. GBP has
both more extensive mathematical guarantees [17] and stronger empirical
performance [6] than general loopy BP. It is also simple to implement, with
all messages between nodes taking the form of usually low-dimensional
Gaussians, which can be represented by small vectors and matrices.

In this article, we give an introduction to Gaussian Belief Propagation as a
strong general purpose algorithmic and representational framework for large
scale distributed inference. Throughout, we present a range of examples across
1D and 2D geometric estimation and image processing in the form of interactive
simulations which allow the reader to experiment with and understand the
properties of GBP. We hope that these simulations emphasize the key properties
of GBP which we summarize below.

Figure 2: INTERACTIVE FIGURE – view at https://gaussianbp.github.io.

In the remainder of the article, we first introduce the GBP algorithm and show
that it is intimately linked to solving a linear system of equations $Ax=b$.
We then explore 4 key practical details for applying GBP to real problems:
extending GBP to handle non-linear relationships and non-Gaussian data
distributions, using local message schedules and lastly how hierarchical
structure can help convergence.

##  2 Technical Introduction

###  2.1 Probabilistic Inference

Inference is the problem of estimating statistical properties of unknown
variables $X$ from known or observed quantities $D$ (the data). For example,
one might be interested in inferring tomorrow’s weather (X) from historic data
(D), or the 3D structure of an environment (X) from a video sequence (D).

Bayesian inference proceeds by first defining a probabilistic model $p(X,D)$
that describes the relationships between data and variables, and then using
the sum and product rules of probability 444 The sum rule is
$p(X)=\sum_{Y}p(X,Y)$ and the product rule is $p(X,Y)=p(Y\rvert X)p(X)$.  to
form the posterior distribution $p(X\rvert D)=\frac{p(X,D)}{p(D)}$. The
posterior summarizes our belief about $X$ after seeing $D$ and can be used for
decision making or other downstream tasks. Given the posterior, we can compute
various properties of $X$, for example:

  1. 1.

The most likely configuration of the variables $X_{\text{MAP}}=\text{arg
max}_{X}p(X\rvert D)$, or

  2. 2.

The marginal posteriors $p(x_{i}\rvert D)=\sum_{X\setminus x_{i}}p(X\rvert
D)$, which summarize our belief about each individual variable given $D$.

These two calculations are known as maximum a posteriori (MAP) inference and
marginal inference respectively. An important difference is that MAP inference
produces a point estimate while marginal inference retains information about
uncertainty.

###  2.2 Factor Graphs

The Hammersley-Clifford theorem tells us that any positive joint distribution
$p(X)$ can be represented as a product of factors $f_{i}$, one per clique,
where a clique is a subset of variables $X_{i}$ in which each variable is
connected to all others:

$p(X)=\prod_{i}f_{i}(X_{i})~{}.$

Figure 3:  INTERACTIVE FIGURE – view at https://gaussianbp.github.io. Hover
over the factors in the equation to highlight them in the factor graph.

Factorized representations can be very convenient as they expose structure in
a model. Factor graphs are the natural visual representation of this
factorization and can be very useful for reasoning about problems [15].

In the factor graphs in this article, circles and squares represent variable
and factor nodes respectively, with edges connecting each factor to the
variables it depends on. An example of a simple factor graph is shown in
Figure 3. By explicitly representing the factors as nodes in the graph, factor
graphs clearly emphasize the conditional independence structure of the problem
\- the lack of a factor directly connecting two variables means they are
conditionally independent 555 Mathematically, two variables $x_{i}$ and
$x_{j}$ are conditionally independent given all other variables $X_{-ij}$ if:
$p(x_{i},x_{j}|X_{-ij})=p(x_{i}|X_{-ij})p(x_{j}|X_{-ij})~{}.$ An equivalent
way condition is:  $p(x_{i}|x_{j},X_{-ij})=p(x_{i}|X_{-ij})~{}.$ Intuitively,
if $X_{-ij}$ causes both $x_{i}$ and $x_{j}$, then if we know $X_{-ij}$ we
don’t need to know about $x_{i}$ to predict $x_{j}$ or about $x_{j}$ to
predict $x_{i}$. Conditional independence is often written in shorthand as:
$x_{i}\bot x_{j}|X_{-ij}$. .

Factor graphs can also be presented as energy based models [33] where each
factor $f_{i}$ defines an energy $E_{i}\geq 0$ associated with a subset of the
variables $X_{i}$ 666 This formalism is closely related to the Boltzmann
distribution in statistical physics which gives the probability of a state $i$
as a function of the energy of the state and the temperature of the system:
$p_{i}=\frac{e^{-E_{i}/kT}}{\sum_{j}e^{-E_{j}/kT}}~{},$ where $k$ is the
Boltzmann constant, $T$ is the temperature of the system and $j$ sums over all
available states. :

$f_{i}(X_{i})\propto e^{-E_{i}(X_{i})}~{}.$

In energy based models, finding the most likely variable configuration is
equivalent to minimizing the negative log probability or the sum of factor
energies:

$X_{\text{MAP}}=\text{arg min}_{X}-\log p(X)=\text{arg
min}_{X}\sum_{i}E_{i}(X_{i})~{}.$

###  2.3 The Belief Propagation Algorithm

Belief propagation (BP) is an algorithm for marginal inference, i.e. it
computes the marginal posterior distribution for each variable from the set of
factors that make up the joint posterior. BP is intimately linked to factor
graphs by the following property: BP can be implemented as iterative message
passing on the posterior factor graph. The algorithm operates by iteratively
updating a node’s locally stored belief by sending and receiving messages from
neighbouring nodes. Each iteration consists of 3 phases:

Figure 4:

Belief Update: The variable node beliefs are updated by taking a product of
the incoming messages from all adjacent factors, each of which represents that
factor’s belief on the receiving node’s variables.
Factor-to-variable message: To send a message to an adjacent variable node, a
factor aggregates messages from all other adjacent variable nodes and
marginalizes over all the other nodes’ variables to produce a message that
expresses the factor’s belief over the receiving node’s variables.
Variable-to-factor message: A variable-to-factor message tells the factor what
the belief of the variable would be if the receiving factor node did not
exist. This is computed by taking the product of the messages the variable
node has received from all other factor nodes.
These 3 operations fully define the algorithm and their equations are
presented below.

Figure 5: INTERACTIVE FIGURE – view at https://gaussianbp.github.io. Figure 6:
INTERACTIVE FIGURE – view at https://gaussianbp.github.io. Green curves
represent variable belief distributions. Left: BP on a tree. Right: BP with
synchronous updates applied to a graph with a loop.

Belief propagation was originally developed for graphs that are trees 777 For
tree graphs, any two nodes are connected by exactly one path.  and the updates
were designed such that the beliefs converge to the exact marginals after one
sweep of messages from a root node to the leaf nodes and back [42]. For models
with arbitrary conditional independence structure, including cycles or
”loops”, loopy BP [39] iteratively applies the same message passing rules to
all nodes. The simplest variant of loopy BP sends messages from all nodes at
every iteration in a synchronous fashion. The videos in Figure 6 illustrate
how BP is applied to trees and graphs with loops.

As BP was originally developed for trees, its application to loopy graphs was
at first empirical [39]. Theoretical grounds for applying the same update
rules to loopy graphs were later developed [55, 54, 53] that explain loopy BP
as an approximate variational inference method in which inference is cast as
an optimization problem. Instead of directly minimizing the factor energies
(as is done in MAP inference), loopy BP minimizes the KL divergence between
the posterior and a variational distribution which we use as a proxy for the
marginals after optimization. Loopy BP can be derived via constrained
minimization of an approximation of the KL divergence known as the Bethe free
energy [55] (see Appendix A for the derivation).

As the Bethe free energy is non-convex, loopy BP is not guaranteed to converge
and even when it does it may converge to the wrong marginals. Empirically,
however BP generally converges to the true marginals although for very loopy
graphs it can fail [39, 53].

Most interesting problems have loopy structures and so for the remainder of
the article we will use BP to refer to loopy BP. So far, although we have
outlined the BP equations, we have not specified the form of the factors,
messages or beliefs. From here, we focus on Gaussian belief propagation which
is a special form of continuous BP for Gaussian models.

###  2.4 Gaussian Models

We are interested in Gaussian models in which all factors and therefore the
joint posterior are univariate / multivariate Gaussian distributions.
Gaussians are a convenient choice for a number of reasons: (1) they accurately
represent the distribution for many real world events [27], (2) they have a
simple analytic form, (3) complex operations can be expressed with simple
formulae and (4) they are closed under marginalization, conditioning and
taking products (up to normalization).

A Gaussian factor or in general any Gaussian distribution can be written in
the exponential form $p(x)\propto e^{-E(x)}$ with a quadratic energy function.
There are two ways to write the quadratic energy which correspond to the two
common parameterizations of multivariate Gaussian distributions: the moments
form888It’s called the moments form as it is parameterized by the first moment
and second central moments of the distribution. and the canonical form. The
key properties of each of these parameterizations are summarized in the table
in Figure 7.

Figure 7: INTERACTIVE FIGURE – view at https://gaussianbp.github.io.

The canonical form is often preferred when performing inference, for two main
reasons. Firstly, taking a product is simple in the canonical form, so it is
easy to form the posterior from the factors. Secondly, the precision matrix is
sparse and relates closely to the structure of the factor graph.

The precision matrix describes direct associations or conditional dependence
between variables. If entry $(i,j)$ of the precision matrix is zero then
equivalently, there is no factor that directly connects $x_{i}$ and $x_{j}$ in
the graph. You can see this in the default preset graph in Figure 8 where
$\Lambda_{13}=\Lambda_{31}=0$ and $x_{1}$ and $x_{3}$ have no factor directly
connecting them.

On the other hand, the covariance matrix describes induced correlations
between variables and is dense as long as the graph is one single connected
component. Unlike the canonical form, the moments form is unable to represent
unconstrained distributions, as can be seen by selecting the unanchored preset
graph in which there is only relative positional information. We encourage the
reader to check out the other preset graphs and edit the graph to explore
Gaussian models and the relationship with the canonical form.

Figure 8:  INTERACTIVE FIGURE – view at https://gaussianbp.github.io. Click on
the factors to turn them on / off and change the factor values by dragging the
underlined numbers. You can choose from preset graphs referenced in the text
or create your own!

###  2.5 From Gaussian Inference to Linear Algebra

The joint distribution corresponding to a factor graph in which all factors
are Gaussian can be represented as a single multivariate Gaussian distribution
(since the energy terms are additive) in the canonical form:

$P(X)\propto\exp(-\frac{1}{2}X^{\top}\Lambda X+\eta^{\top}X)~{}.$

MAP inference corresponds to computing the parameters $X_{\text{MAP}}$ that
maximize the above joint distribution. As usual, we can compute the gradient
of the log-probability (the total energy):

$\nabla_{X}E=\nabla_{X}\log P(X)=-\Lambda X+\eta~{},$

and solve for $\nabla_{X}E=0$. From here, we see that MAP inference in a
Gaussian system reduces simply to solving
$X_{\text{MAP}}=\Lambda^{-1}\eta=\mu$, which as expected is equal to the mean.

Marginal inference computes the per-variable marginal posterior distributions.
In the moments form, the marginal distribution of $x_{i}$ is:

$p(x_{i})=\int
p(X)dX_{-i}\propto\exp\big{(}-\frac{1}{2}(x_{i}-\mu_{i})^{\top}\Sigma_{ii}^{-1}(x_{i}-\mu_{i})\big{)}~{},$

where the mean parameter $\mu_{i}$ is the $i^{th}$ element of the joint mean
vector and the covariance $\Sigma_{ii}$ is entry $(i,i)$ of the joint
covariance matrix. The vector of marginal means for all variables is therefore
the joint mean vector $\mu=\Lambda^{-1}\eta$ = $X_{\text{MAP}}$ and the
marginal variances the diagonal entries of the joint covariance matrix
$\Sigma=\Lambda^{-1}$.

We can therefore summarize inference in Gaussian models as solving the linear
system of equations $Ax=b\Leftrightarrow\Lambda\mu=\eta$. MAP inference solves
for $\mu$ while marginal inference solves for both $\mu$ and the block
diagonal elements of $\Lambda^{-1}$.

###  2.6 Gaussian Belief Propagation

Having introduced Gaussian models, we now discuss Gaussian Belief Propagation
(GBP) a form of BP applied to Gaussian models. Due to the closure properties
of Gaussians, the beliefs and messages are also Gaussians and GBP operates by
storing and passing around information vectors and precision matrices.

Unlike general loopy BP, GBP is guaranteed to compute the exact marginal means
on convergence, although the same is unfortunately not true for the variances
which often converge to the true marginal variances, but are sometimes
overconfident for very loopy graphs [54]. Although GBP does not in general
have convergence guarantees, there some convergence conditions [6, 17, 48] as
well as methods to improve chances of convergence (see chapter 22 in [38]).

The interactive Figure 9 aims to build intuition for GBP by exploring the
effect of individual messages. For easy visualization and interpretation of
the beliefs, we examine 3 spatial estimation problems with increasing
“loopiness”: a chain, a loop and a grid. Click on a variable node to send
messages to its adjacent variables and observe how neighbouring beliefs are
updated. You will see that GBP converges to the true marginals regardless of
the order in which messages are passed 999 Message damping is commonly used to
speed up and improve chances of convergence in very loopy graphs, such as the
grid problem. Message damping both empirically [39] and theoretically [48]
improves convergence without affecting the fixed points of GBP. The idea
behind message damping is to use momentum to reduce chances of oscillation by
replacing the message at time $t$ with a combination of the message at time
$t$ and time $t-1$:
$\tilde{m}_{t}=m_{t}^{\beta}\tilde{m}_{t-1}^{(1-\beta)}~{},$ which is a
weighted sum in log-space:  $\log\tilde{m}_{t}=\beta\,\log
m_{t}+(1-\beta)\,\log\tilde{m}_{t-1}~{}.$ Standard BP is recovered when the
damping parameter $\beta=1$ and $\beta=0$ corresponds to not updating the
message and sending the message from the previous iteration. Message damping
can be applied to both the variable-to-factor messages and factor-to-variable
messages, however we find that applying it just to factor-to-variable messages
is sufficient. For GBP, message damping corresponds to damping the information
vector and precision matrix as a weighted sum:
$\tilde{\eta}_{t}=\beta\,\eta_{t}+(1-\beta)\,\tilde{\eta}_{t-1}\;\;\;\text{and}\;\;\;\tilde{\Lambda}_{t}=\beta\,\Lambda_{t}+(1-\beta)\,\tilde{\Lambda}_{t-1}~{}.$
(1) .

Figure 9:  INTERACTIVE FIGURE – view at https://gaussianbp.github.io.
Send messages by clicking on a variable node or by pressing the random or
synchronous message buttons at the top. Choose from 3 different graphs with
increasing “loopiness”: chain, loop and grid. All problems are anchored by a
variable node in the top left with a strong prior that sets the absolute
position - this is necessary as all other factors only constrain relative
position. All factors are binary constraints (i.e. they connect two variable
nodes) and we do not draw the factor nodes but instead draw an edge connecting
the two variable nodes whose colour represents the factor energy. We employ
message damping for the grid.

##  3 Beyond the standard algorithm

We have introduced Gaussian Belief Propagation in its basic form as a
probabilistic inference algorithm for Gaussian estimation problems. However,
to solve real practical problems with GBP, we often need a number of extra
details and tricks which we discuss in this section.

###  3.1 Non-Gaussian factors

Although requiring all factors to be Gaussian is a convenient assumption, most
interesting problems involve non-linear relationships and / or non-Gaussian
data distributions, both of which result in non-Gaussian factors. GBP can be
extended to handle these problems by linearizing the non-linear relationships
and using covariance scaling to handle non-Gaussian data distributions.

####  3.1.1 Non-linear Relationships

A factor is usually created given some observed data $d$ that we model as
$d\sim h(X)+\epsilon$, where $h$ simulates the data generation process from
the subset of variables $X$ 101010We are overloading $X$ here by using it to
denote a subset of the variables for ease of notation. and
$\epsilon\sim\mathcal{N}(0,\Sigma_{n})$ is Gaussian noise. Rearranging, we see
that the residual is Gaussian distributed
$r=d-h(X)\sim\mathcal{N}(0,\Sigma_{n})$, allowing us to form the factor with
energy:

$E(X)=\frac{1}{2}(h(X)-d)^{\top}\Sigma_{n}^{-1}(h(X)-d)~{}.$

For linear functions $h(X)=\mathtt{J}X+c$, the energy is quadratic in $X$ and
we can rearrange the energy so that the factor is in the Gaussian canonical
form as:

$E(X)=\frac{1}{2}X^{\top}\Lambda X-\eta^{\top}X\;\;\;\;\text{,
where}\;\eta=\mathtt{J}^{\top}\Sigma_{n}^{-1}(d-c)\;\text{and}\;\Lambda=\mathtt{J}^{\top}\Sigma_{n}^{-1}\mathtt{J}~{}.$

If $h$ is non-linear 111111The function $h$ could be any non-linear function,
for example a trained neural network [10] or a Gaussian process [37]., the
energy is no longer quadratic in $X$ meaning the factor is not Gaussian-
distributed. To restore the Gaussian form, it is standard to use a first-order
Taylor expansion about the current estimate $X_{0}$ to approximate the factor
as a Gaussian:

$h(X)\approx h(X_{0})+\mathtt{J}(X-X_{0})~{}.$

Here $\mathtt{J}$ is now the Jacobian matrix and the factor can be written in
the same form as above but with $c=h(X_{0})-\mathtt{J}X_{0}$ [12].

After linearization, the posterior is a Gaussian approximation of the true
posterior and inference is performed by successively solving linearized
versions of the underlying non-linear problem (as in non-linear least squares
optimization).

To see how this linearization works, consider a robot moving in a plane that
measures the 2D distance and angle to a landmark also in the plane, where the
current estimates for the position of the robot and landmark are $r_{0}$ and
$l_{0}$ respectively, and the observed measurement $d=h(r_{0},l_{0})$. In the
interactive Figure 10, we show both the true non-linear factor and the
Gaussian approximated factor with $r$ held constant at $r_{0}$.

Figure 10: INTERACTIVE FIGURE – view at https://gaussianbp.github.io. Move the
two sliders to change the Gaussian noise model. Units of the covariance values
are distance2 and radians2.

The accuracy of the approximate Gaussian factor depends on the linearity of
the function $h$ at the linearization point. As $h$ reasonably is smooth, the
linear approximation is good close to the linearization point $l_{0}$, while
further away, the approximation can degrade. In practice, during optimization
we can avoid this region of poor approximation by relinearizing frequently. As
GBP is local, a just-in-time approach to linearization [43] can be used in
which factors are relinearized individually when the current estimate of the
adjacent variables strays significantly from the linearization point.

####  3.1.2 Non-Gaussian data distributions

Figure 11: INTERACTIVE FIGURE – view at https://gaussianbp.github.io. Move
your mouse over the figure to set the residual and observe how the scaled
covariance $\Sigma$ of the approximate Gaussian changes. The true data
distribution follows a Huber energy and the blue dots indicate the transition
point.

A second cause of non-Gaussian factors is non-Gaussian data distributions. We
usually model observed data as coming from a Gaussian distribution: $d\sim
h(X)+\epsilon$, by choosing $\epsilon$ to be Gaussian noise. Although this is
generally sensible [27], true data distributions often have stronger tails or
are more tightly peaked. In these cases, to retain the Gaussian form for GBP,
we use an approximate Gaussian data distribution via covariance scaling [1,
12]. The covariance of the approximate Gaussian is chosen such that the
quadratic energy matches the true non-Gaussian energy at that residual, as
shown in Figure 11.

Robust data distributions (or M-estimators), such as the Huber energy [25,
26], are a common class of non-Gaussian data distributions which have greater
probability mass in the tails to reduce sensitivity to outliers. The Huber
energy is a continuous function that is quadratic close to the mean and
transitions to a linear function for large residuals to reduce the energy cost
of outliers 121212 The distribution induced by the Huber energy is a Gaussian
distribution close to the mean and a Laplace distribution in the tails. The
probability density function for the Laplace distribution is:
$p(x;\mu,\beta)=\frac{1}{2b}\exp\big{(}\frac{-|x-\mu|}{b}\big{)}~{}.$ :

$E_{\text{huber}}(r)=\begin{cases}\frac{1}{2}r^{\top}\Sigma_{n}^{-1}r,&\text{if}\
\rvert r\rvert<t\\\ A\;+\;B\rvert r\rvert,&\text{otherwise}~{}.\end{cases}$

The approximate quadratic energy for the Huber distribution can be found by
solving $\frac{1}{2}r^{\top}\Sigma_{\text{sc}}^{-1}r=E_{\text{huber}}(r)$ to
give the diagonal scaled covariance:

$\Sigma_{\text{sc}}=\begin{cases}\Sigma_{n},&\text{if}\ \rvert r\rvert<t\\\
\frac{2E_{\text{huber}}(r)}{r^{\top}\Sigma_{n}^{-1}r}\Sigma_{n},&\text{otherwise}~{}.\end{cases}$

Robust energy functions can make GBP much more generally useful - for example,
they are crucial for bundle adjustment [41] and for sharp image denoising (as
we will see later in Figure 14). More generally, our interpretation is that
robust factors can play a similar role to non-linearities in neural networks,
activating or deactivating messages in the graph.

The interactive Figure 12 below gives intuition for the effect of robust
factors for the task of 1D line fitting. The variables we are estimating are
the $y$ values at fixed intervals along the $x$ axis and the blue circles and
lines show the mean and standard deviation of the beliefs. The red squares are
measurements that produce data factors in the graph and there are also
smoothness factors between all adjacent variable nodes encouraging the $y$
values to be close. You can add your own data factors by clicking on the
canvas and a diagram of the factor graph is in the bottom right of the figure.
Press play to run synchronous GBP and observe that a Huber energy can
disregard outliers and retain step discontinuities in the data unlike the
standard squared loss.

Figure 12: INTERACTIVE FIGURE – view at https://gaussianbp.github.io. 1D line
fitting. Create data factors by clicking on the canvas or use a preset data
factor configuration with the outlier, step and random buttons. Note the
variable nodes start at the bottom of the canvas before they have received any
messages.

###  3.2 Local updates and Scheduling

So far, we have assumed that all variable and factor nodes broadcast messages
at each iteration in a synchronous fashion, where all nodes absorb and
broadcast messages in parallel. In fact, this is far from a requirement and as
GBP is entirely local, messages can be sent arbitrarily and asynchronously.

It turns out that the message schedule can have a dramatic effect on the rate
of convergence. For example, swapping synchronous updates for random message
passing tends to improve convergence, while a fixed “round-robin” schedule can
do even better [28]. Better yet, if each message requires some unit of
computation (and therefore energy), it’s possible to prioritize sending
messages that we think will contribute most to the overall convergence of the
system (there is evidence that the brain may apply a similar economical
principle [19]). This is the idea behind residual belief propagation (RBP)
[18] and similar variants [50, 43], which form a message queue according to
the norm of the difference from the previous message.

In Figure 13, we explore message scheduling using the 1D line fitting task
once again. The underlying factor graph is a chain (no loops) and so will
converge after one sweep of messages from left to right and back again. You
can send messages through the graph using the preset schedules (synchronous,
random or sweep) or create your own schedule by clicking on a variable node to
send messages outwards.

Figure 13: INTERACTIVE FIGURE – view at https://gaussianbp.github.io. 1D line
fitting. Experiment with different message passing schedules. Use the
synchronous, random and sweep schedules or click on a variable node to send
messages to its neighbours. Note the variable nodes start at the bottom of the
canvas before they have received any messages.

Playing around with different schedules for surface estimation highlights two
important properties of GBP. First, GBP can converge with an arbitrary message
passing schedule. As a consequence, GBP can readily operate in systems with no
global clock and varying local compute budgets such as on neuromorphic
hardware or between a group of distributed devices [35].

The second property is that GBP can achieve approximate local convergence
without global convergence. Due to the factorized structure of GBP [16],
global inference is achieved by jointly solving many interdependent local
subproblems. There are many instances in which we might only be interested in
local solutions - in these cases, GBP can operate in a just-in-time or
attention-driven fashion, focusing processing on parts of the graph to solve
local subproblems as the task demands. Local message passing can yield
accurate relative local solutions which estimate the marginals up to global
corrections that come from more distant parts of the graph 131313One simple
example is mapping two connected rooms. An accurate local map of one room can
be constructed by focusing processing on the part of the factor graph in that
room. For some applications this may be sufficient while for others it may be
important to build a map with an accurate absolute position which may require
longer range message passing between the parts of the graph corresponding to
each separate room. . This attention-driven scheduling can be very economical
with compute and energy, only sending the most task-critical messages.

In Figure 14 below we explore attention-driven message passing for image
denoising 141414As there are no long-range connections in the image denoising
graph, local message passing can produce the true local marginals as the
effect of more distant parts of the graph is negligible.. Image denoising is
the 2D equivalent of the surface estimation problem from the Figure 13. The
only difference is that previously although variable nodes were at discrete
locations the data factors were at any location, while now the data factors
are at the same discrete locations as the variable nodes with one per node. We
also revisit the use of robust energy functions with GBP via covariance
scaling which is crucial for sharp denoising.

Figure 14: INTERACTIVE FIGURE – view at https://gaussianbp.github.io. GBP for
image denoising [46]. You can choose either all-to-all message passing or an
attention-driven schedule around the mouse. Choosing a squared loss function
results in a blurred denoising while the robust Huber loss gives a sharper
denoising. You can select from a number of preloaded images with salt and
pepper (S&P) or Gaussian noise or use a live feed from your webcam.

###  3.3 Multiscale Learning

Propagating information from one node to another with GBP takes the same
number of iterations as the number of hops between the nodes. For nearby nodes
in a local region, information can be communicated in a small number of
iterations and consensus can be reached quickly, while for distant nodes, a
global consensus can take many more iterations to be established. This is an
inherent property of local algorithms and can be summarized as low frequency
errors decay more slowly than the high frequency errors.

Regular grid structured graphs appear a lot in computer vision (e.g. image
segmentation) and in discretized boundary value problems (e.g. solving for the
temperature profile along a rod). Accelerating convergence in such grid graphs
has been well-studied in the field of Multigrid methods [8]. One simple
approach is to coarsen the grid which transforms low frequency errors into
higher frequency errors that decay faster. After convergence in the coarsened
grid, the solution is used to initialize inference in the original grid which
now has smaller low frequency errors. This is the idea behind coarse-to-fine
optimization which is used in many grid-based problems where it is simple to
build a coarser graph. In one notable work [21], the authors demonstrate much
faster inference for stereo, optical flow and image restoration with
multiscale BP.

Mulitgrid methods can only be applied to graphs with a grid-like structure
where it is possible to build equivalent coarsened representations. In
general, most problems are more unstructured and it is not clear how to build
a coarsened or abstracted representation of the original problem. In the
general case, we see two possible ways to build hierarchy into a model. A
network could be trained to directly predict specific abstractions that form
long range connections when included in the graph. Second, the graph could
contain additional constraints that define a generic hierarchical structure
(much like a neural network) and then the abstractions themselves are also
inferred [22].

##  4 Related Methods

Solving real non-linear problems with GBP is done by iteratively solving
linearized Gaussian versions of the true non-linear problem. This general
pattern of successively solving linearized problems underpins many different
non-linear inference methods. There are efficient libraries [2, 13, 31] for
non-linear inference which use trust region methods like Gauss-Newton or line
search to guide the repeated linear steps 151515Trust region methods
approximate the energy using a model within a trust region. For example, the
Gauss-Newton method uses a quadratic model meaning the factors are
approximated as Gaussians as in GBP. Line search methods choose a descent
direction and then step size at each iteration. In trust region methods, the
most expensive step is solving the linear system, while for line search
methods choosing the direction is the most expensive part..

GBP is just one of many possible algorithms that can be used to solve the
linearized Gaussian model. To place GBP amongst other methods, we present an
overview of a number of related methods for MAP and marginal inference for
Gaussian models in the table below. As a reminder, inference in Gaussian
models is equivalent to solving the linear system $\Lambda\mu=\eta$, for $\mu$
in MAP inference and for $\mu$ and the diagonal elements of $\Lambda^{-1}$ in
marginal inference. You can hover over the circles in the figure to explore
how GBP relates to other methods.

Figure 15:  INTERACTIVE FIGURE – view at https://gaussianbp.github.io. Hover
over the circles to show comments. Green methods are local, by which we mean
that they can be implemented as a distributed message passing algorithm in
which all nodes compute messages on every iteration and each node communicates
only with a small local neighbourhood. That does not mean to say that the
other methods are entirely centralised; many of the methods are somewhat
parallel or have modified versions that are more strongly parallel (e.g.
parallel conjugate gradient). The KL symbol indicates a variational inference
method. Note that this is not a comprehensive list - for example we have
excluded sampling based inference methods. [References: Expectation
propagation [36], Barfoot 2020 [3]]

With so many different inference methods, choosing which method to use can be
a challenge in itself. Judging the speed of each method is complex and depends
on both the sparsity structure of $\Lambda$ and on the implementation on the
available hardware. Our key argument in this article is that we want a general
method that is local, probabilistic and iterative which led us towards
Gaussian Belief Propagation.

Other notable candidate methods that are local, probabilistic and iterative
are Expectation Propagation (EP) [36] and Barfoot’s algorithm [3]. EP is
generally not node-wise parallel and simplifies to GBP in the special case
when it is node-wise parallel, while Barfoot’s algorithm involves extra
communication edges and is yet to be applied to real problems. For these
reasons GBP stands out as the extreme case that maximizes parallelism and
minimizes communication - two principles that are at the core of scalable and
low-power computation.

##  5 Conclusions

We envisage that ML systems of the future will be large scale, heterogeneous
and distributed and as such will require flexible and scalable probabilistic
inference algorithms. In this article, we argued that Gaussian Belief
Propagation is a strong candidate algorithm as it is local, probabilistic,
iterative and asynchronous. Additionally, we showed 1) how GBP is much more
general with a prescription for handling non-linear factors and robust energy
functions, 2) how GBP can operate in an attention-driven fashion and 3) how
hierarchical structure can help convergence. We hope that this visual
introduction will enourage more researchers and practitioners to look into GBP
as an alternative to existing inference algorithms.

We see many exciting directions for future research around GBP and provide a
GBP Library Colab Notebook as a starting point for the interested reader. Some
directions we are most excited about are improving theoretical guarantees,
using learned factors [10, 37, 40], introducing discrete variables, combining
GBP with GNNs [44, 30], incrementally abstracting factor graphs, investigating
numerical precision for messages, using GBP for distributed learning in
overparameterized networks and lastly unifying iterative inference with test-
time self-supervised learning.

##  6 Supplementary Playground

We encourage the interested reader to explore our supplementary GBP
playground. The playground consists of 2 interative diagrams based around 2D
pose graphs. In the first you can construct your own pose graph, set the
initialization and then choose how messages are passed through the graph. The
second is a simulation of a robot exploring a 2D environment with landmarks.

Figure 16: INTERACTIVE FIGURE – view at https://gaussianbp.github.io.

##  7 Acknowledgments

We are grateful to many researchers with whom we have discussed some of the
ideas in this paper, especially from the Dyson Robotics Lab and Robot Vision
Group at Imperial College London. We would particularly like to thank Xiaofan
Mu, Raluca Scona, Riku Murai, Edgar Sucar, Seth Nabarro, Tristan Laidlow,
Nanfeng Liu, Shuaifeng Zhi, Kentara Wada and Stefan Leutenegger.

## References

  * [1] P. Agarwal, G. D. Tipaldi, L. Spinello, C. Stachniss, and W. Burgard.  Robust map optimization using dynamic covariance scaling.  In Proceedings of the IEEE International Conference on Robotics and Automation (ICRA), 2012.
  * [2] S. Agarwal, M. K., and Others.  Ceres solver.  http://ceres-solver.org.
  * [3] T. D. Barfoot.  Fundamental linear algebra problem of gaussian inference.  arXiv preprint arXiv:2010.08022, 2020.
  * [4] P. W. Battaglia, J. B. Hamrick, V. Bapst, A. Sanchez-Gonzalez, V. Zambaldi, M. Malinowski, A. Tacchetti, D. Raposo, A. Santoro, R. Faulkner, et al.  Relational inductive biases, deep learning, and graph networks.  arXiv preprint arXiv:1806.01261, 2018.
  * [5] Y. Bengio.  The consciousness prior.  arXiv preprint arXiv:1709.08568, 2017.
  * [6] D. Bickson.  Gaussian belief propagation: Theory and Application.  PhD thesis, PhD thesis, The Hebrew University of Jerusalem, 2008.
  * [7] C. M. Bishop.  Pattern Recognition and Machine Learning.  Springer-Verlag New York, Inc., 2006.
  * [8] W. L. Briggs, V. E. Henson, and S. F. McCormick.  A multigrid tutorial.  SIAM, 2000.
  * [9] M. M. Bronstein, J. Bruna, Y. LeCun, A. Szlam, and P. Vandergheynst.  Geometric deep learning: going beyond euclidean data.  IEEE Signal Processing Magazine, 34(4):18–42, 2017.
  * [10] J. Czarnowski, T. Laidlow, R. Clark, and A. J. Davison.  Deepfactors: Real-time probabilistic dense monocular slam.  IEEE Robotics and Automation Letters, 5(2):721–728, 2020.
  * [11] A. J. Davison.  Futuremapping: The computational structure of spatial ai systems.  arXiv preprint arXiv:1803.11288, 2018.
  * [12] A. J. Davison and J. Ortiz.  Futuremapping 2: Gaussian belief propagation for spatial ai.  arXiv preprint arXiv:arXiv:1910.14139, 2019.
  * [13] F. Dellaert.  Factor graphs and gtsam.  Technical Report GT-RIM-CP&R-2012-002, Georgia Institute of Technology, 2012.
  * [14] F. Dellaert.  Factor graphs: Exploiting structure in robotics.  Annual Review of Control, Robotics, and Autonomous Systems, 4, 2021\.
  * [15] F. Dellaert and M. Kaess.  Factor graphs for robot perception.  Foundations and Trends in Robotics, 6(1-2):1–139, 2017.
  * [16] P. U. Diehl, J. Martel, J. Buhmann, and M. Cook.  Factorized computation: What the neocortex can tell us about the future of computing.  Frontiers in computational neuroscience, 12:54, 2018.
  * [17] J. Du, S. Ma, Y.-C. Wu, S. Kar, and J. M. Moura.  Convergence analysis of belief propagation on gaussian graphical models.  arXiv preprint arXiv:1801.06430, 2018.
  * [18] G. Elidan, I. McGraw, and D. Koller.  Residual belief propagation: Informed scheduling for asynchronous message passing.  In In Proceedings of the Conference on Uncertainty in Artificial Intelligence (UAI), 2006.
  * [19] T. Evans and N. Burgess.  Coordinated hippocampal-entorhinal replay as structural inference.  In Neural Information Processing Systems (NeurIPS), volume 32, 2019\.
  * [20] W. Fedus, B. Zoph, and N. Shazeer.  Switch transformers: Scaling to trillion parameter models with simple and efficient sparsity.  arXiv preprint arXiv:2101.03961, 2021.
  * [21] P. F. Felzenszwalb and D. P. Huttenlocher.  Efficient belief propagation for early vision.  International journal of computer vision, 70(1):41–54, 2006.
  * [22] D. George, W. Lehrach, K. Kansky, M. Lázaro-Gredilla, C. Laan, B. Marthi, X. Lou, Z. Meng, Y. Liu, H. Wang, et al.  A generative vision model that trains with high data efficiency and breaks text-based captchas.  Science, 358(6368), 2017.
  * [23] Z. Ghahramani.  Probabilistic machine learning and artificial intelligence.  Nature, 521(7553):452–459, 2015.
  * [24] C.-Y. Gui, L. Zheng, B. He, C. Liu, X.-Y. Chen, X.-F. Liao, and H. Jin.  A survey on graph processing accelerators: Challenges and opportunities.  Journal of Computer Science and Technology, 34(2):339–371, 2019\.
  * [25] P. J. Huber.  Robust estimation of a location parameter.  The Annals of Mathematical Statistics, 35(1):pp. 73–101, 1964.
  * [26] P. J. Huber.  Robust Statistics.  Wiley Series in Probability and Statistics. Wiley-Interscience, 1981.
  * [27] E. T. Jaynes.  Probability theory: The logic of science.  Cambridge university press, 2003.
  * [28] D. Koller and N. Friedman.  Probabilistic graphical models: principles and techniques.  MIT press, 2009.
  * [29] F. R. Kschischang, B. J. Frey, and H.-A. Loeliger.  Factor graphs and the sum-product algorithm.  IEEE Transactions on information theory, 47(2):498–519, 2001.
  * [30] J. Kuck, S. Chakraborty, H. Tang, R. Luo, J. Song, A. Sabharwal, and S. Ermon.  Belief propagation neural networks.  In Neural Information Processing Systems (NeurIPS), 2020.
  * [31] R. Kummerle, G. Grisetti, H. Strasdat, K. Konolige, and W. Burgard.  g 2 o: A general framework for graph optimization.  In 2011 IEEE International Conference on Robotics and Automation, pages 3607–3613. IEEE, 2011.
  * [32] D. Lacey.  Graphcore ipu benchmarks.  URL https://www.graphcore.ai/posts/new-graphcore-ipu-benchmarks, 2019\.
  * [33] Y. LeCun, S. Chopra, R. Hadsell, M. Ranzato, and F. Huang.  A tutorial on energy-based learning.  Predicting structured data, 1(0), 2006.
  * [34] R. J. McEliece, D. J. C. MacKay, and J.-F. Cheng.  Turbo decoding as an instance of pearl’s ”belief propagation” algorithm.  IEEE Journal on selected areas in communications, 16(2):140–152, 1998.
  * [35] B. Micusik and G. Evangelidis.  Ego-motion alignment from face detections for collaborative augmented reality.  arXiv preprint arXiv:2010.02153, 2020.
  * [36] T. P. Minka.  Expectation propagation for approximate bayesian inference.  arXiv preprint arXiv:1301.2294, 2013.
  * [37] M. Mukadam, J. Dong, X. Yan, F. Dellaert, and B. Boots.  Continuous-time gaussian process motion planning via probabilistic inference.  The International Journal of Robotics Research, 37(11):1319–1340, 2018.
  * [38] K. P. Murphy.  Machine learning: a probabilistic perspective.  MIT press, 2012.
  * [39] K. P. Murphy, Y. Weiss, and M. I. Jordan.  Loopy belief propagation for approximate inference: An empirical study.  In Proceedings of the Fifteenth Conference on Uncertainty in Artificial Intelligence, 1999.
  * [40] A. Opipari, C. Chen, S. Wang, J. Pavlasek, K. Desingh, and O. C. Jenkins.  Differentiable nonparametric belief propagation.  arXiv preprint arXiv:2101.05948, 2021.
  * [41] J. Ortiz, M. Pupilli, S. Leutenegger, and A. Davison.  Bundle adjustment on a graph processor.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2020.
  * [42] J. Pearl.  Probabilistic reasoning in intelligent systems: networks of plausible inference.  Morgan Kaufmann, 1988.
  * [43] A. Ranganathan, M. Kaess, and F. Dellaert.  Loopy sam.  In Proceedings of the International Joint Conference on Artificial Intelligence (IJCAI), 2007.
  * [44] V. G. Satorras and M. Welling.  Neural enhanced belief propagation on factor graphs.  In International Conference on Artificial Intelligence and Statistics, pages 685–693. PMLR, 2021.
  * [45] F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, and G. Monfardini.  The graph neural network model.  IEEE transactions on neural networks, 20(1):61–80, 2008.
  * [46] R. Scona.  Belief propagation denoising project, 2021.
  * [47] N. Shazeer, A. Mirhoseini, K. Maziarz, A. Davis, Q. Le, G. Hinton, and J. Dean.  Outrageously large neural networks: The sparsely-gated mixture-of-experts layer.  In Proceedings of the International Conference on Learning Representations (ICLR), 2017.
  * [48] Q. Su and Y.-C. Wu.  On convergence conditions of gaussian belief propagation.  IEEE Transactions on Signal Processing, 63(5):1144–1155, 2015.
  * [49] H. Sutter.  Welcome to the jungle.  URL https://herbsutter.com/welcome-to-the-jungle, 2011.
  * [50] C. Sutton and A. McCallum.  Improved dynamic schedules for belief propagation.  In In Proceedings of the Conference on Uncertainty in Artificial Intelligence (UAI), 2007.
  * [51] R. Sutton.  The bitter lesson.  http://www.incompleteideas.net/IncIdeas/BitterLesson.html, 2019\.
  * [52] V. Sze, Y.-H. Chen, T.-J. Yang, and J. S. Emer.  Efficient processing of deep neural networks: A tutorial and survey.  Proceedings of the IEEE, 105(12):2295–2329, 2017.
  * [53] M. J. Wainwright and M. I. Jordan.  Graphical models, exponential families, and variational inference.  Now Publishers Inc, 2008.
  * [54] Y. Weiss and W. T. Freeman.  Correctness of belief propagation in gaussian graphical models of arbitrary topology.  In Neural Information Processing Systems (NeurIPS), 2000.
  * [55] J. S. Yedidia, W. T. Freeman, Y. Weiss, et al.  Generalized belief propagation.  In NIPS, volume 13, pages 689–695, 2000.

## Appendix A: Variational BP Derivation

This standard derivation of the belief propagation equations follows Yedida et
al [55] in showing that the BP update rules follow from constrained
minimization of an approximate free energy known as the Bethe free energy.

We begin by writing the posterior as a product of the factors:

$p(X)=\frac{1}{Z}\prod_{a}f_{a}(X_{a})=\frac{1}{Z}\prod_{a}e^{-E(X_{a})}~{}.$

The goal of variational inference is to find a variational distribution $q(X)$
that approximates the posterior well by minimizing the Kullback-Leibler
divergence between the variational distribution and the posterior. The KL
divergence is a non-negative asymmetric similarity metric that has a minimum
of 0 when $p=q$.

$\displaystyle KL(q\lvert\rvert p)$
$\displaystyle=\sum_{X}q(X)\log\frac{q(X)}{p(X)}$
$\displaystyle=\sum_{X}q(X)\log q(X)-\sum_{X}q(X)\log p(X)$
$\displaystyle=-H_{q}(X)-\mathop{\mathbb{E}}_{q}[\log p(X)]$
$\displaystyle=-H_{q}(X)-\sum_{a}\mathop{\mathbb{E}}_{q}[\log
f_{a}(X_{a})]+\log(Z)$ $\displaystyle=F(p,q)+\log(Z)$

Above, we defined the free energy:

$F(p,q)=-H_{q}(X)-\sum_{a}\mathop{\mathbb{E}}_{q}[\log f_{a}(X_{a})]~{},$

where the first term is the negative of the entropy and the second term is
known as the average energy because $-\log f_{a}(X_{a})=E(X_{a})$. The free
energy has a minimum value of $-\log(Z)$ and by minimizing this free energy we
can also minimizing the KL divergence.

We first consider the form of the free energy for a tree. For tree graphs, the
distribution $q(X)$ can be written in the form:

$q(X)=\prod_{i}b_{i}(x_{i})^{1-d_{i}}\prod_{a}b_{a}(X_{a})~{},$

where the first product is over variables and the second is over the factors.
$b_{i}(x_{i})$ is the marginal distribution over variable $x_{i}$,
$b_{a}(X_{a})$ is the joint marginal distribution over variables $X_{a}$ that
connect to factor $f_{a}$, and $d_{i}$ is the degree of variable node $i$ (the
number of nodes neighbouring node $i$). Plugging this into the expression for
the entropy, we get:

$H_{tree}(X)=-\sum_{i}(1-d_{i})\sum_{x_{i}}b_{i}(x_{i})\log
b_{i}(x_{i})-\sum_{a}\sum_{X_{a}}b_{a}(X_{a})\log b_{a}(X_{a})~{}.$

Similarly, the average energy can be written as:

$-\sum_{i}\mathop{\mathbb{E}}_{q}[\log f_{i}(X_{i})]=-\sum_{a}b_{a}(X_{a})\log
f_{a}(X_{a})$

Putting this together gives the free energy for tree graphs:

$F_{tree}=-\sum_{i}(d_{i}-1)\sum_{x_{i}}b_{i}(x_{i})\log
b_{i}(x_{i})+\sum_{a}\sum_{X_{a}}b_{a}(X_{a})\log\frac{b_{a}(X_{a})}{f_{a}(X_{a})}~{}.$

For general factor graphs with loops, the $F\neq F_{tree}$. The Bethe
approximation is to use the free energy for a tree to approximate the free
energy for arbitrary loopy graphs. The resulting approximate free energy is
known as the Bethe free energy.

Belief propagation can be derived via minimization of the Bethe free energy
subject to two constraints. The first is a marginalization constraint:
$b_{i}(x_{i})=\sum_{X_{a}\setminus i}b_{a}(X_{a})$, and the second is a
normalization constraint: $\sum_{i}b_{i}(x_{i})=1$. With these constraints, we
can form the Lagrangian and then set the derivates with respect to the
parameters to zero:

$L=F_{Bethe}+\sum_{i}\gamma_{i}\bigg{\\{}1-\sum_{x_{i}}b_{i}(x_{i})\bigg{\\}}+\sum_{a}\sum_{i\in
N(a)}\sum_{x_{i}}\lambda_{ai}(x_{i})\bigg{\\{}b_{i}(x_{i})-\sum_{X_{a}\setminus
i}b_{a}(X_{a})\bigg{\\}}$ $\frac{\partial L}{\partial
b_{i}(x_{i})}=0\;\;\;\;\;\;\Rightarrow b_{i}(x_{i})\propto\prod_{a\in
N(i)}\exp\big{(}\lambda_{ai}(x_{i})\big{)}$ $\frac{\partial L}{\partial
b_{a}(X_{a})}=0\;\;\;\;\Rightarrow b_{a}(X_{a})\propto f_{a}(X_{a})\prod_{i\in
N(a)}\exp\big{(}\lambda_{ai}(x_{i})\big{)}$

We now choose the lagrange multiplier to be
$\exp(\lambda_{ai}(x_{i}))=m_{x_{i}\rightarrow f_{a}}(x_{i})=\prod_{c\in
N(i)\setminus a}m_{f_{c}\rightarrow x_{i}}(x_{i})$. This is the familiar
variable-to-factor message equation and substituting this into the above
equations yields the belief propagation fixed point equations (the first of
which the reader will recognize as the belief update rule).

$b_{i}(x_{i})\propto\prod_{a\in N(i)}m_{x_{i}\rightarrow
f_{a}}(x_{i})\propto\prod_{a\in N(i)}m_{f_{a}\rightarrow x_{i}}(x_{i})$
$b_{a}(X_{a})\propto f_{a}(X_{a})\prod_{i\in N(a)}m_{x_{i}\rightarrow
f_{a}}(x_{i})=f_{a}(X_{a})\prod_{i\in N(a)}\prod_{c\in N(i)\setminus
a}m_{f_{c}\rightarrow x_{i}}(x_{i})$

Using the marginalization condition , we can derive an equation for the
messages in terms of other messages and produce the factor-to-variable message
equation:

$\displaystyle m_{f_{a}\rightarrow x_{i}}(x_{i})$
$\displaystyle=\sum_{X_{a}\setminus x_{i}}f_{a}(X_{a})\prod_{j\in
N(a)\setminus i}\;\prod_{b\in N(j)\setminus a}m_{f_{b}\rightarrow
x_{j}}(x_{j})$ $\displaystyle=\sum_{X_{a}\setminus
x_{i}}f_{a}(X_{a})\prod_{j\in N(a)\setminus i}m_{x_{j}\rightarrow
f_{a}}(x_{j})$

This result tells us that the fixed-points of loopy belief propagation are
local stationary points of the Bethe free energy and because the Bethe energy
is bounded from below, BP always has a fixed point.

BP variants have been developed using more accurate or convex approximations
of the free energy [55], however a detailed discussion of the theory behind BP
is beyond the scope of this article and we refer the reader to [53] for a in
depth review.