Paper deep dive
Information Geometry of Message Passing
Mykola Lukashchuk, Kyrylo Yemets, Alex Ledbetter, İsmail Şenöz
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 8/22/2026, 3:18:25 AM
Summary
The paper introduces Natural-Gradient Message Passing (NGMP), a variational inference algorithm that localizes the natural-gradient stationary condition of variational inference onto Forney-style factor graphs. NGMP projects exact belief-propagation log-messages onto the tangent space of the receiving edge's exponential family, preserving representable information better than Variational Message Passing (VMP). Experiments on Poisson smoothing, heteroskedastic regression, and ETTh forecasting demonstrate that NGMP provides superior uncertainty calibration, particularly when edge uncertainty persists.
Entities (10)
Relation Signals (5)
Natural-Gradient Message Passing → derivedfrom → Bethe Free Energy
confidence 95% · We start from the Bethe free energy and constrain a selected edge marginal to an exponential family.
Natural-Gradient Message Passing → operateson → Forney-style Factor Graph
confidence 95% · We show that the natural-gradient stationary condition of variational inference has an edge-local form on a Forney-style factor graph.
Natural-Gradient Message Passing → uses → Natural Gradient
confidence 95% · Each projected message is the natural-gradient projection of the exact belief-propagation log-message
Natural-Gradient Message Passing → comparesfavorablyto → Variational Message Passing
confidence 90% · Compared with variational message passing, NGMP keeps the part of the exact message that the receiving family can represent instead of averaging the factor under the neighboring beliefs.
Natural-Gradient Message Passing → improves → Uncertainty Calibration
confidence 90% · Experiments on Poisson smoothing, heteroskedastic regression, and hourly ETTh forecasting confirm this and show that the gain appears mainly in uncertainty calibration.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We show that the natural-gradient stationary condition of variational inference has an edge-local form on a Forney-style factor graph. We start from the Bethe free energy and constrain a selected edge marginal to an exponential family. At a stationary point, the natural parameter of that edge equals the sum of two projected messages, one from each incident factor. Each projected message is the natural-gradient projection of the exact belief-propagation log-message at the current receiving marginal, or equivalently, the gradient of its expectation in the so-called mean coordinates. We call the resulting scheme natural-gradient message passing (NGMP). The rule is local; each edge may carry its own exponential family, and the message a factor sends depends on the marginal that receives it. Compared with variational message passing, NGMP keeps the part of the exact message that the receiving family can represent instead of averaging the factor under the neighboring beliefs. The two coincide when the uncertainty on the edges entering a non-conjugate factor vanishes, and NGMP is more accurate when that uncertainty persists, for example, along a partially observed latent chain or when parameters are filtered through successive data batches. Experiments on Poisson smoothing, heteroskedastic regression, and hourly ETTh forecasting confirm this and show that the gain appears mainly in uncertainty calibration.
Tags
Links
- Source: https://arxiv.org/abs/2608.15922v1
- Canonical: https://arxiv.org/abs/2608.15922v1
Trouble viewing inline? Open PDF directly →
Full Text
176,043 characters extracted from source content.
Expand or collapse full text
Information Geometry of Message Passing Mykola Lukashchuk m.lukashchuk@tue.nl Affiliation: Eindhoven University of Technology Kyrylo Yemets Affiliation: Lviv Polytechnic National University Alex Ledbetter Affiliation: Eindhoven University of Technology İsmail Şenöz Affiliation: Lazy Dynamics, Utrecht, the Netherlands Abstract We show that the natural-gradient stationary condition of variational inference has an edge-local form on a Forney-style factor graph. We start from the Bethe free energy and constrain a selected edge marginal to an exponential family. At a stationary point, the natural parameter of that edge equals the sum of two projected messages, one from each incident factor. Each projected message is the natural-gradient projection of the exact belief-propagation log-message at the current receiving marginal, or equivalently, the gradient of its expectation in the so-called mean coordinates. We call the resulting scheme natural-gradient message passing (NGMP). The rule is local; each edge may carry its own exponential family, and the message a factor sends depends on the marginal that receives it. Compared with variational message passing, NGMP keeps the part of the exact message that the receiving family can represent instead of averaging the factor under the neighboring beliefs. The two coincide when the uncertainty on the edges entering a non-conjugate factor vanishes, and NGMP is more accurate when that uncertainty persists, for example, along a partially observed latent chain or when parameters are filtered through successive data batches. Experiments on Poisson smoothing, heteroskedastic regression, and hourly ETTh forecasting confirm this and show that the gain appears mainly in uncertainty calibration. 1 Introduction Bayesian inference provides a coherent framework for reasoning under uncertainty, justified by the axioms of consistent belief (18). To apply it, we need a model that specifies a joint density p(,)p( z, y) over observations y and latent quantities z, and Bayes’ rule turns this joint density into a posterior distribution over the unknowns after the observations are revealed p(∣)=p(|)p()∫p(,).p( z y)\;=\; p( y| z)p( z) p( y, z)\,d z. (1) This posterior is the object one would like to use for prediction, filtering, decision making, and model comparison. The only difficulty is evaluating the normalizing constant ∫p(,) p( y, z)\,d z. The normalizing constant is often a high-dimensional integral for which closed-form solutions exist only for special model classes, such as conjugate exponential-family models, linear-Gaussian state-space models, or graphical models with sufficiently simple structures. When a closed-form solution does exist, the obstruction is usually computational: exact probabilistic inference in Bayesian networks is NP-hard (17), and even guaranteed approximate inference is NP-hard in the worst case (19). This is the basic reason why practical Bayesian inference is usually approximate. Variational inference is one of the standard ways to address this, replacing integration with optimization. Instead of trying to integrate the posterior exactly, it chooses a tractable family of candidate distributions and solves an optimization problem inside that family (26; 53; 11). This exchange of integration for optimization is not free: the answer is now limited by the family and the constraints we impose. There are two complementary ways to do Variational Inference. The first is global: choose one approximating distribution qλ()q_λ( z) for the entire latent state and optimize its parameter λ. When qλq_λ is an exponential family, the Fisher geometry gives a clean natural-gradient characterization of variational stationarity; this is the viewpoint developed by 29. The price is that the approximation is specified as one global object. For a large graphical model, the sufficient statistics, Fisher matrix, or required expectations of such a global family may be too large or poorly matched to the factorization that made the model useful in the first place. The second route is local: it starts from a factorization that enables the conversion from a global problem into one of local computations towards the same objective11 1 We drop the relation on y on the right side of the equation because y can always be absorbed into the factors faf_a. p(∣)∝∏a∈fa(a).p( z y)\; \; _a f_a( z_a). (2) Each factor faf_a involves only a small subset of variables a z_a. In the factor-graph representation of a probabilistic model formalized in Section 2.1, these variables are represented by edges, and a belief propagation (BP) message is a local density function passed across such an edge. On trees, products of these messages compute exact edge marginals, while on loopy graphs, the same local view leads to Bethe and structured variational approximations (58). More generally, variational objectives can often be decomposed into local pieces, and their stationary conditions can be interpreted as message-passing algorithms (59; 48). The local view is attractive because it preserves the modularity of the model. A factor contributes a local computation, an edge carries a local belief (message), and different modeling components can be combined without rebuilding a single, monolithic posterior approximation. Locality alone, however, does not make the messages finitely parameterized objects. An exact BP factor-to-edge update through a non-conjugate factor may return a function that is not Gaussian, Gamma, or any other finite-dimensional message supported by the surrounding graph. Variational message passing (VMP) can restore closure by changing local constraints, but this also changes what information enters the outgoing message. This paper asks how the global and local views meet. Can the global natural-gradient stationary condition of variational inference be localized into a rule for messages on a factor graph? The answer developed below is yes. We start from the Bethe free-energy formulation of message passing and impose an exponential-family form constraint on selected edge marginals. fbf_bfcf_c⋮ μl→b _l→ bμk→b _k→ bμb→i _b→ iziz_iBelief Propagation μb→i(zi)=∫fb(b)∏j≠iμj→bdb∖i _b→ i(z_i)\;=\; \!f_b( z_b)\! _j≠ i\! _j→ b\,d z_b i only incoming edges used f~b f_bfcf_c⋮ μl→b _l→ bqkq_kμ~b→i μ_b→ iziz_i(Structured) Variational MP μ~b→i(zi)∝∫f~bμl→bdzl μ_b→ i(z_i)\; \; \! f_b\, _l→ b\,dz_l f~b=exp(qk[logfb]) f_b\;=\; \! (E_q_k[ f_b] ) in-cluster: messages; others: marginals fbf_bfcf_c⋮ μl→b _l→ bμk→b _k→ bmb→im_b→ iziz_i⋮ μm→c _m→ cμn→c _n→ cExpectation Propagation mb→i(zi)=Πℰi[μb→iμc→i]μc→i(zi)m_b→ i(z_i)\;=\; _E_i\! [\, _b→ i\, _c→ i\, ] _c→ i(z_i) uses the message μc→i _c→ i from fcf_c fbf_bfcf_c⋮ μl→b _l→ bμk→b _k→ bqλiq_ 0.75$ _i$μ^b→i μ_b→ iziz_iNatural-Gradient MP (this work) ℓb→i=logμb→i _b→ i= _b→ i ηb→i=iλi[ℓb→i] _b→ i= T_i _i[ _b→ i] μ^b→i(zi)∝expηb→i⊤Ti(zi) μ_b→ i(z_i) \ _b→ i T_i(z_i)\ uses receiving qλiq_ _i Figure 1: Messages produced by belief propagation (BP) rely on the local factor fbf_b and its other incoming messages. Variational message passing (VMP), including its structured variant, first forms a surrogate factor f~b f_b by variational expectation of fbf_b with respect to the current marginals of the out-of-cluster variables, then passes a BP message through this new factor; non-conjugate and projective VMP keep this flow and additionally project the resulting message at the receiving marginal (Section 6). Expectation propagation (EP) instead forms a cavity distribution from the opposing side, restores the exact factor contribution to obtain a tilted marginal, projects that marginal onto the chosen form-constraint family via Πℰi[⋅] _E_i[·], and divides out the cavity to obtain the new site (41). Finally, natural-gradient message passing uses the receiving marginal to define the tangent space for natural-gradient projection of the exact BP log-message. The final message towards ziz_i is obtained as the exponential of the resulting projected natural parameters ηb→i _b→ i applied to the chosen family’s sufficient statistics Ti(zi)T_i(z_i). The idea can be read informally before the formal derivation. Exact BP messages describe the unconstrained local update that a factor would send to an edge if the edge marginal could take an arbitrary functional form. In general, this update does not lie within the finite-dimensional tangent space to the manifold defined by the chosen form constraint (see Figure 3 for a visual intuition). Understanding variational inference as an iterative process moving variable marginals towards stationary points of an approximate Bayesian objective, we therefore project the exact BP update onto the tangent space of the constrained marginal family at the current marginal belief. This yields a natural-gradient message that encodes the component of the BP update that can be represented by the chosen family’s sufficient statistics. The resulting message moves the edge marginal within the constrained family, preserving the locally relevant part of the full BP update while discarding directions that the family cannot encode. This process continues until the BP message encodes only information that the chosen family cannot represent, the difference between the marginal natural parameters and the tangent component is zero, and thus a stationary point has been reached. Figure 1 describes and contrasts different local message-computation strategies based on the information sources they draw from and the message-approximations they make in comparison to BP. Belief propagation is included as a basis, but is understood to not produce conjugate messages in general. What makes natural-gradient message passing (NGMP) unique is its optimal use of the exact BP message information visible to the form-constraint family of the receiving edge. Experiments in Section 7 locate where this matters: whenever the uncertainty on the edges entering a non-conjugate factor cannot be reduced, the projected messages estimate marginals more accurately than VMP and its projective variants at comparable complexity. Hourly ETTh forecasting supports the corresponding predictive consequence: neural gates can preserve competitive point error while their implied predictive scale collapses, whereas a precision-gated NGMP ensemble remains substantially better calibrated. The paper is organized to support different entry points. Sections 2 and 3 reviews Forney-style factor graphs, variational inference, and classical variational message passing as constrained Bethe-free-energy optimization; readers familiar with variational inference and factor graphs may begin with Section 4, which contains our main theoretical contribution. Section 5 gives a conceptual account of how to implement the resulting updates rather than providing a further theoretical contribution, and may be skipped by readers interested mainly in empirical behavior. Section 6 situates the method among existing approaches, while Sections 7 and 8 explains when the method improves upon classical variational message passing and evaluates it in larger models. The final section discusses implications and future directions. Code availability. Code and reproducibility materials for all experiments are available in the accompanying repository. 2 Background 2.1 Forney-style Factor Graphs A Forney-style Factor Graph (FFG) =(,ℰ)G=(V,E) represents a factorized function, f()=∏a∈fa(a),f( z)\;=\; _a f_a( z_a)\,, (3) where a z_a collects the variables that enter factor faf_a. A Forney-style Factor Graph (FFG) represents such a factorization by drawing factors as nodes and variables as edges (21; 36; 48). An edge is incident to a node exactly when the corresponding variable is an argument of that factor. For example, Figure 2 represents the factorization f()= f( z)\;= fa(z1)fb(z2)fc(z1,z2,z3)fd(z4)fe(z3,z4,z5)fg(z5,z6)fh(z6). f_a(z_1)\,f_b(z_2)\,f_c(z_1,z_2,z_3)\,f_d(z_4)f_e(z_3,z_4,z_5)\,f_g(z_5,z_6)\,f_h(z_6). (4) The boxes in the figure are the seven local factors, and the six edges exist because both factors depend on z3z_3, while the edge z6z_6 connects the pair fgf_g and fhf_h. The dashed boxes are not additional factors; they only mark the parts of the graph whose local eliminations produce the displayed messages. faf_afcf_cfbf_bfdf_dfef_efgf_gfhf_h z1z_1z2z_2z3z_3μ3e→ subarrayc _3e\\ → subarray←μ3c subarrayc←\\ _3c subarrayz4z_4z5z_5←μ5e subarrayc←\\ _5e subarrayz6z_6 Figure 2: Forney-style Factor Graph representation of the factorization in equation 4. The basic computational idea is to avoid integrating all variables at once. To form the marginal belief on the middle edge z3z_3, the left subgraph can first be summarized by a function of z3z_3 alone, μ3e(z3)=∫fa(z1)fb(z2)fc(z1,z2,z3)dz1dz2. _3e(z_3)= f_a(z_1)f_b(z_2)f_c(z_1,z_2,z_3)\,dz_1\,dz_2. (5) Similarly, the right subgraph sends a summary toward z3z_3. In the notation of the figure, μ5e(z5) _5e(z_5) =∫fg(z5,z6)fh(z6)dz6, = f_g(z_5,z_6)f_h(z_6)\,dz_6, (6a) μ3c(z3) _3c(z_3) =∫fd(z4)fe(z3,z4,z5)μ5e(z5)dz4dz5, = f_d(z_4)f_e(z_3,z_4,z_5) _5e(z_5)\,dz_4\,dz_5, (6b) q3(z3) q_3(z_3) ∝μ3e(z3)μ3c(z3). _3e(z_3) _3c(z_3). (6c) The updates above are just a repeated use of the distributive law, but the graph makes the required local computations explicit. Belief propagation is the systematic version of this operation: each factor combines its local function with incoming messages from the neighboring edges and sends the resulting function along another edge. We denote the neighboring edges of a node a∈a by ℰ(a)E(a). Conversely, for an edge i∈ℰi , the notation (i)V(i) collects its neighboring nodes. As a notational convention, we index nodes by a,b,ca,b,c and edges by i,ji,j, unless stated otherwise. In this paper, we will frequently refer to a subgraph. We define an edge-induced subgraph by (i)=((i),i)G(i)=(V(i),i), and a node-induced subgraph by (a)=(a,ℰ(a))G(a)=(a,E(a)). Furthermore, we denote a local subgraph by (a,i)=((i),ℰ(a))G(a,i)=(V(i),E(a)), which collects all local nodes and edges around i and a, respectively. The FFG formalism of Forney allows edges of degree at most two, max|(i)|=2 |V(i)|=2. We use the terminated FFG formalism of 48, in which every edge has degree two, |(i)|=2|V(i)|=2, by terminating each half-edge with a factor proportional to 11. 2.2 Variational Inference For the remainder of this section, we keep the model equation 3 but momentarily forget that it factorizes. We treat f()f( z) as a single unnormalized density on the joint state =(i)i∈ℰ z=( z_i)_i , so that the Bayesian posterior of interest is p()=f()Z,Z=∫f(),p( z)\;=\; f( z)Z, Z\;=\; f( z)\,d z, (7) where Z is the model evidence. As remarked in the introduction, evaluating Z and the marginals of p is intractable in all but the most favorable cases. Variational inference (26; 53; 24) sidesteps this by replacing integration with optimization: pick a tractable family of candidate densities Q and select the q∈q closest to p in some divergence. The canonical choice is the (reverse) Kullback–Leibler divergence. KL[q∥p]=∫q()logq()p()d⩾ 0,D_KL\! [\,q\,\|\,p\, ]\;=\; q( z) q( z)p( z)\,d z\;\ \;0, (8) with equality iff q=pq=p almost everywhere. Substituting equation 7 and rearranging, KL[q∥p]=∫q()logq()f()d⏟=:ℱ[q]+logZ.D_KL\! [\,q\,\|\,p\, ]\;=\; q( z) q( z)f( z)\,d z_ =:\,F[q]\;+\; Z. (9) The functional ℱ[q]F[q] is the variational free energy. Because logZ Z is constant in q, minimizing ℱ[q]F[q] over Q is equivalent to minimizing KL[q∥p].D_KL[\,q\,\|\,p\,]. The further rewriting makes the role of ℱF transparent ℱ[q]=q[−logf()]⏟average energy−ℍ[q]⏟entropy of q=KL[q∥p]−logZ.F[q]\;=\; E_q[- f( z)]_average energy\;\;-\! H\! [\,q\, ]_entropy of q\!=D_KL\! [\,q\,\|\,p\, ]\;-\; Z. (10) The energy–entropy split is the form that appears in statistical physics and underlies the maximum-entropy principle of 25; the KL form shows that ℱF is an upper bound on −logZ- Z, tight when q=pq=p holds almost everywhere. Equivalently, −ℱ-F is the Evidence Lower Bound (ELBO) used throughout the modern variational inference literature (24; 33). Restriction to an exponential family. The minimization equation 9 is still over the infinite-dimensional set of densities. The standard tractability move is to restrict Q to a regular minimal exponential family (14; 53) qλ()=h()exp⊤T()−A(),∈Λ⊆ℝd,q_λ( z)\;=\;h( z) \! \ λ T( z)-A( λ) \, λ∈ ^d, (11) with sufficient statistic T()T( z), base measure h()h( z), and log-partition function A()=log∫h()exp(⊤T())d.A( λ)\;=\; h( z) \! ( λ T( z) )\,d z. (12) Two further objects that we will need throughout the paper follow from the direct differentiation of equation 12: ()=∇λA()=qλ[T()], μ( λ)\;=\; _λA( λ)\;=\;E_q_λ\! [T( z) ], (13a) G()=∇λ2A()=Covqλ[T()]≻0. G( λ)\;=\;∇^2_λA( λ)\;=\;Cov_q_λ\! [T( z) ] 0. (13b) The map ↦() λ μ( λ) is a global diffeomorphism: bijective and smooth, with the smooth inverse ↦() μ λ( μ). Thus ∈Λ λ∈ and ∈ℳ μ are two coordinate systems for the same exponential-family distribution, namely natural parameters and mean parameters. Here ℳM is the interior of the convex hull of the support of T, and the Jacobian of the coordinate map is precisely the Fisher information G()G( λ) in natural coordinates (5). The Legendre dual of A is A∗()=supλ(⊤−A())=qλ(μ)[logqλ(μ)()]−qλ(μ)[logh()],A^*( μ)\;=\; _λ\; ( λ μ-A( λ) )\;=\;E_q_λ(μ)\! [ q_λ(μ)( z) ]-E_q_λ(μ)[ h( z)], (14) so that −ℍ[qλ]-H\! [\,q_λ\, ] is, up to the base-measure constant, the dual A∗()A^*( μ) (53). Substituting equation 14 into the free energy, equation 10, gives the free energy in natural-parameter coordinates ℱ()=A∗(())−qλ[logf()]+const,F( λ)\;=\;A^*( μ( λ))\;-\;E_q_λ[ f( z)]\;+\;const, (15) which the variational problem now minimizes over the finite-dimensional set Λ . The natural gradient and Khan stationary condition. The Fisher information equation 13b equips the exponential family with the Fisher geometry: local displacements in λ are measured by their effect on the corresponding distribution, rather than by the Euclidean distance between parameter vectors. Informally, the natural gradient is the steepest ascent direction for a scalar F:Λ→ℝF: measured in this geometry. More precisely, it is the Riesz representation of the differential of F under the Fisher metric (4; 5) ∇~λF:=G()−1∇λF. ∇_λF\;:=\;G( λ)^-1 _λF. (16) A clean identity simplifies this enormously. Because =∇A() μ=∇ A( λ) has Jacobian G()G( λ), the chain rule gives ∇λF=G()∇μF _λF=G( λ)\, _μF for any F that is smooth on the manifold, hence ∇~λF=∇μF. ∇_λF\;=\; _μF. (17) The natural gradient in natural coordinates equals the ordinary gradient in mean coordinates. This duality is the cornerstone of 4 and underlies the natural gradient variational inference literature (24; 27). Applied to the variational free energy equation 15, there are two elementary derivative facts ∇λA∗(()) _λA^*( μ( λ)) =G(), =G( λ) λ, (18a) ∇λqλ[logf] _λE_q_λ[ f] =Covqλ[T(),logf()] =Cov_q_λ\! [T( z), f( z) ] (18b) that combine into the stationary condition ∇λℱ()= 0⟺=G()−1∇λqλ[logf()]=∇~λqλ[logf()]=∇μqλ[logf()]. _λF( λ)\;=\;0\; \; λ\;=\;G( λ)^-1 _λE_q_λ[ f( z)]\;=\; ∇_λE_q_λ[ f( z)]\;=\; _μE_q_λ[ f( z)]. (19) This is the Khan stationary condition (30; 27; 28; 29): at any minimizer of the variational free energy over an exponential family, the natural parameter λ equals the natural gradient of the expected log-model. Equation equation 19 is therefore both an information-geometric characterization of the variational fixed point and a recipe for finding it by mirror-descent or natural-gradient iteration (30; 28). 3 Variational Message Passing Section 2.2 treated f()f( z) as a structureless joint and committed the variational density to a single global exponential family. We now do the opposite: we keep the variational objective equation 9, but reinstate the factorization of f through the FFG G. The functional form of the approximation is then not prescribed as one global density qλq_λ, as in equation 11. It is induced by the admissible belief space and by any additional constraints imposed on that space (48). This is the point we need for the next section: the variational view makes the admissible belief space itself a knob. Bethe inference and the local polytope. The Bethe construction (58; 59) replaces the joint q by local beliefs: a node belief qa(a)q_a( z_a) for every factor a∈a and an edge belief qi(zi)q_i(z_i) for every edge i∈ℰi . On a TFFG, where every edge has degree two, these beliefs define the Bethe approximation q()=∏a∈qa(a)∏i∈ℰqi(zi),q( z)\;=\; _a q_a( z_a) _i q_i(z_i), (20) with admissible beliefs constrained by the local polytope (53) ℒ()=(qa,qi):∫qa(a)da= 1,∀a∈,∫qa(a)da∖i=qi(zi),∀a∈,i∈ℰ(a).L(G)\;=\; \\,(q_a,q_i)\;:\; aligned & q_a( z_a)\,d z_a\;=\;1,&&∀\,a ,\\ & q_a( z_a)\,d z_a i\;=\;q_i(z_i), &&∀\,a ,\ i (a) aligned \. (21) Substituting equation 20 into the variational free energy gives the Bethe free energy (59; 48) ℱℬ[qa,qi]=∑a∈∫qa(a)logqa(a)fa(a)da+∑i∈ℰ∫qi(zi)log1qi(zi)dzi.F_B [\q_a\,\q_i\ ]\;=\; _a q_a( z_a) q_a( z_a)f_a( z_a)\,d z_a\;+\; _i q_i(z_i) 1q_i(z_i)\,dz_i. (22) Thus Bethe inference is the constrained variational problem minqa,qi∈ℒ()ℱℬ[qa,qi]. _\q_a,q_i\ (G)F_B [\q_a\,\q_i\ ]. We will use the standard Lagrangian notation for this problem. With scalar multipliers ψa,ψi _a, _i for normalization and function-valued multipliers λai(zi) _ai(z_i) for marginalization, ℒ= \;= ℱℬ[qa,qi]+∑a∈ψa(∫qa(a)da−1)+∑i∈ℰψi(∫qi(zi)dzi−1) _B [\q_a\,\q_i\ ]+ _a _a\! ( q_a( z_a)\,d z_a-1 )+ _i _i\! ( q_i(z_i)\,dz_i-1 ) (23) + + ∑a∈∑i∈ℰ(a)∫λai(zi)(qi(zi)−∫qa(a)da∖i)dzi. _a _i (a) _ai(z_i)\! (q_i(z_i)- q_a( z_a)\,d z_a i )dz_i. Known stationary consequences. By Theorem 1 of 48, stationarity of equation 23 over the unmodified local polytope gives the ordinary sum-product equations, recovering the Bethe stationary conditions of 58; 59. Writing the marginalization multipliers as messages, μi→a(zi):=exp(λai(zi)), _i→ a(z_i)\;:=\; \! ( _ai(z_i) ), (24) the stationary beliefs have the BP product form qa∗(a)∝fa(a)∏i∈ℰ(a)μi→a(zi),q_a^*( z_a)\; \;f_a( z_a) _i (a) _i→ a(z_i), (25) and, on a degree-two edge (i)=b,cV(i)=\b,c\, qi∗(zi)∝μb→i(zi)μc→i(zi),μb→i(zi):=exp(λbi(zi)).q_i^*(z_i)\; \; _b→ i(z_i)\, _c→ i(z_i), _b→ i(z_i)\;:=\; \! ( _bi(z_i) ). (26) Closing the stationary system with the marginalization constraints recovers the BP update μa→j(k+1)(zj)=∫fa(a)∏i∈ℰ(a)∖jμi→a(k)(zi)da∖j. _a→ j^(k+1)(z_j)\;=\; f_a( z_a) _i (a) j _i→ a^(k)(z_i)\,d z_a j. (27) The important point for us is that BP is what one obtains when the admissible belief space is exactly ℒ()L(G): no finite-dimensional form for qaq_a or qiq_i has been imposed. A known way to tune the admissible belief space is to add factorization constraints on node beliefs. Structured VMP splits the incident edges of a factor a into a cluster partition (a)⊆(ℰ(a))C(a) (E(a)) and requires (20; 48) qa(a)=∏n∈(a)qan(an),q_a( z_a)\;=\; _n (a)q_a^n( z_a^n), (28) where an z_a^n collects the variables of cluster n. By Theorem 2 of 48, the corresponding stationary update is qan,∗(an)∝f~an(an)∏i∈nμi→a(zi),f~an(an):=exp(∏m≠nqam[logfa(a)]).q_a^n,*( z_a^n)\; \; f_a^\,n( z_a^n) _i∈ n _i→ a(z_i), f_a^\,n( z_a^n)\;:=\; \! (\,E_ _m≠ nq_a^m\! [ f_a( z_a) ] ). (29) The coarsest partition recovers BP, while the finest partition gives mean-field VMP; intermediate partitions give the structured hierarchy described by 48. Why constraints matter. The computational effect is visible in the Normal–Gamma observation model x∼(μ,σ2),τ∼(a,b),y∣x,τ∼(x,τ−1).x (μ,σ^2), τ (a,b), y x,τ (x,τ^-1). Exact BP through the likelihood factor fy(x,τ)=(y∣x,τ−1)f_y(x,τ)=N(y x,τ^-1) does not preserve Gaussian and Gamma messages simultaneously: a Gamma belief on τ yields a heavy-tailed, non-Gaussian message to x, while a Gaussian belief on x yields a non-Gamma message to τ. The mean-field split qy(x,τ)=qyx(x)qyτ(τ)q_y(x,τ)=q_y^x(x)q_y^τ(τ) instead replaces fyf_y by tilted factors of the form equation 29. Since logfy(x,τ)=12logτ−12τ(y−x)2+const f_y(x,τ)= 12 τ- 12τ(y-x)^2+const, these factors are Gaussian in x and Gamma in τ. Factorization therefore restores finite closure by averaging out the coupling; Section 5.3 shows the uncertainty correction that this averaging omits. From node factorizations to edge form constraints. Structured VMP illustrates the new degree of freedom exposed by the variational formulation: we can change the feasible set and derive the corresponding message updates from stationarity. A factorization constraint changes the dependence structure of a node belief, but it does not by itself introduce a natural parameter carried by an edge. The next section studies a different modification of the admissible belief space. On selected edges, we impose an exponential-family form constraint qi∈ℰi,qi=qλi.q_i _i, q_i=q_ _i. This is the step that introduces a finite coordinate i λ_i into the otherwise functional Bethe Lagrangian. The exact Lagrange multipliers remain functions, but stationarity with respect to i λ_i sees only the tangent directions of the receiving edge family. The result is a finite message: the Fisher-metric projection of the local log-message onto that edge family. With no additional node factorization, this local log-message is the exact BP log-message; under structured VMP, it is the tilted log-message induced by equation 29. 4 Natural-Gradient Message Passing We now localize the global natural-gradient stationarity condition equation 19. In the global formulation, one chooses a single exponential-family density qλ()q_λ( z) over all latent variables and obtains =∇μqλ[logf()]. λ= _μE_q_λ[ f( z)]. The factorization of f is invisible to this identity except through the global expectation. The Bethe formulation provides a different knob: instead of choosing one family for the whole joint density, we can change the admissible belief space edge by edge. The question of this section is whether imposing an exponential-family form constraint on a single edge belief recovers a local version of Khan’s identity. Edge-constrained belief space. Fix an edge i∈ℰi with incident factors (i)=b,cV(i)=\b,c\, and equip it with a regular minimal exponential family ℰi=qλi:i∈ΛiE_i=\q_ _i\,:\, λ_i∈ _i\. We absorb the carrier into the edge reference measure, so qλi(zi)=expi⊤Ti(zi)−Ai(i),i∈Λi⊆ℝdi,q_ _i(z_i)\;=\; \! \ λ_i T_i(z_i)-A_i( λ_i) \, λ_i∈ _i ^d_i, (30) with mean parameter i=∇λiAi(i) μ_i= _ _iA_i( λ_i) and Fisher information Gi(i)=∇λi2Ai(i)G_i( λ_i)=∇^2_ _iA_i( λ_i). Compared with the local polytope equation 21, the new admissible belief space replaces the free edge density qiq_i with qλiq_ _i: ℒℰi()=(qa,qjj≠i,i):∫qa(a)da=1,∀a∈,∫qa(a)da∖j=qj(zj),∀a∈,j∈ℰ(a)∖i,∫qa(a)da∖i=qλi(zi),∀a∈(i).L_E_i(G)= \\,(\q_a\,\q_j\_j≠ i, λ_i)\;:\; aligned & q_a( z_a)\,d z_a=1,&&∀ a ,\\ & q_a( z_a)\,d z_a j=q_j(z_j),&&∀ a ,\ j (a) \i\,\\ & q_a( z_a)\,d z_a i=q_ _i(z_i),&&∀ a (i) aligned \. (31) The factor beliefs and exact marginalization multipliers remain functional; only the edge marginal has acquired the finite coordinate i λ_i. Substituting equation 30 into the Bethe Lagrangian gives the only new edge-coordinate term ℒi(i,λaia∈(i))=−Ai∗(i(i))+∑a∈(i)qλi[λai(zi)].L_i ( λ_i,\ _ai\_a (i) )\;=\;-\,A_i^*\! ( μ_i( λ_i) )\;+\; _a (i)E_q_ _i\! [ _ai(z_i) ]. (32) The first term is the edge entropy written as the negative Legendre dual, and the edge normalization multiplier drops out because qλiq_ _i is normalized. Differentiating equation 32 gives the new marginal stationarity condition 0=−Gi(i)i+∑a∈(i)Covqλi[Ti,λai].0=-G_i( λ_i) λ_i+ _a (i)Cov_q_ _i\! [T_i, _ai ]. (33) Theorem 1 (Natural-gradient message passing on a form-constrained edge). Consider the Bethe free energy equation 22 over the edge-constrained belief space ℒℰi()L_E_i(G) in equation 31. For a real-valued function r(zi)r(z_i), define the tangent update at the receiving edge marginal by iλi[r]:=Gi(i)−1Covqλi[Ti,r]=∇μiqλi[r]. T_i _i[r]:=G_i( λ_i)^-1Cov_q_ _i\! [T_i,r ]= _ _i\,E_q_ _i[r]. (34) The operator iλi T_i _i maps a function-valued update on the edge to the finite tangent, equivalently natural-parameter directions visible at the current receiving marginal qλiq_ _i. Let λaj _aj denote the marginalization multipliers of the constrained Bethe Lagrangian. For each incident factor a∈(i)a (i), define the exact BP log-message μa→i(zi) _a→ i(z_i) :=∫fa(a)∏j∈ℰ(a)∖iexp(λaj(zj))da∖i, :=\; f_a( z_a) _j (a) i \! ( _aj(z_j) )\,d z_a i, (35) ℓa→i(zi) _a→ i(z_i) :=logμa→i(zi), :=\; _a→ i(z_i), and its projected natural parameter a→i:=iλi[ℓa→i]=∇μiqλi[ℓa→i(zi)]. η_a→ i:= T_i _i\! [ _a→ i ]= _ _iE_q_ _i\! [ _a→ i(z_i) ]. (36) At any stationary point of the constrained Bethe Lagrangian on the edge (i)=b,cV(i)=\b,c\, the edge natural parameter satisfies i=b→i+c→i. \; λ_i\;=\; η_b→ i\;+\; η_c→ i\;.\; (37) For each a∈(i)a (i), the corresponding projected message is μ^a→i(zi)∝expa→i⊤Ti(zi). μ_a→ i(z_i)\; \; \! \ η_a→ i T_i(z_i) \. (38) Proof sketch. The only new variational variable on edge i is i λ_i, so differentiating equation 32 gives equation 33. Factor-side stationarity still has the BP product form. Closing the marginalization constraint between qaq_a and qλiq_ _i gives, for each a∈(i)a (i), the identity logqλi(zi)=λai(zi)+ℓa→i(zi)+const. q_ _i(z_i)= _ai(z_i)+ _a→ i(z_i)+const. Applying iλi T_i _i to this identity and combining the two incident factors with equation 33 yields equation 37. The full algebra is given in Appendix A. Appendix B gives a smaller tutorial derivation of the same mechanism on a one-variable, two-factor graph. Readers who are less familiar with constrained Bethe Lagrangians may find it useful to read that appendix before the full proof. From the stationary condition to a finite message. Theorem 1 should be read as a projected stationarity statement. It is the localized form of the Khan stationary condition equation 19: the global logf() f( z) is replaced edge by edge by the incoming BP log-message ℓa→i(zi) _a→ i(z_i), the only function of the surrounding graph visible from edge i. The exact BP log-message need not lie in the tangent space of the receiving exponential family. At the current marginal qλiq_ _i, it can be decomposed into the part visible through the sufficient statistics and an orthogonal residual, ℓa→i=a→i⊤Ti+ℓa→i⟂+const,iλi[ℓa→i⟂]=0. _a→ i= η_a→ i T_i+ _a→ i +const, T_i _i[ _a→ i ]=0. Here, the projection is written in log-message coordinates: a log-potential r represents the infinitesimal multiplicative update qλiϵ∝qλiexpϵrq_ _i^ε q_ _i \ε r\, whose density tangent is qλi(r−qλi[r])q_ _i(r-E_q_ _i[r]). In these coordinates, the tangent–normal decomposition and the induced projected message are shown schematically in Figure 3. Pℰ=qλ:λ∈ΛE=\q_λ:λ∈ \TqλℰT_q_λEprojected messageμ^a→i∝eiλi[ℓa→i]⊤Ti μ_a→ i e T_i _i[ _a→ i] T_iℓ⟂ _ exact BP messageμa→i=eℓa→i _a→ i=e _a→ iqλq_λ Figure 3: Geometric view of the projected message. The exact BP message is μa→i=expℓa→i _a→ i= \ _a→ i\, a tangent in the ambient probability space P. The decomposition is drawn in log-message coordinates: ℓa→i _a→ i is split into the tangent component a→i⊤Ti η_a→ i T_i and the orthogonal residual ℓa→i⟂ _a→ i at the current marginal qλiq_ _i. The actual finite message sent along the constrained edge is obtained by exponentiation of the tangent component, μ^a→i∝expa→i⊤Ti μ_a→ i \ η_a→ i T_i\. The residual may remain in the functional Lagrange multipliers, but it is invisible to the finite-dimensional edge marginal. Making the ambient space picture fully rigorous would require a non-parametric density-manifold formulation; see 55 for tangent and cotangent spaces of probability densities and Fisher–Rao/Wasserstein information metrics. We use the ambient manifold here informally to explain what is the finite-dimensional projection. The exact functional Lagrange multipliers can contain this residual, but the finite edge marginal responds only to the tangent component defined in equation 36. Natural-gradient message passing. We call natural-gradient message passing the finite message passing scheme that uses this tangent component as the natural parameter of the message propagated to the rest of the graph. That is, instead of sending the exact BP message μa→i _a→ i when it is not representable in the chosen family, NGMP sends μ^a→i(zi)∝expa→i⊤Ti(zi). μ_a→ i(z_i) \ η_a→ i T_i(z_i)\. Equivalently, NGMP sets the orthogonal residual ℓa→i⟂ _a→ i to zero for communication purposes and keeps only the actual finite-dimensional change visible to the receiving family. Thus Theorem 1 motivates NGMP by identifying the projected fixed-point condition that these finite messages must satisfy. The extension from one constrained edge to an arbitrary subset ⊆ℰS is edge-local. Imposing qi∈ℰiq_i _i on every i∈i gives one equation of the form equation 37 per constrained edge, while the exact BP log-messages equation 35 still couple those equations through the surrounding graph. Unconstrained edges continue to obey the classical BP update of Section 3. If the edge form constraint is the only additional constraint touching factor a, then ℓa→i _a→ i in equation 36 is the exact BP log-message. If a structured-VMP factorization constraint is also imposed on a, the same gradient rule applies with ℓa→i _a→ i replaced by the tilted-factor log-message induced by equation 29. In that mean-field setting one projection step is the non-conjugate VMP update of 32; see Section 6. The message depends on the marginal it is sent to. A striking feature of equation 36 is now immediate: a→i η_a→ i, and hence μ^a→i μ_a→ i, depends on the receiving edge’s own marginal qλiq_ _i. This is the formal version of the NGMP information flow previewed in Figure 1: the factor still forms the exact BP log-message from its local factor and incoming messages, but the finite message sent onward is obtained by projecting that log-message at the current receiving marginal. In BP and structured VMP, the outgoing message is computed from the factor and its other incoming messages; in EP, the outgoing site update is built from a cavity/tilted construction. NGMP is different because the receiving marginal itself is the projection point. Operationally, the update is therefore a Fisher-metric fixed-point iteration that repeatedly re-projects the exact BP log-message at the current edge marginal. The next section gives the same replacement a surrogate-model interpretation: the projected message is implemented by replacing a non-representable local contribution with a conjugate surrogate whose natural parameter is the tangent projection. 5 Surrogate Models Theorem 1 gives a projected message, not a new global posterior family. This section explains how such messages can be implemented by replacing non-representable local contributions with conjugate surrogate leaves. The first part treats non-conjugate unary leaves. This is essentially conjugate-computation variational inference in FFG language: replacing a difficult local likelihood with a conjugate pseudo-observation creates an auxiliary graph on which ordinary BP can be run. Figure 4 shows the complete tree both as a Gaussian surrogate model and as NGMP on the original graph. Later in Section 5.2, We treat factors with multiple constrained interfaces, such as (y∣x,τ−1)N(y x,τ^-1). This is where the local FFG view matters more. Rather than inventing one global Normal–Gamma variational family for the joint (x,τ)(x,τ) block, each interface receives its own projected conjugate message. Figure 5 then shows how replication turns these per-interface updates into a loop. 5.1 Surrogate Leaves Consider the Gaussian state-space model with non-conjugate Poisson observations used by 30. Writing (y∣r)P(y r) for a Poisson likelihood with rate r, p(z0:K,y1:K)=(z0∣m0,v0)∏k=1K(zk∣zk−1,σ2)∏k=1K(yk∣expzk).p(z_0:K,y_1:K)\;=\;N(z_0 m_0,v_0) _k=1^KN(z_k z_k-1,σ^2) _k=1^KP(y_k z_k). (39) Its FFG is a tree: the Gaussian transition factors form a chain, and each Poisson likelihood is a unary leaf attached to one state edge. The transition factors are conjugate to Gaussian messages, but the observation leaf is not. The exact BP message from the observation leaf to the edge zkz_k is the likelihood itself, μyk→k(zk) _y_k→ k(z_k) ∝(yk∣expzk), \;P(y_k z_k), (40) ℓk(zk) _k(z_k) :=logμyk→k(zk)=ykzk−expzk+const. :=\; _y_k→ k(z_k)\;=\;y_kz_k- z_k+const. This exact log-message is not quadratic in zkz_k, so it cannot be passed unchanged inside a Gaussian message-passing scheme. The surrogate-leaf implementation is to keep the exact log-message as the local object being projected, but communicate only its Gaussian tangent component at the receiving edge marginal. Figure 4 places this unary replacement back into the full model for two observation slices. The right panel keeps the original Poisson tree and labels the projected NGMP messages; the left panel represents the same messages as Gaussian pseudo-observation factors. Because every Poisson leaf touches only one latent edge, both views remain trees. The iteration is required because each projected site must be recomputed at its updated Gaussian marginal, not because of a cycle in the FFG. (a) Surrogate-model viewNp(z0)p(z_0)N==N==N⋯·sz0z_0z1z_1z2z_2z3z_3⋆N ⋆N y~1(t) y_1^(t)Λ1(t) _1^(t)y~2(t) y_2^(t)Λ2(t) _2^(t)(b) Message passing viewNp(z0)p(z_0)N==N==N⋯·sz0z_0z1z_1z2z_2z3z_3PPμ^y1→z1(t) μ_y_1→ z_1^(t)μ^y2→z2(t) μ_y_2→ z_2^(t)q1(t−1)q_1^(t-1)q2(t−1)q_2^(t-1)y1y_1y2y_2 Figure 4: Two equivalent implementations of the projected updates on two displayed observation slices of equation 39; both factor graphs are trees and the dots indicate continuation to later states. (a) Surrogate-model view. Each Poisson observation leaf is represented by a Gaussian pseudo-observation ⋆=(y~k(t)∣zk,(Λk(t))−1)N =N( y_k^(t) z_k,( _k^(t))^-1), so ordinary Gaussian BP runs on the entire auxiliary tree. (b) NGMP message view. The original Poisson factors remain visible; each blue label is the projected observation message and the adjacent red marginal defines its tangent space. Inserting a Gaussian factor with the blue message’s natural parameters produces the corresponding leaf in panel (a). Applying this conversion independently at every k gives equation 45. The same idea generalizes per interface: a multi-interface factor sends a separate projected message through each constrained interface, as shown in Figure 5. Constrain the edge marginal to the univariate Gaussian family qk(zk)=(zk∣mk,vk)q_k(z_k)=N(z_k m_k,v_k) with sufficient statistics Tk(zk)=(zk,zk2)⊤T_k(z_k)=(z_k,z_k^2) . It is convenient to use mean coordinates (mk,sk)(m_k,s_k), where sk=qk[zk2]=mk2+vks_k=E_q_k[z_k^2]=m_k^2+v_k. Under qkq_k, qk[ℓk]=ykmk−ρk+const,ρk:=qk[expzk]=exp(mk+vk2).E_q_k[ _k]\;=\;y_km_k- _k+const, _k\;:=\;E_q_k[ z_k]\;=\; \! (m_k+ v_k2 ). (41) Since vk=sk−mk2v_k=s_k-m_k^2, differentiating equation 41 with respect to the Gaussian mean parameters gives the projected natural parameter ηyk→k⋆=∇(mk,sk)qk[ℓk]=[yk+(mk−1)ρk−ρk/2]. _y_k→ k \;=\; _(m_k,s_k)E_q_k[ _k]\;=\; bmatrixy_k+(m_k-1) _k\\[2.0pt] - _k/2 bmatrix. (42) This is equation 36 in the unary case: the non-quadratic log-message is replaced by the Gaussian message whose natural parameter is the mean-coordinate gradient of its expectation. Writing the projected Gaussian message in canonical form, μ^yk→k⋆(zk)∝expξk⋆zk−12τk⋆zk2,ξk⋆=yk+(mk−1)ρk,τk⋆=ρk, μ_y_k→ k (z_k)\; \; \! \ _k z_k- 12 _k z_k^2 \, _k =y_k+(m_k-1) _k, _k = _k, (43) and completing the square yields the pseudo-observation form μ^yk→k⋆(zk)∝(y~k⋆∣zk,(Λk⋆)−1),Λk⋆=ρk,y~k⋆=mk+yk−ρkρk. μ_y_k→ k (z_k)\; \;N( y_k z_k,( _k )^-1), _k = _k, y_k =m_k+ y_k- _k _k. (44) After every Poisson leaf has been rewritten this way, the auxiliary model at an outer iteration t is a normal-normal chain, p~(t)(z0:K∣y~1:K(t))∝(z0∣m0,v0)∏k=1K(zk∣zk−1,σ2)∏k=1K(y~k(t)∣zk,(Λk(t))−1). p^(t)(z_0:K y_1:K^(t))\; \;N(z_0 m_0,v_0) _k=1^KN(z_k z_k-1,σ^2) _k=1^KN( y_k^(t) z_k,( _k^(t))^-1). (45) On this surrogate graph, ordinary Gaussian BP is sufficient. The outer loop carries the non-conjugacy: run Gaussian BP, read the updated Gaussian edge marginals, recompute the projected Poisson leaves, and repeat. This is the FFG version of conjugate-computation variational inference. The old ingredient is the conversion of a non-conjugate local term into a conjugate surrogate site; the NGMP theorem identifies that site’s natural parameter as the projected BP log-message at the receiving marginal. Since the Poisson factor has only one latent interface, the two views in Figure 4 are particularly simple. We next retain this local conversion but let one factor touch two shared latent interfaces. 5.2 From Leaves to Multi-Interface Factors The leaf case hides the main reason for the edge-local formulation. Consider the simplest repeated factor with two latent interfaces: N conditionally i.i.d. Normal observations sharing one unknown mean and one unknown precision, x x ∼(m0,v0),τ∼(a0,b0), (m_0,v_0), τ (a_0,b_0), (46) yn|x,τ y_n x,τ ∼(x,τ−1),n=1,…,N. (x,τ^-1), n=1,…,N. The observation factor fn(x,τ)=(yn∣x,τ−1)f_n(x,τ)=N(y_n x,τ^-1) touches both the Gaussian mean interface x and the Gamma precision interface τ. As discussed in Section 3, exact BP through this factor does not preserve both families simultaneously. A single joint Normal–Gamma surrogate family for all local quantities would be awkward, and it is unnecessary for message passing. The FFG view only requires each constrained interface to receive a finite message in its own family. Despite having only two unknowns, the FFG is already loopy with two observations. In Figure 5, the two Normal likelihoods give two paths between the equality constraints for x and τ, closing the diamond cycle. (a) Variational message passingNp(x)p(x)Gp(τ)p(τ)====xτ 1y_1y2y_2q(t−1)(x)q^(t-1)(x)μ~y1→τ(t) μ_y_1→τ^(t)(b) Natural-gradient message passingNp(x)p(x)Gp(τ)p(τ)====xτ 1y_1y2y_2μx∖1(t) _x 1^(t)μ^y1→τ(t) μ_y_1→τ^(t)q(t−1)(τ)q^(t-1)(τ)μτ∖2(t) _τ 2^(t)μ^y2→x(t) μ_y_2→ x^(t)q(t−1)(x)q^(t-1)(x) Figure 5: The two-observation instance of equation 46 and the message from the first observation toward τ. The Normal prior for x and Gamma prior for τ are at the top; the two (yi∣x,τ−1)N(y_i x,τ^-1) factors close the diamond loop. Solid edges carry x and dashed edges carry τ. (a) VMP uses the previous global marginal q(t−1)(x)∝p(x)∏n=1Nμ~yn→x(t−1)(x)q^(t-1)(x) p(x) _n=1^N μ_y_n→ x^(t-1)(x) to compute μ~y1→τ(t) μ_y_1→τ^(t). (b) NGMP instead uses the cavity μx∖1(t) _x 1^(t) and projects the resulting exact BP log-message at qτ(t−1)q_τ^(t-1). The blue labels enumerate the four messages around the loop; their subscripts give the directions. The red labels are the receiving marginals at which the two blue projected messages are evaluated. The new y1→τy_1\!→\!τ message changes the cavity entering the lower Normal, whose message toward x changes the cavity returning to the upper Normal on the next sweep. The message toward x is analogous. Thus, VMP iterates between two global marginals, whereas NGMP keeps the cavity-specific information required by the joint-belief Bethe fixed point. With N observations, the same pattern couples 2N2N projected likelihood messages, so even exact local projections require loopy message passing on this replicated FFG. 5.3 VMP as the Classical Multi-Interface Surrogate Classical VMP restores finite messages by replacing the coupled observation factor with tilted factors obtained by averaging the log-likelihood over the other interface. For fn(x,τ)=(yn∣x,τ−1)f_n(x,τ)=N(y_n x,τ^-1), f~nx(x)=exp(q(τ)[logfn(x,τ)]),f~nτ(τ)=exp(q(x)[logfn(x,τ)]). f_n^\,x(x)= \! (E_q(τ)[ f_n(x,τ)] ), f_n^\,τ(τ)= \! (E_q(x)[ f_n(x,τ)] ). (47) Averaging the local log-factor over the other interface gives the conjugate forms: logf~nx(x)=−12[τ](yn−x)2+const,logf~nτ(τ)=12logτ−12τ[(yn−x)2]+const. f_n^\,x(x)=- 12E[τ](y_n-x)^2+const, f_n^\,τ(τ)= 12 τ- 12 [(y_n-x)^2]+const. (48) Thus, each observation contributes a Gaussian site on x and a Gamma site on τ. On the precision interface, the 12logτ 12 τ term increases the Gamma shape by 12 12 per observation. For an incoming Gaussian cavity with variance V, however, the exact BP log-message contains the determinant correction −12log(1+Vτ)- 12 (1+Vτ), which VMP omits. This omission produces the familiar precision overconfidence illustrated in (10, Figure 10.4); Section 7 measures its size and when it matters. 5.4 NGMP as the Improved Local Surrogate NGMP keeps the same surrogate-graph construction as VMP–one conjugate surrogate leaf per constrained interface–but changes how the leaf parameters are chosen. For the message from fnf_n to x, form the exact BP log-message by integrating out the other interface against the incoming message it sends to the factor: ℓn→x(x)=log∫(yn∣x,τ−1)μτ∖n(τ)dτ. _n→ x(x)= (y_n x,τ^-1)\, _τ n(τ)\,dτ. (49) For the message from fnf_n to τ, ℓn→τ(τ)=log∫(yn∣x,τ−1)μx∖n(x)dx. _n→τ(τ)= (y_n x,τ^-1)\, _x n(x)\,dx. (50) Each exact log-message is then projected at its own receiving marginal, n→i=∇μiqi[ℓn→i],i∈x,τ. η_n→ i= _ _iE_q_i\! [ _n→ i ], i∈\x,τ\. (51) Then, with μ^n→i∝expn→i⊤Ti μ_n→ i \ η_n→ i T_i\ (remind Equation 38), the x projection is a Gaussian leaf μ^n→x(x)∝expξnx−12Λnx2, μ_n→ x(x)\; \; \! \ _nx- 12 _nx^2 \, (52) and the τ projection is a Gamma leaf μ^n→τ(τ)∝expΔanlogτ−Δbnτ. μ_n→τ(τ)\; \; \! \ a_n τ- b_nτ \. (53) The blue labels in Figure 5b show this dependency for two observations. The projected y1→τy_1\!→\!τ site changes the τ-cavity entering the second Normal factor; its projected message toward x then changes the x-cavity entering the first factor. The receiving marginals shown in red set the tangent spaces of the two projections. The message therefore depends indirectly on its own previous value through the loopy cavity update, rather than only through a global marginal average. The point of the local-surrogate construction is that this NGMP computation also has a surrogate-model interpretation, exactly as the unary Poisson update has the two equivalent views in Figure 4. Once the two projected messages in equation 52 and equation 53 have been computed, freeze their natural parameters and insert one conjugate leaf on each receiving interface. Figure 6 shows this conversion for one observation factor. The original factor remains the object used to compute the two exact cavity log-messages, but ordinary BP is run on the auxiliary graph in which its two projected contributions are represented separately. What distinguishes NGMP from VMP is therefore not the final Gaussian and Gamma leaf families, but their source: NGMP projects the exact cavity log-messages at the receiving marginals, whereas VMP uses the mean-field tilted factors in equation 47. (a) Original multi-interface factorff==ggxx(y∣x,τ−1)N(y x,τ^-1)yyhh==kτ(b) Per-interface surrogate leavesff==ggxx⋆N y~⋆ y Λ⋆ h==kτ⋆G a⋆a b⋆b Figure 6: Surrogate-model interpretation of the multi-interface NGMP computation in Figure 5b. (a) The original observation factor (y∣x,τ−1)N(y x,τ^-1) touches the x and τ equality constraints and produces one projected message for each interface. (b) With those message parameters frozen, the factor’s contributions are represented by a Gaussian leaf ⋆N on x and a Gamma leaf ⋆G on τ; their observed ports denote the projected natural parameters. Neighboring factors f,g,h,kf,g,h,k keep their ordinary BP rules. This is the direct multi-interface analogue of Figure 4a, except that one original factor yields one surrogate leaf per constrained interface. This gives an outer fixed-point map. Let λ collect the natural parameters of the constrained edge marginals, and let ^() G( λ) denote the conjugate surrogate graph obtained by projecting all required exact log-messages at those marginals. During the inner BP sweep, this surrogate graph is frozen. Running ordinary BP on ^() G( λ) returns new constrained-edge natural parameters, which define Φ():=natural parameters returned by BP on ^(). ( λ):=natural parameters returned by BP on G( λ). (54) For a single constrained degree-two edge, Φi() _i( λ) is the sum of the two projected incoming message parameters in equation 37; for a larger graph, Φ stacks these local edge updates after the conjugate BP sweep has propagated their consequences through the surrogate model. The NGMP iteration is therefore (t+1)=Φ((t)), λ^(t+1)= ( λ^(t)), (55) possibly with damping or momentum (Appendix C explains how these stabilize the iteration in practice). A fixed point =Φ() λ= ( λ) is exactly the surrogate-implementation form of the projected stationarity equations in Theorem 1. This is the new contribution relative to the classical surrogate story. NGMP does not claim that the original Normal precision factor has become conjugate, nor does it require a global Normal–Gamma joint approximation. It preserves the edge-local constrained-Bethe stationarity of Theorem 1: each constrained edge chooses its own family, receives its own projected message, and then participates in ordinary BP on the current conjugate surrogate graph. 6 Related Work Variational message passing and local constraints. VMP realizes coordinate-ascent variational inference as local expected-log factor updates (56); structured VMP retains selected clusters, as in equation 29 (20; 48). Factor-graph fragments and Extended VMP change how these local quantities are evaluated but retain this information flow (54; 2). When the model is conjugate, the tilted message of equation 29 already lies in the edge family and no projection is needed; this is the setting in which we use the plain name VMP. Non-conjugate VMP. 32 were the first to bring the Fisher information of the receiving marginal into VMP. When the tilted message μ~a→i=expℓ~a→i μ_a→ i= \ _a→ i\ leaves the edge family, non-conjugate VMP (NCVMP) sends the message whose natural parameter is the Fisher-metric gradient of the expected log-factor at the current marginal, Gi(i)−1∇λiqλi[ℓ~a→i]G_i( λ_i)^-1 _ _iE_q_ _i[ _a→ i]. In the notation of equation 34, this is iλi[ℓ~a→i] T_i _i[ _a→ i]: one tangent projection of the tilted VMP log-message, applied once per update. NCVMP is therefore a hybrid. Its information flow is that of VMP, because out-of-cluster variables enter only through their marginals in f~a f_a, and its projection point is the receiving marginal qλiq_ _i, as in NGMP. Its use of natural gradients is limited to this single step. Conjugate-computation variational inference. 30 generalized this step to CVI, which applies mirror descent to a global ELBO over a fixed-form, possibly mean-field-factorized exponential family. The duality in equation 17 and the stationary condition equation 19 turn the mean-gradient of an expected non-conjugate term into an auxiliary conjugate natural parameter. Khan and Lin call these additive contributions to a global or mean-field coordinate “messages.” Their derivation does not, however, start from an FFG local polytope or derive a distinct outgoing message for every interface of a multi-interface factor. Later work extends the same conjugate-computation principle (27; 28; 34). Khan and Lin’s Gaussian pseudo-observation/Kalman construction also supplies the computation used by the Poisson model in equation 39; our claim concerns its factor-to-edge semantics. FFG-local CVI marginal optimization. Akbayrak et al. make CVI local to an FFG edge by collecting its scheduled messages and optimizing the complete marginal (3). We call this schedule instantiated with the CVI objective of 30 Projective VMP (PVMP). Given the VMP or structured-VMP messages on edge i, write their normalized product as mi(zi):=∏a∈(i)ma→i(zi)m_i(z_i):= _a (i)m_a→ i(z_i). Then qi∗=argminqi∈iKL[qi∥mi].q_i^*= _q_i _iD_KL\! [\,q_i\, \|\,m_i\, ]. (56) The optimizer qi∗q_i^* returns as the next VMP marginal without EP-like division by an opposing message. PVMP runs the natural-gradient iteration on the edge to convergence with the other marginals fixed, whereas NCVMP takes one step of it and moves on; the two share stationary points on the edge but differ in cost per update, which is the budget-matched comparison of Appendix D. Both take VMP messages as input and add a projection at the receiving marginal, so both are hybrids in the above sense. Neither identifies one factor’s natural-gradient contribution with its outgoing message. Q-conjugacy makes this objective analytic for the Poisson–Gaussian model in equation 39 (38). NGMP recovers NCVMP and PVMP under mean-field constraints. NGMP accepts factorization constraints in addition to the edge form constraint. If a mean-field or structured-VMP constraint is imposed on factor a, the log-message ℓa→i _a→ i in equation 36 is the tilted log-message ℓ~a→i _a→ i induced by equation 29, and one projection step of NGMP is exactly the NCVMP update of 32, with the same numerical result. Iterating that projection on the edge marginal to convergence instead of sending the one-step message gives PVMP. Without the factorization constraint, NGMP projects the exact BP log-message, so the factor’s other neighbors enter through their messages rather than their marginals; this is the generalization made in Theorem 1. We follow this distinction in naming the baselines of Section 7 and Section 8: VMP when the model is conjugate and no projection is needed, NCVMP when a single projection step is taken, and PVMP when the projection is iterated to convergence. Expectation propagation and local divergence projections. EP removes a site to form a cavity, restores the factor to form a tilted distribution, performs an inclusive-KL moment projection, and divides by the cavity to recover the site (41). Minka’s divergence framework relates such local objectives more broadly (42). Thus, EP projects a tilted marginal, NCVMP takes one projection step of the tilted VMP message, PVMP optimizes a complete edge marginal, and NGMP projects a factor-to-edge BP log-message. Gaussian surrogate inference. Gaussian surrogate likelihoods are well established: variational Gaussian sites yield a posterior precision equal to the prior precision plus local site precisions (46), consistent with Gaussian Markov structure (47). Such surrogates recover classical filtering and smoothing computations (8; 40). The pseudo-observation and Kalman computation are therefore not contributions; the derived factor-to-edge semantics is. From global gradients to factor-to-edge messages. The formal distinction is already visible in the preceding equations. Khan’s condition equation 19 uses the global log-model; BP instead defines one functional message per factor interface in equation 35. Theorem 1 proves that its projection equation 36 is the corresponding factor-to-edge natural parameter. Hence, an arbitrary-degree factor yields a separate message at each interface, with its own family and receiving marginal, without postulating a global Gaussian, Normal–Gamma, or other joint family. This per-interface statement is stronger than a global additive-gradient decomposition, a one-step NCVMP message, or a complete-marginal PVMP update. 7 Comparing VMP to NGMP: When Does NGMP Matter? VMP and NGMP address the same local approximation problem: both keep the graph and chosen marginal families fixed, but they differ in what information enters each update. VMP computes the update from the expected log-factor under the current neighboring beliefs, whereas NGMP projects the exact cavity log-message. The ablations below hold everything else fixed and ask when this difference matters. We name the baselines as in Section 6. In the first ablation, the model is conjugate, so the baseline is plain VMP with closed-form updates. In the two non-conjugate ablations, the tilted VMP messages leave the edge families, so the baseline is PVMP, which iterates the projective update of equation 56 on every edge to convergence (3); the single-step variant of this projection is NCVMP (32) and serves as the budget-matched control in Appendix D. The computational complexity favors NGMP: it takes a single natural-gradient step per edge update, whereas PVMP runs an inner optimization scheme on every edge until convergence and may therefore evaluate up to a hundred gradients per update. The relevant distinction is edge uncertainty and how often an approximate update is composed. When the neighboring beliefs concentrate, variational expectations become evaluations, and the competing updates agree. When uncertainty persists, the local discrepancy can accumulate either while smoothing information along a latent-state chain or while filtering shared parameters through successive data batches. We isolate these three regimes below: vanishing edge uncertainty, error accumulation along a Poisson smoothing chain, and error accumulation along an online heteroscedastic filtering chain. The broader comparison on real data is deferred to Section 8. Every result is averaged over 20 instances or masks, with scripts provided in the when_ngmp_helps directory of the accompanying code repository. Vanishing edge uncertainty. We first return to the normal mean–precision model of equation 46, with x∼(0,25)x (0,25), τ∼Gamma(2,1)τ (2,1), and N observations. For VMP, this model is conjugate, so its coordinate updates are in closed-form; NGMP’s tangent projection has no closed form here and is computed numerically with an unscented approximation (Section B.3). Both constrained edges receive every observation, so their uncertainty vanishes with N. Figure 7 compares VMP and NGMP with numerically exact marginals, obtained by integrating out x and applying one-dimensional quadrature over τ. For the precision marginal, KL[pexact∥q]D_KL[p_exact\|q] decays approximately as N−2N^-2 for VMP and N−3N^-3 for NGMP: the accuracy ratio grows from 1.7×1.7× at N=4N=4 to ≈700×≈700× at N=512N=512, where matching NGMP’s divergence would require VMP to observe ≈15,000≈15,000 samples. The mean marginals are already nearly identical. The difference is therefore an uncertainty correction: VMP omits the determinant term (Section 5.3) in the exact message toward τ, whereas NGMP retains its tangent component. In this isolated model, the difference has no practical consequence: at N=512N=512, the entire remaining VMP divergence amounts to a 0.2%0.2\% deficit of the posterior precision variance, too small to affect any prediction. The remaining ablations show when it starts to matter: first when the error is repeated while smoothing through the unobserved stretches of a chain, and then when state-dependent approximate batch updates are composed during filtering. (a) Posterior over x (b) Posterior over τ Figure 7: Marginal inclusive divergence in the isolated Normal mean–precision model: (a) KL[p(x∣y)∥q(x)]D_KL[p(x y)\|q(x)] and (b) KL[p(τ∣y)∥q(τ)]D_KL[p(τ y)\|q(τ)] (mean ± 95% CI over 20 instances; each instance draws a fresh true (x,τ)(x,τ) from the priors and generates its own stream of samples). For τ, VMP decays as N−2N^-2 and NGMP as N−3N^-3; for x, NGMP is a roughly constant 2.3×2.3× closer. Error accumulation along a Poisson smoothing chain. Next, we use the Poisson state-space model of equation 39 on the monthly sunspot series (50), with σ2=0.1σ^2=0.1. Unlike in the first ablation, every state edge of the chain is non-conjugate, so the baseline is PVMP: it computes its marginals with the projective update of equation 56 (configuration details in Appendix D). Both inference methods are nevertheless closed-form in this model: NGMP’s projection is a single analytic natural-gradient step per edge, and the gradients iterated by PVMP’s inner optimizer are likewise available analytically (Section B.3). At 55–20%20\% holdout, most missing months neighbor observed ones and the methods tie within their confidence intervals (Table 1). At 50%50\%, long unobserved stretches keep the state edges uncertain and the methods separate. Figure 8 identifies the failure mode. PVMP repeats the mean-field transition update through each missing month, pinning its posterior variance near σ2/2=0.05σ^2/2=0.05 independently of distance from the data. NGMP instead propagates the cavity uncertainty, so the variance grows from 0.090.09 beside an observation to 0.270.27 in the deepest gaps. Accordingly, PVMP’s held-out negative log-likelihood rises from 4.84.8 at distance one to 33.533.5 at distances five to eight, whereas NGMP remains near 5.75.7. Both methods reach their Bethe plateaus within the 20-sweep budget — NGMP by sweep three to four — so the gap is a property of the fixed points, not of early stopping (traces in Figure 15 of Appendix D). Figure 8: Sunspot results at 50%50\% holdout. Top: a five-year window around the deepest held-out stretch of a representative mask, one method per panel on identical axes; lines are the plug-in predictive rates exp([zk]+Var(zk)/2) (E[z_k]+Var(z_k)/2), shaded bands push the pointwise 95%95\% credible intervals of zkz_k through the exponential link, and diamonds are held-out counts. PVMP’s band keeps a fixed width inside the gap and its rate drifts confidently away from the held-out counts, while NGMP’s band widens toward the gap center and keeps covering them. Bottom: held-out negative log-likelihood (log scale) and posterior state variance against distance to the nearest observation over all 20 masks. Table 1: Sunspot state-space model: held-out predictive metrics (mean ± 95% CI over 20 random masks). PVMP and NGMP agree for short gaps and separate when the 50%50\% mask creates long uncertain stretches. The budget-matched NCVMP control is reported in Table 5 of Appendix D. NLL RMSE Held out PVMP NGMP PVMP NGMP 5% 4.534 ± 0.087 4.532 ± 0.085 13.918 ± 0.385 13.935 ± 0.378 10% 4.569 ± 0.058 4.566 ± 0.058 13.972 ± 0.260 13.995 ± 0.260 20% 4.672 ± 0.039 4.661 ± 0.037 14.429 ± 0.214 14.446 ± 0.208 50% 5.411 ± 0.162 4.930 ± 0.035 16.839 ± 0.471 15.789 ± 0.162 Error accumulation along a heteroscedastic filtering chain. Finally, consider heteroscedastic regression: the observation noise is not constant but depends on the input location, so the model must learn two functions from the same data: the mean and the noise level. The oscillatory mean is −(x+12)sin(3πx)-(x+ 12) (3π x), whose true noise standard deviation is 0.45|x+0.5|0.45\, x+0.5 , so the data are nearly noise-free around x=−0.5x=-0.5 and increasingly noisy toward the ends of the input range. The model uses separate random Fourier feature maps ϕ(⋅)φ(·) for the mean and ψ(⋅)ψ(·) for the log-precision, so s_o ∼(ψ(xo)⊤w,ctop−1), \! (ψ(x_o) w,c_top^-1 ), (57a) yo y_o ∼(ϕ(xo)⊤v,e−so). \! (φ(x_o) v,e^-s_o ). (57b) The latent score sos_o is the log-precision of observation o: the observation noise variance at input xox_o is e−soe^-s_o, and sos_o depends on the input through ψ(xo)⊤wψ(x_o) w, so the first line of equation 57 is the learned noise map over input space. Both conditional means are dot products between a fixed feature vector and a weight vector; Figures 9 and 12 draw each such conditional as a softdot node22 2 The softdot (soft dot product) factor is f∗(z∣b,ϕ,τ)=(z∣b⊤ϕ,τ−1),f_*\! (z b,φ,τ )=N\! (z b φ,\,τ^-1 ), (58) a Gaussian over z whose mean is the dot product of the coefficients b and the features ϕφ and whose precision τ arrives on a separate edge; as τ→∞τ→∞ the factor degenerates to a deterministic dot product, hence the name. Figure 9(b) opens the node graphically., reusing the notation of 39. The 400 observations arrive in ten batches of 40, and the model is refit after each batch. With θ=(w,v)θ=(w,v) and qb(θ)=qb(w)qb(v)q_b(θ)=q_b(w)q_b(v) denoting the Gaussian weight marginals after batch b, each filtering step multiplies the previous posterior by the approximate message of the new batch, qb(θ)∝qb−1(θ)μ^ℬb(θ;qb−1),b=1,…,10.q_b(θ) q_b-1(θ)\, μ_B_b(θ;q_b-1), b=1,…,10. (59) If each batch message were an exact Gaussian likelihood, this recursion would be ordinary Kalman-style filtering, and the ten sequential updates would return the same posterior as one fit to all 400 observations at once. Approximate messages break this equivalence: writing ℬU_B for the update induced by a batch ℬB, in general ℬ10∘⋯∘ℬ1≠⋃b=110ℬbU_B_10 ·s _B_1 _ _b=1^10B_b, because each μ^ℬb μ_B_b depends on the posterior it is applied to. The experiment measures how far each method falls short of this ideal: the batching penalty in Table 2 is the increase in held-out NLL of the sequential fit over the full-batch fit, and a method with a small penalty can be trusted to filter online. Figure 9 shows the corresponding factor graph and online schedule. Table 2: Sequential heteroscedastic model over 20 paired seeds (mean ± 95% CI). The same 400 observations are fitted jointly or filtered in ten batches; the batching penalty is the resulting increase in held-out NLL. full batch sequential NLL RMSE NLL RMSE batching penalty PVMP 0.276±0.0310.276± 0.031 0.588±0.0370.588± 0.037 0.512±0.0440.512± 0.044 0.583±0.0390.583± 0.039 0.235±0.0300.235± 0.030 NGMP 0.235±0.0310.235± 0.031 0.591±0.0380.591± 0.038 0.244±0.0310.244± 0.031 0.607±0.0420.607± 0.042 0.009±0.0100.009± 0.010 softdotexpN_\! ==qb−1(w)q_b-1(w)==qb−1(v)q_b-1(v)==softdotexpN_\! ==p(w)p(w)p(v)p(v)softdotexpN_\! ====softdotexpN_\! ====q1(w)q_1(w)q1(v)q_1(v)q2(w)q_2(w)⋯·sq2(v)q_2(v)⋯·sq10(w)q_10(w)q10(v)q_10(v)ℬ1B_1: 40 observations ℬ2B_2: 40 observations ℬ10B_10: 40 observations (c) Filtering schedule over the shared feature weightssoftdotexpN_\! sos_o==qb−1(w)q_b-1(w)w⋯·s==qb−1(v)q_b-1(v)v⋯·sctopc_topψ(xo)ψ(x_o)ϕ(xo)φ(x_o)yoy_oo∈ℬbo _b, |ℬb|=40|B_b|=40(a) Batch box, branch ooexp ⊤ϕb^\! φτ=esτ=e^s softdot expN_\! ssbbϕ ⊤ϕb^\! φτ=esτ=e^s softdot expN_\! ssbbϕ (b) Opened expN_\! Figure 9: Factor graph and online schedule for the sequential heteroscedastic model. Panel (a) shows branch o of the batch box; each softdot node is the soft dot-product factor equation 58. The upper softdot is the log-precision path of equation 57, with b=wb=w, ϕ=ψ(xo)φ=ψ(x_o), τ=ctopτ=c_top, and z=soz=s_o; the observation path below it is the same factor with the uncertain precision τ=esoτ=e^s_o. The lateral stubs with ellipses at the equality nodes continue to the other observations o′∈ℬb∖oo _b \o\, so the same w and v are shared by all 40 branches. Panel (b) opens the double-bordered expN_\! composite and, inside it, the softdot box; the inner box is equation 58 drawn graphically. The deterministic dot node forms the dot product b⊤ϕb φ, the Gaussian N receives it as its mean, and the precision of that Gaussian is the exponentiated incoming score τ=esτ=e^s. In panel (a) the ports are b=vb=v, ϕ=ϕ(xo)φ=φ(x_o), s=sos=s_o, and z=yoz=y_o. Panel (c) repeats the compact batch box for ℬ1,ℬ2,…,ℬ10B_1,B_2,…,B_10 beneath the two shared-weight chains. Batch 1 starts from the model priors p(w)p(v)p(w)p(v); after batch b, the Gaussian marginals qb(w)qb(v)q_b(w)q_b(v) are frozen and copied left-to-right as the next priors. A perpendicular stop bar marks that no message from a later batch is sent back into an earlier batch; the bar specifies the online schedule and is not an additional model factor. With exact state-independent Gaussian messages, this posterior-as-prior multiplication is ordinary associative Kalman-style filtering. Solid edges carry weights, features, and observations; dashed edges carry precisions; dash-dotted edges carry the latent log-precision scores. We use 128 Matérn-3/23/2 random Fourier features plus an intercept for the mean path and 32 RBF random Fourier features plus an intercept for the log-precision path. The mean-path coefficients have a prior standard deviation 2.02.0, while the log-precision RFF coefficients retain a standard deviation 1.61.6. Each of the 400 likelihood factors has an uncertain local log-precision sos_o, but all sites contribute to the same shared noise weights w. Both methods share the graph and marginal families. The observation factor of equation 57 is non-conjugate with respect to its log-precision edge, so the baseline is again PVMP: it computes its expected-log-factor sites there with the projective update of equation 56, whereas NGMP projects the exact cavity messages at the receiving marginals, with the required expectations computed by respective tangent projections (Section B.3). Both methods run under the same 240240-sweep budget; Bethe traces and final-sweep residual changes are given in Appendix D. With all 400 observations fitted together, the held-out RMSE is effectively tied, while NGMP has a small NLL advantage (Table 2). When the same observations are processed in ten batches of 40, the RMSE remains similar, but PVMP’s NLL rises by 0.235±0.0300.235± 0.030 nats. NGMP’s 0.009±0.0100.009± 0.010-nat penalty is consistent with batching invariance. Thus, the predictive distributions, not merely an internal covariance diagnostic, establish PVMP’s batch-partition dependence in this model. The mechanism is visible in the posterior over the noise weights w. After the sequential fit, logdetΣw _w is −72.34-72.34 under PVMP and −8.01-8.01 under NGMP; across the 3333 coordinates, PVMP’s geometric-mean posterior variance per coordinate is therefore roughly 7.0×7.0× smaller, because it has accumulated the local under-dispersion of its 400400 sites into a shared overconfident posterior. The collapse happens early: after the first 4040-point batch, PVMP is already more concentrated (logdetΣw=−27.75 _w=-27.75) than NGMP becomes with all 400 points, and a prior that certain cannot be corrected by later batches (Figure 17 in Appendix D). NGMP instead concentrates as evidence accumulates and ends near its full-batch level (−6.12-6.12 versus −8.01-8.01). Figure 10 shows the two sequential fits. Both methods resolve the oscillatory mean, but NGMP’s sequential predictive density stays close to its full-batch fit while PVMP’s NLL degrades. The learned variance functions in Figure 11 localize the difference: batching distorts PVMP’s posterior-mean variance into a sharp peak near x=2x=2 while its pointwise credible interval remains narrow, whereas NGMP retains substantially more uncertainty about the variance, especially where data are sparse. The separately plotted mean-weight contribution distinguishes this epistemic component from the learned aleatoric variance. The full-batch fits add no visible information at this scale and are deferred to Figures 18 and 19 in Appendix D. (a) PVMP, sequential batches (b) NGMP, sequential batches Figure 10: Posterior predictive bands after the ten-step filtering chain of Figure 9, on one representative seed with identical axes. The orange dashed curve is PVMP and the blue solid curve is NGMP; the matching shaded regions are pointwise 95% posterior-predictive bands. Gray dash-dotted curves show the true mean, and gray points are the 400 training observations. Both methods resolve the oscillatory mean; under posterior reuse, PVMP’s noise-weight covariance collapses and its NLL rises, whereas NGMP remains close to its full-batch NLL. The full-batch fits of both methods are visually indistinguishable from panel (b); Figure 18 in Appendix D shows all four fits. (a) NGMP, sequential batches (b) PVMP, full batch (c) PVMP, sequential batches Figure 11: Posterior predictive variance on the same representative seed as Figure 10. The blue solid curve in panel (a) is NGMP and the orange long-dashed curves in panels (b)–(c) are PVMP; panel (b) uses all observations at once, panels (a) and (c) use the ten sequential batches. Each curve is the posterior mean total predictive variance, S¯(x)=Vepi(x)+Valea(x) S(x)=V_epi(x)+V_alea(x), where Vepi(x)=ϕ(x)⊤Σvϕ(x)V_epi(x)=φ(x) _vφ(x) is epistemic variance from the uncertain mean weights and Valea(x)=q[e−s(x)]V_alea(x)=E_q[e^-s(x)] is learned aleatoric variance. Shading gives a pointwise 95% posterior credible interval for the total. The purple dotted curve shows Vepi(x)V_epi(x) alone, so its vertical gap to the orange or blue total is Valea(x)V_alea(x). The gray dash-dotted curve is the benchmark’s true aleatoric variance [0.45(x+0.5)]2[0.45(x+0.5)]^2. These pointwise credible intervals describe posterior uncertainty for one fit and are distinct from the across-seed confidence intervals in Table 2. All panels use identical logarithmic axes. NGMP’s full-batch fit is visually indistinguishable from panel (a); Figure 19 in Appendix D shows all four fits. Practical rule. Across these ablations, NGMP is indistinguishable from VMP and PVMP when the relevant edge uncertainty is reducible, and it improves calibration or prediction when uncertain updates are repeatedly composed. Hence, when both updates are available for the same constrained interface, NGMP is the safer default. However, the undamped NGMP fixed-point iteration can oscillate and, in practice, may not converge, so every NGMP run in this paper damps its updates in natural coordinates, and the deeper hierarchies of Section 8 additionally use vector-transport momentum (Appendix C). Damping rescales the same single natural-gradient step, and momentum reuses the previous one, so neither adds gradient evaluations. VMP and PVMP remain adequate when neighboring beliefs are already concentrated; persistent uncertainty during smoothing or state-dependent approximate messages reused during filtering are signals that their local error may accumulate. The filtering result is model-specific: exact Gaussian message multiplication remains associative, and NGMP is empirically almost invariant to batching here. Closing PVMP’s cost gap by granting the projective update of equation 56 a single inner step per edge, which is NCVMP, does not help: this control underfits, still collapses along the same concentration path, and yields uniformly worse results at the same runtime, so when PVMP is used, its inner projections should be run to convergence (Tables 6 and 17 in Appendix D). 8 Experiments This section evaluates NGMP in two larger modeling settings. The point of both is that probabilistic hierarchies assembled from the factors of Section 7 learn non-linear dependencies well. The first setting is supervised regression on six UCI data sets, where the depth-three heteroscedastic hierarchy of Figure 12 is compared with several variational Bayesian learning frameworks. This is a framework-level comparison rather than another NGMP–PVMP ablation: direct comparisons with PVMP are given in Section 7, while applying PVMP’s inner manifold optimization at every non-conjugate edge is computationally prohibitive for these data sets and the depth-three hierarchy. Dataset dimensions and an exact model-capacity accounting are reported in Tables 7 and 8. The second composes frozen pretrained forecasters into a probabilistic ensemble, following the setup of 39; the Bayesian model wraps existing predictors without retraining them. Homoscedastic baseline variants never win either comparison, so the main text reports the heteroscedastic variants and defers the homoscedastic ones to Appendix E. 8.1 Regression We evaluate whether the local natural-gradient construction remains useful in a larger supervised-learning setting on six UCI regression data sets: Concrete, Energy, Boston Housing, Power Plant, Wine Quality Red, and Yacht. We use the depth-three instance of the heteroscedastic hierarchy in Figure 12, which generalizes the sequential heteroscedastic model of Section 7. Its first linear-Gaussian path models the predictive mean, and two uncertain precision layers produce input-dependent observation noise. Every layer has its own fixed feature realization and length scale. softdotexpN_\! ⋯·sexpN_\! expN_\! sL,os_L,os1,os_1,octopc_topyoy_oψL(xo) _L(x_o)level LψL−1(xo) _L-1(x_o)level L−1L-1ψ1(xo) _1(x_o)level 11ϕ(xo)φ(x_o)observation==NwLw_L⋯·s==NwL−1w_L-1⋯·s⋯·s==Nw1w_1⋯·s==Nv⋯·s Figure 12: One observation branch of the arbitrary-depth heteroscedastic hierarchy used in the regression experiment. Every level has its own Gaussian weights wℓw_ , shared across observation branches through the equality nodes, and its own clamped feature map ψℓ(xo) _ (x_o). The top-level softdot is the soft dot-product factor equation 58. Each level turns its incoming score into a precision and emits the score sℓ,os_ ,o read by the level below, and the observation level turns the last score into the noise precision of yoy_o. The double-bordered expN_\! composites are opened in Figure 9(b). The experiments in this section use the depth-three instance. Edge styles are as in Figure 9: solid for weights, features, and observations; dashed for precisions; dash-dotted for scores. Each of the feature maps (ϕφ for the mean layer, ψℓ _ for the precision layers) contains 1,000 random Fourier features, sampled once per train–test split and then held fixed; posterior uncertainty is therefore over the linear weights and latent precisions, not over the random features. The Gaussian weight blocks admit closed-form updates, the non-conjugate exponential links on the precision paths are updated with NGMP, and the sweeps are run until convergence. Feature construction, initialization, and optimization settings are given in Appendix E. Prediction integrates the Gaussian weight posteriors: the predictive variance is the sum of the mean-weight epistemic variance and the expected aleatoric variance propagated through the two precision layers. Protocol. We use 20 deterministic repeated-holdout splits with 90% of each data set for training and 10% for testing. Features and targets are standardized using training statistics only. All methods use the same saved split indices, so their results are paired at the level of both splits and held-out examples. The data sets range from 308 to 9,568 observations and from 4 to 13 input features; the per-data-set counts are given in Table 7. We compare NGMP with Bayes by Backprop (B) (12), diagonal and full-covariance deterministic variational inference (dDVI and DVI) (57), Bayesian predictive coding (BPC) (52), and IVON (49). B, dDVI, and DVI exist in a heteroscedastic variant, shown in Table 3, and a homoscedastic variant with fixed observation noise, marked by the prefix ho (hoBBB, hodDVI, hoDVI); BPC and IVON are available only in homoscedastic form (hoBPC, hoIVON). Appendix E explains the difference between the two variants and reports the full tables. The data sets come from the UCI Machine Learning Repository (6), and the repeated-holdout convention follows the one established by 22. The negative log likelihood (NLL) includes the Jacobian required to return from standardized targets to the original target units; RMSE is also reported in original units. Table 3 gives means and approximate 95% confidence-interval half-widths, 1.961.96 times the standard error across the 20 paired splits; the RMSE comparison is Table 10 in Appendix E. Results. NGMP has the lowest mean NLL on Concrete, Energy, Power, and Yacht. On Wine, its interval includes the best point estimate, while on Boston, the two DVI variants have the strongest NLL. The RMSE comparison (Table 10) is more mixed: NGMP is best on Energy and Power, is confidence-interval compatible with the best results on Yacht and Wine, and does not improve on DVI for Concrete or Boston. This distinction is expected for a probabilistic model: RMSE assesses only the predictive mean, whereas NLL also rewards the quality of the learned heteroscedastic uncertainty. In particular, the large NLL gains on Energy and Yacht, together with competitive RMSE, indicate that the precision hierarchy improves the predictive distribution rather than merely shifting its mean. Table 3: Test negative log likelihood (NLL) in original target units; lower is better. Values are means ± approximate 95% confidence-interval half-widths over 20 paired splits. Bold values include the best point estimate within their confidence interval. B, dDVI, and DVI are the heteroscedastic variants; IVON is homoscedastic and uses K=20K=20 posterior network draws for its predictive mixture. The remaining homoscedastic variants are in Table 9. Method Concrete Energy Boston Power Wine Yacht B 3.3328±0.02153.3328± 0.0215 2.4229±0.03502.4229± 0.0350 2.6555±0.03852.6555± 0.0385 2.7785±0.01432.7785± 0.0143 0.9438±0.0232 0.9438± 0.0232 2.6565±0.04512.6565± 0.0451 dDVI 3.0477±0.04323.0477± 0.0432 1.1920±0.4625 1.1920± 0.4625 2.4646±0.0921 2.4646± 0.0921 2.8238±0.01832.8238± 0.0183 0.9440±0.0314 0.9440± 0.0314 0.4641±0.10000.4641± 0.1000 DVI 3.0406±0.04383.0406± 0.0438 1.0094±0.2973 1.0094± 0.2973 2.4489±0.0917 2.4489± 0.0917 2.8249±0.02002.8249± 0.0200 0.9448±0.0307 0.9448± 0.0307 0.4543±0.10960.4543± 0.1096 IVON (K=20K=20) 3.7966±0.00823.7966± 0.0082 3.2693±0.04123.2693± 0.0412 3.3284±0.08553.3284± 0.0855 3.7839±0.00133.7839± 0.0013 1.0430±0.01561.0430± 0.0156 3.6584±0.01013.6584± 0.0101 NGMP 2.9845±0.0433 2.9845± 0.0433 0.9783±0.0655 0.9783± 0.0655 2.6720±0.07522.6720± 0.0752 2.7610±0.0226 2.7610± 0.0226 0.9571±0.0293 0.9571± 0.0293 0.2782±0.1437 0.2782± 0.1437 8.2 Ensemble Forecasting We follow the ensemble-forecasting protocol of 39 and replace only its PVMP inference with NGMP. Interested readers can find details of the ETTh tasks, frozen expert bank, and training and evaluation splits in Section E.2. The controlled edge-uncertainty experiments in Section 7 predict a larger separation in NLL than in point error when uncertain reliability updates are composed. We test that prediction against four neural gates and a precision-gated probabilistic ensemble. The neural baselines are an affine gate and a two-layer ReLU gate trained with Adam (31), and an affine gate and a two-layer ReLU gate trained with IVON (49). Here K=1000K=1000 means that the IVON posterior predictive is approximated with 1,000 gate-weight draws. Each draw produces one Gaussian ensemble component; the reported NLL uses the equally weighted mixture density, and the point forecast averages the component means. The precise component construction is given in Section E.2. The Adam and IVON rows are independently trained neural benchmarks; no causal optimizer effect is inferred from their contrast. For readability, the main table retains the stronger IVON neural baselines; all Adam results remain in the appendix. The gating model is the precision-gated ensemble (PGE) of 39: a depth-one instance of the precision hierarchy in which the frozen experts form the lower level and their predictions enter a precision-weighted mixture, while one softdot with an exponential link maps the shared 65-dimensional context to each expert’s precision. The PVMP row reports the numbers from that paper, computed on the identical test split; the NGMP row runs this paper’s inference on the same model. The affine IVON gate is the closest neural comparison to the precision-gated model: both consume the same 65-dimensional context and frozen seven-expert bank, and both map the context to seven reliability scores. Their scale semantics differ. IVON optimizes the softmax-weighted expert squared error; NGMP, on the other hand, performs inference in an explicit probabilistic precision model whose priors and likelihood determine both relative weights and absolute predictive scale. Table 4: ETTh forecasting in standardized OT units. Entries are test-set point estimates ± approximate 95% normal confidence-interval half-widths over test predictions; RMSE intervals use the delta method. Lower is better. Among the methods shown, the lowest point estimate is bold; the second-lowest is also bold when its paired 95% normal confidence interval relative to the lowest, computed on identical test predictions, includes zero. For IVON, K=1000K=1000 is the number of gate-weight samples drawn from the fitted variational posterior to approximate the posterior predictive. The predictive-mixture construction and complete experimental setup are described in Section E.2; the full results, including Adam-trained gates, are reported in Table 11. ETTh1 Method 96 192 336 720 RMSE MoE – IVON, affine gate (K=1000K=1000) 0.3832±0.00900.3832± 0.0090 0.3695±0.00910.3695± 0.0091 0.3672±0.00840.3672± 0.0084 0.5428±0.01270.5428± 0.0127 MoE – IVON, ReLU gate (K=1000K=1000) 0.3832±0.00900.3832± 0.0090 0.3697±0.00910.3697± 0.0091 0.3677±0.00840.3677± 0.0084 0.5334±0.01250.5334± 0.0125 PGE – PVMP (39) 0.3583±0.00830.3583± 0.0083 0.3386±0.00850.3386± 0.0085 0.3105±0.00730.3105± 0.0073 0.3347±0.00760.3347± 0.0076 PGE – NGMP (this work) 0.3554±0.00840.3554± 0.0084 0.3337±0.00860.3337± 0.0086 0.3115±0.00740.3115± 0.0074 0.3300±0.00800.3300± 0.0080 NLL MoE – IVON, affine gate (K=1000K=1000) 10.2777±0.736610.2777± 0.7366 ×104±×1043.5515\!×\!10^4± 1.1067\!×\!10^4 429.0098±143.4073429.0098± 143.4073 168.4529±26.4000168.4529± 26.4000 MoE – IVON, ReLU gate (K=1000K=1000) 2216.8905±676.27762216.8905± 676.2776 ×108±×1072.2344\!×\!10^8± 9.2280\!×\!10^7 58.6475±9.698058.6475± 9.6980 94.2422±16.263894.2422± 16.2638 PGE – PVMP (39) 0.4120±0.01730.4120± 0.0173 0.3701±0.01710.3701± 0.0171 0.3141±0.01350.3141± 0.0135 0.3763±0.01400.3763± 0.0140 PGE – NGMP (this work) 0.3888±0.02100.3888± 0.0210 0.3378±0.02000.3378± 0.0200 0.2877±0.01610.2877± 0.0161 0.3571±0.01530.3571± 0.0153 ETTh2 Method 96 192 336 720 RMSE MoE – IVON, affine gate (K=1000K=1000) 0.5821±0.01340.5821± 0.0134 0.5198±0.01290.5198± 0.0129 0.5869±0.01420.5869± 0.0142 0.8249±0.01810.8249± 0.0181 MoE – IVON, ReLU gate (K=1000K=1000) 0.5821±0.01340.5821± 0.0134 0.5655±0.01390.5655± 0.0139 0.5874±0.01420.5874± 0.0142 0.6202±0.01290.6202± 0.0129 PGE – PVMP (39) 0.5882±0.01200.5882± 0.0120 0.5798±0.01270.5798± 0.0127 0.5940±0.01280.5940± 0.0128 0.5669±0.01310.5669± 0.0131 PGE – NGMP (this work) 0.5631±0.01240.5631± 0.0124 0.5408±0.01280.5408± 0.0128 0.5463±0.01250.5463± 0.0125 0.6101±0.01430.6101± 0.0143 NLL MoE – IVON, affine gate (K=1000K=1000) ×1021±×10215.9904\!×\!10^21± 6.4398\!×\!10^21 304.6987±52.2098304.6987± 52.2098 ×105±×1053.7880\!×\!10^5± 1.0683\!×\!10^5 ×105±×1043.2209\!×\!10^5± 7.4793\!×\!10^4 MoE – IVON, ReLU gate (K=1000K=1000) ×1018±×10181.9494\!×\!10^18± 1.3836\!×\!10^18 ×105±×1041.1309\!×\!10^5± 6.8836\!×\!10^4 ×106±×1052.0843\!×\!10^6± 5.8045\!×\!10^5 ×106±×1051.1263\!×\!10^6± 2.1811\!×\!10^5 PGE – PVMP (39) 0.9342±0.03060.9342± 0.0306 0.9237±0.03280.9237± 0.0328 0.9612±0.03340.9612± 0.0334 0.8699±0.02990.8699± 0.0299 PGE – NGMP (this work) 0.9397±0.03860.9397± 0.0386 0.8602±0.03650.8602± 0.0365 0.8702±0.03480.8702± 0.0348 0.9774±0.03570.9774± 0.0357 Results. Across the benchmark, RMSE alone often makes the neural gates appear competitive, whereas their NLL is substantially worse than that of both precision-gated models. This larger separation in NLL than in RMSE supports the prediction from the controlled experiments. Against the PVMP results of 39, NGMP has lower NLL at all four ETTh1 horizons and at ETTh2 horizons 192 and 336. It has lower RMSE at ETTh1 horizons 96, 192, and 720 and at ETTh2 horizons 96, 192, and 336. Normal confidence intervals are shown with every point estimate. Full results, including the Adam gates omitted from the compact main table, are reported in Table 11. Figure 13: ETTh1 at H=192H=192: NGMP (solid blue) and PVMP (hatched orange), with 95% predictive intervals, over an illustrative one-week segment. The black solid line is the observed target. This panel illustrates the predictions only; the full-test RMSE and NLL comparison is reported in Table 4. 9 Discussion and Conclusion This paper establishes an edge-local counterpart of the natural-gradient stationarity condition for fixed-form variational inference. Each finite message is the Fisher projection of an exact BP log-message at its receiving marginal, which keeps BP, structured VMP, EP, and NGMP in one local variational language while preserving their different information flows. The experiments locate the practical consequence: when uncertainty persists on the inputs to a non-conjugate factor and approximate updates are repeatedly composed, retaining its representable component improves calibration and prediction. The ETTh study exhibits the same effect in a larger model, where competitive point forecasts alone conceal a collapse of predictive scale. Compiling factor graphs to surrogate programs. The surrogate interpretation also suggests a systems-level research direction. Figure 4 shows two representations of the same projected computation, not two approximations whose posteriors happen to agree: projected messages on the original graph and conjugate leaves on the surrogate graph implement the same update. Figure 6 extends this identity to one surrogate leaf per constrained interface, and Equations 54 and 55 collect the resulting computation in the fixed-point map. RxInfer.jl combines the expressive GraphPPL.jl model representation with ReactiveMP.jl, where messages and marginals are reactive streams whose subscribers trigger local updates (43; 7). This execution model is valuable for streaming and selective recomputation, but it also leaves subscription dispatch and graph traversal in the runtime. A complementary GraphPPL.jl backend could instead traverse the factor graph once, generate the required per-interface projection kernels and conjugate surrogate leaves, fix the inner BP schedule, and emit a typed numerical implementation of the message passing. This lowering would not change the asymptotic arithmetic complexity or the NGMP fixed point; it would move graph orchestration out of the hot loop and expose independent projections and batched linear algebra to compilation by, for example, Reactant.jl, whose MLIR tracing system is being extended with control-flow-aware Julia static analysis (37), or JAX (13). Specialized libraries such as Dynamax already expose pure JAX inference kernels for state-space model classes (35); the opportunity here is to retain GraphPPL’s general modeling interface while obtaining a similarly static execution path. We performed this lowering by hand for the models in this paper. Automating it would remove substantial implementation work and turn the surrogate view into a practical compiler target. Continuous-state active-inference message passing. A second prospect concerns planning. The AIF-MP construction of 44, building on the EFE-as-variational-inference formulation of 45, uses channel-reparameterized observation and dynamics factors but restricts its implementation to discrete state spaces with exact factor evaluations. For continuous states, the exact messages from these modified factors will generally not remain in a tractable finite family. NGMP suggests a principled relaxation: retain the AIF-MP factors and their entropy corrections, impose exponential-family form constraints on the continuous state edges, and project each factor-to-edge log-message according to equation 36. The required local expectations could be analytic or approximated with the sigma-point and quadrature strategies of Section B.3. Unlike a mean-field factor update, this projection preserves the portion of cavity uncertainty visible to the chosen state family, which is precisely the information needed when planning depends on uncertainty about future states. Deriving the full channel-augmented constrained-Bethe scheme and testing its convergence and control performance remain future work; the present result identifies NGMP as a concrete route from discrete exact AIF-MP to approximate continuous-state inference. Together, these prospects suggest a separation between modeling and execution: using factor graphs as the expressive modeling interface and surrogate programs as their compiled numerical representation. NGMP supplies the mathematical translation between the two while retaining the uncertainty needed for downstream problems such as continuous-state planning. Automating this translation and establishing when its continuous-state projections are sufficiently accurate for control are the central next steps. Acknowledgements We gratefully acknowledge financial support by the Dutch Ministry of Economic Affairs (PPS funding), by the Dutch Research Council (NWO) and by hearing aid manufacturer GN Hearing, under contracts TKI-HTSM/21.0161/2112P09 (project: Auto-AR) and KICH3.LTP.20.006 (Project: ROBUST). We are grateful for insightful discussions with Bert de Vries and other BIASlab members. References Absil et al. (2008) P.-A. Absil, R. Mahony, and R. Sepulchre Optimization algorithms on matrix manifolds. Princeton University Press, Princeton, N.J. ; Woodstock (en). Note: OCLC: ocn174129993 External Links: ISBN 978-0-691-13298-3 Cited by: §C.1. Akbayrak et al. (2021) S. Akbayrak, I. Bocharov, and B. De Vries Extended Variational Message Passing for Automated Approximate Bayesian Inference. Entropy 23 (7), p. 815 (en). External Links: ISSN 1099-4300, Link, Document Cited by: §6. Akbayrak et al. (2022) S. Akbayrak, İ. Şenöz, A. Sarı, and B. de Vries Probabilistic programming with stochastic variational message passing. International Journal of Approximate Reasoning 148, p. 235–252 (en). External Links: ISSN 0888613X, Link, Document Cited by: §6, §7. Amari (1998) S. Amari Natural Gradient Works Efficiently in Learning. Neural Computation 10 (2), p. 251–276. External Links: ISSN 0899-7667, Link, Document Cited by: §2.2, §2.2. Amari (2016) S. Amari Information Geometry and Its Applications. Applied Mathematical Sciences, Vol. 194, Springer Japan, Tokyo (en). External Links: ISBN 978-4-431-55977-1 978-4-431-55978-8, Link, Document Cited by: §2.2, §2.2. Asuncion et al. (2007) A. Asuncion D. Newman et al. UCI machine learning repository. Irvine, CA, USA. Cited by: §8.1. Bagaev and De Vries (2023) D. Bagaev and B. De Vries Reactive Message Passing for Scalable Bayesian Inference. Scientific Programming 2023, p. 1–26 (en). External Links: ISSN 1875-919X, 1058-9244, Link, Document Cited by: §9. Bell (1994) B. M. Bell The iterated kalman smoother as a gauss–newton method. SIAM Journal on Optimization 4 (3), p. 626–636. Note: tex.eprint: https://doi.org/10.1137/0804035 External Links: Link, Document Cited by: §6. Bergmann (2022) R. Bergmann Manopt.jl: Optimization on Manifolds in Julia. Journal of Open Source Software 7 (70), p. 3866. External Links: Document Cited by: §C.1. Bishop (2006) C. M. Bishop Pattern Recognition and Machine Learning (Information Science and Statistics). Springer-Verlag. External Links: ISBN 0-387-31073-8, Link Cited by: §5.3. Blei et al. (2017) D. M. Blei, A. Kucukelbir, and J. D. McAuliffe Variational Inference: A Review for Statisticians. Journal of the American Statistical Association 112 (518), p. 859–877. External Links: ISSN 0162-1459, Link, Document Cited by: §1. Blundell et al. (2015) C. Blundell, J. Cornebise, K. Kavukcuoglu, and D. Wierstra Weight uncertainty in neural network. In Proceedings of the 32nd international conference on machine learning, F. Bach and D. Blei (Eds.), Proceedings of machine learning research, Vol. 37, Lille, France, p. 1613–1622. External Links: Link Cited by: §8.1. Bradbury et al. (2018) J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang JAX: composable transformations of Python+NumPy programs. External Links: Link Cited by: §9. Brown (1986) L. D. Brown Fundamentals of statistical exponential families with applications in statistical decision theory. SPIE (en). External Links: Link, Document Cited by: §2.2. Caticha (2004) A. Caticha Relative Entropy and Inductive Inference. AIP Conference Proceedings 707, p. 75–96. Note: arXiv: physics/0311093 External Links: ISSN 0094243X, Link, Document Cited by: Appendix B. Caticha (2011) A. Caticha Entropic Inference. p. 20–29 (en). Note: arXiv:1011.0723 [cond-mat, physics:physics, stat] External Links: Link, Document Cited by: Appendix B. Cooper (1990) G. F. Cooper The computational complexity of probabilistic inference using bayesian belief networks. Artificial Intelligence 42 (2), p. 393–405. External Links: ISSN 0004-3702, Link, Document Cited by: §1. Cox (2001) R. T. Cox Algebra of Probable Inference. Johns Hopkins University Press (en). Note: Read_Status: To Read Read_Status_Date: 2025-06-17T14:44:14.287Z External Links: ISBN 978-0-8018-6982-2, Link, Document Cited by: §1. Dagum and Luby (1993) P. Dagum and M. Luby Approximating probabilistic inference in Bayesian belief networks is NP-hard. Artificial Intelligence 60 (1), p. 141–153. External Links: ISSN 0004-3702, Link, Document Cited by: §1. Dauwels (2007) J. Dauwels On Variational Message Passing on Factor Graphs. In IEEE International Symposium on Information Theory, Nice, France, p. 2546–2550. External Links: Link, Document Cited by: §3, §6. Forney (2001) G. Forney Codes on graphs: normal realizations. IEEE Transactions on Information Theory 47 (2), p. 520–548. External Links: ISSN 0018-9448, Link, Document Cited by: §2.1. Hernandez-Lobato and Adams (2015) J. M. Hernandez-Lobato and R. Adams Probabilistic backpropagation for scalable learning of bayesian neural networks. In Proceedings of the 32nd international conference on machine learning, F. Bach and D. Blei (Eds.), Proceedings of machine learning research, Vol. 37, Lille, France, p. 1861–1869. External Links: Link Cited by: §8.1. Hochreiter and Schmidhuber (1997) S. Hochreiter and J. Schmidhuber Long Short-Term Memory. Neural Comput. 9 (8), p. 1735–1780. External Links: ISSN 0899-7667, Link, Document Cited by: §E.2. Hoffman et al. (2013) M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley Stochastic Variational Inference. Journal of Machine Learning Research 14 (40), p. 1303–1347. External Links: Link Cited by: §2.2, §2.2, §2.2. Jaynes (2003) E.T. Jaynes Probability Theory: The Logic of Science. Cambridge University Press. External Links: Link Cited by: §2.2. M. I. Jordan and T. J. Sejnowski (Eds.) (2001) M. I. Jordan and T. J. Sejnowski (Eds.) Graphical models: foundations of neural computation. Computational neuroscience, MIT Press, Cambridge, Mass. External Links: ISBN 978-0-262-60042-2 Cited by: §1, §2.2. Khan and Nielsen (2018) M. E. Khan and D. Nielsen Fast yet Simple Natural-Gradient Descent for Variational Inference in Complex Models. arXiv:1807.04489 [cs, math, stat]. Note: arXiv: 1807.04489 External Links: Link Cited by: §2.2, §2.2, §6. Khan and Rue (2023) M. E. Khan and H. Rue The Bayesian learning rule. Journal of Machine Learning Research 24 (281), p. 1–46. Cited by: §2.2, §6. Khan (2025) M. E. Khan Information Geometry of Variational Bayes. Information Geometry 8 (S1), p. 275–289 (en). External Links: ISSN 2511-2481, 2511-249X, Link, Document Cited by: §1, §2.2. Khan and Lin (2017) M. Khan and W. Lin Conjugate-Computation Variational Inference : Converting Variational Inference in Non-Conjugate Models to Inferences in Conjugate Models. In Proceedings of the 20th International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research, Vol. 54, p. 878–887. External Links: Link Cited by: §B.2, Appendix B, §2.2, §5.1, §6, §6. Kingma and Ba (2015) D. Kingma and J. Ba Adam: A Method for Stochastic Optimization. In International Conference on Learning Representations (ICLR), San Diega, CA, USA. Cited by: §8.2. Knowles and Minka (2011) D. A. Knowles and T. P. Minka Non-conjugate Variational Message Passing for Multinomial and Binary Regression. In Advances in Neural Information Processing Systems, Vol. 24, Granada Spain (en). External Links: Link Cited by: Appendix D, §4, §6, §6, §7. Kucukelbir et al. (2017) A. Kucukelbir, D. Tran, R. Ranganath, A. Gelman, and D. M. Blei Automatic Differentiation Variational Inference. Journal of Machine Learning Research 18 (1), p. 430–474. External Links: Link Cited by: §2.2. Lin et al. (2018) W. Lin, N. Hubacher, and M. E. Khan Variational message passing with structured inference networks. In International conference on learning representations, External Links: Link Cited by: §B.2, §6. Linderman et al. (2025) S. W. Linderman, P. Chang, G. Harper-Donnelly, A. Kara, X. Li, G. Duran-Martin, and K. Murphy Dynamax: A Python package for probabilistic state space modeling with JAX. Journal of Open Source Software 10 (108), p. 7069. External Links: Link, Document Cited by: §9. Loeliger et al. (2007) H. Loeliger, J. Dauwels, J. Hu, S. Korl, L. Ping, and F. R. Kschischang The Factor Graph Approach to Model-Based Signal Processing. Proceedings of the IEEE 95 (6), p. 1295–1322. External Links: ISSN 0018-9219, Document Cited by: §2.1. Lounes et al. (2026) G. Lounes, R. Gerzaguet, and M. Gautier Control-flow aware MLIR tracing. In Proceedings of the 36th international workshop on rapid system prototyping, Rsp ’25, New York, NY, USA, p. 56–62. Note: Number of pages: 7 External Links: ISBN 979-8-4007-2224-0, Link, Document Cited by: §9. Lukashchuk et al. (2024) M. Lukashchuk, I. Senöz, and B. de Vries Q-conjugate message passing for efficient bayesian inference. In International conference on probabilistic graphical models, p. 295–311. Cited by: §6. Lukashchuk et al. (2026) M. Lukashchuk, K. Yemets, W. M. Kouw, D. Bagaev, İ. Şenöz, J. Beck, and B. de Vries Composing non-conjugate factor graphs with closed-form variational inference. arXiv preprint arXiv:2605.29467. Cited by: §E.2, Table 11, Table 11, Table 11, Table 11, Table 11, Table 11, §7, §8.2, §8.2, §8.2, Table 4, Table 4, Table 4, Table 4, §8. Mangion et al. (2011) A. Z. Mangion, K. Yuan, V. Kadirkamanathan, M. Niranjan, and G. Sanguinetti Online variational inference for state-space models with point-process observations. Neural computation 23 (8), p. 1967–1999. Cited by: §6. Minka (2001) T. P. Minka Expectation propagation for approximate Bayesian inference. In Proceedings of the seventeenth conference on uncertainty in artificial intelligence, UAI’01, Seattle, Washington, p. 362–369. Note: Number of pages: 8 tex.address: San Francisco, CA, USA External Links: ISBN 1-55860-800-1 Cited by: Figure 1, Figure 1, §6. Minka (2005) T. Minka Divergence Measures and Message Passing. Technical report Microsoft Research. Cited by: §6. Nuijten et al. (2024) W. W. L. Nuijten, D. Bagaev, and B. de Vries GraphPPL.jl: A Probabilistic Programming Language for Graphical Models. Entropy 26 (11). External Links: ISSN 1099-4300, Link, Document Cited by: §9. Nuijten et al. (2026a) W. W. L. Nuijten, M. Lukashchuk, T. van de Laar, and B. de Vries What Type of Inference is Active Inference?. External Links: Document Cited by: §9. Nuijten et al. (2026b) W. W. L. Nuijten, T. van de Laar, and B. de Vries Expected free energy-based planning as variational inference. Transactions on Machine Learning Research. External Links: ISSN 2835-8856 Cited by: §9. Opper and Archambeau (2009) M. Opper and C. Archambeau The variational gaussian approximation revisited. Neural Computation 21 (3), p. 786–792. Note: tex.eprint: https://direct.mit.edu/neco/article-pdf/21/3/786/820223/neco.2008.08-07-592.pdf External Links: ISSN 0899-7667, Link, Document Cited by: §6. Rue and Held (2005) H. Rue and L. Held Gaussian Markov random fields: theory and applications. Chapman and Hall/CRC. Cited by: §6. Senöz et al. (2021) I. Senöz, T. van de Laar, D. Bagaev, and B. de Vries Variational Message Passing and Local Constraint Manipulation in Factor Graphs. Entropy 23 (7), p. 807 (en). External Links: ISSN 1099-4300, Link, Document Cited by: Appendix B, §1, §2.1, §2.1, §3, §3, §3, §3, §3, §3, §6. Shen et al. (2024) Y. Shen, N. Daheim, B. Cong, P. Nickl, G. M. Marconi, B. C. E. M. Raoul, R. Yokota, I. Gurevych, D. Cremers, M. E. Khan, and T. Möllenhoff Variational learning is effective for large deep networks. In Proceedings of the 41st international conference on machine learning, R. Salakhutdinov, Z. Kolter, K. Heller, A. Weller, N. Oliver, J. Scarlett, and F. Berkenkamp (Eds.), Proceedings of machine learning research, Vol. 235, p. 44665–44686. External Links: Link Cited by: §8.1, §8.2. SILSO World Data Center (1945) SILSO World Data Center The international sunspot number. International Sunspot Number Monthly Bulletin and online catalogue. Cited by: §7. Tan (2025) L. S. L. Tan Analytic natural gradient updates for Cholesky factor in Gaussian variational approximation. Journal of the Royal Statistical Society Series B: Statistical Methodology 87 (4), p. 930–956. Note: tex.eprint: https://academic.oup.com/jrsssb/article-pdf/87/4/930/61670827/qkaf001.pdf External Links: ISSN 1369-7412, Link, Document Cited by: §C.1. Tschantz et al. (2025) A. Tschantz, M. Koudahl, H. Linander, L. Da Costa, C. Heins, J. Beck, and C. Buckley Bayesian predictive coding. arXiv preprint arXiv:2503.24016. Cited by: §8.1. Wainwright and Jordan (2008) M. J. Wainwright and M. I. Jordan Graphical Models, Exponential Families, and Variational Inference. Foundations and Trends® in Machine Learning 1 (1–2), p. 1–305 (English). External Links: ISSN 1935-8237, 1935-8245, Link, Document Cited by: §1, §2.2, §2.2, §2.2, §3. Wand (2017) M. P. Wand Fast Approximate Inference for Arbitrarily Large Semiparametric Regression Models via Message Passing. Journal of the American Statistical Association 112 (517), p. 137–168 (en). External Links: ISSN 0162-1459, 1537-274X, Link, Document Cited by: §6. Wang and Li (2020) Y. Wang and W. Li Information Newton’s flow: second-order optimization method in probability space. arXiv preprint arXiv:2001.04341. Cited by: Figure 3, Figure 3. Winn and Bishop (2005) J. Winn and C. M. Bishop Variational Message Passing. Journal of Machine Learning Research 6 (23), p. 661–694. External Links: ISSN 1533-7928, Link Cited by: §6. Wu et al. (2019) A. Wu, S. Nowozin, T. Meeds, R. E. Turner, J. M. Hernadez-Lobato, and A. L. Gaunt Deterministic variational inference for robust bayesian neural networks. In International conference on learning representations, External Links: Link Cited by: §8.1. Yedidia et al. (2001) J. S. Yedidia, W. T. Freeman, and Y. Weiss Bethe free energy, Kikuchi approximations, and belief propagation algorithms. Advances in neural information processing systems 13, p. 24 (en). Cited by: §1, §3, §3. Yedidia et al. (2005) J. S. Yedidia, W.T. Freeman, and Y. Weiss Constructing free-energy approximations and generalized belief propagation algorithms. IEEE Transactions on Information Theory 51 (7), p. 2282–2312. External Links: ISSN 0018-9448, Link, Document Cited by: §1, §3, §3, §3. Zeng et al. (2023) A. Zeng, M. Chen, L. Zhang, and Q. Xu Are transformers effective for time series forecasting?. In Proceedings of the AAAI conference on artificial intelligence, Vol. 37, p. 11121–11128. Cited by: §E.2. Zhou et al. (2021) H. Zhou, S. Zhang, J. Peng, S. Zhang, J. Li, H. Xiong, and W. Zhang Informer: Beyond Efficient Transformer for Long Sequence Time-Series Forecasting. In The Thirty-Fifth AAAI Conference on Artificial Intelligence, AAAI 2021, Virtual Conference, Vol. 35, p. 11106–11115. Cited by: §E.2. Appendix A Proof of Theorem 1 Proof of Theorem 1. We suppress additive constants, since all message multipliers are defined only up to normalization gauge. First consider the factor-side variation. The edge form constraint changes the edge variable but does not constrain the adjacent factor beliefs. Therefore the same variation with respect to qa(a)q_a( z_a) as in the Bethe Lagrangian gives qa∗(a)∝fa(a)∏j∈ℰ(a)exp(λaj(zj)).q_a^*( z_a)\; \;f_a( z_a) _j (a) \! ( _aj(z_j) ). (60) Now consider the edge-side variation. The only terms depending on i λ_i are collected in equation 32. The exponential-family identities are ∇λiAi∗(i(i))=Gi(i)i,∇λiqλi[r]=Covqλi[Ti,r]. _ _iA_i^*\! ( μ_i( λ_i) )=G_i( λ_i) λ_i, _ _iE_q_ _i[r]=Cov_q_ _i\! [T_i,r ]. (61) Applying equation 61 to equation 32 gives 0=−Gi(i)i+∑a∈(i)Covqλi[Ti,λai],0=-G_i( λ_i) λ_i+ _a (i)Cov_q_ _i [T_i, _ai ], which is exactly the edge stationarity condition equation 33. Multiplying by Gi(i)−1G_i( λ_i)^-1 gives the projected form i=∑a∈(i)iλi[λai]. λ_i= _a (i) T_i _i\! [ _ai ]. It remains to close the multipliers through the marginalization constraints. Substitute equation 60 into the constraint qλi(zi)=∫qa(a)da∖i.q_ _i(z_i)= q_a( z_a)\,d z_a i. For each a∈(i)a (i) this gives qλi(zi)∝exp(λai(zi))∫fa(a)∏j∈ℰ(a)∖iexp(λaj(zj))da∖i.q_ _i(z_i) \! ( _ai(z_i) ) f_a( z_a) _j (a) i \! ( _aj(z_j) )d z_a i. Using the definition equation 35, the same identity in log-coordinates is logqλi(zi)=λai(zi)+ℓa→i(zi)+const. q_ _i(z_i)= _ai(z_i)+ _a→ i(z_i)+const. (62) Apply the projection iλi T_i _i to both sides of equation 62. Constants vanish under covariance with TiT_i, and iλi[logqλi]=i T_i _i[ q_ _i]= λ_i because logqλi q_ _i is affine in TiT_i. Therefore i=iλi[λai]+iλi[ℓa→i]=iλi[λai]+ηa→i. λ_i= T_i _i\! [ _ai ]+ T_i _i\! [ _a→ i ]= T_i _i\! [ _ai ]+ _a→ i. (63) Finally use that (i)=b,cV(i)=\b,c\. From equation 63, iλi[λbi]=i−ηb→i,iλi[λci]=i−ηc→i. T_i _i\! [ _bi ]= λ_i- _b→ i, T_i _i\! [ _ci ]= λ_i- _c→ i. Substituting these two identities into the projected form of equation 33 gives i=(i−ηb→i)+(i−ηc→i), λ_i= ( λ_i- _b→ i )+ ( λ_i- _c→ i ), and hence i=ηb→i+ηc→i. λ_i= _b→ i+ _c→ i. This is equation 37. ∎ Appendix B Supplementary Derivation for Information-Geometric Message Passing This appendix gives a deliberately small example of inductive inference. The derivations below are not meant as a new algorithmic contribution; the same stationary equations are obtained from the conjugate-computation view of 30. The purpose here is narrower and more pedagogical. We use a graph with one latent variable and two unary factors, so there is no factorization approximation hiding in the notation. In this setting, the global conjugate-computation view and the local constrained-Bethe view coincide completely, which makes it a useful place to clarify language. The main point is to clarify which object is meant by a message once a marginal is constrained to a finite exponential family. The exact Lagrange multipliers of a constrained free-energy problem are functional objects, whereas the messages used by conjugate message-passing algorithms are finite exponential-family surrogates. This appendix shows how the latter can be computed from the former without claiming that the underlying non-conjugate factor has become conjugate. We follow the same variational logic as the constrained free-energy derivations of 15; 16 and the constrained Bethe-free-energy construction of 48. A model induces a positive unnormalized density; a free-energy functional ranks candidate posterior beliefs; and constraints restrict which beliefs are admissible. The calculus of variations then tells us what the stationary condition requires. Here, the graph has only one latent variable, so all of the bookkeeping is visible on the page. B.1 Exact Constrained Model Let the latent variable z be connected to two unary factors, f(z)=fc(z)fn(z),f(z)=f_c(z)\,f_n(z), (64) where fcf_c is conjugate to the chosen family ℰE, and fnf_n is not. We constrain the marginal belief on z to qλ∈ℰq_λ , with qλ(z)=exp⊤T(z)−A(),=qλ[T(z)].q_λ(z)= \ λ T(z)-A( λ)\, μ=E_q_λ[T(z)]. (65) Here and below, densities are written with respect to the carrier reference measure dν(z)=h(z)dzdν(z)=h(z)\,dz. With respect to the Lebesgue measure, the density in equation 65 is multiplied by h(z)h(z). The two factor beliefs qcq_c and qnq_n are left unrestricted. For this product graph, the local Bethe objective is ℱ(qc,qn,qλ)=∑a∈c,n∫qa(z)logqa(z)fa(z)ν(z)−∫qλ(z)logqλ(z)ν(z).F(q_c,q_n,q_λ)= _a∈\c,n\ q_a(z) q_a(z)f_a(z)\,dν(z)- q_λ(z) q_λ(z)\,dν(z). (66) The local-polytope constraints are qc(z)=qλ(z),qn(z)=qλ(z).q_c(z)=q_λ(z), q_n(z)=q_λ(z). (67) Introducing functional multipliers ρc _c and ρn _n for these constraints gives the Lagrangian ℒ =ℱ(qc,qn,qλ)+∑a∈c,n∫ρa(z)qλ(z)−qa(z)ν(z) =F(q_c,q_n,q_λ)+ _a∈\c,n\ _a(z)\q_λ(z)-q_a(z)\\,dν(z) +∑a∈c,nτa(∫qa(z)dν(z)−1)+τ(∫qλ(z)dν(z)−1). + _a∈\c,n\ _a ( q_a(z)\,dν(z)-1 )+τ ( q_λ(z)\,dν(z)-1 ). (68) The last normalization term is redundant when qλq_λ is parameterized as a normalized density, but keeping it visible makes the constrained variational problem explicit. Taking the variational derivative with respect to qaq_a gives δℒδqa(z)=logqa(z)fa(z)+1−ρa(z)+τa. δ q_a(z)= q_a(z)f_a(z)+1- _a(z)+ _a. (69) Setting equation 69 to zero yields qa(z)∝fa(z)expρa(z),a∈c,n.q_a(z)\; \;f_a(z) \ _a(z)\, a∈\c,n\. (70) Therefore, the exact marginalization constraints can be satisfied by ρa(z)=logqλ(z)−logfa(z)+const. _a(z)= q_λ(z)- f_a(z)+const. (71) This equation is a useful sanity check. Exact marginalization is not the obstruction: the functional multipliers can always absorb whatever is needed to make qa=qλq_a=q_λ. The obstruction is finite representation. If faf_a is non-conjugate, then equation 71 contains the negative non-conjugate residual of logfa f_a. Thus, the exact Lagrange multiplier is generally not an exponential-family message. The stationarity of the constrained marginal itself adds only a projected condition. Variation with respect to the finite parameter λ, holding the multipliers fixed, gives 0=−G()+Covqλ[T,ρc]+Covqλ[T,ρn].0=-G( λ) λ+Cov_q_λ\! [T, _c ]+Cov_q_λ\! [T, _n ]. (72) Premultiplying by G()−1G( λ)^-1 gives =Πλ[ρc]+Πλ[ρn], λ= ^λ[ _c]+ ^λ[ _n], (73) where Πλ[r]=∇μqλ[r(z)] ^λ[r]= _μE_q_λ[r(z)]. Thus , the stationary condition identifies only the components of the multipliers visible through the sufficient statistics T. The orthogonal functional components remain free to enforce the exact constraints. Substituting equation 71 into equation 73 and using Πλ[logqλ]= ^λ[ q_λ]= λ yields =Πλ[logfc]+Πλ[logfn]. λ= ^λ[ f_c]+ ^λ[ f_n]. (74) This is a stationary condition on the constrained marginal; it is not the density equality qλ(z)∝fc(z)fn(z)q_λ(z) f_c(z)f_n(z). This also fixes the terminology. The model factor faf_a is not a Lagrange multiplier, and the exact multiplier ρa _a is not the finite message passed by the algorithm. The ordinary exact BP log-message from a unary factor is ℓa(z)=logfa(z)+const. _a(z)= f_a(z)+const. (75) Consequently, equation 74 is the unary-factor specialization of equation 36 and equation 37: each mean-coordinate gradient is the natural parameter of its outgoing message. For a non-conjugate factor, this exact log-message is itself non-conjugate. The exact constrained Lagrange problem remains well-defined because the multiplier in equation 71 can cancel the non-conjugate part inside the unrestricted factor belief. A finite message-passing algorithm cannot communicate that arbitrary residual unless we replace the factor with a surrogate. B.2 Projected Surrogate Update The point of the surrogate construction is narrower than that of exact inference in a new model. For an arbitrary reference value 0 λ_0, the surrogate below does not claim to have the same marginal as the original constrained problem. Rather, it provides a finite-dimensional update map whose fixed points coincide with the stationary points of the original constrained problem. The objective to prove is therefore the fixed-point equivalence =Φ()⟺∇μℱexact()=0, λ= ( λ) _μF_exact( λ)=0, (76) not equality between the one-step surrogate marginal and the original constrained optimum for every choice of 0 λ_0. The finite-dimensional update is obtained by replacing the original product graph, locally and temporarily, with a conjugate surrogate product graph. Continue with the carrier-measure convention used above: write all densities with respect to dν(z)=h(z)dzdν(z)=h(z)\,dz. Thus, the exponential-family marginal has a log density ⊤T(z)−A() λ T(z)-A( λ) in this reference measure. Restoring Lebesgue densities simply multiplies the final marginal by h(z)h(z). At a reference marginal qλ0q_ _0, define the projected natural contribution of factor a by a(0)=Πλ0[ℓa]=∇μqλ[ℓa(z)]|λ=λ0. η_a( λ_0)= _0[ _a]= _μ\,E_q_λ[ _a(z)] |_λ= _0. (77) The surrogate graph has the same topology as the original two-factor graph, but its unary factors are conjugate sites m^a(z,0)∝expa(0)⊤T(z),a∈c,n. m_a(z; λ_0) \ η_a( λ_0) T(z)\, a∈\c,n\. (78) When faf_a is conjugate, this projection returns the usual natural-parameter contribution of the factor. When faf_a is non-conjugate, it returns only the component of ℓa _a visible to the sufficient statistics T under the reference marginal qλ0q_ _0. Now solve the exact constrained Bethe problem for this new surrogate graph. Introduce unrestricted surrogate factor beliefs q^c q_c and q^n q_n, keep the edge belief constrained to qλ∈ℰq_λ , and define ℱ^λ0(q^c,q^n,qλ)=∑a∈c,n∫q^a(z)logq^a(z)m^a(z,0)ν(z)−∫qλ(z)logqλ(z)ν(z). F_ _0( q_c, q_n,q_λ)= _a∈\c,n\ q_a(z) q_a(z) m_a(z; λ_0)\,dν(z)- q_λ(z) q_λ(z)\,dν(z). (79) The local constraints are q^c(z)=qλ(z),q^n(z)=qλ(z). q_c(z)=q_λ(z), q_n(z)=q_λ(z). (80) With functional multipliers ρ^c ρ_c and ρ^n ρ_n, the surrogate Lagrangian is ℒ^λ0 L_ _0 =ℱ^λ0+∑a∈c,n∫ρ^a(z)qλ(z)−q^a(z)ν(z) = F_ _0+ _a∈\c,n\ ρ_a(z)\q_λ(z)- q_a(z)\\,dν(z) +∑a∈c,nτ^a(∫q^a(z)dν(z)−1)+τ^(∫qλ(z)dν(z)−1). + _a∈\c,n\ τ_a ( q_a(z)\,dν(z)-1 )+ τ ( q_λ(z)\,dν(z)-1 ). (81) Taking the variational derivative with respect to q^a q_a gives δℒ^λ0δq^a(z)=logq^a(z)m^a(z,0)+1−ρ^a(z)+τ^a. δ L_ _0δ q_a(z)= q_a(z) m_a(z; λ_0)+1- ρ_a(z)+ τ_a. (82) Setting equation 82 to zero yields q^a(z)∝m^a(z,0)expρ^a(z). q_a(z)\; \; m_a(z; λ_0) \ ρ_a(z)\. (83) Enforcing the marginal constraints therefore identifies the exact surrogate multipliers as ρ^a(z)=logqλ(z)−logm^a(z,0)+const. ρ_a(z)= q_λ(z)- m_a(z; λ_0)+const. (84) Unlike equation 71, these multipliers are finite exponential-family objects. Indeed, using equation 78, ρ^a(z)=(−a(0))⊤T(z)+const. ρ_a(z)= ( λ- η_a( λ_0) ) T(z)+const. (85) It remains to take stationarity with respect to the constrained marginal parameter λ. Using ∇λA∗(())=G() _λA^*( μ( λ))=G( λ) λ and ∇λqλ[r]=Covqλ[T,r] _λE_q_λ[r]=Cov_q_λ[T,r] for fixed r, the derivative of equation 81 with respect to λ is 0=−G()+∑a∈c,nCovqλ[T,ρ^a].0=-G( λ) λ+ _a∈\c,n\Cov_q_λ\! [T, ρ_a ]. (86) Premultiplying by G()−1G( λ)^-1 gives =Πλ[ρ^c]+Πλ[ρ^n]. λ= ^λ[ ρ_c]+ ^λ[ ρ_n]. (87) Substituting equation 84, and using Πλ[logqλ]= ^λ[ q_λ]= λ together with Πλ[logm^a(⋅,0)]=a(0) ^λ[ m_a(·; λ_0)]= η_a( λ_0), yields ⋆(0)=c(0)+n(0). λ ( λ_0)= η_c( λ_0)+ η_n( λ_0). (88) Thus the exact solution of the surrogate constrained problem has natural parameter ⋆(0) λ ( λ_0) and density qλ⋆(λ0)(z)=exp(c(0)+n(0))⊤T(z)−A(c(0)+n(0))w.r.t. ν.q_λ ( _0)(z)= \! \ ( η_c( λ_0)+ η_n( λ_0) ) T(z)-A ( η_c( λ_0)+ η_n( λ_0) ) \ .r.t. ν. (89) Equivalently, with respect to Lebesgue measure the same density is multiplied by h(z)h(z). Equation equation 88 is the inner exact message-passing solution on the surrogate graph. There is no additional Lagrange multiplier enforcing 0= λ_0= λ in equation 81. The parameter 0 λ_0 defines the surrogate factors and is held fixed while the inner constrained problem is solved. In other words, the inner problem defines an update map Φ(0)≔⋆(0)=c(0)+n(0). ( λ_0) λ ( λ_0)= η_c( λ_0)+ η_n( λ_0). (90) The projected-message fixed point is the self-consistency condition that the reference marginal used to build the surrogate is also the marginal returned by that surrogate, that is, Φ()= ( λ)= λ: =c()+n()=Πλ[logfc]+Πλ[logfn]. λ= η_c( λ)+ η_n( λ)= ^λ[ f_c]+ ^λ[ f_n]. (91) One could introduce a separate copy of the reference parameter and a constraint forcing that copy to equal λ, but this would either be redundant or change the problem. It is redundant if the surrogate sites are already frozen before the inner optimization. It changes the stationarity equations if one differentiates through the dependence a(0)=Πλ0[ℓa] η_a( λ_0)= _0[ _a] because then extra ∂a/∂0∂ η_a/∂ λ_0 terms enter. The message-passing algorithm uses the first interpretation: build a local surrogate at the current marginal, solve the surrogate problem exactly, and seek a fixed point of the resulting map. This fixed-point condition is exactly the first-order condition of the original constrained problem, not a new exactness claim for every intermediate surrogate. Indeed, after enforcing the exact constraints qc=qn=qλq_c=q_n=q_λ, the original constrained objective reduces, up to constants, to ℱexact()=A∗(())−qλ[logfc(z)+logfn(z)].F_exact( λ)=A^*( μ( λ))-E_q_λ[ f_c(z)+ f_n(z)]. (92) Its derivative in mean coordinates is ∇μℱexact=−Πλ[logfc]−Πλ[logfn]=−Φ(). _μF_exact= λ- ^λ[ f_c]- ^λ[ f_n]= λ- ( λ). (93) Thus, the outer iteration t+1=Φ(t) λ^t+1= ( λ^t) is a fixed-point method, equivalently a unit natural-gradient step, for solving the exact constrained stationarity condition. A converged fixed point is therefore a stationary marginal of the original constrained problem. This does not imply that the intermediate surrogate marginals are exact marginals of the original model, nor that the fixed point is the global minimizer when the reduced objective is non-convex. At such a fixed point, the surrogate multipliers are explicit ρ^c(z)=n()⊤T(z)+const,ρ^n(z)=c()⊤T(z)+const. ρ_c(z)= η_n( λ) T(z)+const, ρ_n(z)= η_c( λ) T(z)+const. (94) This is the payoff of the surrogate construction. At a fixed point, the surrogate graph has the same constrained marginal as the stationary point of the original constrained problem, but its half-edge multipliers are now finite-dimensional exponential-family objects. Equivalently, the exponentiated multipliers are computable messages in the chosen sufficient statistics rather than arbitrary functional objects. This matters in message passing because a local solution is not enough: a factor must also propagate a message to the rest of the graph. The exact constrained problem can satisfy marginalization by using the functional ratios in equation 71, but those ratios are not generally messages that a finite-dimensional algorithm can pass. The surrogate replaces that functional object with a conjugate site whose exact Lagrangian multipliers lie in the same finite family as the constrained marginal. In Gaussian-chain examples with non-conjugate observation factors, this is the mechanism that turns a non-Gaussian likelihood contribution into a Normal site message while preserving the fixed-point marginal condition. Thus equation 94 is the exact marginalization condition for the surrogate graph, not a claim that the original non-conjugate factor has become conjugate. The exact Lagrange multipliers of the original constrained problem remain the functional ratios in equation 71. The finite message in equation 78 is a different object: it is the conjugate surrogate factor whose constrained Bethe Lagrangian has finite-dimensional multipliers and an exact exponential-family marginal. This is also the convention used by conjugate-computation methods. In that literature, the word “message” refers to the conjugate surrogate obtained from a mean-parameter or natural-gradient projection (30; 34). The point of the local information-geometric derivation is to obtain the same surrogate message from the edge-local stationary condition while keeping separate the exact functional multipliers that solve the constrained Lagrangian in the background. B.3 Computing the Tangent Projection in Practice The projection Πλ[ℓ]=∇μqλ[ℓ(z)] ^λ[ ]= _μE_q_λ[ (z)] of equation 77 reduces, via ∇λqλ[ℓ]=Covqλ[T,ℓ] _λE_q_λ[ ]=Cov_q_λ[T, ] and preconditioning by G()−1G( λ)^-1, to covariances of the sufficient statistics with the log-message under the current marginal. Whether the update is closed-form therefore depends only on whether these expectations are available analytically for the factor and family at hand. The ablations of Section 7 cover the three cases that occur in practice. Closed form: Poisson state-space model. Toward a Gaussian state edge, the observation log-factor is ℓ(z)=yz−ez (z)=yz-e^z up to constants, and Gaussian expectations of z, z2z^2, and eze^z—hence all required covariances—are analytic. NGMP’s projection is a single closed-form natural-gradient step per edge, and the PVMP baseline iterates the same closed-form gradients inside its inner manifold optimizer, so both inference methods are closed-form in this model. The budget-matched NCVMP control of Appendix D uses the same closed-form gradient but takes a single inner step. Unscented approximation: Normal mean–precision model. The exact cavity log-message toward the precision edge is not a linear combination of the Gamma sufficient statistics (logτ,τ)( τ,τ), and its Gamma expectations have no closed form. The implementation evaluates the tangent projection with a deterministic unscented (sigma-point) approximation of the required expectations under the current marginal. VMP, by contrast, needs no projection in this model: the mean-field coordinate updates are conjugate and closed-form. Quadrature: sequential heteroskedastic model. The cavity messages through the likelihood factor of equation 57 have no closed-form expectations under the joint (μ,s)(μ,s) cluster marginal; the implementation evaluates the tangent projection by one-dimensional numerical quadrature with 3232 nodes over the log-precision marginal, while the PVMP baseline again iterates closed-form gradients of its mean-field sites, and the NCVMP control again takes a single such step. In every case, the projection consumes only expectations under the current marginal, so the accuracy of the numerical variants is governed by the chosen sigma-point or quadrature rule rather than by sampling noise. Appendix C Damping and Momentum The surrogate construction of Section 5 defines the outer fixed-point map Φ() ( λ) in equation 54: build projected surrogate leaves at the current constrained marginals, freeze the surrogate graph, and run the conjugate BP sweep. In mean coordinates, the exact reduced constrained-Bethe objective has the same stationarity form as the global exponential-family free energy: ∇μFexact()=−Φ(). _μF_exact( λ)= λ- ( λ). (95) Thus, the undamped update (t+1)=Φ((t)) λ^(t+1)= ( λ^(t)) is a unit natural-gradient or fixed-point step. Damping changes only the step length in natural coordinates, (t+1)=(1−α)(t)+αΦ((t)),0<α≤1. λ^(t+1)=(1-α) λ^(t)+α\, ( λ^(t)), 0<α≤ 1. (96) Equivalently, for an individual surrogate message with natural parameter ⋆ η computed by projection, the damped message parameter is (t)=(1−α)(t−1)+α⋆. η^(t)=(1-α) η^(t-1)+α\, η . (97) For the Gaussian Poisson leaf, this means averaging the canonical coordinates (ξ,τ)(ξ,τ) before converting back to pseudo-observation parameters (y~,Λ)( y, ). In message-function form μ^(t)∝(μ^(t−1))1−α(μ^⋆)α. μ^(t)\; \; ( μ^(t-1) )^1-α ( μ )^α. (98) Whether damping is needed depends on the data. On the sunspot masks of Section 7, the undamped iteration (α=1α=1), natural damping with α=0.25α=0.25, and the damped heavy-ball setting used for the sunspot results (α=0.5α=0.5, β=0.2β=0.2; equation 99 below) all converge to the same fixed point by iterations 5, 14, and 8 respectively: the counts are moderate, and the initialization log(yk+1) (y_k+1) is close to the fixed point, so damping only costs iterations. The fixed-point map stops contracting when the latent log-rate makes long excursions to extreme values; for example, a stretch of zero counts at a very negative log-rate, where the natural-gradient message reacts strongly to the marginal variance. Figure 14 uses data simulated from the Poisson state-space model equation 39 itself (z0=0z_0=0, σ2=0.1σ^2=0.1 as in the sunspot experiments, every count observed): the latent random walk makes larger excursions and produces longer stretches of zero counts as the chain length N grows. At N=100N=100 , every setting converges. From N=250N=250 on, the undamped map settles into a period-two oscillation on a growing share of seeds (3, 9, and 9 of 20 for N=250N=250, 500500, 10001000) and at N=1000N=1000 diverges on three of them; damping with α=0.25α=0.25 converges on all 80 fits at the price of more iterations; the heavy-ball setting is the fastest whenever it converges, but its momentum step can produce a message with negative precision, which happens on 3 of 20 seeds at N=500N=500 and 9 of 20 at N=1000N=1000 (see Section C.1). These are the conditions under which damping, and in harder cases momentum, are necessary, and why they are part of the practical outer solver rather than a cosmetic post-processing step. Figure 14: Convergence of the NGMP fixed-point iteration on the Poisson state-space model equation 39 for series simulated from the same model (σ2=0.1σ^2=0.1, all counts observed) of length N=100N=100, 250250 (top row), 500500, and 10001000 (bottom row), 20 seeds each, 200 variational iterations. BFE on the vertical axis is the Bethe free energy of the model per observation minus the value reached at convergence by the α=0.25α=0.25 setting on the same seed, so every seed converges to zero (log scale); curves are the geometric mean across seeds with a 95% confidence band. Orange: undamped fixed-point iteration (α=1α=1). Blue: messages damped in natural coordinates with α=0.25α=0.25. Green dashed: the damped heavy-ball setting used for the sunspot experiments (α=0.5α=0.5, β=0.2β=0.2); its momentum step left the natural domain and the fit stopped on 3 of 20 seeds at N=500N=500 and 9 of 20 at N=1000N=1000 (the undamped iteration failed on one seed at N=1000N=1000), and the curves average over the seeds that reached each iteration. Longer chains make larger excursions of the latent log-rate; the undamped map then stops converging while the damped map still does. Momentum can be applied to the same natural-parameter residual Φ()−=−∇μFexact() ( λ)- λ=- _μF_exact( λ). A simple heavy-ball form is (t+1)=β(t)+α(Φ((t))−(t)),(t+1)=(t)+(t+1),0≤β<1. aligned v^(t+1)&=β v^(t)+α ( ( λ^(t))- λ^(t) ),\\ λ^(t+1)&= λ^(t)+ v^(t+1), aligned 0≤β<1. (99) Other accelerated fixed-point variants differ in how they choose the point at which Φ is evaluated, but the coordinate principle is the same: damping and momentum operate on natural/message parameters, not on derived pseudo-observation parameters. When the iteration converges, the residual Φ()− ( λ)- λ is zero, so damping and momentum change the convergence path but not the stationary equations of Theorem 1. C.1 Vector-transport momentum The heavy-ball rule equation 99 adds the previous momentum vector to a residual computed at a different iterate, as if both lived in the same vector space. An exponential family is a smooth manifold on which the Fisher information defines a Riemannian metric. Therefore, the Euclidean rule generalizes directly: the momentum vector (t) v^(t) is a tangent vector at the iterate where it was formed, and before it is combined with the new residual, it is carried to the current iterate by a vector transport ⋅←⋅T_·←·. (t+1) v^(t+1) =β(t)←(t−1)(t)+α(Φ((t))−(t)), =β\,T_ λ^(t)← λ^(t-1) v^(t)+α ( ( λ^(t))- λ^(t) ), (100) (t+1) λ^(t+1) =(t)+(t+1). = λ^(t)+ v^(t+1). A full account of retractions and vector transports is outside the scope of this paper; we refer to 1 for the theory and to the momentum gradient rule of Manopt.jl (9) for the algorithmic template that equation 100 follows. The transport is approximated elementwise: the previous momentum vector is rescaled by the square root of the ratio of the diagonal Fisher metrics at the previous and current natural parameters. The resulting step is capped at a maximum Euclidean norm; the settings used in the regression experiments are reported in Section E.1. With β>0β>0, the update is no longer a convex combination of natural parameters and can leave the natural domain, for example, a Gamma message with a non-positive rate, so on such edges, we use damping alone (β=0β=0). The motivation for momentum, beyond the stability argument of Figure 14, comes from 51, who shows that natural-gradient updates for Gaussian variational approximations benefit substantially from momentum. The same holds here: in the depth-three hierarchies of Section 8.1, damping alone stalls within the sweep budget, and the vector-transport momentum update is what converges. Appendix D Convergence Diagnostics and Additional Results for the Comparison Study This appendix collects the convergence diagnostics behind the fixed-point claims of Section 7, together with the full predictive comparison for the sequential heteroscedastic study. In both non-conjugate ablations PVMP evaluates its own mean-field variational objective, so its trace is a true Bethe free energy on the fitted graph. NGMP rebuilds its local Gaussian surrogates between outer sweeps, so its trace scores the node energies through the current surrogate marginals: it is a convergence diagnostic whose level is not directly comparable to PVMP’s, because the two curves approximate different Bethe functionals. What the traces establish is that every comparison in Section 7 probes fixed points rather than truncated optimization. Poisson state-space model. Figure 15 shows the Bethe free energy per observed count on the sunspot series at all four holdout fractions. Both methods reach their plateaus well within the 20-sweep budget: NGMP by sweep three to four, PVMP shortly after (at 50%50\% holdout in roughly ten). The held-out separation at 50%50\% in Table 1 is therefore a property of the fixed points, not of stopping early. One configuration detail of the projective update of equation 56 matters for this model. With the default gradient-norm bound of its inner manifold optimizer, a marginal that drifts far into the tail becomes trapped: the projection returns its own input, neither damping nor additional sweeps escape, and rare masks blow up catastrophically — mask-average held-out negative log-likelihood up to 6666, driven by a single month assigned a rate of ≈104≈10^4 where the exact tilted update at the same state returns ≈log114≈ 114. Widening the bound from 11 to 100100 frees the trap and is used throughout Section 7. NGMP’s projection on these edges is closed-form and has no inner optimizer to tune. Figure 15: Bethe free energy per observed count on the sunspot series (mean ± 95% CI across 20 masks, log scale) at 5%5\% and 10%10\% holdout (top row) and 20%20\% and 50%50\% holdout (bottom row). Both methods plateau within the 20-sweep budget: NGMP by sweep three to four, PVMP shortly after. Each PVMP trace evaluates its own variational objective; the NGMP trace is a surrogate diagnostic, since the local Gaussian surrogates are rebuilt between outer sweeps. The dotted trace is the budget-matched NCVMP control of Table 5. Sequential heteroscedastic filtering. Figure 16 shows the corresponding traces for the filtering study. At the end of the shared 240240-sweep budget, PVMP’s Bethe free energy decreases by 5.3×10−45.3× 10^-4 nats per observation in the final full-batch sweep and by 3.8×10−43.8× 10^-4 in the first sequential batch. The final-sweep changes in batches 22–1010 are below 7×10−87× 10^-8: the collapsed prior of Figure 17 leaves batches 22–1010 essentially converged at initialization, so additional variational iteration cannot change the filtering outcome. NGMP’s surrogate diagnostic visibly plateaus; its final-sweep changes range from 5.6×10−65.6× 10^-6 to 3.3×10−53.3× 10^-5 nats per observation across the full and sequential fits. Figure 18 shows the posterior predictive bands of all four fits on identical axes. The two full-batch fits and NGMP’s sequential fit are visually indistinguishable from one another; only PVMP’s sequential fit differs, which is why the main text shows the sequential pair alone. Figure 19 repeats the comparison for the predictive-variance decomposition: NGMP’s full-batch and sequential fits coincide visually, so the main text omits the full-batch NGMP panel. Figure 16: Bethe free energy per observation for the sequential heteroscedastic model (mean ± 95% CI over 20 seeds). Left: the full-batch fit against sweep (first two sweeps omitted for scale). Right: the final-sweep value of each filtering update against batch index. PVMP evaluates its own mean-field variational objective; the NGMP trace scores the same node energies through moment-matched joint (μ,s)(μ,s) cluster marginals, so the two curves approximate different Bethe functionals and their levels are not directly comparable. The dotted trace is the budget-matched NCVMP control of Table 6. Figure 17: Noise-weight concentration logdetΣw _w along the filtering chain (mean ± 95% CI over 20 seeds; dotted lines show full-batch values; lower means more certain). PVMP’s sequential endpoint coincides with its full-batch value, and the budget-matched NCVMP control (Table 6) traces the identical curve; NGMP concentrates as evidence accumulates and ends near its full-batch level. (a) PVMP, full batch (b) PVMP, sequential batches (c) NGMP, full batch (d) NGMP, sequential batches Figure 18: Posterior predictive bands of all four fits on one representative seed with identical axes. Orange dashed curves are PVMP and blue solid curves are NGMP; the matching shaded regions are pointwise 95% posterior-predictive bands. Gray dash-dotted curves show the true mean, and gray points are the 400 training observations. Panels (a) and (c) use all observations at once, whereas panels (b) and (d) use the ten-step filtering chain of Figure 9. Panels (a), (c), and (d) are visually indistinguishable; only PVMP’s sequential fit in panel (b) collapses. (a) PVMP, full batch (b) PVMP, sequential batches (c) NGMP, full batch (d) NGMP, sequential batches Figure 19: Posterior predictive variance of all four fits on the same representative seed as Figure 18. Panels (a) and (c) use all observations at once, whereas panels (b) and (d) use ten sequential batches. The orange long-dashed curves in panels (a)–(b) are PVMP and the blue solid curves in panels (c)–(d) are NGMP. Each is the posterior mean total predictive variance, S¯(x)=Vepi(x)+Valea(x) S(x)=V_epi(x)+V_alea(x), where Vepi(x)=ϕ(x)⊤Σvϕ(x)V_epi(x)=φ(x) _vφ(x) is epistemic variance from the uncertain mean weights and Valea(x)=q[e−s(x)]V_alea(x)=E_q[e^-s(x)] is learned aleatoric variance. Shading gives a pointwise 95% posterior credible interval for the total. The purple dotted curve shows Vepi(x)V_epi(x) alone, so its vertical gap to the orange or blue total is Valea(x)V_alea(x). The gray dash-dotted curve is the benchmark’s true aleatoric variance [0.45(x+0.5)]2[0.45(x+0.5)]^2. All panels use identical logarithmic axes. Budget-matched control: NCVMP. The runtime of the two methods is dominated by different quantities: NGMP takes one natural-gradient step per edge update, while PVMP’s projective update runs an inner manifold optimizer on every edge — up to 100100 gradient evaluations per update in our configuration. The control reported here matches the budgets from PVMP’s side: same graph, constraints, initialization, and outer sweep schedule, but every ProjectedTo call is limited to a single inner iteration, so each edge update costs one gradient evaluation, as in NGMP. A single natural-gradient step of the projective update at the current marginal is the non-conjugate VMP update of 32, so this control is NCVMP (Section 6). NCVMP degrades predictive likelihood at the matched runtime. On the sunspot model (Table 5) its held-out negative log-likelihood is 12.512.5–16.116.1 against 4.54.5–5.45.4 for both converged methods, at every holdout fraction. On the sequential heteroscedastic model (Table 6), its full-batch NLL is 0.9130.913 against 0.2760.276 for PVMP and 0.2350.235 for NGMP, while its sequential NLL is 0.9160.916 against 0.5120.512 and 0.2440.244, respectively; its noise-weight posterior nevertheless collapses along converged PVMP’s exact concentration path (Figure 17), so the overconcentration of Section 7 is not an artifact of over-optimizing the projections, weaker per-edge optimization neither repairs the calibration nor preserves the fit. The NCVMP Bethe traces are the dotted curves in Figures 15 and 16. Table 5: Sunspot state-space model with the budget-matched control (mean ± 95% CI over 20 random masks): the two converged methods of Table 1 alongside NCVMP, that is, PVMP with single-step projections. NLL RMSE Held out PVMP NCVMP NGMP PVMP NCVMP NGMP 5% 4.534 ± 0.087 12.492 ± 0.546 4.532 ± 0.085 13.918 ± 0.385 26.441 ± 0.951 13.935 ± 0.378 10% 4.569 ± 0.058 12.998 ± 0.453 4.566 ± 0.058 13.972 ± 0.260 27.137 ± 0.768 13.995 ± 0.260 20% 4.672 ± 0.039 13.547 ± 0.410 4.661 ± 0.037 14.429 ± 0.214 27.616 ± 0.549 14.446 ± 0.208 50% 5.411 ± 0.162 16.116 ± 0.391 4.930 ± 0.035 16.839 ± 0.471 30.729 ± 0.507 15.789 ± 0.162 Table 6: Sequential heteroscedastic model with the budget-matched NCVMP control (mean ± 95% CI over 20 paired seeds), extending Table 2. full batch sequential NLL RMSE NLL RMSE batching penalty PVMP 0.276±0.0310.276± 0.031 0.588±0.0370.588± 0.037 0.512±0.0440.512± 0.044 0.583±0.0390.583± 0.039 0.235±0.0300.235± 0.030 NCVMP 0.913±0.0150.913± 0.015 0.577±0.0350.577± 0.035 0.916±0.0170.916± 0.017 0.576±0.0350.576± 0.035 0.003±0.0050.003± 0.005 NGMP 0.235±0.0310.235± 0.031 0.591±0.0380.591± 0.038 0.244±0.0310.244± 0.031 0.607±0.0420.607± 0.042 0.009±0.0100.009± 0.010 Appendix E Experiment details and additional results E.1 Regression details Table 7 reports the dimensions before the repeated 90/10 train–test splits. The feature count is the number of raw input variables before standardization and random-feature construction. Table 7: UCI regression data sets used in Section 8.1. The sample count is the complete data set before each repeated holdout. Data set Samples Input features Concrete 1,030 8 Energy 768 8 Boston Housing 506 13 Power Plant 9,568 4 Wine Quality Red 1,599 11 Yacht 308 6 Table 8 separates two numerical notions of capacity. P is the number of predictive weights and biases, whereas U is the number of additional independent scalars required to represent global posterior uncertainty beyond its mean. Both exclude fixed random features, optimizer state, and observation-local messages. Here d is the raw feature count from Table 7. Table 8: Numerical model and uncertainty capacity in the UCI comparison. NGMP uses the hierarchy in Figure 12 with L=2L=2: one mean-weight vector and two precision-weight vectors. Each vector has 1,000 random-feature coefficients, d direct-input coefficients, and one intercept. All ranges are exact over the six data sets. Method Architecture Model coefficients P Posterior / propagation Uncertainty parameters U Noise hoBBB d→1d\!→\!50\!→\!50\!→\!1 2,851–3,301 Diag.-G 2,851–3,301 fixed B d→2d\!→\!50\!→\!50\!→\!2 2,902–3,352 Diag.-G 2,902–3,352 input hodDVI d→1d\!→\!50\!→\!1 301–751 Diag.-G / 50 301–751 fixed dDVI d→2d\!→\!50\!→\!2 352–802 Diag.-G / 50 352–802 input hoDVI d→1d\!→\!50\!→\!1 301–751 Diag.-G / 1,275 301–751 fixed DVI d→2d\!→\!50\!→\!2 352–802 Diag.-G / 1,275 352–802 input hoBPC d→1d\!→\!50\!→\!50\!→\!1 2,851–3,301 MN–W 5,221–5,311 global hoIVON d→1d\!→\!50\!→\!50\!→\!1 2,851–3,301 Diag.-G 2,851–3,301 fixed NGMP Hierarchy, L=2L=2 Figure 12 3,015–3,042 Full-G / layer 1,516,545–1,543,815 2-level input Diag.-G and Full-G denote diagonal- and full-covariance Gaussian weight beliefs; MN–W denotes Matrix-Normal–Wishart. Exactly, U=PU=P for Diag.-G, U=((d+1)(d+2))/2+5,206U=((d+1)(d+2))/2+5,206 for MN–W, and U=3(d+1,001)(d+1,002)/2U=3(d+1,001)(d+1,002)/2 for NGMP. For DVI/dDVI, the number after the slash is the propagated hidden second-moment state per observation: 50 marginal variances for dDVI versus 1,275 unique covariance entries for DVI. Each of the feature maps (ϕφ for the mean layer, ψℓ _ for the precision layers) contains 1,000 random Fourier features approximating a Matérn-3/23/2 kernel, concatenated with the standardized input and an intercept. The frequencies are divided equally among three spectral scales (0.5,1,2)(0.5,1,2), giving a multiscale representation. The layer length scales decrease geometrically as (1.5,1.5/2,0.75)(1.5,1.5/ 2,0.75) from the mean layer to the deepest precision layer. The feature maps are sampled once per train–test split and then held fixed. Gaussian weight priors are used at all layers. The first precision-layer intercept is initialized from the homoscedastic residual variance, clipped to [10−3,103][10^-3,10^3], and the upper carrier precision is fixed to τtop=25 _top=25. We use at most 60 sweeps, stop when the relative change of the concatenated weight means falls below 10−510^-5, and stabilize linear solves with diagonal jitter 10−810^-8. Natural parameters are updated with the vector-transport momentum rule of Section C.1, with step parameter α=0.6α=0.6, momentum parameter β=0.8β=0.8, and maximum step 0.50.5. The ho prefix in Tables 9 and 10 marks a homoscedastic baseline. B, DVI, and IVON use fixed unit observation variance in standardized target units; BPC instead learns one global output covariance through its final Matrix-Normal–Wishart layer. The heteroscedastic B and DVI variants add a clipped input-dependent log-variance head. The hoIVON predictive is evaluated as a mixture over K=20K=20 posterior network draws. Table 9: Full test negative log likelihood (NLL) in original target units for all baseline variants; lower is better. Values are means ± approximate 95% confidence-interval half-widths over 20 paired splits. Bold values include the best point estimate within their confidence interval. Method Concrete Energy Boston Power Wine Yacht hoBBB 3.8613±0.00593.8613± 0.0059 3.2892±0.00433.2892± 0.0043 3.2488±0.01573.2488± 0.0157 3.7960±0.00193.7960± 0.0019 1.0336±0.01731.0336± 0.0173 3.6925±0.00593.6925± 0.0059 B 3.3328±0.02153.3328± 0.0215 2.4229±0.03502.4229± 0.0350 2.6555±0.03852.6555± 0.0385 2.7785±0.01432.7785± 0.0143 0.9438±0.0232 0.9438± 0.0232 2.6565±0.04512.6565± 0.0451 hodDVI 3.8602±0.00663.8602± 0.0066 3.2970±0.00273.2970± 0.0027 3.2571±0.01113.2571± 0.0111 3.7918±0.00113.7918± 0.0011 1.0360±0.01641.0360± 0.0164 3.7236±0.00973.7236± 0.0097 dDVI 3.0477±0.04323.0477± 0.0432 1.1920±0.4625 1.1920± 0.4625 2.4646±0.0921 2.4646± 0.0921 2.8238±0.01832.8238± 0.0183 0.9440±0.0314 0.9440± 0.0314 0.4641±0.10000.4641± 0.1000 hoDVI 3.8602±0.00663.8602± 0.0066 3.2971±0.00273.2971± 0.0027 3.2572±0.01113.2572± 0.0111 3.7918±0.00113.7918± 0.0011 1.0360±0.01641.0360± 0.0164 3.7238±0.00973.7238± 0.0097 DVI 3.0406±0.04383.0406± 0.0438 1.0094±0.2973 1.0094± 0.2973 2.4489±0.0917 2.4489± 0.0917 2.8249±0.02002.8249± 0.0200 0.9448±0.0307 0.9448± 0.0307 0.4543±0.10960.4543± 0.1096 hoBPC 3.9874±0.02023.9874± 0.0202 2.6357±0.05262.6357± 0.0526 2.9819±0.08622.9819± 0.0862 3.5396±0.01593.5396± 0.0159 1.0545±0.02561.0545± 0.0256 3.0251±0.07123.0251± 0.0712 hoIVON (K=20K=20) 3.7966±0.00823.7966± 0.0082 3.2693±0.04123.2693± 0.0412 3.3284±0.08553.3284± 0.0855 3.7839±0.00133.7839± 0.0013 1.0430±0.01561.0430± 0.0156 3.6584±0.01013.6584± 0.0101 NGMP 2.9845±0.0433 2.9845± 0.0433 0.9783±0.0655 0.9783± 0.0655 2.6720±0.07522.6720± 0.0752 2.7610±0.0226 2.7610± 0.0226 0.9571±0.0293 0.9571± 0.0293 0.2782±0.1437 0.2782± 0.1437 Table 10: Test root mean squared error (RMSE) in original target units for all baseline variants; lower is better. Values are means ± approximate 95% confidence-interval half-widths over 20 paired splits. Bold values include the best point estimate within their confidence interval. Method Concrete Energy Boston Power Wine Yacht hoBBB 7.2822±0.20727.2822± 0.2072 2.8175±0.14592.8175± 0.1459 3.8392±0.4243 3.8392± 0.4243 4.2831±0.07034.2831± 0.0703 0.6469±0.0185 0.6469± 0.0185 3.4875±0.40633.4875± 0.4063 B 6.7086±0.31176.7086± 0.3117 2.7537±0.15482.7537± 0.1548 3.7959±0.4387 3.7959± 0.4387 4.0251±0.0737 4.0251± 0.0737 0.6470±0.0176 0.6470± 0.0176 3.0494±0.41503.0494± 0.4150 hodDVI 6.7179±0.27646.7179± 0.2764 2.2486±0.11592.2486± 0.1159 3.4801±0.3374 3.4801± 0.3374 4.3009±0.07194.3009± 0.0719 0.6484±0.0190 0.6484± 0.0190 2.1779±0.23012.1779± 0.2301 dDVI 5.6591±0.3352 5.6591± 0.3352 2.1088±0.30882.1088± 0.3088 3.4741±0.4119 3.4741± 0.4119 4.1321±0.07634.1321± 0.0763 0.6463±0.0198 0.6463± 0.0198 0.8491±0.1387 0.8491± 0.1387 hoDVI 6.7179±0.27666.7179± 0.2766 2.2493±0.11612.2493± 0.1161 3.4802±0.3375 3.4802± 0.3375 4.3009±0.07204.3009± 0.0720 0.6484±0.0190 0.6484± 0.0190 2.1787±0.23012.1787± 0.2301 DVI 5.6654±0.3249 5.6654± 0.3249 1.9434±0.33531.9434± 0.3353 3.4291±0.4013 3.4291± 0.4013 4.1428±0.08014.1428± 0.0801 0.6469±0.0197 0.6469± 0.0197 0.8528±0.1533 0.8528± 0.1533 hoBPC 12.9488±0.267112.9488± 0.2671 3.3192±0.20533.3192± 0.2053 4.7579±0.50244.7579± 0.5024 8.2997±0.11628.2997± 0.1162 0.6971±0.01820.6971± 0.0182 4.8077±0.47864.8077± 0.4786 hoIVON (K=20K=20) 5.6951±0.3808 5.6951± 0.3808 2.2252±0.76922.2252± 0.7692 5.0444±1.14875.0444± 1.1487 3.9703±0.0843 3.9703± 0.0843 0.6627±0.01610.6627± 0.0161 2.6593±0.36562.6593± 0.3656 NGMP 6.1986±0.38716.1986± 0.3871 0.7198±0.0534 0.7198± 0.0534 3.8201±0.3952 3.8201± 0.3952 3.9541±0.0729 3.9541± 0.0729 0.6582±0.0189 0.6582± 0.0189 0.9295±0.1548 0.9295± 0.1548 E.2 Ensemble forecasting details We consider univariate long-horizon forecasting on the ETTh1 and ETTh2 electricity-transformer data sets (61), whose observations are hourly. From the preceding 96-hour multivariate context, the task is to predict one standardized oil-temperature target (OT) H∈96,192,336,720H∈\96,192,336,720\ hours ahead. Following 39, the frozen expert bank consists of five independently trained neural forecasters: CNN, DLinear, NLinear (60), LSTM (23), and NConv. Constant forecasts at the 10th and 90th training-set quantiles complete the frozen seven-expert bank. The gate is fitted on the original validation partition and evaluated on the original test partition; no expert or VAE representation is retrained. For the IVON gates, K=1000K=1000 denotes 1,000 draws of the gate weights from the fitted variational posterior. For logits z(k)z^(k) from draw k, the Gaussian component mean is μ(k)=∑jsoftmax(z(k))jfjμ^(k)= _jsoftmax(z^(k))_jf_j, where fjf_j is the jjth frozen expert prediction, and its variance is (σ2)(k)=(∑jexpzj(k))−1(σ^2)^(k)=( _j z_j^(k))^-1. The reported posterior predictive is the equally weighted mixture of these 1,000 components. Its NLL is evaluated with log-sum-exp over the component densities, and its point prediction is the average of the component means. Table 11 reports the full ETTh results, including the Adam-trained gates omitted from the main table. Table 11: Full descriptive ETTh1 and ETTh2 RMSE and NLL results. Intervals are 95% descriptive confidence-interval half-widths computed per test origin; RMSE intervals use the delta method. Counts are 3446, 3426, 3398, and 3321 at horizons 96, 192, 336, and 720. The PGE–PVMP rows are taken from 39, computed on the same test split. ETTh1 Method 96 192 336 720 RMSE MoE – Adam, affine gate 0.3836±0.00900.3836± 0.0090 0.3699±0.00910.3699± 0.0091 0.4152±0.00990.4152± 0.0099 0.4639±0.01110.4639± 0.0111 MoE – Adam, ReLU gate 0.3955±0.00970.3955± 0.0097 0.3699±0.00910.3699± 0.0091 0.4801±0.01230.4801± 0.0123 0.4642±0.01120.4642± 0.0112 MoE – IVON, affine gate (K=1000K=1000) 0.3832±0.00900.3832± 0.0090 0.3695±0.00910.3695± 0.0091 0.3672±0.00840.3672± 0.0084 0.5428±0.01270.5428± 0.0127 MoE – IVON, ReLU gate (K=1000K=1000) 0.3832±0.00900.3832± 0.0090 0.3697±0.00910.3697± 0.0091 0.3677±0.00840.3677± 0.0084 0.5334±0.01250.5334± 0.0125 PGE – PVMP (39) 0.3583±0.00830.3583± 0.0083 0.3386±0.00850.3386± 0.0085 0.3105±0.00730.3105± 0.0073 0.3347±0.00760.3347± 0.0076 PGE – NGMP (this work) 0.3554±0.00840.3554± 0.0084 0.3337±0.00860.3337± 0.0086 0.3115±0.00740.3115± 0.0074 0.3300±0.00800.3300± 0.0080 NLL MoE – Adam, affine gate ×1025±×10252.0553\!×\!10^25± 1.5536\!×\!10^25 ×1022±×10221.1050\!×\!10^22± 1.0354\!×\!10^22 ×1022±×10211.1159\!×\!10^22± 6.0082\!×\!10^21 ×1011±×10119.0836\!×\!10^11± 6.4103\!×\!10^11 MoE – Adam, ReLU gate ×1024±×10242.2759\!×\!10^24± 1.1906\!×\!10^24 ×1019±×10181.2482\!×\!10^19± 8.2942\!×\!10^18 ×1014±×10142.4524\!×\!10^14± 1.7209\!×\!10^14 ×1017±×10171.5279\!×\!10^17± 1.0479\!×\!10^17 MoE – IVON, affine gate (K=1000K=1000) 10.2777±0.736610.2777± 0.7366 ×104±×1043.5515\!×\!10^4± 1.1067\!×\!10^4 429.0098±143.4073429.0098± 143.4073 168.4529±26.4000168.4529± 26.4000 MoE – IVON, ReLU gate (K=1000K=1000) 2216.8905±676.27762216.8905± 676.2776 ×108±×1072.2344\!×\!10^8± 9.2280\!×\!10^7 58.6475±9.698058.6475± 9.6980 94.2422±16.263894.2422± 16.2638 PGE – PVMP (39) 0.4120±0.01730.4120± 0.0173 0.3701±0.01710.3701± 0.0171 0.3141±0.01350.3141± 0.0135 0.3763±0.01400.3763± 0.0140 PGE – NGMP (this work) 0.3888±0.02100.3888± 0.0210 0.3378±0.02000.3378± 0.0200 0.2877±0.01610.2877± 0.0161 0.3571±0.01530.3571± 0.0153 ETTh2 Method 96 192 336 720 RMSE MoE – Adam, affine gate 0.5821±0.01340.5821± 0.0134 0.5655±0.01390.5655± 0.0139 0.5190±0.01210.5190± 0.0121 0.6201±0.01290.6201± 0.0129 MoE – Adam, ReLU gate 0.5592±0.01280.5592± 0.0128 0.5655±0.01390.5655± 0.0139 0.5190±0.01210.5190± 0.0121 0.6201±0.01290.6201± 0.0129 MoE – IVON, affine gate (K=1000K=1000) 0.5821±0.01340.5821± 0.0134 0.5198±0.01290.5198± 0.0129 0.5869±0.01420.5869± 0.0142 0.8249±0.01810.8249± 0.0181 MoE – IVON, ReLU gate (K=1000K=1000) 0.5821±0.01340.5821± 0.0134 0.5655±0.01390.5655± 0.0139 0.5874±0.01420.5874± 0.0142 0.6202±0.01290.6202± 0.0129 PGE – PVMP (39) 0.5882±0.01200.5882± 0.0120 0.5798±0.01270.5798± 0.0127 0.5940±0.01280.5940± 0.0128 0.5669±0.01310.5669± 0.0131 PGE – NGMP (this work) 0.5631±0.01240.5631± 0.0124 0.5408±0.01280.5408± 0.0128 0.5463±0.01250.5463± 0.0125 0.6101±0.01430.6101± 0.0143 NLL MoE – Adam, affine gate ×1023±×10233.0584\!×\!10^23± 2.2858\!×\!10^23 ×1023±×10232.8850\!×\!10^23± 2.5057\!×\!10^23 ×1019±×10195.8244\!×\!10^19± 3.7745\!×\!10^19 ×1019±×10192.0804\!×\!10^19± 1.3583\!×\!10^19 MoE – Adam, ReLU gate ×1031±×10301.1333\!×\!10^31± 9.7716\!×\!10^30 ×1024±×10247.0543\!×\!10^24± 6.0491\!×\!10^24 ×1023±×10231.4496\!×\!10^23± 1.2198\!×\!10^23 ×1038±×10384.6866\!×\!10^38± 4.7607\!×\!10^38 MoE – IVON, affine gate (K=1000K=1000) ×1021±×10215.9904\!×\!10^21± 6.4398\!×\!10^21 304.6987±52.2098304.6987± 52.2098 ×105±×1053.7880\!×\!10^5± 1.0683\!×\!10^5 ×105±×1043.2209\!×\!10^5± 7.4793\!×\!10^4 MoE – IVON, ReLU gate (K=1000K=1000) ×1018±×10181.9494\!×\!10^18± 1.3836\!×\!10^18 ×105±×1041.1309\!×\!10^5± 6.8836\!×\!10^4 ×106±×1052.0843\!×\!10^6± 5.8045\!×\!10^5 ×106±×1051.1263\!×\!10^6± 2.1811\!×\!10^5 PGE – PVMP (39) 0.9342±0.03060.9342± 0.0306 0.9237±0.03280.9237± 0.0328 0.9612±0.03340.9612± 0.0334 0.8699±0.02990.8699± 0.0299 PGE – NGMP (this work) 0.9397±0.03860.9397± 0.0386 0.8602±0.03650.8602± 0.0365 0.8702±0.03480.8702± 0.0348 0.9774±0.03570.9774± 0.0357