Gaussian Ensemble Belief Propagation for Efficient Inference in High-Dimensional, Black-box Systems
Authors: Dan MacKinlay
Russell Tsuchida
Daniel Edward Pagendam
Petra Kuhnert
Venue: ICLR 2025 Poster (2025)
Keywords: factor graph
data assimilation
ensemble kalman filter
gaussian belief propagation
geospatial
probabilistic graphical model
hierarchical model
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)
TL;DR: If we carefully combine the structures of the Ensemble Kalman Gfilter and Gaussian Belief Propagation, the result is an efficient method for inference in hierarchical spatial systems

Abstract: Efficient inference in high-dimensional models is a central challenge in machine learning.
We introduce the Gaussian Ensemble Belief Propagation (GEnBP) algorithm, which combines the strengths of the Ensemble Kalman Filter (EnKF) and Gaussian Belief Propagation (GaBP) to address this challenge.
GEnBP updates ensembles of prior samples into posterior samples by passing low-rank local messages over the edges of a graphical model, enabling efficient handling of high-dimensional states, parameters, and complex, noisy, black-box generative processes.
By utilizing local message passing within a graphical model structure, GEnBP effectively manages complex dependency structures and remains computationally efficient even when the ensemble size is much smaller than the inference dimension --- a common scenario in spatiotemporal modeling, image processing, and physical model inversion.
We demonstrate that GEnBP can be applied to various problem structures, including data assimilation, system identification, and hierarchical models, and show through experiments that it outperforms existing belief propagation methods in terms of accuracy and computational efficiency.

Supporting code is available at https://github.com/danmackinlay/GEnBP}{github.com/danmackinlay/GEnBP

Published as a conference paper at ICLR 2025
GAUSSIAN ENSEMBLE BELIEF PROPAGATION FOR EF- FICIENT INFERENCE IN HIGH-DIMENSIONAL, BLACK- BOX SYSTEMS
Dan MacKinlay* CSIRO's Data6l
Dan het..aaacknnaay@laaall.....
Russell Tsuchida Monash University!
Dan Pagendam CSIRO's Data6l
Petra Kuhnert CSIRO's Data61
ABSTRACT
Efficient inference in high-dimensional models is a central challenge in machine learning. We introduce the Gaussian Ensemble Belief Propagation (GEnBP) al- gorithm, which combines the strengths of the Ensemble Kalman Filter (EnKF) and Gaussian Belief Propagation (GaBP) address this challenge. GEnBP up- dates ensembles of prior samples into posterior samples by passing low-rank local messages over the edges graphical model, enabling efficient handling high- dimensional states, parameters, and complex, noisy, black-box generative processes. By utilizing local message passing within a graphical model structure, GEnBP effectively manages complex dependency structures and remains computationally efficient even when the ensemble size is much smaller than the inference dimension common scenario in spatiotemporal modeling. image processing, and physical model inversion. We demonstrate that GEnBP can be applied to various problem structures, including data assimilation, system identification, and hierarchical mod- els, and show through experiments that it outperforms existing belief propagation methods in terms accuracy and computational efficiency.
Supporting code is available aigthhbbccoddammakkiaaaytttiii
poitzzoor sample
sutpless
Soound bnth
gicund tnth
10 20 300
40 50 50
10
20
310
40
60
80
(a) GaBP MSF=11993 log (8)=
228
(b) GEnBP MSF=0.161 log 1og.pgg
49.6
Figure 1: Prior and posterior samples for latent q in the ID system identification problem (Section 4.1 The GEnBP ensemble has size N=64, and N=64 samples drawn from the GaBP posterior. The GEnBP prior comprises samples; the GaBP prior is sampled from a Gaussian density with the same moments.
pittps://aat http://///anacckiiinng name Work completed while CSIRO's Data61

---

Published as a conference paper at ICLR 2025
INTRODUCTION
We combine the Ensemble Kalman Filter (EnKF) and Gaussian Belief Propagation (GaBP) perform inference in high- high-dimensional hierarchical systems. While both methods are well-established, their combination seems novel and empirically outperforms existing approaches in key problems.
GaBP (Yedidia et al. 2005) is a specific Gaussian message-passing algorithm designed for inference in graphical models. It performs bensity-based inference by leveraging a joint density over the model's state space. This density is approximated as product of Gaussian factors, defining graphical model (Koller Friedman 2009) over the graph G. Sum-product messages passed along G's edges help infer target marginals (Pearl 2008). dlessage-passing methods have been successfully applied various applications, including Bayesian hierarchical models (Wand) 2017), error- cerre-corg codes ( (Forney] 2001, Kschischang et al. (001), and Simultaneous Localization and Mapping (SLAM) tasks DDellaert & Kaess| 2017) 7 Ortiz et al. (021).
Graphical models possess various useful properties, such as permitting inference distributed over many computational nodes (Vchtari et al. 2201)). domain adaptation (Bareinboim & Pearl 2013). a natural means for estimating treatment effects (Pearl et al.] 20166) and online, streaming updates (Dellaert & Kaess 2017 Eustice et al. 2006). They perform well with low-dimensional variables but become computationally expensive in high-dimensional scenarios.
The EnKF (Evensen 2003) is widely used for inference in state-spaut space models. It constructs a posterior sample of hidden state from prior samples by moment-matching the observation--condiiooal distribution. Although the state update relies on Gaussian approximation, EnKF never explicitly evaluates the Gaussian density. This approach works well in high-dimensional settings, such as climate models (Houtekamer & Zhang 2016), where computing densities can be prohibitively expensive. However, the standard EnKF is limited to state filtering and doesn't extend more general inference tasks. In recent times, various techniques have extended EnKF to the system- identification setting, where not only states but latent parameters of dynamics are jointly inferred (e.g. Fearnhead & Künsch [2018] Chen et al.) (002).. EnKF is prized for its ability to handle correlated spatial fields efficiently, as with other spatial Gaussian approximations such as the Integrated Nested Laplace Approximation (Rue et al., 2009) and the related Stochastic Partial Differential Equation approach (Lindgren et al. 2011) 1p. Unlike those methods, it can exploit existing simulators to encode complex physical dynamics, and does not require the Jacobian of the simulator. Thus far, there are few methods which extend the attractive simulator-based high-dimensional inference of the EnKF to the more general setting of graphical models.
GEnBP borrows strength from both EnKF and GaBP, achieving the EnKF's efficiency in high- dimensional data processing and GaBP's capability to handle complex graphical model structures (Table GEnBP is capable of performing inference in models with noisy and moderately non-linear observation processes, unknown process parameters, and deeply nested dependencies among latent variables, scaling to millions of dimensions in variables and observations. The key insight is that, despite the name, Gaussian Belief Propagation as commonly used is not truly generic in the class of possible Gaussian approximations (Section 2.2.3). Like GaBP and EnKF, GEnBP uses Gaussian approximations. However, it relies on empirical samples for statistics over non-linear nodes, similar to EnKF, rather than the linearisation at the mode used in GaBP. By using ensembles to represent Gaussian distributions, GEnBP avoids the computational burden of full covariance matrices, enabling efficient belief propagation with potentially superior accuracy.
This approach is highly relevant for many practical problems, including physical and geospatial systems like computational fluid dynamics, geophysical model inversion, and weather prediction. Such systems typically feature high- nigh-dimensional representations with spatial correlations induced by known physical dynamics that create low-rank structure, which empirical covariance approximations can effectively exploit.
Main Contributions
1. We introduce Gaussian Ensemble Belief Propagation (GEnBP), a novel message-paste passing method for inference in high-dimensional graphical models.
2

---

Published as a conference paper at ICLR 2025
2. GEnBP leverages EnKF-like ensemble statistics instead of the traditional linearization used in GaBP, improving the handling of non-linear relationships and reducing computational complexity.
3. We demonstrate that GEnBP scales higher dimensions than traditional GaBP while achieving comparable or better accuracy in practically important physical problems, such as spatiotemporal modeling and physical model inversion.
4. GEnBP is enabled by our development of:
(a) rank-efficient techniques for propagating high-dimensional Gaussian beliefs efficiently, and
(b) computationally efficient methods for converting between Gaussian beliefs and en- semble approximations without high- nigh-dimensional matrix operations,
which are, to our knowledge, novel contributions in this context.
PRELIMINARIES
We introduce our notation and essential concepts. See Appendix Cor belief propagation background and Appendix for the Ensemble Kalman filter.
2.1 MODELS AS GENERATIVE PROCESSES AND DENSITIES
Random variables are denoted in sans-serif (e.g. x) and their values in bold (e.g. x); we assume all variables have densities, p(æ:).
The model state is partitioned into variables linked by a structural equation model M (Wright 1934).
Specifically, M comprises J generating equations (P,}/=n where each P, is defined as"
Pj F
Dr, R + +R De; x X
Each P, relates an input set I an output 0j.
The ancestral variables. defined as ot =
, Gj have no inputs. We classify variables as:
*
Evidence 8: Observed variables.
Latent L: Unobserved or nuisance variables.
Query 2: Variables whose widencc-condiiionn distribution is of interest.
We assume that 2 d. This structure allows us to represent the joint density of variables as product of conditional densities,
p(x) II p(æo, (..g).
(2)
Our main objective is specifically compute the posterior distribution p( (( & p( æy ===='s ) i.e. to update our beliefs about the ancestral variables by assimil- ating observations of the evidence variables while accounting for the latent variables.
2.2 BELIEF PROPAGATION IN GRAPHICAL MODELS
Graphical models associate graph G with the factorization of the model density (see Equation 3) Here, BP refers to loopy belief propagation in factor graphs (Kschischang et al. 2001 additional details are in Appendix C
The factor graph is constructed by representing each conditional probability p(ze, x as a factor potential S over N₁ 03 U Sg:
p(x) II (x.)))
(3)
g is bipartite, with factor nodes for each f, connected to eachh x in its neighborhood. MYz.
Faah P, may be stochastie; we omit the noise terms for brevity.
3

---

Published as a conference paper at ICLR 2025
BP approximates the marginal (or belief) of query node
bg (axx) x / p(z))da dayp-
Proposition 1 (Belief Propagation on Factor Graphs). By iteratively propagating the messages,
{ (1.N3)
12a'\Im,ma m/
dæ. NgVK
iEMgie
mxp +f5 = II
(6)
SE V
BP approximates the marginals
m. xx 20 (poo day
bg (xx)
(7)
SE.14
Proof: See Yedidia al. 2000).
Though its analysis is complex (Yedidia al. 2005
Weiss & Freeman, 2001). BP achieves state-of-
the-art performance in many applications (Davison & Ortiz 2019).
2.2.1 THE POSTERIOR GRAPH
In practice, we are interested in the evidence-conditional posterior graph G which conditions on X 2x To construct G*, we modify the factors f3 for all j No by replacing them with conditioned factors,
pl Xg x
(8)
Observed variables and their edges are removed. The target marginal is then
p (x
X = x bg- (#x),
with bg- (Ixr) the belief over ancestral variables. See Figure in the Appendix.
2.2.2 NECESSARY OPERATIONS FOR BELIEF PROPAGATION
To compute bg. (I we perform conditioning, marginalization, and multiplication on densities:
Definition 1. Consider state space x: = [ af with density fæx; 0). The following opera-
tions suffice for BP in an observation-condiionn graph ((Eq. 8)):
Conditioning:
{(x 0), a (a.;; 02) == / (æzg 1xx = =x
(9)
Marginalization: J (ax; 0) + (xi; 0x) == / (x;; 0) ddaas
(10)
Multiplication: (x; 0) .f (æx; 0 + f (x: 0¹) = (x:; 0) (zz; 0x)
(11)
2.2.3 GAUSSIAN DISTRIBUTIONS IN BELIEF PROPAGATION
When variable relationships linear with additive Gaussian noise, and ancestral distributions Gaussian, each factor is Gaussian. In these cases, all operations Definition | have closed-form
solutions (Bickson 2009 Yedidia et al. 2000), used in Gaussian Belief Propagation (GaBP)
(Dellaert & Kaess 2017).
We represent Gaussian densities in two forms:
Definition 2 (Gaussian Density Forms).
1 m)
T m)
Moments Form: @M (x; m, K) = 2
exp
- (x
(12)
Canonical Form: (c(2; =
exp
1
r.
Pa
+
n
x
5"h
Tp-
(13)
2x
2
4

---

Published as a conference paper at ICLR 2025
Here, and K are the mean and covariance, while P K and n Pm are the precision and information vectors. Multiplication is simplest in canonical form:
(c(;; n. P) Pc (2x; n', P') de (æ; + n', P. P*)
(14)
For nonlinear simulators P, the joint covariance is non-Gaussian; standard GaBP methods (Eustice et al. 2006, Ranganathan et al. 2007] Dellaert & Kaess 2017) then use the d-method (Dorfman
1938
to approximate covariances via a first-order Taylor expansion:
Cov(g(0)) (0) Cov (,)
where Jg is the Jacobian of g at
The accuracy of the 8-method is challenging to analyze for nonlinear P. Additionally, GaBP scales unfavorably with the dimensionality of nodes, incurring a memory cost (D22) and a time cost of O(D)) whenever D x D covariance matrix is inverted.
2.3
ENSEMBLE KALMAN FILTERING
The Ensemble Kalman Filter (EnKF) reduces the computational burden large 1) by representing prior distributions with ensembles. An ensemble matrix of N samples, X [x] ,..., (N)], with x (2) ~ ΦMs (m, K). We define the ensemble mean = XA and deviation X = -XB, where A = and B 1
By overloading the mean and variance operations to describe the empirical moments of ensembles, define:
EX X
1 N-1 Vary N
V,
Cov Cov Y)
N-1 XYT
(15)
Here, a diagonal matrix known as the nugger term, typically set I Setting or > 0 useful for numerical stability and to encode model uncertainty. These diagonal terms are usually treated as algorithm hyperparameters and also appear in GaBP.
The statistics of ensemble X define an implied Gaussian density:
Φa (z;; X, |arv[X]).
Prior model ensembles are sampled via ancestral sampling using the generative model described in Equation 25
Two of the Belief Propagation (BP) operations from Definition are applicable the EnKF,
Proposition 2. Partition x follows the Gaussian distribution:
x
x
such that X
[X] X Assume the ensemble X
akk xt
X X
Vary Xx Cov Xe)
Cov(Xe, Xx)
X OM
(16)
Vary (Xe)
In ensemble form, conditioning (Eq.
9 is performed as:
X, a + X Cov(X,, Xx) Vary (xx) - X2))
(17)
The computational cost of solving Equation 17 is O(N3 + 22 Dx.). Marginalization (Eq. simply trancation, i.e., X
10 is
Proof: See Appendix D
While the EnKF reduces computational costs and potentially improves approximation accuracy over GaBP for nonlinear relationships, it does not generalise to other model structures without the multiplication operation (Eq. II needed for BP.
In the following sections, we extend the EnKF-like approach to the BP setting by defining the necessary operations to handle general graphical models.
5

---

Published as a conference paper at ICLR 2025
Table Relations in Gaussian Ensemble Belief Propagation
Empirical statistics
Generative
Density-basee
Sample Condition
Operations
Propagate
Directed
Factor
X2
D2
Graph type
X3
Decomposition
P X2)
f(11, T2, 23)
Empirical moments m, K
Canonical parameters 72, P
Node Parameters
Ensemble recovery
GAUSSIAN ENSEMBLE BELIEF PROPAGATION
GEnBP consists of: (1) sampling from the generative prior; (2) converting ensemble statistics to canonical form; (3) performing BP with low-rank representations; and (4) recovering ensemble samples match updated beliefs (see Figure T. Although conceptually simple, the implementation details are delicate (see Appendix H).
We consider an approximation @M (mn,,, K) optimal if vn = m and the Frobenius norm |K - K|| is ceriiizzed-equivalllee if the Maximum Mean Discrepancy (MMD) with second-order polyno- mial kernel is minimized (Sriperumbudur et al. 2010 Example 3).
3.1 RANK-EFFICIENT GAUSSIAN PARAMETERISATIONS
A key component of GEnBP is the use of Diagonal Matrix with Low-rank perturbation (DLR) representations. We represent a symmetric positive-definite matrix K E DDD in DLR form as
K V +
with L R DxN, diagonal V, and s E (-1,+t +1} (see Appendix In the EnKF, the empirical covariance from ensemble samples is naturally DLR; here, N denotes the ensemble size (or rank of the DLR factors) and D the variable dimension.
To initialize BP, set each factor fj to the empirical distribution of the ensemble at that factor, Xy:
fj (M; Var,,rX,).
(18)
Among possible DLR approximations, the ensemble-based approach is most favorable; see Appendix for cost comparison.
Proposition 3. Suppose x ~ &M (a:; m, K) has a DLR covariance = LLT + V, with DxN and V diagonal. Then, the canonical form parameters n and P can be computed efficiently:
n Pm, P-K-t=U-RRi
where U is diagonal and R R DxN Both the conversion canonical form and the recovery of moments have time complexity O(N3 + D2
Proof: Standard application of the Woodbury identity, See Appendix F3

---

Published as a conference paper at ICLR 2025
3.2 BELIEF PROPAGATION WITH DLR REPRESENTATIONS
Standard GaBP requires expensive full-matrix operations for density multiplication. In contrast, GEnBP leverages the DLR structure for efficiency. For two DLR matrices,
K V+sLL'' K' +sLLL
their sum
K+ (V
L LLLL LTT
Thus, density multiplication (ef. Equation 14 can be performed efficiently while maintaining the DLR representation. Other operations such as vector products, rank reduction, and marginalisation also efficient (see Appendix
3.3 ENSEMBLE CONFORMATION
After message propagation, update ensemble samples to reflect the new beliefs (ensemble conformation) by finding an affine transformation the prior ensemble that best matches the posterior mean and covariance. Specifically, we seek affine transformation
Tp.T : X μB + XT,
with μ m (the posterior mean) and T E B N chosen to minimize
= arg min Var XT*) - (L.. V)
(19)
Since identifiable only up to a unitary transform, we have:
Proposition Any symmetric positive semi-definite matrix = T satisfying
(N (N-1 -
(V
defines an affine transformation that minimizes the loss in Equation
T9 Such can be computed
with memory complexity O(M² N2) and time complexity O(N3 DN2 + DM?).
Proof: See Appendix
H3
For a variable with K neighbors, in the worst case M KNN and the ensemble recovery cost is O3 + DN2K2)
3.4 COMPUTATIONAL COMPLEXITY
The overall cost of GEnBP depends the graph structure and ensemble size N. In many applications N < D, yielding significant savings compared to GaBP, which scales poorly with D. Table 2 summarizes the computational complexities of GaBP and GEnBP. Notably, GEnBP avoids the O(D)) costs of full covariance operations, making well-suited for high-dimensional problems.
GEnBP scales O(D) with D (versus O(D3) for GaBP) but O(K3) with node degree (versus (1))) In practice, large factors can be decomposed via Forney factorization (Forney) 20011 de Vries Friston 2017) reduce K without altering marginal distributions. If the R scaling is prohibitive, one may convert to full covariance matrix and use GaBP updates (see Appendix [[); however, in our experiments this was unnecessary.
EXPERIMENTS
In this section, we compare GEnBP against an alternative belief propagation method, GaBP, and, for reference, global Laplace approximation (Mackay) 1992). We use synthetic benchmarks designed to assess performance in bighh-iimnssonal,, nonlinear dynamical systems. In both cases, the graph structure is randomized system identification task (Appendix A. where static parameter influences a noisily observed, nonlinear dynamical system. Although GEnBP applies beyond system identification, this benchmark is chosen for its simplicity and popularity. More sophisticated graphical
7

---

Published as a conference paper at ICLR 2025
Table 2: Computational complexities for Gaussian Belief Propagation (GaBP) and Gaussian Ensemble Belief Propagation (GEnBP), where D is the node dimension, K the node degree, and N the ensemble size. Time complexities are measured in floating-polt multiplications.
Operation
GaBP
GEnBP
Time Complexity
Simulation
0(1)
O(N)
Error propagation
O(D*)
Jacobian calculation
O(D)
Covariance matrix
O(D2)
O(ND)
noddor-to-node message
O(D*)
O(K NO 3 K2)
Node-to-factor message
O(D*)
O(1) O(K + DN2K2)
Ensemble recovery
Canonical- -omments conversion
O(D3)
O(N3 N2D)
Space Complexity
Covariance matrix
O
O(ND) O(NDK)
Precision matrix
(
model problems discussed in Appendix C and Appendix B.3 demonstrates GEnBP applied to few-shot domain adaptation a system dynamics emulator on unobserved states.
For t ..., T we define
x = (x:
y, = Py(xx)
(20)
where the evidence is 8 {y: y! V and the query is 2
The GaBP implementation is modified from Ortiz et al. (2021) to allow arbitrary prior covariance. Hyperparameters not directly comparable between methods; we address this by choosing favorable values for each algorithm. We measure performance by the mean squared error (MSE) of the posterior mean estimate relative to the ground truth qo and by the log-likelihood log by- (go: q) go- In gradient-based methods such as Laplace or GaBP, the estimated posterior covariance may not be positive definite, rendering the log likeellloood undefined; in such cases, it is omitted. Additional experiments and details are provided in Appendix B
4.1 ID TRANSPORT MODEL
The transport problem (Appendix A.2) is a simple nonlinear dynamical system defined over dimensional spatial extent, chosen for ease of visualization. States are subject to both transport and diffusion, with the transport term introducing nonlinearity. Observations consist of subsampled state vectors perturbed by additive Gaussian noise. Figure [ shows sample prior and posterior distributions produced by GaBP and GEnBP, demonstrating that GEnBP recovers the posterior substantially better. Further analysis of this example is provided in Appendix A.2.T where it is also compared against global Laplace approximation.
4.2 2D FLUID DYNAMICS MODEL
In the computational fluid dynamics (CFD) problem (Appendix A.3), states governed by dis- cretized Havier-Stoker equations over a 2D spatial domain. The parameter of interest a static, latent forcing field q. In this setting, GaBP fails for D 1024, and the Laplace approximation is unavailable because the estimated posterior covariance is far from positive definite. In our 2D incompressible model, both the state and forcing fields are scalar fields over d d domain, which we stack into vectors so that D Dx, d². For comparison with classic sampler that does use belief propagation. we implement Langevin Monte Carlo (LMC) algorithm (Roberts & Tweedie, 1996) as a baseline, running it until it achieves MSE comparable to that of GEnBP (2000) iterations after 1000-iteration burn-in).
Figure shows results as the dimension increases. The Langevin Monte Carlo algorithm achieves the best posterior RMSE, but roughly 10s times the computational cost GEnBP. GEnBP attains
IIf q known, estimating (x.). from (y.) would be filtering problem, solvable by EnKE
8

---

Published as a conference paper at ICLR 2025
Execution Time (8)
Mean-Sope -eaan-Squared Error
Log Likelihaod
le+03
2
+
GFnBP
+
Langevin
le+01
+
-1
le+00
0
10a
10a
10
10a
10a
10
10a
10a
10*
b)
D
D
Figure 2: Influence of dimension D in the 2D fluid dynamics model. Error bars indicate empirical 50% intervals from n 10 runs. In grey-shaded regions, GaBP ran out of memory.
Execution Time (s) 4
Mean- -quared Error
Log Likelihood f
x x10
GaBP GEnBP
let02
Langevin
1c+01
1c-03
10-1
10
10-1
a)
Figure 3: Influence viscosity v on 32 x 32 2D fluid dynamics model. Error bars indicate empirical 50% intervals from 40 simulations.
superior MSE and posterior likelihood compared to GaBP while scaling to much higher Dæ. GaBP experiments are truncated because trials with o 1024 failed due resource exhaustion.
We caution that absolute runtimes should be interpreted carefully. In particular, the hime-dimensional Jacobian calculation GaBP is not effectively parallelized by PyTorch library, penalizing GaBP by a constant factor. The asymptotic ((D) scaling of these Jacobian calculations (Margossian 2019) will eventually favor GEnBP.
Figure illustrates the effect of varying viscosity y from laminar to turbulent regimes, producing diverse nonlinear behaviors (see Figure (). GEnBP strietly dominates GaBP in speed. generally outperforms in terms of MSE, although there are ranges around 1 where performance is similar or slightly worse. Regarding the posterior log- log-likelihood at the ground truth, GEnBP superior in the low-z (turbulent) regime, while GaBP performs better in the high-v (laminar) regime. Langevin Monte Carlo is not competitive in terms of posterior log- logllikelhoood, although its performance might be improved through hyperparameter tuning.
5 CONCLUSIONS
Gaussian Ensemble Belief Propagation (GEnBP) advances inference in probabilistic graphical models by combining the strengths of the Ensemble Kalman Filter (EnKF) and Gaussian Belief Propagation (GaBP). It scales to higher factor dimensions than GaBP and achieves superior accuracy in complex, high- nigh-dimensional problems. By employing ensemble approximations, GEnBP accommodates the larger and more intricate factors common in real-world applications. Unlike EnKF, GEnBP handles complex dependencies without requiring gradients and often outperforms GaBP both accuracy and speed.

---

Published as a conference paper at ICLR 2025
While GEnBP introduces additional tuning parameters and requires more model executions during sampling, these costs are often offset by avoiding the high-dimensional Jacobian calculations needed in GaBP. Its effectiveness is best in systems where the dynamics can be well-approximated in low-dimensional subspace (i.e., via DLR approximations). As with GaBP and EnKF, GEnBP relies on Gaussian approximations and is thus constrained to unimodal distributions; convergence analysis remains an open question. Moreover, comparing overall computational cost is complex and depends on the graph structure. Our empirical results demonstrate that there are problems where GEnBP is significantly more efficient than GaBP, although the full range of such problems has not been characterized.
Despite these limitations, GEnBP's scalability, flexibility, and case of use make it promising tool for applications such as geospatial prediction and high-dimensional data assimilation. Future work will focus on integrating practical improvements from both GaBP and EnKR-such covariance loc- alization, covariance inflation, adaptive ensemble selection, message damping, incremental updating, robust covariance scaling, and graph truncation- enhance stability and convergence at lower cost. Exploring connections with particle-based and score-based methods may further expand GEnBP's capabilities.
ACKNOWLEDGEMENTTS
We are indebted Laurence Davies, Rui Tong, Cheng Soon Ong and Hadi Afshar for helpful discussions and feedback this work.
REFERENCES
Elias Bareinboim and Judea Pearl. A General Algorithm for Deciding Transportability of Ex- perimental Results. Journal of Causal Inference, 1(1):107-134, May 2013. doi: 10.1515/ jci-2012-0004.
Danny Bickson. Gaussian Belief Propagation: Theory and Application. PhD thesis, July 2009.
Yuming Chen, Daniel Sanz-Alonso, and Rebecca Willett. Reduced-order autodifferentiable ensemble Kalman filters. Inverse Problems, 39(12):12400) October 2023. doi:
Andrew J. Davison and Joseph Ortiz. FutureMapping 2: Gaussian Belief Propagation for Spatial AI. arXiv:1910.14: 14939999.199911 [cs], October 2019.
Bert de Vries and Karl J. Friston. A Factor Graph Description of Deep Temporal Active Inference. Frontiers in Computational Neuroscience, 11, 2017.
Frank Dellaert and Michael Kaess. Factor Graphs for Robot Perception. Foundations and Trends@ in Robotics, 6(1-2):1-139 August 2017. doi: 10.1561/2300000043.
R. A. Dorfman. A note the 5-method for finding variance formulae. Biometric Bulletin, 1938.
Araad Doucet. A Note on Efficient Conditional Simulation of Gaussian Distributions. Technical report, University of British Columbia, 2010.
Ryan M. Eustice, Hanumant Singh, and John J. Leonard. Exactly Sparse Delayed-State Filters for View-Based SLAM. JEEE Transactions on Robotics, 2206):1100-1114 December 2006. doi: 10.1 10.1TOO 00.100990022000888666666 2006.886400.8
Geir Evensen. The Ensemble Kalman Filter: Theoretical formulation and practical implementation. Ocean Dynamics, 53(4):343-367 November 2003. doi: 10. 100077/10033-0000000----
Geir Evensen. Data Assimilation - The Ensemble Kalman Filter. Springer, Berlin; Heidelberg, 2009. ISBN 978-3-642-0371131 3-422-337113 doi:
Paul Fearnhead and Hans R. Künsch. Particle Filters and Data Assimilation. Annual Review of Stat-
isties and Its Application, 5(1):421-449 2018. doi: 10.
10

---

Published as a conference paper at ICLR 2025
Joel H. Ferziger, Milovan Perié, and Robert L. Street. Computational Methods for Fluid Dynamics. Springer, Cham, 4th ed. 2020 edition edition, August 2019. ISBN 978-3-319-9969 99691-2.
G.D. Forney. Codes on graphs: Normal realizations. IEEE Transactions on Information Theory, 47 (2):520-548, February 2001. doi: 1.1109918.91057.
Brendan J. Frey, Frank R. Kschischang, Hans-Andrea Loeliger, and Niclas Wiberg, Factor graphs and algorithms. In Proceedings the Annual Allerton Conference on Communication Control and Computing, volume 35, pp. 666-680 Citeseer, 1997.
R. Furrer and T. Bengtsson. Estimation of high-dimensional prior and posterior covariance matrices in Kalman filter variants. Journal of Multivariate Analysis, 98(2):27 February 2007. doi: 1.1016//.jmv..20066
Nathan Halko, Per-Gunnar Martinsson, and Joel A. Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions, December 2010.
L. Houtekamer and Fuqing Zhang. Review of the Ensemble Kalman Filter for Atmospheric Data Assimilation. Monthly Weather Review, 144(12):4489 4532, December 2016. doi: 10.1 MWR-D-15-0440.1.
Daphne Koller and Nir Friedman. Probabilistic Graphical Models Principles and Techniques. MIT Press, Cambridge, MA, 2009, ISBN 0-262-01319-3.
Frank R Kschischang, Brendan J Frey, and Ands-AAndea Locliger. Factor graphs and the sum- product algorithm. IEEE Transactions Information Theory, 47(2):498-519 February 2001. doi: 10.1
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier Neural Operator for Parametric Partial Differential Equations. arXiv:2010.08895 [cs, math), October 2020.
Finn Lindgren, Håvard Rue, and Johan Lindström. An explicit link between Gaussian fields and Gaussian Markov random fields: The stochastic partial differential equation approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 73(4):423 498, September 2011 1.
doi: 10. 1461-9/146779888.. 00177-9988..201.0000 2011.00777.x
David J. C. Mackay. A Practical Bayesian Framework for Backpropagation Networks. Neural Computation, 4(3):448-472, May 1992. doi: 10.1162/neco. (elp///cco..99..4... 1992.4.3.448.
Dan MacKinlay. The Ensemble Kalman Update is an Empirical Matheron Update, February 2025.
Charles C. Margossian. A Review of automatic differentiation and its efficient implementation. WIREs Data Mining and Knowledge Discovery, 944):e3305 July 2019. doi: 10. 10.100//WDMM.300
Kevin P. Murphy, Yair Weiss, and Michael I. Jordan. Loopy belief propagation for approximate inference: An empirical study. Proceedings of the Fifteenth Conference on Uncertainty in Artificial Intelligence, UAI'99, pp. 467-475, San Francisco, CA, USA, July 1999. Morgan Kaufmann Publishers Inc. ISBN 978-1-55860-6144
Joseph Ortiz, Talfan Evans, and Andrew J. Davison. A visual introduction to Gaussian Belief Propagation. 21X.:22107.02308 [cs]. July 2021.
Judea Pearl. Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference. The Morgan Kaufmann Series in Representation and Reasoning. Kaufmann, San Francisco, Calif, rev. 2. print, 12. [dr.] edition, 2008. ISBN 978- 978-1-55860-47977
Judea Pearl, Madelyn Glymour, and Nicholas Jewell. Causal Inference in Statistics: A Primer. Wiley, March 2016. ISBN 118-1-119-188841
Ananth Ranganathan, Michael Kaess, and Frank Dellaert. Loopy SAM. In Proceedings of the 20th International Joint Conference on Artifical Intelligence, IJCAI'07, pp. 2191-2196, San Francisco, CA, USA, January 2007. Morgan Kaufmann Publishers Inc.

---

Published as a conference paper at ICLR 2025
Gareth O. Roberts and Richard L. Tweedie. Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, 2(4):341-363 December 1996.
Håvard Rue, Sara Martino, and Nicolas Chopin. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 71(2):11-332, 2009. doi: 10.111 1467-9666-9988.
14,71181 200.9988820000.... 000000.
Bharath K. Sriperumbudur, Arthur Gretton, Kenji Fukumizu, Bernhard Schölkopf, and Gert R. G. Lanckriet. Hilbert Space Embeddings and Metries on Probability Measures. Journal of Machine Learning Research, 11:1517-1561, April 2010.
Aki Vehtari, Andrew Gelman, Tuomas Sivula, Pasi Jylänki, Dustin Tran, Swupnil Sahai, Paul Bloms- tedt, John P. Cunningham, David Schiminovich, and Christian Robert. Expectation propagation as way of life: A framework for Bayesian inference partitioned data. arXii:1112.448 [stat], November 2019.
Martin J. Wainwright and Michael I. Jordan. Graphical Models, Exponential Families, and Variational Inference, volume 1 of Foundations and Trends Machine Learning. Now Publishers, 2008. doi: 10.1561/2200000000
M. P. Wand. Fast Approximate Inference for Arbitrarily Large Semiparametric Regression Models via Message Passing. Journal of the American Statistical Association, 112(517):137-668, January
2017. doi: 10. 10.1080/01621455 (10000662149.20166 2016.. 1197833.
Yair Weiss and William T. Freeman. Correctness of Belief Propagation in Gaussian Graphical Models of Arbitrary Topology. Neural Computation, 13(10)22133-2200 October 2001. doi: 10.1 (0.116220899660077554177777
James Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Mare Peter Deisenroth. Pathwise Conditioning of Gaussian Processes. Journal of Machine Learning Research, 22(105):1-47, 2021.
Sewall Wright. The Method of Path Coefficients. The Annals of Mathematical Statistics, 5(3): 161-215, September 1934. doi:
Jonathan Yedidia, William T. Freeman, and Yair Weiss. Generalized belief propagation. In Proceedings of the 13th International Conference on Neural Information Processing Systems, NIPS'00, pp. 668 674, Cambridge, MA, USA, January 2000. MIT Press.
J.S. Yedidia, W.T. Freeman, and Weiss. Constructing Free-Energy Approximations and Generalized Belief Propagation Algorithms. IEEE Transactions on Information Theory, 51(7):2282-2311 July 2005. doi: 10. 10.1109/TT.22005555000
12

---

Published as a conference paper at ICLR 2025
A BENCHMARK PROBLEMS
To demonstrate our method's utility, we consider example problems of the system identification type, where latent parameter must be inferred from its influence the dynamics of hidden Markov model. In this context, the problem is closely related to but more challenging than than-state filtering.
Unless otherwise specified, the hyperparameters for both GaBP and GEnBP are set to 72 0.01 and の² 0.001. GEnBP has additional hyperparameters 17² = 0.1 and ensemble size N = 64. We cap the number pessage-propagatiio descent iterations 150 and relinearise or re-simulate after every 10 steps.
A.1 SYSTEM IDENTIFICATION
In the system identification problem, our goal is to estimate the fime-invariant latent system parameter q. This parameter influences all states in the hidden Markov model, where X represents the unobserved initial state and subsequent states X1, X evolve over time. Each state x; is associated with an observation y, as depicted in Figure
q
X
X1
Y1
Y2
Figure 4: Generative model for the system identification problem. The latent parameter q influences all states, and observed states are shaded.
A.2 ONE-DIMENSIONAL TRANSPORT PROBLEM
We consider a one-dimensional state-space model where the state at time t, x E Rd, evolves according
Pr Xt X1+1
Sz (C + 7)4) E.,
(1-))()))))))
(21)
Py Yr
De(xx) + n
(22)
where:
X is the state vector at time
q E Rd is a fixed parameter vector to be estimated.
Y E [0, 1] is a decay parameter.
C' denotes convolution of with kernel C' (modeling diffusion).
Sg is a circular shift operator shifting the vector by positions (modeling advection).
6 NN(0, (Ta)) is process noise.
13

---

Published as a conference paper at ICLR 2025
y, is the observation vector at time
Dr is a downsampling operator selecting every 8-th element.
77 -N o² Im is observation noise.
The latent parameter q is drawn from a prior that generates smooth periodic functions. Specifically, each element is defined
A exp x cos d -n)) k ==0,,....,dd d 1,
2
(23)
where
# Uniform(0. 2x] is random phase.
K ~ x² df smoothness d d-1/2) controls the concentration (i.e. smoothness) of the
function.
A is a scaling constant.
O Operator Definitions
Convolution (C' + z): Applies blur kernel C to the vector z, modeling diffusion.
Circular Shift (S²): Shifts z by k positions to the right, wrapping around.
Downsampling (Dy): Selects every &-th element from z, reducing spatial resolution.
Intuitive Interpretation The state update in Equation 21 models transport process where the
state x, evolves by:
1. Decaying toward background field: 7xx + (1 - 4.
2. Diffusing via convolution with C.
3. Being advected via the circular shift Sk.
4. Incorporating process noise Ct-
The observation equation Equation 22 represents measurements of the state at reduced resolution with additive noise 7n Our inference goal is estimate q from the sequence of observations {y,} and the initial state Xg. We use T = 10 timesteps for this problem.
A.2.1 RESULTS
pigtt the influence of the state dimension Dg 2 on various inference quality measures. GEnBP estimates substantially more accurate than GaBP, indicated by lower MSE and more stable posterior likelihoods. In this relatively low-dimensional problem, global Laplace approximation is also available. We observe that while the Laplace approximation is faster than both GaBP and GEnBP, its MSE performance is intermediate, and its posterior likelihood is less staecayyspeciallyy D ø increases beyond 100. In contrast, GEnBP's posterior log likelihood degrades more slowly and remains defined.
A.3 NAVIER-STOKES SYSTEM
The Navier- Stocke equation is classic problem fluid modeling (Ferziger al. 2019). Here we consider 2D incompressible flow solved with spectral method using PyTorch implementation from Li et al. (2020). Defining the vorticity w
&n and streamfunction & via
Br By
u Oy'
Oxb
dr
the Navier-Stokes equations are:
Poisson equation:
2 w.
2. Vorticity equation:
dw
0w
Pw
14

---

Published as a conference paper at ICLR 2025
Execution Time (s)
Mean-Squared Error
Log Likelihood f
-10
10a
+ GEaBP
GaRP
+
I Laplace
10a
-3
10
200
400
200
400
200
400
a.)
D
b)
D
22
Figure 5: Influence of dimension D in the transport example. Error bars are empirical 90% intervals from 40 runs. Missing log likelihood values indicate undefined values from non-posiinn non-positive definite covariance estimates.
X2
Y1
y2
ys
(a) 10
(b) 10
X2
(c) v/ 10a
Figure 6: Three Navier-Stokes simulations with At = 0.2 and varying viscosity The simulation is on a 64 x 64 grid; observations are x 8 and corrupted by white Gaussian noise with standard deviation 0.1 normalized scale.
15

---

Published as a conference paper at ICLR 2025
The equations are discretized onto d X d grid. At each timestep, inject additive Gaussian white noise V4 ~ (N,, (a and a static forcing term q to the velocity field. The state x obtained by stacking the 2D streamfunction & into vector. Observations are given by
yr (x)nsample,xx)) 77,
with 77 ~ 00, 1.
The prior state and forcing are sampled from a discrete periodic Gaussian random field with power spectral density
PSD(K) o² (4 + +
(24)
B EXPERIMENTAL DETAILS AND ADDITIONAL RESULTS
B. HARDWARE CONFIGURATION
Experiments were conducted on a Dell PowerEdge C6525 Server with AMD EPYC 7543 32-Core Processors rumning at 2.8 GHz (3.7 GHz turbo) with 256 MB cache. Float precision is set to 64 bits, and memory usage capped at 32 GB. Execution time limited to 119 minutes.
B.2 ENSEMBLE SIZE
In the main text, leave the choice of ensemble size N open. In practice, ensembles the order of N z 102 seem sufficient. The computational cost of GEnBP scales O(DN2 + N (see Section 3.4), but the extra precision gained by increasing is not unbounded. Figure shows an experiment where N is varied from 16 to 512 in steps of 16. The performance improves rapidly up to N = 64 with diminishing returns thereafter. We use N 64 as our default.
Mean-Squared Error
Exccution Time
x10 3
x10
Log Likelihood
10a
1.6
1.4
10a
10
100
200
100
200
100
200
Figure 7: Performance of GEnBP versus ensemble size N for the 1D transport system identification problem (Appendix A.2). Confidence intervals are empirical 95% intervals based on n = 80 runs.
B.3 EMULATION USING GENBP
In Section we demonstrated GEnBP in system identification problems, which are pedagogically useful but may not showcase the full range belief propagation methods. GEnBP can, however, be applied to more complex problems.
Here, we devise proof-of-concept experiment in few-shot domain adaptation emulator for unobserved states. The graphical model structure is shown in Figure I0
neural network
The interpretation is as follows: We have a forward model P X
1, u + Xr, a known observation operator H. X, + yı- and a neural network Q u trained to predict the state given the unknown parameter u. In the system identification problem, we stop here. However, in few-shot domain adaptation, introduce Bayesian neural network N 8 w trained on related (possibly synthetic) data. Our goal is to adapt the neural network to a new target problem by estimating the posterior distribution p(wly) using GEnBP.
with weights previously
16

---

Published as a conference paper at ICLR 2025
In our experiment, use the CFD model from Section as the ground truth. For the neural network, we employ the Fourier Neural Operator (FNO) architecture (Li al. 2020) because it can produce high-fidelity results with fewer than 10t parameters. We simulate 10 time steps on 32 x 32 grid with viscosity v 0.001 and A 0.1. GEnBP is applied to the first 5 time steps for domain adaptation, with observations downsampled by factor of 5 and corrupted with zero-mean Gaussian noise (o 0.01). An ensemble size of N 30 is used. The FNO employs 12 Fourier modes, 32 hidden channels, and 4 layers, with Tensor Fourier Neural Operator variant that incorporates 5% rank factor, totaling 1 1,590,030 weights. Figure 8 shows that the domain-adapp emulator closely tracks the target density over subsequent timesteps, while the unadapted emulator diverges.
Domain adapation of the latent Navier Stokes emulator
Adapted emulator
Original emulator
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
Figure 8: Performance of GEnBP in the emulator domain adaptation problem. The graph shows the Maximum Mean Discrepancy (MMD) between the observation-condin state distribution, the unadapted emulator, and the domain-adapter emulator, computed using isotropic Gaussian kernel of width 1/VD, over 50 randomized emulation problems.
C MINIMUM VIABLE INTRODUCTION TO PROBABILISTIC GRAPHICAL MODELS
The field of probabilistic graphical models is mature and extensive (see, e.g., Koller Friedman (2009); Wainwright & Jordan (2008)). ). For our purposes, it suffices to understand that this framework enables interence in complex hierarchical systems with partial observations, For example, Figure illustrates a simplified weather model where an oceanic physics simulator and an atmospheric physics simulator provide forward predictions and interact through physical coupling. Their outputs are observed by different sensors (thermal IR satellite imagery for the oceans and phased radar for the atmosphere). These systems are high-dimensional, noisy, and governed by nonlinear PDEs. Our primary target here is the one-step-ahead prediction for the atmospherie model. though extensions to jointly learn observation and physics parameters are straightforward.
Figure 10 shows model closer the system identification problem described in Section In this model, we extend the system identification task (recovering u) by also recovering an emulator (e.g., a neural network) that predicts approximate state X; from the previous emulated state x The goal is to estimate a posterior over neural network weights w to improve predictions when only noisy observations y, are available.
17

---

Published as a conference paper at ICLR 2025
Plased ratle observations
T4
Y1
Y2
Y3
Y5
Y6
y7
X5
Themaal IR
t1
t5
Figure complex graphical model structure from a weather prediction problem. includes an occanic physies simulator, atmospheric physics simulator, phased radar observations, and thermal IR satellite imagery. The model high-dimensional, noisy, and governed by nonlinear partial differential equations.
X1
X2
W
XD
X1
X2
Xt
u
y1
y2
y
Figure 10: Directed graph of the generative model Memu an emulation problem, showing temporal dependencies (l. = 1,2,.... 7). Observed variables y, shaded, and the global parameter w is highlighted.
Graphical model formalisms provide general framework for estimating states and parameters in such systems. The key idea is that if establish of inference rules for arbitrary models, can extend them to complex models like that in Figure
C.I DIRECTED GRAPII AND FACTOR GRAPIIS
Starting from the structural equations defining a model, by sampling from the ancestral variable distributions and iteratively applying the generative model,
X61
(x(xx))
(25)
x
(x(x)))
we obtain samples from the joint prior distribution.
For valid structural equation model, the following must hold:
18

---

Published as a conference paper at ICLR 2025
1. For all j. O n = 0 (i.e., self-loops).
2. If E j, then for k: j. v appears only in Kk.
This defines directed acyelic graph (DAG) over the variables, where nodes represent variables and edges represent generating equations. Each generating equation j induces density p(zo, (. so that
( II p( (.,.))
(26)
Under mild conditions (local Markov, faithfulness, and consisteney; see Koller & Friedman 2009 this representation is equivalent the joint distribution.
As running example, consider:
PiO
X1
Pi2(x1)
X2
Pia(xi)
X3
(27)
X4
P234 (X2, X3)
X5
(x5(x2)
X6
Pa5e (x4. xs)
This model is diagrammed in Figure Iaa and its density factorizes as
x (11)))22 )P13 (2 [z )233) (x, [22 23' 233 (ws (2g (xo (x1, 25).
(28)
The directed generative model M is shown in Figure 11a
The directed generative model is intuitive but not always convenient for inference, The factor graph representation (Frey et al. 1997 Kschischang et al. 2001) discards directional arrows and uses generic factors:
p x22 (x1 32).
(29)
Thus, the running example density Equation
28
becomes
p(z) (2))))
22a) 3234 (22, 03, (24) (25 25) (1:4, 4t,, me).
(30)
The corresponding factor graph is shown in Figure 1lb Factor graphs are bipartite, with nodes for both variables and factors, and they allow simple, local approximate belief-updating rules (introduced Section C.2).
We now introduce the conditional graph transformation. In practice, we are interested in the evidence- conditional posterior p Xg = x,g), i.e..
p =xx')
px =x's)=
(31)
plæs = =x²))
1p 22, 48 x'') dx
(32)
J p D5 = x²x)) dx dx
This conditional distribution can be represented as factor graph, as shown in Figure 11c. where the observed nodes are fixed and the corresponding factors are replaced by conditional factors (denoted fis and fise).
C.2 FACTOR GRAPH GAUSSIAN BELIEF PROPAGATION
Loopy belief propagation Murphy et al. inference in factor graphs. Although its convergence properties have been extensively studied (Wain wright & Jordan 2008; Yedidia et al.) 2005). we follow industrial practice (e.g. Dellaert & Kaess 2017) and treat sufficiently accurate for purposes.
1999 is an efficient algorithm for approximate marginal
19

---

Published as a conference paper at ICLR 2025
L
2
22
375
23
37G
(a) Prior directed graph of generative model M: p(æ) = Pi (xi (2 ( )p1 (2) P2s (ws) (2 P234 (20 (22, 40s (0s) / (Do a24, wt))
F1
f12
D22
f25
375
f1
f13
S234
f456
33
24
E6
(b) Prior factor graph 9 for M:
p(z) = fi f1 fis f as) f23 (22, mz, a. f25 ats) (a56 ms, 26).
L
1
31
f12
T2
f25
f1
f13
f231
S156
E3
24
(c) Posterior factor graph G* for M after observing Is and æe, where p(1,, 22, Da, æa (2s, æt) x (mi))faa(a (ia))i (x aas).fff f234 (ar2, æs, ar fās (ist)))sss
Figure 11: Graphical model variants for the example model in Figure 11 (a) prior generative graph, (b) prior factor graph, and (c) posterior factor graph; with query nodes 2 {1), latent nodes L (2, 3, 4), and observed nodes & (3.
Proposition 1 (Belief Propagation on Factor Graphs). By iteratively propagating the messages,
/ (5, (I.))) II mx-t dx
Mf
dar.
(5)
VEWg/V
Mx++11
mx
(6)
SELAe V
BP approximates the marginals as
my. xx plle) day
bg (xr))
1I SE
The belief propagation algorithm for generic factor graphs, using the messages defined in Equations
(Eq. 55. (Eq. and (Eq. [). summarized in Algorithm [ollowwig the terminology of Ortiz et al. (222
20

---

Published as a conference paper at ICLR 2025
Algorithm Loopy Low-rank Belief Propagation over Factor Graph G
Require: Factor graph 9 with variable nodes {xx} and factor nodes {ff}
Require: Initial messages (x,: + fÿ} and {/r - X1:
Ensure: Approximate marginal beliefs b(x²) for all X5: E E
1: Initialize message queues Qx--1, and (2 as empty.
2: while not converged do
3: for each factor f3 E 9 do
4: for each variable X E M2 do
5: Send factor-to-varee variable message (cf. (Eq. 5
6:
end for end for
7:
8:
for each variable X E G do
9:
for each factor E Kx do
10:
Send variablc-to-factor message (ef. (Eq, 6)
11:
end for
12: Update belief b(xr))) using (Eq. 7
end for
14: Check convergence criteria.
end while
16:
return b(x+) for all X E G.
C.3 GAUSSIAN FACTOR UPDATES
The factor graph formalism is especially convenient for Gaussian models because the Gaussian density is closed under multiplication, conditioning, and marginalisation. Although we introduced these ideas in Section 2.2.3 we detail them here.
Working with Gaussians in canonical form advantageous (Eustice et al. 2006). We denote by (x:; m. K) the Gaussian density with mean m and covariance K.
Proposition 5. Partition the random vector as
m& me
KE
Kit
Xj
ΦMs
(33)
ace
Kek
Kee
Define K == (KT++Ki, take the following forms.
Then, for Gaussian nodes and factors, the operations Definition
Conditioning:
dM KK);
2t + das xg: me KecKil (x mz), Kee KerKit Kae)
(34)
Marginalisation:
dn aoj mj, mj,Kis) + dx axxi max Kurk)
(35)
Multiplication (moments):
O'ar (z5; mly K3)) Φss m.j; + ans (1r; R (B::KKK,,, m, + K3J -mm;) R)
(36)
Multiplication (canonical): 6'c ( n P'r). de (o(p:,n,,P,)) + s0 (zjin' + n,, P/J + P3})
(37)
Proof: See Bickson 2009)
Note that Equation 37 is Equation 14 in the main text.
21

---

Published as a conference paper at ICLR 2025
C.4 GAUSSIAN APPROXIMATION OF NON-GAUSSIAN FACTOR POTENTIALS
For simplicity, assume that all factor potentials are bivariate with P(x1). (If not, the variates can be stacked and relabeled accordingly.)
A classic approach to approximating the factor potential generated by a nonlinear process is the propagation-of-erve also known as the 5-method (Dorfman, 1938). The 8-method uses first-order Taylor expansion to approximate how a nonlinear function transforms the moments (e.g.. the mean and covariance) a random variable. This is analogous to a Laplace approximation specifically used for error propagation in moments. In its simplest form, it approximates
Ex(x)) of [Ex1)
(38)
which is justified by Taylor expansion.
To estimate the joint covariance of the transformed variable, define
[(((tt) E
[
[xx) = ((xx) P(x)]
E
P(x))
(39)
The difference
Ejensen (x1. = Ex(())) ( - f (EX1)
(40)
is called the Jensen gap and represents one source of сπтог in the lincarisation.
Since the function P typically nonlinear and intractable for exact moment calculation, the ô-method further approximates ff using first-order Taylor expansion about the mean:
f(x) (ExI) Va' (x Ex1
(41)
The error in this Taylor approximation is defined
TTyloo (1 S. = Exxx) -
(42)
In summary, the overall error in approximating the transformed covariance using the 8-method includes the Jensen gap, Taylor error, any additional error due to non-Gaussiayity of the inputs. Although the d-method is simple and computationally efficient, it is not always optimal in terms of accuracy; alternative methods (such as higher-order approximations or sigma-point approaches, ensemble approximations) might sometimes yield better results.
D MINIMUM VIABLE INTRODUCTION TO THE ENSEMBLE KALMAN FILTER
The field of Ensemble Kalman methods mature and extensive; see, e.g., Evensen 2009) for comprehensive treatment. Here provide brief introduction sufficient for our purposes.
X
XI
X2
y1
Y2
Figure 12: Generative model for state filtering problem, where hidden states x, evolve over time and produce observations y, for t 1,2.....
22

---

Published as a conference paper at ICLR 2025
The central idea is as follows: Given the state-space model in Figure 12 where the joint prior density factors as
T
p( F1:T) = II ( p(ye
filtering seeks the conditional density p(Zo.rly1.T). By induction, if we know an estimate of 1a-o::r-1/9t-r-- then incorporating the new observation yr only requires multiplying by p(ap/er-1) (Vrlarr) and normalizing.
Under the assumptions of linearity and Gaussianity, let
= Φx (2r-1; Kr-1)-
Linearity implies that the state transition can be written as
Tr Pr-1(0r-1) Fr-1 + dr -1 + ET:
(43)
with cT ~ (0, Qr). Consequently, the joint density for ær-1 and oT (conditioned on y1.r-1) is
[
KT-1 Fr-IKr-
Pr-1Kr-1F7 1 + QT 1)
Kr-1FT-1
)
Qa (1
mT-1
FT
dr-1 1r-1mr-1
The marginal for æт is then
p 91.T-1) = @m (M mr, Kr),
(44)
with
Qr. Next, conditioning on the observation yr (assuming the linear observation model yr = HTæT - Cr + ET with gT ~ (O, RT) yields the standard Kalman filter update:
mr Fr--mr-1 + dr-1, K2 = Fr
Ky = Kr - Kr (HrRTH) + RT) Hr
(45)
my = m + Kr (HrKrH? + Rt) -1 ((T (yr-Hrru Hymm eT)
(46)
In practice, state transitions observation models may be nonlinear. The Extended Kalman Filter (EKF) approximates these models by linearizing via Taylor expansion (propagation of errors). For example, one sets
d
Fy-1 2
V P
and similarly for the observation model, while choosing diagonal covariance matrices Qr and Rr. However, the EnKF takes a different approach.
Rather than computing updates in closed form, the Ensemble Kalman Filter (EnKF) approximates the filtering distribution using an ensemble of Monte Carlo samples. Let
= [1x(1)
[N(N)
be a matrix of N samples drawn from p(æ; m, K). The empirical ensemble statistics are defined as follows:
1 X XAN, where Ax N
1
(47)
X. = X XBN, with By = 1IxN,
ÉX X,
(48)
VaryX = N XX +V,
(49)
Cov(X, N
T
(50)
23

---

Published as a conference paper at ICLR 2025
where V is diagonal "nugget" term (often set to I to ensure numerical stability.
The EnKF propagates the ensemble forward using the model:
Xr-1 XT
yITT-1 = Xr VI-T-1
(51)
P 11.7-1)
Thus, the EnKF approximates the joint density by the empirical distribution defined by the ensemble without ever explicitly evaluating the Gaussian density.
A key step is the ensemble update (often referred to as the Matheron update).
Proposition 2. Partition x [x] follows the Gaussian distribution:
x
such that X
[ xT]. Assume the ensemble X
Vary Xx Cov(Xx, Xc)
(ox((xxxx))
X~OM
(16)
Et
X
Vary (Xe)
In ensemble form, conditioning (Eq. 9 is performed as:
x + Xx + (Xv((x, Xx) Xx) Vary (x -Xa)
(17)
The computational cost of solving Equation 17 O(V3 + N Dxx). Marginalizaaion (Eq. 10 is
simply truncation, i.e., X + X3.
Proof: The equality Equation 17 follows from 6 with the substitution of Equation I6 which can be computed with a cost of ((N 0 Dx, using the Woodbury formula to efficiently solve the linear
system involving Vary (Xx) MacKinlay 2025). Marginalization is simply performed by truncating
the ensemble the relevant components.
We can use Equation 17 to calculate XTlæıT by using y/ as the observations ø in the update. We have thus obtained sampled from the filtered distribution without ever evaluating its density. This update then used to condition the ensemble on new observations. For example, if the joint ensemble of states and observations is
XT Yr
XT VI-T-1
H ((Xx(31x11
one uses Equation 17 with the observed yr in place of æद, to obtain the updated ensemble Xr .... In this way, samples from the filtered distribution are obtained without explicit density evaluations.
D.I GAUSSIAN APPROXIMATION OF NON-GAUSSIAN FACTOR POTENTIALS
Unlike the GaBP, which approximates non-Gaussian factors using the 8-method (i.e. via Taylor expansion and linearisation). the EnKF avoids this entirely by using empirical samples to estimate the moments of the joint distribution. Although this stochastic approximation introduces additional (aleatoric) noise, empirical studies (e.g. Furrer & Bengtsson (2007)) suggest that the EnKF often produces accurate estimates than the Gaussian approximations used in GaBP. In particular, by sampling from the full joint distribution, the EnKF bypasses the Jensen gap and Taylor errors inherent the 6-method.
Heuristically, this means that while the GaBP may produce an approximate joint variance of the form Var
GaRP
xi
J(m1) Var(x;) Var(x1) J J(m1) Var(x1) J (m1))
Var(x1)
(52)
X2
with J(m:) the Jacobian of P at the mean, the EnKF estimates the joint variance empirically as X(N)
FaKF
X1
X() ((((())
Var
Var
P(X(N))
(53)
X2
where the latter includes nugget term to maintain numerical stability. Although the relative quality of these approximations is complex to characterize theoretically, empirical evidence supports the superior performance of the EnKF approach in many high-dimensional problems.
24

---

Published as a conference paper at ICLR 2025
E MATHERON UPDATES FOR GAUSSIAN VARIATES
The pathwise update for Gaussian processes, known as the Matheron update (see, e.g., Wilson al. (202]); Doucet (2010), provides method to generate samples from a conditional distribution without directly computing its density. Suppose
Kyy Kwy
Kyw Kww
M
(54)
W
Mmw
Then the conditional moments are given by
E = w0 = my + Kwy Kwwl (w0" Kwo(e0" mw).
(55)
Varly w = w"1 Kyy Kwy K yyw.
(56)
Proposition 6. For
w as in Equation 54 the mapping
Y Kwy Kow (w" -w)
(57)
produces a random variable y" that has the same distribution as conditioned on w w=, i.e.
y" w ww==w.
Proof: Taking the expectation of Equation 57 gives
E(y"] = My + Kwy Kww Kw - -ww)
which matches Equation 55 For the variance, since wn is fixed and w is Gaussian, show that
Var 1y| = Kyy Kwy Kyw
which is Equation 56 Thus, the mapping in Equation [57 indeed produces samples with the desired
conditional moments.
Note that this update does not require explicit evaluation of the density of y and forms the basis for pathwise (sample-based) updates in Gaussian process models.
F DIAGONAL MATRICES WITH LOW-RAND PERTURBATION
Suppose
K= V + +LLL where Ls x N matrix, = diag(v) (with v a D-vector), and
E 1-1,1}
specifies the sign of the low-rank update. We are primarily interested in such matrices when D; in that case we call them Diagonal Matrix with Low-rank perturbation because their computational properties are favourable. (In contrast, matrices without such exploitable structure we refer as dense.)
Throughout this section we assume that the matrices in question are positive definite and that all dimensions are conformable. We refer to the matrix L the component of the DLR matrix, and diagonal matrix as the nuggel term (by analogy with classic kriging). Note that if 2 D the name is misleading since the matrix is not truly low-rank; the identities below still hold but the computational advantages are lost.
F.1 MULTIPLICATION BY DENSE MATRIX
The product of an arbitrary dense matrix X with DLR matrix can be computed efficiently by grouping operations. particular,
K K=V+sLL
then
KX = VX + (L
This requires only ODD dim(X)) operations for the diagonal term and O(DN dim(X)) for the low-rank term, and it can be computed without forming the full dense matrix K. (In general the product is not itself DLR matrix.)
25

---

Published as a conference paper at ICLR 2025
F.2 ADDITION
The sum of two DLR matrices with the same sign is also DLR matrix. Suppose
v= = +
Then
K+K's + V + LL +LL"))
(58)
= = + + [T. Ml LT
(59)
Thus the sum a DLR matrix with nugget + V and component L L'].
F.3 DLR INVERSES OF DLR MATRICES
Inverses of DLR matrices can themselves be written DLR form and computed efficiently.
Proposition 7. Let
K V + TL
be a DLR matrix with positive sign
1). Then its inverse is
-- RR
(60)
where
R R= V= chol (1 (1+L'V-L)""))
Here chol(A) denotes (conventional) Cholesky factor satisfying chol(A) chol(A)T A.
Proof: The Woodbury identity gives
-- V-=------((---V-1) '
Defining R as above yields Equation 60
Proposition 8. Let
P= U RR
be a DLR matrix with negative sign (s 11. Then its inverse is given by
- +IL
(61)
where
L= U-U-R choll (1o0((1- -I-RII-URUUT (I-RJU-IR)")
Proof: Applying the Woodbury identity for a negative low-rank update, we have
-U-RRT) -
R RR--RR R U
Defining as above, we obtain Equation 61
The computational cost of either inversion is O(N2D + N roughly O(N3) for constructing the Cholesky factor and (DD for the necessary matrix multiplications. The space cost is OND). For given K V + s L.I. the term
8 S II VV
conventionally referred the capacitance of the matrix.
26

---

Published as a conference paper at ICLR 2025
F.4 EXACT INVERSION WHEN THE COMPONENT HIGH- RANK
Suppose instead that the low-rank component is not low-rank at all (i.e. the D x N component satisfies N >> D). In this the cost of computing the low-rank mnvrsion-ddmmnnaeed by O(N)) Cholesky factorisation-exeee exceeds the cost of naively inverting the dense x matrix (which is O(D"). In such circumstances it can be preferable to compute the full dense inverse and then recover a DLR representation via spectral decomposition.
For example, if one computes p-1 (with as in Proposition 8). then since
p-1
p- = U LLI
one may obtain a low-rank factor by forming the spectral decomposition
p-1- =
and setting
L QA
(In the special case where is a constant diagonal, more efficient methods may be available.) A similar strategy applies for inverting DLR matrix of the form
== = V + LL
by computing its dense inverse K- and then recovering the low-rank representation from the identity
K-1= -RR
F.5 REDUCING COMPONENT RANK
We can use the singular value decomposition (SVD) to reduce the rank of the component in DLR matrix in way that minimizes the Frobenius norm error. Suppose
K V + I.I..
where is x matrix and V = diag(v) v). Let the thin SVD of L be
YSZ
with Y DXM S E R MXM (diagonal with nonnegative entries), and Z MXM, that = and 7 zT = IM. Then
LL = Y
Any zero singular values can be removed exactly. Moreover, if we retain only the largest singular values (setting the others to zero), we obtain the best Frobenius-norm approximation of LL of rank r. (A conventional thin SVD of L costs ODDM²) ); if M is large one may consider randomized methods.)
G GAUSSIANS WITH DLR PARAMETERS
We recall the common parameterisations of the multivariate Gaussian used in section Section 2.2.3 In moments form, we write
M=, K),
where E[x] and K = Var(x). In canonical form the same density is written
x ~ Nc P),
where
K m, P P=
(Up to normalizing constants the exponent in the density may be written either as - m) T K- K-1 m) or equivalently as -f Pæ n T x.)
A given Gaussian ensemble is naturally associated with a moments-form Gaussian,
X,2 NM(m K),
27

---

Published as a conference paper at ICLR 2025
where write
7 X]
K = (X) = - -
introducing parameter o² to ensure invertibility of the covariance if needed, more generally substituting an arbitrary diagonal matrix for I²I.
By applying Proposition 7 to the DLR form of the covariance, may compute the canonical parameters cheaply. That is, writing
K-1= Var,2 (X) = =P= 1
we have
P m.
Conversely, given n and P = - RR one may recover m and via
K = p- 0 2 RRT,
m=P--n= -
RR
using the alternative low-rank inverse formula (see Proposition 8
H
GENBP DETAILS
H.1 GENBP ALGORITHM
Here we detail the GEnBP algorithm, specifying the matrix operations needed to perform the Gaussian updates while maintaining the DLR forms for the matrix parameters.
Algorithm 2: GEnBP
Require: Graph 9, generative processes {Pr observations æs, and an ancestral ensemble
sample X
Ensure: An intivaation-oonitionnnl sample X
1: while not converged do
2: Sample an ensemble from the generative processes P, on G (using (Eq. 3
Convert & into the conditional graph g" by conditioning the observed factors (using
(Eq. () and (Eq. [7p.
Convert variables and factors to DLR canonical form. {Section 3.1
5: Propagate DLR messages on g". {Section H.2 Algorithm
Conform the ancestral nodes to the belief: update the ensemble sample via transforma-
tion T so that T(X)) bg- (x.r). {Section 3.3
end while
8: return Approximate posterior sample X
H.2 GENBP PACTOR-TO-VARIAN MESSAGE
We now outline the steps for GEnBP f; + factor-to- voriiiale message when there is a single incoming message X₁K + fy- The outgoing message Gaussian with canonical parameters (n', Pt)) where P/ U/- R is in DLR form. Extending the algorithm multiple incoming messages straightforward but verbose.
Algorithm 3: GEnBP /j + X₁ Message (Single Incoming)
Require: Factor fy with canonical parameters:
Information vector: n n;
T&
Precision matrix in DLR form:
U RRT. where
28

---

Published as a conference paper at ICLR 2025
diag
Ut [
Re R&
with Rs E R DxNN and Rş D
Require: Incoming message X + fy with canonical parameters:
Information vector: níg RR Pr,
Precision matrix in DLR form: PA U R T, where U/ diag(u))) and
R2g E IR XN'
Ensure: Outgoing message f - X with canonical parameters:
Information vector: nt/.
Precision matrix in DLR form: P = U R
Combine the information vectors for xg: n2 Tok: + n'y-
2: Combine the precision diagonals: üg f- uk + ukg
3: Concatenate the low-rank components: R -- RR
RA..
4: Form the joint information vector: 7111 <-
n2
[i
Uu
5: Form the joint precision diagonal: Un + diag
R. 0
6: Form the joint low-rank component: RuT <-
R
7: Store DLR form joint precision matrix Pn < Un RyBB.
Invert to obtain the DLR joint covariance: K₁ <- Pni Vn + Lu
9: Compute the joint mean: mi +- Kn V Vnnu nn
10: Extract the marginal mean for Xy: my <- (1
11: Extract the marginal covariance for Xy:
Diagonal: Ve Vn [1:,; 1:Dz).
Low-rank: L, < Lm [1:D,;: :].
Store DL.Rform K2 +- Vz + I.g.T
12: Reduce the rank of Kg to via SVD:
Compute SVD: L2 = ASB
Retain the top N components: 1/2 Ap.1in)S
Store DLR covariance: K, + Ve + L,L,
13: Compute the DLR marginal precision: P/P K, NVR-RR,
14: Compute the marginal information vector: n/ Pime U.me R&R m
15: return n/ and Py U R.R.T
H.3 ENSEMBLE RECOVERY
After belief propagation, suppose the belief for query variable is 6(x2) Nzr m, K with
DLR covariance Kg = V + وو To convert this belief into ensemble samples, we choose an affine transformation T that maps the prior ensemble X to updated ensemble X', T(X2) so that the empirical distribution the transformed ensemble approximates Nas (x m. K under metric d. (For brevity we suppress the subscript Q.)
We restrict T to the family of affine transformations
T T X + JB + XT.
29

---

Published as a conference paper at ICLR 2025
with parameters e R² and T E R RNxN, Minimizing the distance between the empirical moments and the target moments yields
= EX,
(62)
Var (XT) - -
-K|
T arg min T
(63)
arg min
Var, (T(X)) - (LL v)
(64)
Since the cost depends only on TTT, the minimizers are nonunique. We therefore optimize over M TT in the space of N N positive semidefinite matrices M
M arg min L(M),
(65)
MEM.
where
L(M) = N -
XMXT - -LLL + -v
(66)
with V = V - The gradient is given by
VML 2
x
XMXT
LLT V 7/x
N-
N
An unconstrained solution M' satisfies
XXMXXX= (N- (X + *TVX).
This system can be solved at O(N3) cost (using, e.g., pivoted LDL decomposition) and then projected onto the PSD cone. Finally, extract T via the SVD of M = USU by setting T 2 With careful ordering, the total cost is O(N + DN2 + DM²).
I ALTERNATIVE LINEARISATIONS AND LOW-RANK DECOMPOSITIONS
In response to reader's 's question, we investigate whether GEnBP merely low-rank decomposition of the GaBP algorithm. We argue that it is not; rather, the relationship between the two is as illustrated in Figure 13]
Although both GEnBP and GaBP use Gaussian approximations, they do not employ the same approximations. This is why, as shown Section GEnBP can surpass GaBP in accuracy as well as speed. We noted in Appendix C.3 that GaBP not only selects particular Gaussian density family (which it shares with GEnBP) but also uses specifically the 5-method for approximating nonlinear factors, opposed to the sample approximation of GEnBP. This apparently minor difference has significant implications, which we discuss below.
Consider the joint density of factor x2 P(x1) highlight these differences. For simplicity, suppose that the input and output dimensions are equal, i.e. D1 D2 = D. Two alternatives for covariance estimation have been discussed in this work.
GaBP uses propagation of errors:
Var X1
J(m1) Var X
X1 X2
Var
(67)
Var Xa T
J(m1) Varxi J JT (m1)
where (mm1) denotes the Jacobian of P evaluated mi. In contrast, GEnBP uses an empirical estimate:
(1)
((N)
Var
Xi
Var
(((((1)
P(2(N)))
(68)
30

---

Published as a conference paper at ICLR 2025
Gaussian Ensemble
Ensemble Kalman Filter
Belief propagation
Belief Propagation
(GEnBP)
Easembling
Ensembling
Linear Gaussian
Bellef propagation
Linear Gaussian
Kalman Filter
Belief Propagation
5 S-methad
S-merhed
Nonlinear Gaussian Belief
Extended Kalman Filter
Belief propagation
Propagation
(GaBP)
Figure 13: Relationship between GEnBP and GaBP.
Here the empirical covariance Var is computed using N samples from the prior, with an inflation term I ensure invertibility.
As noted in Section 3.4 the cost of computing the empirical covariance in Equation 68 is higher than that of the propagation-of main components to this cost:
estimate in Equation 67 both in time and space. There are two
Firstly, the cost of generating the ensemble (in GEnBP) versus computing the Jacobians (in GaBP). In GEnBP, generating N samples from the prior incurs cost that scales O(N). In contrast, computing the Jacobian in GaBP depends on the function but is generally O(D).
Secondly, the cost of performing the necessary matrix multiplications. The matrix product in Equation involves three 7) matrices and hence costs O(D3). In Equation 68, the empirical covariance is computed at cost of O(DN). (It would be O(D2N) if we were to calculate the full covariance matrix; however, as argued in this paper, it suffices to compute the deviance matrices.)
For simplicity, suppose we are in regime where N < D and the node degree K is small. In this setting, the belief propagation steps in GEnBP are generally cheaper. One might then ask whether GaBP could also benefit from low-rank belief updates in the SAEM setting.
Suppose we wished to construct a third option low-rank GaPB (LRGaBP)-whnn uses a low-rank decomposition of the covariance matrix to perform approximate GaBP in the hope of achieving efficiency similar to GEnBP. First, one would compute a rank N decomposition of the prior covariance, i.e
Varx LLT L LE DxN
perhaps by eigendecomposition. Naively, this operation is D3 O(D² log + + N3) using the method of (Halko et 2010, 6.1).
31

---

Published as a conference paper at ICLR 2025
The joint variance arising from the propagation-of-erios (or d-method (Dorfman 1938) can then be written as
Varx1
J(m1) Varx:
Xi
Var
~
X2
Varxi J(m1)
J(m1) Var x1 J (m1))
LLT
I( J(m1)Li J [mm))
LL J(m1)
J(mi) I. [(m.)
(69)
This is indeed low-rank decomposition and is amenable to the same techniques the other low-rank methods outlined in Appendix (F) However, is not computationally competitive. Since J ((m1) is a D x D matrix, the product J (m, L costs O(D²N), in addition to the (D) cost of computing the Jacobian.
We summarise the costs of this hypothetical LRGaBP step in Table
Notably, although the precise relationship between these algorithms depends on the problem structure, generally expect GEnBP to be more efficient than LRGaBP for digh-dimensional problems at fixed N; in many cases GEnBP exhibits a lower exponent D. Furthermore, since LRGaBP is an approximation to GaBP, its accuracy inherently bounded by that of GaBP, whereas being--eing a different approximation to the target estimand not subject to this limitation.
Nonetheless, the LRGaBP algorithm has not been pursued further given its computational inefficiency and limited accuracy gains relative to GEnBP.
Table 3: Computational Costs for Gaussian Belief Propagation (GaBP), Ensemble Belief Propagation (GEnBP), and the hypothetical Low-Rank Gaussian Belief Propagation (LRGaBP) for node dimension D and ensemble size/component rank N.
Operation Time
GaBP
GEnBP
LRGaBP
Simulation
O(1) O(D*) ((D)
O(N)
O(1)
Error propagation
O(D'N)
Jacobian calculation
O(D)
Covariance SVD
ONN + N
O(D*log N + +N3 + N
Space
Covariance Matrix
O(D2)
O(ND) O(ND)
O(ND + D O(ND)
Precision Matrix
O(D2)
32


Reviews:
[{"review_id": "X7gpU2QHHR", "reviewer": "Reviewer_dsRA", "summary": "The authors introduce a novel method named Gaussian Ensemble Belief Propagation (GEnBP). They combine EnKF and GaBP together, in order to improve the efficiency in high-dimensional and non-linear cases. The method uses low-rank message-passing algorithms.", "questions": "Please see weaknesses.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The authors try to use the strength of GaBP to solve (or I would like to say mitigate) the bad performance of EnKP in high-dim and non-linear problems. In the given examples, the new method works well. The computational complexity reduces largely for the dimensions of the data. It seems that the paper demonstrates superiority over existing methods in various benchmarks (but the figures are hard to read). I believe this method is potential for real-world applications for example in geospatial modeling.", "weaknesses": "1. Most important question for me: the scalability of the method to extremely large models, like weather simulations, remains to be fully demonstrated. Computational complexity still depends heavily on the degree of nodes in the graphical model, which means you must restrict the scale of the model. However, models in climate research always be extremely huge. How to choose N, how to choose K? If too small, can the model catches the nonlinearity well?\n2. EnKF is somehow popular in climate or geophysical models, because of its computational complexity. But GaBP, in my impression, is a little bit outdated. For the nonlinear cases, we may use variational inference, or particle filter. Have you compared with them?\n3. The innovation: GaBP mostly for low-dim problem. I am not sure it can help to improve the EnKF itself. More importantly, please see Bickson et al (2008) for the similar idea.\n4. Even for small model, for example, Laplace seems better after consideration of trade-off. And Fig 2(c) seems not complete.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "7eZp5Yuqg5", "reviewer": "Reviewer_HdeM", "summary": "The paper introduces the GEnBP algorithm, a method that combines Ensemble Kalman Filter (EnKF) and Gaussian Belief Propagation (GaBP) to efficiently perform inference in high-dimensional systems. GEnBP updates ensembles of prior samples into posterior samples using a message-passing strategy on a graphical model, making it suitable for applications with complex dependencies and high-dimensional variables, such as in geospatial and physical models.", "questions": "1. Figure 2 is not cited in the text, which may confuse readers\n2. Consider put Figure 1 in the experimental section\n3. The paper highlights potential applications in fields like geospatial prediction and physical model inversion, which may involve real-time data. How does GEnBP handle streaming data, and what would be required to extend it to real-time applications with time-varying dependencies?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well-structured, making it easy to follow the ideas presented. The proposed method can be considered to be novel in its approach to addressing the computational challenges of GaBP in high-dimensional systems.", "weaknesses": "While GEnBP demonstrates strong performance in high-dimensional contexts, its complexity scales disproportionately with the node degree in graphical models. This presents a significant limitation in highly connected graphs, where methods such as Forney factorization have been proposed by the authors as partial solutions. However, the effectiveness of these methods remains uncertain. Can some simple numerical results or insightful analyses be added to showcase this potential solution?\n\nAdditionally, a notable drawback is the absence of comparative analysis with other related works, both those employing and not employing EnKF, for efficient inference in high-dimensional, black-box systems. Key related works include:\n\n  - Chen Y, Sanz-Alonso D, Willett R. \"Autodifferentiable Ensemble Kalman Filters,\" *SIAM Journal on Mathematics of Data Science*, 2022; 4(2):801-833.\n\n  - Chen Y, Sanz-Alonso D, Willett R. \"Reduced-order Autodifferentiable Ensemble Kalman Filters,\" *Inverse Problems*, 2023; 39(12):124001.\n\n  - Lin Z, Sun Y, Yin F, Thiéry A. \"Ensemble Kalman Filtering-Aided Variational Inference for Gaussian Process State-Space Models.\" *arXiv preprint arXiv:2312.05910*, 2023.\n\n  - Girin L, Leglaive S, Bie X, Diard J, Hueber T, Alameda-Pineda X. \"Dynamical variational autoencoders: A comprehensive review.\" *arXiv preprint arXiv:2008.12595*, 2020.\n\nIncluding a discussion or comparison on system identification and data assimilation from these (and other existing) works would enhance the comprehensiveness of the proposed method's positioning.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "pjmX2lzl3h", "reviewer": "Reviewer_Mt1V", "summary": "The paper presents a generalisation of the ensemble Kalman filter to make inference on arbitrary graphical models using belief propagation  (BP) with ensembles. By relying on ensembles rather than keeping track of the full moments, the method is able to scale to problems where the state dimension of the nodes are high dimensional, e.g. weather forecasting. The authors propose several numerical tricks that make this possible while retaining scalability. This is tested against a standard Gaussian BP on some examples, showcasing its computational efficiency and performance improvements, offered by its improved handling of nonlinear conditional dependence between nodes.", "questions": "- It would be great if the authors could please clarify why we need a particle representation of the beliefs when we can just make computations with the DLR representation? My understanding is that in EnKF, it is used to propagate the particles by the nonlinear dynamics model and update them using Matheron's rule, but what does this correspond to in the general BP framework (in terms of steps (5)-(7))?\n- Do we need tricks like covariance localisation and inflation in this general BP setting? These are almost always used to make EnKF work robustly in high dimensions.\n- In Figure 2, why are some results missing for GaBP and Laplace method in the log-likelihood comparison? Was it producing NaNs due to overconfident predictions?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The paper is generally well-written and well-motivated. The methodology presented, while seemingly a small step beyond standard BP, is not at all trivial to execute, requiring various tricks to make it work in practice. In particular, I find the efficient computation of the factor-to-variable message and the ensemble conformation step to be an interesting contribution that is essential to retain the ensemble-based representation of the beliefs. The experiments showcase interesting settings where this technique could be applied, such as latent forcing identification in fluid models. In addition, it shows clear benefits of the approach over vanilla BP, with orders of magnitude increase in computation speed and performance gains.", "weaknesses": "While the paper is clearly written for the most part, there are some parts that I find need more explanation. See my questions below for details. I also find that the experiments can be improved with further baseline comparisons. Is GaBP really the only baseline that can solve the system identification problems in the experiments? One may probably also consider methods like the integrated nested Laplace approximation (INLA), or particle MC methods (e.g. the particle marginal Metropolis-Hastings algorithm).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "tU1WcFHr9D", "reviewer": "Reviewer_QpY9", "summary": "Efficient inference is challenging in graphical models with high-dimensional variables. \nThe authors introduce an inference algorithm based on Gaussian Ensemble Belief Propagation (GEnBP) and the Ensemble Kalman Filter (EnKF), combining strengths of both algorithms to arrive at an algorithm which is competitive in accuracy, in scalability with respect to the dimensionality of the variables, and which can be applied to a wide range of inference problem structures.", "questions": "Regarding Section 4.1:\n- Figure 2 is never referenced. Although I assume the first paragraph describes observations made from it? \n- Why is the log likelihood only shown for GaBP in the case of the lowest dimensionality setting used? \n- I struggle to reach the following interpretation based on Figure 2: “We see that the Laplace approximation … but its posterior likelihood, while similar to GaBP, is even less stable, and both are inferior to GEnBP”:\n- - For dimensionalities greater than around 100, it seems GEnBP is underperforming Laplace in terms of log likelihood (while the results provides no comparison to GaBP)? \n- - Why is GaBP and Laplace interpreted as less stable than GEnBP based on this figure? \n\nSection 4.2:\n- Figure 3 is from the text implied to contain results using Langevin MC, but it does not --- only Figure 4 does. \n- The results show (and the text says) that GaBP performs better in terms of log likelihood for high viscosity settings. What do you suspect is the explanation for this? \n\nOverall, I think the paper is good and makes valuable contributions. However, I think the manuscript needs an overhaul to fix the language errors, missing punctuation, and figure references before it is ready for publication. If this is done and the experiments section is clarified I intend to raise my evaluation score.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 4, "strengths": "The authors present a novel methodology with well motivated approximations and algorithmic structure, where important connections and details are highlighted. The methodology contains contributions and explanations that are valuable also on their own. Graphical illustrations and notation is used in an effective manner to improve clarity, and the text is well structured.", "weaknesses": "Although I find the text to be well structured in terms of its content, the text contains many language errors, forgotten punctuation, and incorrect figure references. I have listed some errors I encountered below which the authors may find helpful for improvement of the manuscript.  \n\nList errors/suggestions style and language:\n- Line 22: no spaces around “–”. \n- Line 48: why “seems” novel and not “is” novel, if that is the case? \n- Line 68: “table” is not capitalised.  This is also the case in many other places. Note that it should be capitalised since “Table 1” is a name.   \n- Line 71: “section” not capitalised.  This is also the case in many other places. \n- Line 72: missing punctuation after “approximations”.\n- Line 81-95: In my opinion, excessive use of bolds. Is it needed more than say a couple of times here, or even at all?\n- Line 83: How about “a novel message-passing method for inference in high-dimensional graphical models”?\n- Line 91: “Rank” is (incorrectly) capitalised even though the list is written as a sentence. \n- Line 93: Same as above but for “Compute”. \n- Line 99: How about “We will now introduce our notation and essential concepts”?\n- Line 100: “appendix” is not capitalized. This is also the case in many other places. \n- Line 118: How about “We will assume that queries are always ancestral variables, i.e. that …”?\n- Line 134: “equation” is not capitalized.  This is also the case in many other places. \n- Line 146: Missing word 'with'. “Edges connect each factor node j with each variable node…”? \n- Line 173: Unnecessary abbreviations, “we’ve” and “we’re”.\n- Line 180: observed is misspelt.\n- Line 184: “figure” not capitalised.  This is also the case in many other places. \n- Line 191: Missing punctuation before “The following”.\n- Line 202: “definition” not capitalised. This is also the case in many other places. \n- Line 227: Double punctuation. \n- Line 261: Missing punctuation before “Marginalization”. \n- Line 294: Missing punctuation before “Throughout”. \n- Line 312: Missing punctuation before “Here”. \n- Line 314: Misspelt equivalently. \n- Line 323: Incorrect grammar, “See appendix I for a comparison with a naive attempt to do without the ensemble”.\n- Line 325: Missing punctuation before “Then”. \n- Line 336: Missing punctuation after “high dimensions”. \n- Line 353: “algorithm” not capitalised. This is also the case in many other places. \n- Line 376: Incorrect grammar: “For a variable with K neighbors, in the worst case M = KN , when the cost becomes…” and missing punctuation. \n- Line 381: How about “In many practical applications it holds that N ≪ D, resulting in significant computational savings compared to GaBP which scales poorly with D.”. \n- Line 385: How about “better suited”?\n- Line 404: Incorrect grammar “While GEnBP scales favorably with D, but unfavourably with respect to the node degree K, scaling as O(K 3 ) for some operations.”\n\n\nThe presentation of the experiments is unclear at times, please see Questions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)