Deep Gaussian Processes on Directed Acyclic Graphs

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

Deep Gaussian Processes on Directed Acyclic Graphs

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

Deep Gaussian Processes on Directed Acyclic Graphs Federico L. Perlino1

Oliver Hamelijnck2

Adam M. Johansen1

Theodoros Damoulas3,1,2

Department of Statistics, University of Warwick 2 AI Team, Unilink Software Ltd 3 Department of Computer Science, University of Warwick

arXiv:2607.09645v1 [stat.ML] 10 Jul 2026

Abstract Many real-world processes can be represented as compositions of functions along a directed acyclic graph (DAG). In causal modelling, these correspond to the underlying mechanisms; in engineering, to multiple fidelity levels; and in gene-regulatory networks, to transcription factors. These functions are partially observed across the DAG, with noisy and heterogeneously sampled measurements, posing significant challenges for reconstruction, uncertainty propagation, and inference. To tackle these challenges, we place priors over functions and naturally arrive at Deep Gaussian Processes over DAGs. We theoretically study their prior-collapse behaviour, and the effect of graph topology and intermediate observations on the preservation of information. We obtain almost-sure lower bounds on the asymptotic frequency of depths at which the distinction between inputs is preserved, identify broad kernel classes for which these hold, and prove an observation by Dunlop et al. (2018) on the role of input connections. We offer a structured variational approximation that retains graph dependencies, preserves compositional uncertainty, and captures the explaining-away behaviour of colliders. Finally, we empirically validate our theoretical results and our methodology, and model a latent-collider DAG, a protein signalling network, and a multi-fidelity heavy-ion collision emulation task, attaining state-of-the-art performance while recovering low-fidelity contributions and yielding interpretability over the simulator hierarchy.

Introduction

Various phenomena across the sciences, and beyond, can be represented as compositions of interdependent latent functions along a Directed Acyclic Graph (DAG). In probabilistic modelling, the DAG encodes conditional dependencies between quantities of interest, although the interpretation of these dependencies varies by setting. In causal models, it represents mechanistic relations (Spirtes et al., 2001; Pearl, 2009); in multi-fidelity modelling across physics and engineering (Fernández-Godino, 2023), it encodes dependencies between information sources of different fidelity (Kennedy and O’Hagan, 2000; Perdikaris et al., 2017; Ji et al., 2024); and in systems biology it is used to describe regulatory links between transcription factors in generegulatory networks (Friedman, 2004). The graph itself may be elicited from expert knowledge (Del Sagrado and Moral, 2003; O’Hagan et al., 2006), derived from mechanistic constraints and the natural directionality of the phenomenon (Boudali and Bechta Dugan, 2005), or inferred from data through structure discovery (Heckerman et al., 1995). More broadly, DAGs are often formulated by scientists as explicit representations of their hypotheses (Greenland et al., 1999). In such DAG settings, we rarely have complete, noise-free observations at every node. Data may be available only at a subset of nodes, at varying sample sizes and resolutions, and are often affected by missingness, measurement error, or model discrepancy (Little and Rubin, 2019; Carroll et al., 2006; Kennedy and O’Hagan, 2001). Inference then becomes a coupled inverse problem, in which DAG-dependent latent functions at different nodes must be jointly recovered from indirect, heterogeneous evidence. Furthermore, when an observed downstream quantity can be explained by several upstream functions, evidence for one explanation changes the posterior plausibility of the others, which is the classical explaining-away effect (Pearl, 1988; Lauritzen, 1996). At the same time, latent quantities that are never directly observed are typically weakly or non-identifiable (Raue et al., 2009; Allman et al., 2009), and hence collapsing onto any single explanation misrepresents the uncertainty in the system (Ustyuzhaninov et al., 2020).

A natural modelling response to these challenges is to place Gaussian Process priors (e.g., Rasmussen and Williams, 2006) over the latent functions composing the DAG, endowing each node with a principled probabilistic representation. Deep Gaussian Processes (DGPs) are hierarchical compositions of GP mappings (Lawrence and Moore, 2007; Damianou and Lawrence, 2013). In the standard chain setting, such compositions have been shown to improve contraction rates for compositional targets (Finocchio and Schmidt-Hieber, 2023; Giordano et al., 2022), yet deep GP priors may also collapse (Neal, 1995; Duvenaud et al., 2014; Dunlop et al., 2018; Tong and Choi, 2021), failing to preserve information with depth and motivating input (or skip) connections as a practical remedy (Neal, 1995; Duvenaud et al., 2014). DGPs therefore provide a natural, but delicate, language for compositional probabilistic modelling. We take this viewpoint to define DGPs on DAGs (DAG-DGPs), where the inductive bias of the system is directly reflected in the architecture. Compared with a standard chain DGP, a DAG-DGP exposes two modelling choices that are central in scientific applications. First, a node may have several parents, so the kernel at that node must specify how parent contributions are fused. Second, observations may be available at arbitrary internal nodes, so the model must propagate uncertainty through the graph while using intermediate measurements to anchor internal representations. This viewpoint is related to specialised multi-fidelity and information-fusion models, where lower-fidelity or intermediate simulator outputs are propagated as uncertain inputs to improve high-fidelity emulation (Perdikaris et al., 2017; Cutajar et al., 2019), and to recent graphical multi-fidelity emulators that organise such dependencies over directed trees (Ji et al., 2024). These are designed for specific DAG topologies with nested data (Le Gratiet and Garnier, 2014), leaving joint latent-function inference and uncertainty composition across a general DAG largely open. DAGs carry a rich structure, which we exploit by modelling them directly rather than using layerisation, a non-trivial graph-drawing problem utilising dummy nodes to preserve dependencies (Sugiyama et al., 1981; Harrigan and Healy, 2006). Posterior inference over the DAG composition of functions is challenging as it requires marginalisation over intermediate latent functions. Even in simpler chain DGP settings this difficulty has motivated a large literature on approximate inference, including variational, expectation-propagation, sampling, and doubly stochastic methods (Damianou and Lawrence, 2013; Hensman and Lawrence, 2014; Dai et al., 2016; Bui et al., 2016; Salimbeni and Deisenroth, 2017; Havasi et al., 2018; Salimbeni et al., 2019). In a DAG the challenge goes beyond layer-wise uncertainty propagation: observed colliders and descendants induce posterior dependence between a priori independent branches, precisely what is needed to represent explaining-away. Structured approximations for DGPs have shown that richer posterior dependence is important for calibrated uncertainty and compositional ambiguity (Ustyuzhaninov et al., 2020; Lindinger et al., 2020; Ober and Aitchison, 2021), yet existing constructions target only chain architectures or cross-layer dependence. This motivates tractable, structured variational inference that preserves posterior dependencies across the DAG. We make the following contributions: • DAG-DGP framework. We formulate the first unified DGP framework on DAGs, with noisy observations at arbitrary nodes, recovering chain DGPs, multi-fidelity DGPs, and graphical multi-fidelity emulators as special cases. • Structured variational inference for DAG-DGPs. We offer a structured variational approximation for DAG-DGPs retaining compositional uncertainty (Ustyuzhaninov et al., 2020) and explaining-away, recovering (Ustyuzhaninov et al., 2020, Sec. 4.1) and our own extension of Salimbeni and Deisenroth (2017) to DAGs as special cases. • Theoretical results on DAG-DGPs. We offer lower bounds on the frequency of DAG depths for avoiding prior-collapse and characterise how graph structure, kernel families, and intermediate observations affect information propagation. We provide explicit lower bounds for bounded-curvature exponential family-observation models. • Theoretical results on standard DGPs (i.e. chains). We prove non-collapse for the inputconnected chain DGP setting considered by (Dunlop et al., 2018) and offer the first theoretical account of standard DGPs with intermediate observations. • Empirical validation and scientific applications. We validate our theory, and demonstrate explaining-away in collider DAGs, while assessing performance, recovery, and scalability on proteinsignalling (Sachs et al., 2005) and multi-fidelity heavy-ion collision (Ji et al., 2024). 2

(a)

(b) ···

. . .

U11

(d) ···

U11

Figure 1: DAGs discussed in Secs. 2, 3, and 4 (top row) and the corresponding block sparsity patterns of the structured precision matrix Λ (bottom row): chain DGP, disjoint routes (Thm. 4), V-structure with branching, and a three-layered DAG. Circles denote latent nodes and squares nodes with partial observations. With observations placed at the terminal nodes, each DAG coincides with its moralised ancestral graph, so every inducing block contributes to Λ. Dashed grey edges denote moralisation (co-parents joined within each V-structure) and dotted black edges denote chordal fill-ins added to obtain the chordal completion H. Colours identify maximal cliques of H, with black outlines marking separator blocks.

Deep Gaussian Processes: from Chains to Directed Acyclic Graphs

We begin by reviewing a chain DGP that organises latent variables along a total order of L layers, i.e., F0 (x) = x,

fℓ ∼ GP(0, Kℓ ),

Fℓ (x) = fℓ (Fℓ−1 (x)),

ℓ = 1, . . . , L,

(1)

where x ∈ X is an input case (e.g., Fig 1(a)). Each latent layer receives a single latent input. The kernel Kℓ is therefore defined on the state space of layer ℓ − 1, and no node has multiple parents. The data enter separately from this latent recursion. In the standard supervised formulation, the dataset is D = (X, Y), and observations are linked to the final layer through an observational distribution p(Y | FL ). Intermediate observations with their own distributions have been included in specialised models, notably multi-fidelity DGPs (Cutajar et al., 2019), but their placement is then tied to the fidelity ordering. Despite being flexible for modelling hierarchical structures, the construction in Eq. (1) can suffer from prior collapse (Neal, 1995; Duvenaud et al., 2014), a problem we address in Section 4. Scalable DGP inference typically augments each layer with inducing variables and optimises a variational ELBO using Monte Carlo propagation through the latent GP conditionals (Hensman and Lawrence, 2014; Salimbeni and Deisenroth, 2017). Using Markov chain Monte Carlo can provide fully Bayesian inference, but typically at substantially higher computational cost (Havasi et al., 2018; Sauer et al., 2023). Motivated by the need to model compositional functions whose dependencies are specified by a given DAG, we develop a DAG-DGP architecture that propagates information along the graph-induced partial order, accommodates multi-parent dependencies through nodewise fusion rules, and incorporates heterogeneous observations within a single compositional model. Setup. Let G = (V, E) be a DAG (with vertices V and edges E) with roots R, non-roots U = V \ R. Denote with Pa(w) the parent set of w ∈ U. For each node w, let Fw ∈ Rn×dw collect the latent values (i) of the n observations, with row Fw ∈ Rdw . Roots r ∈ R are supplied with design matrices Xr ∈ Rn×dr , encoding deterministic inputs such as covariates, spatial locations, or time indices, and we set Fr = Xr . For notational simplicity, observations are indexed over the same set as the corresponding latent values, so a non-root node w may have response Yw ∈ Rn×dw along with an observation mask Ow ⊆ [n] × [dw ]. Only

entries in Ow enter the observational distribution, while the remaining entries are treated as missing. We let D = {Xr }r∈R , {(Yw , Ow )}w∈U denote the full dataset. DAG-DGP prior. Each non-root latent node takes as input the latent values of its parents and therefore operates Q on the product of the parent latent spaces. Concretely, for each w ∈ U we introduce a function fw : p∈Pa(w) Rdp → Rdw with prior fw ∼ GP(0, Kw ), where Kw is defined on the corresponding product space. The DAG-DGP prior is thus given by the recursion Fw(i) = fw {Fp(i) }p∈Pa(w) , w ∈ U, i ∈ [n]. (2) An additive Gaussian innovation at each node can equivalently be absorbed into Kw (Salimbeni and Deisenroth, 2017) and is omitted throughout. Since the nodewise GP modules are mutually independent a priori, the joint prior over the latent nodes factorises over the DAG into a product of Gaussian conditionals, one per non-root node given its parents; these factors are precisely those appearing in the posterior of Eq. (4). Fusion kernels. In DAGs nodes may have multiple parents (see, e.g., Fig. 1), so the way their contributions are combined is a modelling choice that shapes the inductive bias of local latent evaluations. We encode this choice through a node-wise fusion rule that directly captures parent contributions and goes beyond the typical concatenations employed in Gaussian process networks (GPNs) (Friedman and Nachman, 2000; Giudice dw ×dw et -valued positive-semidefinite on Q al., 2023;dpKiroriwal et al., 2025). The kernel Kw at node(p)w isdpR dp R . Assigning to each parent p ∈ Pa(w) a kernel K : R × R → Rdw ×dw capturing its isolated w p∈Pa(w) contribution, and setting (p) Kw := Φw Kw , (3) p∈Pa(w) where the fusion rule Φw returns a valid kernel on the full parent space is a convenient way to achieve this. The fusion rule determines whether parent effects enter independently, interact, or gate one another, and can vary from node to node to reflect heterogeneous domain knowledge across the graph. Natural instances include additive fusion, which sums the per-parent contributions and treats them as independent effects (Duvenaud et al., 2011); product fusion, which multiplies per-parent contributions and thereby allows the effect of each parent to be modulated by the others; and ANOVA-type fusion, which supplements the additive main effects with explicit pairwise interaction terms (Álvarez et al., 2012). More specialised fusion mechanisms can encode domain-specific structure, as in multifidelity models (Perdikaris et al., 2017; Cutajar et al., 2019; Ji et al., 2024). Intermediate observations. Heterogeneous observations are incorporated locally at the nodes where they are available. For each non-root node w ∈ U with observations, we specify a nodewise conditional distribution pw (Yw | Fw , Ow ) for the observed entries given the latent evaluations matrix Fw . For nodes without observations, the corresponding factor is set to one as they carry no information. These nodewise factors anchor internal latent nodes wherever data exist. Thus, the posterior is Y p {Fw }w∈U | D ∝ p0 Fw | {Fp }p∈Pa(w) pw (Yw | Fw , Ow ) . (4) w∈U

Structured Variational Inference for DAG-DGPs

Marginalising over the latent hierarchies makes the posterior in Eq. (4) computationally challenging, as in standard DGPs but over a more complex object. We introduce two doubly stochastic variational families for DAG-DGPs. The first, DAG-VI, is a mean-field inducing posterior: a scalable DAG adaptation of Salimbeni and Deisenroth (2017). The second, DAG-SVI, retains more posterior dependencies and is more expressive. The first family is a special case of the second. For each non-root node w, consider Mw inducing locations Zw in the input space of each latent, and the corresponding inducing values Uw = fw (Zw ); write F and U for the collections of latent and inducing values. Following Salimbeni and Deisenroth (2017), we define: Y q(F, U) := q(U) p0 Fw | {Fp }p∈Pa(w) , Uw . (5) w∈U

To complete the specification, we choose q(U) so as to retain the posterior couplings most directly informed by the evidence. Since only nodes w with Ow = ̸ ∅ contribute observational terms, we restrict attention to their ancestral graph. The exact latent-state posterior is Markov with respect to its moralized graph (e.g., Lauritzen, 1996), in which co-parents of a common child become adjacent. For example, in the collider in Fig. 1(c), conditioning on an observed descendant of the child induces posterior dependence between the co-parents, a phenomenon known as explaining away. This moralized graph captures the posterior couplings closest to the observational distributions, but it need not be decomposable. We therefore pass to a chordal completion H, which admits a clique-separator Gaussian representation and sparse Cholesky elimination (Lauritzen, 1996; Rue and Held, 2005). Let C1 , . . . , Ck denote the maximal cliques of H, with separators Si := Ci ∩(C1 ∪· · ·∪Ci−1 ). Decomposability of H then yields the clique factorization (Green and Thomas, 2013, Eq. (2)) qH (U) =

k Y

q UCi \Si | USi ,

(6)

i=1

in which each factor is a conditional Gaussian on the clique residual given its separator. Equivalently, qH is jointly Gaussian N (m, Σ) with precision Λ = Σ−1 supported on H; the maximal cliques in Fig. 1 are the dense blocks of Λ. Inducing blocks lying in a common clique are freely correlated, so explaining-away between co-parents and the couplings induced along ancestral paths to observed nodes are preserved; blocks that share no clique are factorised out, dropping in particular dependencies between disjoint branches of the DAG with no common observed descendant. In the standard chain, H is a path, so DAG-SVI recovers the block-tridiagonal precision family of Ustyuzhaninov et al. (2020, Sec. 4.1). The resulting ELBO for our model is, where qH (Fw ) is the marginal induced by qH (U), ! X Y LH = EqH (Fw ) [log pw (Yw | Fw , Ow )] − KL qH (U) p0 (Uw ) , (7) w∈U

w: Ow ̸=∅

and the expectation is estimated with mini-batching by ancestral Monte Carlo over the latent DAG in topological order, while the KL term is analytic. Imposing the stronger restriction ̸ w Q Λvw = 0 for all v = gives a mean-field approximation over the inducing outputs (DAG-VI), qVI (U) = w∈U qw (Uw ), recovering the DAG-DGP adaptation of Salimbeni and Deisenroth (2017). This restriction still propagates uncertainty through the DAG, but removes posterior dependence between distinct latent nodes, and therefore cannot represent explaining-away (see App. E.6) or compositional uncertainty, i.e., posterior uncertainty over the unobserved latent functions (Ustyuzhaninov et al., 2020). Scaling up to larger DAGs. The expectation in Eq. (7) is Branching-tree depth estimated by marginal ancestral sampling (see Prop. (9)). Since qH is in canonical form, the sampler requires marginal and conditional moments, equivalently selected applications of Λ−1 . Let J = |U|, M = maxw dim(Uw ), and K = SB for S Monte Carlo samples and minibatch size B. For large DAGs, DAG-VI is the most scalable option, achieving O(JM 3 +KJM 2 ) complexity at the expense of expressivity. DAG-SVI trades higher cost for more expressivity. A dense implementation forms Σ = Λ−1 , costing O((JM )3 + K(JM )2 ) per ELBO evaluation, and is preferable when JM is moderate or H is close to dense. For larger sparse Total inducing dimension (20 per latent node) DAGs, the block sparsity of Λ (Fig. 1) can instead be exploited directly: sparse Cholesky provides the Gaussian solves required Figure 2: Wall-clock time per ELBO evaluaby the sampler, while Takahashi selected-inversion recursions tion on increasingly deep branching trees. provide the diagonal covariance blocks needed for the analytic KL term (Takahashi et al., 1973; Erisman and Tinney, 1975). In favourable sparse regimes, the base sparse cost is O(Jc2 M 3 + KJM 2 ), where c is the largest number of inducing blocks in any clique of H, up to the conditioning-update factors accounted for in App. E.5. Fig. 2 reports the resulting wall-clock behaviour on a branching-tree benchmark. 234 5

102

DAG-SVI (dense ¤

Wall-clock time (s, log scale)

DAG-VI

¡1

)

DAG-SVI (sparse ¤ ¡1 )

101

100

10−1

10−2

2500

5000

7500

10000

12500

15000

17500

20000

0.4 0.2

lower bound 1 ¡ (1 ¡ p" ) s empirical ½^"

0.0 0

Separating nodes per slice s

10−1 10−2 10−3 10−4 10−5

Effect of outdegree

Refresh from Gaussian intermediate observations 1.0

0.6

Empirical survival probability

0.6

Expected squared contrast C`

Persistent fraction ½^"

0.8

Effect of indegree under product fusion 100

k=1 k=2 2

k=4 k=8 4

L=3 L=6 L=9

0.5

0.8

0.4

0.6

0.3

0.4

0.2

0.1 0.0

Chain empirical Chain theorem Disjoint-route DAG empirical Disjoint-route DAG theorem

c Predictive ¦ `0 ! `0 + h (M`0 + h > ")

Separating nodes prevent prior collapse 1.0

0.0 1

Depth `

Branching factor b

Distance from observed antichain h

Figure 3: Empirical validation of the main theoretical results in Sec. 4. See also App. G.2.

On Prior and Posterior Non-collapse in DAG-DGPs

A central pathology of DGP priors is the loss of separation between distinct inputs under repeated composition (Duvenaud et al., 2014; Dunlop et al., 2018; Tong and Choi, 2021), hindering information preservation across depth. Following recent usage (Meng and Zhang, 2024), we call this prior collapse. For DAG-DGPs this raises a richer question: given two cases a = ̸ b, when does a difference at the roots, or refreshed by internal observations, remain visible downstream? Following Dunlop et al. (2018), we measure distinguishability at (a) (b) node w via the two-case contrast ∆w := Fw − Fw , and say w carries an ε-contrast for (a, b) if ∥∆w ∥2 > ε. Input connections in chain DGPs (Neal, 1995; Duvenaud et al., 2011) can be viewed as deterministic refreshes of this contrast; general DAGs admit further such mechanisms via topology, intermediate observations, and root-dependent fusion kernels.

4.1

Repeated separating nodes prevent prior collapse

Unlike chains, DAGs lack a unique notion of layer. We therefore consider progressive antichain decompositions Fh−1 V = ℓ=0 Aℓ , where each Aℓ is an antichain (no two nodes connected by a directed path) and every node in Aℓ reaches some node in Aℓ+1 . Each Aℓ acts as a depth slice, recovering single-layer slices in chains. Every finite DAG admits such a decomposition (App. B.1); asymptotic statements are read along increasing-depth sequences or truncations. We track the largest contrast on an antichain, Mℓ := maxw∈Aℓ ∥∆w ∥2 , since a single non-collapsed node suffices to distinguish the two cases at depth ℓ. Prior collapse for (a, b) means Mℓ → 0 a.s. For a non-root node w, let Γw (a, b) be the conditional covariance of ∆w given parent evaluations (a) (b) {Fp , Fp }p∈Pa(w) . We call w v⋆ -separating for (a, b) if Γw (a, b) ⪰ v⋆ Idw a.s.; this requires the parents to expose a contrast visible to the kernel at w. A separating node injects at least v⋆ of conditional variance into the difference, giving a uniformly positive chance of counteracting collapse. Theorem 1 (Repeated separating nodes prevent prior collapse). Under the DAG-DGP prior, suppose that every Aℓ in a progressive antichain sequence contains at least s ≥ 1 v⋆ -separating nodes for the pair (a, b). m−1 Then, 1 X ε ∀ε > 0 : lim inf 1 {Mℓ > ε} ≥ 1 − (1 − pε )s a.s. pε := 2ΦN − √ . m→∞ m v⋆ ℓ=0 Hence non-trivial contrasts occur on a positive fraction of depths, ruling out prior collapse (Fig. 3(a)). In a chain each antichain has a single latent node, and the input connection induces separation; for the squared-exponential input-connected chain of (Dunlop et al., 2018, Remark 5(3)), Corollary 2 shows this node is separating at every depth, so the bound applies with s = 1. Anchors need not be placed everywhere: if separating nodes occur infinitely often, permanent collapse is ruled out, and if they occur on a positive fraction of depths the bound is multiplied by that fraction (App. B.4). Separation is easy to certify if the kernel retains non-degenerate dependence on root coordinates. Additive and ANOVA root-only components give transparent sufficient conditions, as their contribution cannot be cancelled by other parents; this covers fusion kernels in multifidelity models (Perdikaris et al., 2017; Cutajar et al., 2017; Ji et al., 2024). Perfectly observed internal nodes also act as separating coordinates after conditioning (App. B.7).

4.2

Effect of the DAG topology

Theorem 1 identifies anchors against collapse; we next isolate topological effects, starting with indegree. Consider a local layered radial block (Fig. 1(d) shows a layered DAG), where parents of Aℓ lie in Aℓ−1 and each node uses a radial kernel on the concatenated parent state. The expected contrast then admits an explicit recursion (full assumptions in App. C). For clarity we display product fusion of squared-exponential parent kernels; App. C extends the recursion to general monotonic Laplace-radial kernels, including rational-quadratic kernels and layer-dependent dimensions or hyperparameters. Let kℓ = maxw∈Aℓ |Pa(w)| and Cℓ := maxw∈Aℓ E[∥∆w ∥22 ]. For product squared-exponential fusion with variance τ 2 , length-scale λ, and output dimension d, App. C gives h −dkℓ /2 i Cℓ ≤ 2dτ 2 1 − 1 + Cℓ−1 /(dλ2 ) . For kℓ = 1 this recovers the chain recurrence of Tong and Choi (2021); for kℓ > 1, the kernel sees contrast accumulated across incoming edges, broadening the recurrence (Fig. 3(b)). Near zero the slope is dτ 2 kℓ /λ2 , yielding the sufficient contraction criterion dτ 2 k̄/λ2 < 1 with k̄ := supℓ kℓ . Indegree thus enters the contraction threshold linearly, making the chain collapse mechanism less transferable to high-indegree regions under product fusion. Additive root-retaining fusion is insensitive to indegree, provided its root-retaining component is separating (Prop. 5, App. C.2). Outdegree provides a complementary mechanism: large outdegree yields many parallel descendants, hence many conditionally independent chances for a contrast to survive. In the same layered radial regime, consider a b-ary branching subgraph through successive antichains, a topology arising in various probabilistic models (Jordan and Jacobs, 1994; Liu et al., 2024). If the seed contrast exceeds t > 0 with positive probability and the one-step exceedance probability satisfies pt > 1/b, then with positive probability the maximum contrast stays above t at every depth (Prop. 6, App. C.3; Fig. 3(c)).

4.3

Intermediate observations as stochastic skip connections

Skip connections preserve information by reintroducing input-dependent variation at later depths; intermediate observations have an analogous posterior effect. For two cases a = ̸ b and an observed internal node u, after (a) (b) (a) (b) observing Yu , Yu the posterior may assign substantial mass to ∥Fu − Fu ∥2 > ε, making the observation a local source of separation at u. Being noisy, it does not guarantee a genuine latent contrast, but updates the posterior mass on the latent ε-contrast event at u. To affect later antichains, this contrast must then propagate along directed routes. Let Πℓ0 be the filtering posterior after assimilating observations on the strict ancestors of Aℓ0 and on Aℓ0 itself, Hℓ0 the sigma-field of strict ancestral states, and Πℓ0 →ℓ1 the forward predictive law from Aℓ0 to Aℓ1 without assimilating later observations. An admissible route has off-route parents fixed at the source antichain, so its contrast can be tracked in isolation; an ε-retaining route with factor ρ carries an ε-contrast from source to target with probability at least ρ; routes that share no nodes between source and target propagate independently (definitions in App. D). Theorem 2 (Intermediate observations act as stochastic skip connections). Consider two input cases a = ̸ b, a threshold ε > 0, and antichain levels ℓ1 ≥ ℓ0 . Assume that s distinct observed source nodes u1 , . . . , us ∈ Aℓ0 can reach s distinct target nodes v1 , . . . , vs ∈ Aℓ1 through pairwise interior-disjoint admissible routes γj : uj ⇝ vj . (a) (b) Assume further that each γj is ε-retaining with factor ρj ∈ [0, 1]. Define Ej := ∥Fuj − Fuj ∥2 > ε , and s qj := Πℓ0 (Ej | Hℓ0 ). Then Y Πℓ0 →ℓ1 (Mℓ1 > ε) ≥ EΠℓ0 1 − (1 − ρj qj ) . (8) j=1

The bound has a direct interpretation: qj is the posterior probability (after observations up to level ℓ0 ) that source uj carries a genuine latent ε-contrast, while ρj measures how likely route γj propagates it downstream. Thus ρj qj is the contribution of route j, and the product in (8) is the probability that no route succeeds. In a chain (Fig. 1(a)) with one observed internal layer and one partially observed downstream route, this reduces to Πℓ0 →ℓ1 (Mℓ1 > ε) ≥ EΠℓ0 [ρq]. In DAGs—e.g. the disjoint-route DAG of Fig. 1(b)—several observed sources combine via the probability that at least one propagated contrast reaches the target antichain; Fig. 3(d) compares the two regimes empirically.

1.899

×165

Prior VI posterior SVI posterior

1 1.574 1.249

1 1

Fw2 (x2* )

0.3292

DAG-SVI

1 ×119

0.275

0.049

0.600

0.1646

Absolute error

0.925

DAG-VI

w3 value

0.0000

RMSE = 0.106 NLL = -0.234

1 3

Fw1 (x1* )

RMSE = 0.039 NLL = -0.604

1 1

Figure 4: Latent-collider experiment. Left: ground-truth w3 , input marked. Centre: parent posterior under DAG-VI and DAG-SVI. Right: w3 and per-method errors, with root mean square error (RMSE) and negative log likelihood (NLL) vs. observations. (i)

(i)

(a)

For Gaussian observations Yu = Fu + ξu , ξu ∼ N (0, σu2 ), set Γ = Γu (a, b) and DY = Yu Γ > 0, then conditionally on Hℓ0 , Γ 2σu2 Γ Fu(a) − Fu(b) | Yu(a) , Yu(b) , Hℓ0 ∼ N D , . Y Γ + 2σu2 Γ + 2σu2

(b)

− Yu . If (9)

The shrinkage factor Γ/(Γ + 2σu2 ) interpolates the source strength between observation-driven (σu2 ≪ Γ) and noise-dominated (σu2 ≫ Γ) regimes, so the qj of Thm. 2 reduce to explicit Gaussian tails, empirically validated in Fig. 3(d). App. D.4 extends the analysis to bounded-curvature one-parameter exponential families, including Bernoulli and binomial nodes.

Compositional Uncertainty and Explaining Away in Latentcolliders

We consider a synthetic latent-collider experiment Table 1: Test performance on both the real-world (COLLIDER) in which independent GP parents, w1 Sachs and the synthetic collider dataset. and w2 , feed a partially observed child w3 through an additive RBF fusion kernel. Fig. 4 (centre) visualises Dataset Metric DAG-VI DAG-SVI the posterior over the two parent latent nodes at a RMSE ↓ 0.646 0.642 SACHS Interpolation CRPS ↓ 0.355 0.359 fixed input, illustrating two posterior phenomena that RMSE ↓ 0.737 0.642 DAG-SVI captures and that DAG-VI cannot represent. SACHS Extrapolation CRPS ↓ 0.325 0.299 In terms of variability, the DAG-SVI samples cover a RMSE ↓ 0.114 0.055 broad range of latent parent representations, preservCOLLIDER CRPS ↓ 0.075 0.060 ing compositional uncertainty; DAG-VI, restricted by construction to factorised marginals latent evaluations, instead concentrates onto a near-degenerate point. In terms of dependencies, the DAG-SVI is negatively correlated, encoding the reciprocal compensation through which the two parents explain the child, signature of explaining away; mean-field DAG-VI produces uncorrelated samples by design and cannot represent this coupling. App. G.3 provides a geometric visualization of this behaviour, and the right panel of Fig. 4(right) shows that coupling translates into improved reconstruction of the partially observed child.

Protein Signalling Network from the Sachs Flow Cytometry Dataset

We evaluate the proposed DAG-DGP framework on the realworld protein signalling network dataset of Sachs et al. (2005), a widely used benchmark in the causal discovery literature (Mooij et al., 2020). Since our focus is observational modelling, we consider the cd3cd28 icam2 subset, which corresponds to a specific intervention regime. The dataset contains 902 observations over 11 variables, which we log-transform prior to modelling. We use the DAG structure provided by bnlearn (Scutari, 2010) (Fig. 5) and construct a random 80/20 train– test split. It is well known that the Sachs data were not generated under ideal interventions and may therefore contain latent confounding effects. To account for this, we explicitly model a confounder through a 2D latent variable layer (Salimbeni et al., 2019). At test time, inference requires the posterior distribution over these latent variables; consequently, we assume that both PKA and Raf are fully observed in both the training and test sets. Following Lindinger et al. (2020), we consider an interpolation and extrapolation task detailed, along with the training procedure, in App. G.4. Results, averaged over all nodes, are reported in Tab. 1 and are evaluated in log space. Both DAG-VI and DAG-SVI successfully capture the joint distribution induced by the DAG structure in the interpolation task, and DAG-SVI improves prediction, uncertainty quantification, and coverage (see Tab. 3) in extrapolation.

Unobserved Observed Input

PKC Plcg

PIP3

Jnk

PKA

P38 PIP2

Mek Erk Akt

Raf

Figure 5: DAG for the Sachs protein signalling network. Diamond nodes denote exogenous/root inputs. Blue nodes correspond to observed variables on which deep GP priors are placed. The dotted node denotes a latent confounder, modelled using a 2D latent variable layer (Salimbeni et al., 2019). The first latent dimension acts as input to PKA, and the second to Raf.

Deploying DAG-DGPs for Multi-Fidelity Heavy-ion Emulation YL1 YH L1 X

H L2

YL2

Figure 6: DAG-DGP for the heavy-ion emulation task.

We evaluate DAG-DGPs on the heavy-ion collision real dataset of Ji et al. (2024), a graphical multi-fidelity emulation problem with a shared nine-dimensional input and a scalar pion-yield-ratio output. The elicited simulator graph has two lower-fidelity nodes, L1 and L2 , feeding the high-fidelity node H. Instead of being sequential, these lower fidelities are complementary: L1 uses simplified linearized conformal hydrodynamics followed by Cooper–Frye conversion, whereas L2 uses 1 + 1D ideal QCD (quantum chromodynamics) hydrodynamics but omits this conversion stage. This makes it meaningful to study how each approximation contributes to explaining the high-fidelity response, and whether one provides more informative support for H than the other. Accordingly, we fit a DAG-DGP over the elicited simulator graph in Figure 6.

We report three evaluations. On the published split of Ji et al. (2024), with 200 observations at each lower fidelity, 25 observations at H, and the original 75 high-fidelity test points, DAG-VI improves on the graphical multi-fidelity Gaussian process (GMGP) family and their benchmark comparison; DAG-SVI gives the best performance on root-mean-square error (RMSE), normalised RMSE (N-RMSE; see Ji et al. (2024)), and continuous ranked probability score (CRPS); see Tab. 2. On the same test task, for our best model (DAG-SVI) we further evaluate how the two parents of H contribute to its posterior via Shapley values (Fig 7): both lower-fidelity branches contribute substantially for a large fraction of high-fidelity test points, while L2 provides the dominant contribution overall. This is 9

consistent with the simulator construction, where L2 preserves a more realistic QCD hydrodynamic evolution, whereas L1 retains the conversion stage but uses a simpler hydrodynamic approximation. The remaining protocols isolate two different effects. The high-fidelity-scarce protocol performs repeated 5fold cross-validation over the 25 observations at H, testing cross-fidelity transfer when high-fidelity supervision is most limited; DAG-SVI again gives the strongest H-level predictions. The full-hierarchy protocol uses 200, 200, and 100 observations at L1 , L2 , H, respectively, and performs 10-fold cross-validation with held-out data at all fidelities to evaluate joint prediction. In this data-richer setting, DAG-SVI improves both joint and marginal high-fidelity metrics over DAG-VI, while both models remain competitive. Results and training procedures are in App. G.5.

High-fidelity GP KO-path (Kennedy and O’Hagan, 2000) KO-misspecified (Kennedy and O’Hagan, 2000) NARGP (Perdikaris et al., 2017) r-GMGP (Ji et al., 2024) d-GMGP (Ji et al., 2024) DAG-VI† DAG-SVI†

RMSE /10−2 ↓

N-RMSE ↑

CRPS /10−2 ↓

5.49 3.48 3.95 3.66 2.92 2.17

0.46 0.66 0.71 0.73 0.72 0.79

3.54 1.99 2.30 2.13 1.64 1.34

1.0

Local normalized Shapley share

Model

2.12 ± 0.02 0.796 ± 0.002 1.18 ± 0.03 2.03 ± 0.01 0.804 ± 0.001 1.16 ± 0.02

Table 2: Predictive performance on the published heavy-ion split, evaluated on the 75-point highfidelity test set. The top block reproduces results from (Ji et al., 2024). † Our methods are reported over five runs as mean ± std.

0.8

0.6

0.4

0.2

0.0 L1

Predictive Mean

Posterior Variance

Figure 7: Normalized Shapley shares for L1 and L2 on the published heavy-ion test set under DAG-SVI. Left: predictive mean of H. Right: posterior variance in H.

Discussion

8.1

DAG-DGPs as a General Framework

DAG-DGPs recover and extend several Gaussian-process architectures by restricting three components: the DAG topology, the nodewise observation pattern, and the fusion rule. This yields two concrete benefits. First, it lets us flexibly model noisy, heterogeneously observed data organised by general DAG topologies, including graphical multi-fidelity problems and Bayesian network datasets. Second, whenever these restrictions recover an existing architecture, the resulting special case inherits our theoretical analysis and structured variational inference scheme. DGP models. Standard DGPs (Lawrence and Moore, 2007; Damianou and Lawrence, 2013) are recovered by choosing a chain DAG, using the trivial single-parent fusion rule, and observing only the terminal layer. Input-connected DGPs (Duvenaud et al., 2014) add the deterministic input as a parent of every latent layer. Multi-fidelity DGPs (Cutajar et al., 2019) arise by interpreting the chain as an ordered fidelity hierarchy, observing the corresponding fidelity nodes, and using a multi-fidelity fusion rule. d-GMGPs (Ji et al., 2024) further restrict this multi-fidelity topology to a directed in-tree with nested designs, where each simulator node is connected to the shared deterministic input and its lower-fidelity parents. Stochastic deep Gaussian processes over graphs (Li et al., 2020) target input–output transformations between signals on a fixed graph. Although their modelling aim differs from ours, they can be recovered as layered DAG-DGPs by unrolling the fixed graph over depth, as we show in Prop. 10 of App. F. GPN-based models. GPN-based models (Friedman and Nachman, 2000; Giudice et al., 2023; Kiroriwal et al., 2025) are obtained by choosing the DAG to be a process network, e.g. a multi-stage system where subprocess outputs feed downstream stages. Classical GPNs correspond to the fully observed case, where all measurements are available. When full observability is relaxed, as recently proposed in Bayesian optimisation

(Kiroriwal et al., 2025), subprocess GPs are conditioned on deterministic stage inputs, yielding an inputconnected DAG; the RBF kernel specified in that model over the concatenated parent/input space then corresponds to a product fusion rule in our framework.

8.2

Open Challenges and Future Directions

We introduced a unifying modelling framework for composition of functions over Directed Acyclic Graphs, motivated by the need to represent such inductive biases of domain-knowledge systems in science and engineering within probabilistic machine learning. Modelling such systems naturally calls for a probabilistic treatment, in which uncertainty over latent functions and their graph-induced dependencies is retained. When the DAG is given a causal interpretation, the framework opens the door for structural causal modelling and causal representation learning (Pearl, 2009; Spirtes et al., 2001; Peters et al., 2017; Schölkopf et al., 2021). Several directions remain open. Beyond the collapse phenomenon addressed in our theoretical analysis, an intriguing avenue is to explore posterior contraction rates of DAG-DGPs, extending recent advances developed for chain DGPs (Finocchio and Schmidt-Hieber, 2023). We assumed a well-specified DAG; in practice, domain knowledge specifies the graph only imperfectly, raising the question of how to robustify DAG-DGPs against DAG misspecification or perform joint inference over DAGs and composing latent functions, drawing inspiration from e.g. Branchini et al. (2023); Chickering (2002); Zheng et al. (2018); Witty et al. (2020); Aglietti et al. (2020); Giudice et al. (2024). In terms of robustness to likelihood or prior misspecification, our variational framework could be easily extended towards Generalised Variational Inference (GVI) (Knoblauch et al., 2022). Our structured and mean-field doubly stochastic VI schemes trade off posterior dependencies and uncertainty quantification against computational efficiency; alternative trade-offs could be explored via different graph-layering strategies for the approximate posterior Harrigan and Healy (2006), or by extending sampling methodologies or recent hybrid optimization sampling schemes developed for chain DGPs (Havasi et al., 2018; Sauer et al., 2023; Latz et al., 2025) to the full DAG setting. Finally, scaling DAG-DGPs to much larger graphs remains an open challenge, with potential directions including asynchronous distributed training, message passing, state-space formulations, and back-propagation of evidence.

Acknowledgments We are especially grateful to Yi Ji, Simon Mak, Derek Soeder, J.-F. Paquet, and Steffen A. Bass for making available the code and data for the heavy-ion collision experiment in Ji et al. (2024). This work was supported by United Kingdom Research and Innovation (UKRI) via grant number EP/Y014650/1, as part of the ERC Synergy project OCEAN.

References Virginia Aglietti, Theodoros Damoulas, Mauricio A. Álvarez, and Javier González. Multi-task causal learning with Gaussian processes. In Advances in Neural Information Processing Systems, volume 33, pages 6293–6304, 2020. URL https://proceedings.neurips.cc/paper/2020/hash/ 45c166d697d65080d54501403b433256-Abstract.html. Elizabeth S. Allman, Catherine Matias, and John A. Rhodes. Identifiability of parameters in latent structure models with many observed variables. The Annals of Statistics, 37(6A):3099–3132, 2009. doi: 10.1214/ 09-AOS689. Mauricio A. Álvarez, Lorenzo Rosasco, and Neil D. Lawrence. Kernels for vector-valued functions: A review. Foundations and Trends in Machine Learning, 4(3):195–266, 2012. doi: 10.1561/2200000036. Stefan Arnborg, Derek G. Corneil, and Andrzej Proskurowski. Complexity of finding embeddings in a k-tree. SIAM Journal on Algebraic Discrete Methods, 8(2):277–284, 1987. doi: 10.1137/0608024. Krishna B. Athreya and Peter E. Ney. Branching Processes, volume 196 of Die Grundlehren der mathematischen Wissenschaften. Springer, Berlin, Heidelberg, 1972. ISBN 9783642653711. doi: 10.1007/978-3-642-65371-1.

Kazuoki Azuma. Weighted sums of certain dependent random variables. Tohoku Mathematical Journal, 19(3): 357–367, 1967. doi: 10.2748/tmj/1178243286. Jørgen Bang-Jensen and Gregory Z. Gutin. Digraphs: Theory, Algorithms and Applications. Springer Monographs in Mathematics. Springer, London, second edition, 2009. ISBN 9781848009981. doi: 10.1007/ 978-1-84800-998-1. Anne Berry, Jean R. S. Blair, Pinar Heggernes, and Barry W. Peyton. Maximum cardinality search for computing minimal triangulations of graphs. Algorithmica, 39(4):287–298, 2004. doi: 10.1007/s00453-004-1100-0. Christopher M. Bishop. Pattern Recognition and Machine Learning. Information Science and Statistics. Springer, New York, NY, 2006. ISBN 9780387310732. Hichem Boudali and Joanne Bechta Dugan. A discrete-time Bayesian network reliability modeling and analysis framework. Reliability Engineering & System Safety, 87(3):337–349, 2005. doi: 10.1016/j.ress.2004.06.004. Nicola Branchini, Virginia Aglietti, Neil Dhir, and Theodoros Damoulas. Causal entropy optimization. In Francisco Ruiz, Jennifer Dy, and Jan-Willem van de Meent, editors, Proceedings of The 26th International Conference on Artificial Intelligence and Statistics, volume 206 of Proceedings of Machine Learning Research, pages 8586–8605. PMLR, 25–27 Apr 2023. URL https://proceedings.mlr.press/v206/branchini23a. html. Thang D. Bui, Daniel Hernández-Lobato, José Miguel Hernández-Lobato, Yingzhen Li, and Richard E. Turner. Deep Gaussian processes for regression using approximate expectation propagation. In Proceedings of the 33rd International Conference on Machine Learning, volume 48 of Proceedings of Machine Learning Research, pages 1472–1481. PMLR, 2016. URL https://proceedings.mlr.press/v48/bui16.html. Raymond J. Carroll, David Ruppert, Leonard A. Stefanski, and Ciprian M. Crainiceanu. Measurement Error in Nonlinear Models: A Modern Perspective, volume 105 of Monographs on Statistics and Applied Probability. Chapman & Hall/CRC, Boca Raton, FL, second edition, 2006. ISBN 9781420010138. doi: 10.1201/9781420010138. David Maxwell Chickering. Optimal structure identification with greedy search. Journal of Machine Learning Research, 3(Nov):507–554, 2002. URL https://jmlr.org/papers/v3/chickering02b.html. Kurt Cutajar, Edwin V. Bonilla, Pietro Michiardi, and Maurizio Filippone. Random feature expansions for deep Gaussian processes. In Doina Precup and Yee Whye Teh, editors, Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 884–893. PMLR, 2017. URL https://proceedings.mlr.press/v70/cutajar17a.html. Kurt Cutajar, Mark Pullin, Andreas Damianou, Neil D. Lawrence, and Javier González. Deep Gaussian processes for multi-fidelity modeling. arXiv preprint arXiv:1903.07320, 2019. doi: 10.48550/arXiv.1903.07320. URL https://arxiv.org/abs/1903.07320. Zhenwen Dai, Andreas C. Damianou, Javier González, and Neil D. Lawrence. Variational auto-encoded deep Gaussian processes. In Proceedings of the 4th International Conference on Learning Representations, 2016. doi: 10.48550/arXiv.1511.06455. URL https://openreview.net/forum?id=TqCgctDVcEz. ICLR 2016, Conference Track. Andreas Damianou and Neil D. Lawrence. Deep Gaussian processes. In Carlos M. Carvalho and Pradeep Ravikumar, editors, Proceedings of the Sixteenth International Conference on Artificial Intelligence and Statistics, volume 31 of Proceedings of Machine Learning Research, pages 207–215, Scottsdale, Arizona, USA, 2013. PMLR. URL https://proceedings.mlr.press/v31/damianou13a.html. A. Philip Dawid. Conditional independence in statistical theory. Journal of the Royal Statistical Society: Series B (Methodological), 41(1):1–15, 1979. doi: 10.1111/j.2517-6161.1979.tb01052.x. José Del Sagrado and Serafín Moral. Qualitative combination of Bayesian networks. International Journal of Intelligent Systems, 18(2):237–249, 2003. doi: 10.1002/int.10086. 12

Matthew M. Dunlop, Mark A. Girolami, Andrew M. Stuart, and Aretha L. Teckentrup. How deep are deep Gaussian processes? Journal of Machine Learning Research, 19(54):1–46, 2018. URL https: //jmlr.org/papers/v19/18-015.html. David Duvenaud, Oren Rippel, Ryan Adams, and Zoubin Ghahramani. Avoiding pathologies in very deep networks. In Samuel Kaski and Jukka Corander, editors, Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics, volume 33 of Proceedings of Machine Learning Research, pages 202–210, Reykjavik, Iceland, 22–25 Apr 2014. PMLR. URL https://proceedings.mlr.press/v33/ duvenaud14.html. David K. Duvenaud, Hannes Nickisch, and Carl Edward Rasmussen. Additive Gaussian processes. In Advances in Neural Information Processing Systems, volume 24, pages 226– 234. Curran Associates, Inc., 2011. URL https://proceedings.neurips.cc/paper/2011/hash/ 4c5bde74a8f110656874902f07378009-Abstract.html. Albert M. Erisman and William F. Tinney. On computing certain elements of the inverse of a sparse matrix. Communications of the ACM, 18(3):177–179, March 1975. doi: 10.1145/360680.360704. M. Giselle Fernández-Godino. Review of multi-fidelity models. Advances in Computational Science and Engineering, 1(4):351–400, 2023. doi: 10.3934/acse.2023015. Gianluca Finocchio and Johannes Schmidt-Hieber. Posterior contraction for deep Gaussian process priors. Journal of Machine Learning Research, 24(66):1–49, 2023. URL https://jmlr.org/papers/v24/21-0556. html. Nir Friedman. Inferring cellular networks using probabilistic graphical models. Science, 303(5659):799–805, 2004. doi: 10.1126/science.1094068. Nir Friedman and Iftach Nachman. Gaussian process networks. In Proceedings of the Sixteenth Conference on Uncertainty in Artificial Intelligence, pages 211–219. Morgan Kaufmann, 2000. Matteo Giordano, Kolyan Ray, and Johannes Schmidt-Hieber. On the inability of Gaussian process regression to optimally learn compositional functions. In Advances in Neural Information Processing Systems, volume 35, pages 22341–22353. Curran Associates, Inc., 2022. URL https://proceedings.neurips.cc/paper_files/ paper/2022/hash/8c420176b45e923cf99dee1d7356a763-Abstract-Conference.html. Enrico Giudice, Jack Kuipers, and Giusi Moffa. A Bayesian take on Gaussian process networks. In Advances in Neural Information Processing Systems, volume 36, pages 56602– 56614, 2023. URL https://proceedings.neurips.cc/paper_files/paper/2023/hash/ b146e7c87685fa208bd95ce4b08e330c-Abstract-Conference.html. Enrico Giudice, Jack Kuipers, and Giusi Moffa. Bayesian causal inference with Gaussian process networks, 2024. URL https://arxiv.org/abs/2402.00623. arXiv:2402.00623. Peter J. Green and Alun Thomas. Sampling decomposable graphs using a Markov chain on junction trees. Biometrika, 100(1):91–110, 2013. doi: 10.1093/biomet/ass052. Sander Greenland, Judea Pearl, and James M. Robins. Causal diagrams for epidemiologic research. Epidemiology, 10(1):37–48, 1999. Martin Harrigan and Patrick Healy. On layering directed acyclic graphs. In Michael Jünger, Stephen Kobourov, and Petra Mutzel, editors, Graph Drawing, volume 5191 of Dagstuhl Seminar Proceedings (DagSemProc), Dagstuhl, Germany, 2006. Schloss Dagstuhl – Leibniz-Zentrum für Informatik. doi: 10.4230/DagSemProc. 05191.6. URL https://drops.dagstuhl.de/entities/document/10.4230/DagSemProc.05191.6. Marton Havasi, José Miguel Hernández-Lobato, and Juan José Murillo-Fuentes. Inference in deep Gaussian processes using stochastic gradient Hamiltonian Monte Carlo. In Advances in Neural Information Processing Systems, volume 31, pages 7506–7516. Curran Associates, Inc., 2018. URL https://proceedings.neurips. cc/paper/2018/hash/4172f3101212a2009c74b547b6ddf935-Abstract.html. 13

David Heckerman, Dan Geiger, and David M. Chickering. Learning Bayesian networks: The combination of knowledge and statistical data. Machine Learning, 20(3):197–243, 1995. doi: 10.1007/BF00994016. James Hensman and Neil D. Lawrence. Nested variational compression in deep Gaussian processes, 2014. URL https://arxiv.org/abs/1412.1370. arXiv:1412.1370. Wassily Hoeffding. Probability inequalities for sums of bounded random variables. Journal of the American Statistical Association, 58(301):13–30, 1963. doi: 10.1080/01621459.1963.10500830. Yi Ji, Simon Mak, Derek Soeder, Jean-François Paquet, and Steffen A. Bass. A graphical multi-fidelity Gaussian process model, with application to emulation of heavy-ion collisions. Technometrics, 66(2):267–281, 2024. doi: 10.1080/00401706.2023.2281940. Michael I. Jordan and Robert A. Jacobs. Hierarchical mixtures of experts and the EM algorithm. Neural Computation, 6(2):181–214, 1994. doi: 10.1162/neco.1994.6.2.181. Marc C. Kennedy and Anthony O’Hagan. Predicting the output from a complex computer code when fast approximations are available. Biometrika, 87(1):1–13, 2000. doi: 10.1093/biomet/87.1.1. Marc C. Kennedy and Anthony O’Hagan. Bayesian calibration of computer models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(3):425–464, 2001. doi: 10.1111/1467-9868.00294. Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations, 2015. URL https://openreview.net/forum?id=8gmWwjFyLj. Saksham Kiroriwal, Julius Pfrommer, and Jürgen Beyerer. Bayesian optimization using partially observable Gaussian process network. In NeurIPS 2025 Workshop MLxOR: Mathematical Foundations and Operational Integration of Machine Learning for Uncertainty-Aware Decision-Making, 2025. Jeremias Knoblauch, Jack Jewson, and Theodoros Damoulas. An optimization-centric view on Bayes’ rule: Reviewing and generalizing variational inference. Journal of Machine Learning Research, 23(132):1–109, 2022. URL https://jmlr.org/papers/v23/19-1047.html. Jonas Latz, Aretha L. Teckentrup, and Simon Urbainczyk. Sparse techniques for regression in deep Gaussian processes, 2025. URL https://arxiv.org/abs/2505.11355. arXiv:2505.11355. Steffen L. Lauritzen. Graphical Models, volume 17 of Oxford Statistical Science Series. Clarendon Press, Oxford, 1996. ISBN 9780198522195. Neil D. Lawrence and Andrew J. Moore. Hierarchical Gaussian process latent variable models. In Proceedings of the 24th International Conference on Machine Learning, pages 481–488, New York, NY, 2007. Association for Computing Machinery. doi: 10.1145/1273496.1273557. URL https://dl.acm.org/doi/10.1145/ 1273496.1273557. Loic Le Gratiet and Josselin Garnier. Recursive co-kriging model for design of computer experiments with multiple levels of fidelity. International Journal for Uncertainty Quantification, 4(5):365–386, 2014. doi: 10.1615/Int.J.UncertaintyQuantification.2014006914. URL https://doi.org/10.1615/Int.J. UncertaintyQuantification.2014006914. Naiqi Li, Wenjie Li, Jifeng Sun, Yinghua Gao, Yong Jiang, and Shu-Tao Xia. Stochastic deep Gaussian processes over graphs. In Advances in Neural Information Processing Systems, volume 33, pages 5875–5886, 2020. URL https://proceedings.neurips.cc/paper_files/paper/2020/file/ 415e1af7ea95f89f4e375162b21ae38c-Paper.pdf. Jakob Lindinger, David Reeb, Christoph Lippert, and Barbara Rakitsch. Beyond the mean-field: Structured deep Gaussian processes improve the predictive uncertainties. In Advances in Neural Information Processing Systems, volume 33, pages 8498–8509. Curran Associates, Inc., 2020. URL https://proceedings.neurips. cc/paper/2020/hash/60a70bb05b08d6cd95deb3bdb750dce8-Abstract.html.

Roderick J. A. Little and Donald B. Rubin. Statistical Analysis with Missing Data. John Wiley & Sons, Hoboken, NJ, third edition, 2019. ISBN 9781119482260. doi: 10.1002/9781119482260. Yuhao Liu, Marzieh Ajirak, and Petar M. Djurić. Gaussian process-gated hierarchical mixtures of experts. IEEE Transactions on Pattern Analysis and Machine Intelligence, 46(9):6443–6453, 2024. doi: 10.1109/ TPAMI.2024.3381936. Qiuxian Meng and Yongyou Zhang. Amortized variational inference for deep Gaussian processes. arXiv preprint arXiv:2409.12301, 2024. doi: 10.48550/arXiv.2409.12301. L. Mirsky. A dual of Dilworth’s decomposition theorem. The American Mathematical Monthly, 78(8):876–877, 1971. doi: 10.1080/00029890.1971.11992886. Joris M. Mooij, Sara Magliacane, and Tom Claassen. Joint causal inference from multiple contexts. Journal of Machine Learning Research, 21(99):1–108, 2020. URL https://jmlr.org/papers/v21/17-123.html. Radford M. Neal. Bayesian Learning for Neural Networks. PhD thesis, University of Toronto, Toronto, Canada, 1995. URL https://glizen.com/radfordneal/ftp/thesis.pdf. Sebastian W. Ober and Laurence Aitchison. Global inducing point variational posteriors for Bayesian neural networks and deep Gaussian processes. In Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 8248–8259. PMLR, 18–24 Jul 2021. URL https://proceedings.mlr.press/v139/ober21a.html. Anthony O’Hagan, Caitlin E. Buck, Alireza Daneshkhah, J. Richard Eiser, Paul H. Garthwaite, David J. Jenkinson, Jeremy E. Oakley, and Tim Rakow. Uncertain Judgements: Eliciting Experts’ Probabilities. John Wiley & Sons, Chichester, UK, 2006. ISBN 9780470029992. doi: 10.1002/0470033312. Judea Pearl. Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference. Morgan Kaufmann, San Francisco, CA, 1988. ISBN 9780080514895. doi: 10.1016/C2009-0-27609-4. Judea Pearl. Causality: Models, Reasoning, and Inference. Cambridge University Press, Cambridge, second edition, 2009. ISBN 9780521895606. doi: 10.1017/CBO9780511803161. Paris Perdikaris, Maziar Raissi, Andreas Damianou, Neil D. Lawrence, and George Em Karniadakis. Nonlinear information fusion algorithms for data-efficient multi-fidelity modelling. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2198):20160751, 2017. doi: 10.1098/rspa.2016.0751. Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of Causal Inference: Foundations and Learning Algorithms. Adaptive Computation and Machine Learning. The MIT Press, 2017. ISBN 9780262037310. URL https://mitpress.mit.edu/9780262037310/elements-of-causal-inference/. Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning. Adaptive Computation and Machine Learning. The MIT Press, Cambridge, MA, 2006. ISBN 9780262182539. doi: 10.7551/mitpress/3206.001.0001. Andreas Raue, Clemens Kreutz, Thomas Maiwald, Julie Bachmann, Marcel Schilling, Ursula Klingmüller, and Jens Timmer. Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics, 25(15):1923–1929, 2009. doi: 10.1093/bioinformatics/ btp358. Walter Rudin. Principles of Mathematical Analysis. McGraw–Hill, New York, third edition, 1976. Håvard Rue and Leonhard Held. Gaussian Markov Random Fields: Theory and Applications, volume 104 of Monographs on Statistics and Applied Probability. Chapman & Hall/CRC, Boca Raton, FL, 2005. ISBN 9780203492024. doi: 10.1201/9780203492024. Karen Sachs, Omar Perez, Dana Pe’er, Douglas A. Lauffenburger, and Garry P. Nolan. Causal proteinsignaling networks derived from multiparameter single-cell data. Science, 308(5721):523–529, 2005. doi: 10.1126/science.1105809. Erratum in: Science. 2005 Aug 19;309(5738):1187. 15

Hugh Salimbeni and Marc Deisenroth. Doubly stochastic variational inference for deep Gaussian processes. In Advances in Neural Information Processing Systems, volume 30, pages 4588– 4599. Curran Associates, Inc., 2017. URL https://proceedings.neurips.cc/paper/2017/hash/ 8208974663db80265e9bfe7b222dcb18-Abstract.html. Hugh Salimbeni, Vincent Dutordoir, James Hensman, and Marc Deisenroth. Deep Gaussian processes with importance-weighted variational inference. In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 5589–5598. PMLR, 2019. URL https://proceedings.mlr.press/v97/ salimbeni19a.html. Annie Sauer, Robert B. Gramacy, and David Higdon. Active learning for deep Gaussian process surrogates. Technometrics, 65(1):4–18, 2023. doi: 10.1080/00401706.2021.2008505. I. J. Schoenberg. Metric spaces and completely monotone functions. Annals of Mathematics, 39(4):811–841, 1938. doi: 10.2307/1968466. Bernhard Schölkopf, Francesco Locatello, Stefan Bauer, Nan Rosemary Ke, Nal Kalchbrenner, Anirudh Goyal, and Yoshua Bengio. Toward causal representation learning. Proceedings of the IEEE, 109(5):612–634, 2021. doi: 10.1109/JPROC.2021.3058954. URL https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber= 9363924. Marco Scutari. Learning bayesian networks with the bnlearn R package. Journal of Statistical Software, 35 (3):1–22, 2010. doi: 10.18637/jss.v035.i03. Peter Spirtes, Clark Glymour, and Richard Scheines. Causation, Prediction, and Search. The MIT Press, second edition, 2001. ISBN 9780262194402. URL https://mitpress.mit.edu/9780262194402/ causation-prediction-and-search/. Richard P. Stanley. Enumerative Combinatorics, Volume 1. Number 49 in Cambridge Studies in Advanced Mathematics. Cambridge University Press, Cambridge, second edition, 2012. ISBN 9781107602625. doi: 10.1017/CBO9781139058520. URL https://www.cambridge.org/core/books/ enumerative-combinatorics/3155CDE1D973D49F873BDE2EAF8D7651. Kozo Sugiyama, Shojiro Tagawa, and Mitsuhiko Toda. Methods for visual understanding of hierarchical system structures. IEEE Transactions on Systems, Man, and Cybernetics, 11(2):109–125, February 1981. doi: 10.1109/TSMC.1981.4308636. URL https://doi.org/10.1109/TSMC.1981.4308636. K. Takahashi, J. Fagan, and M.-S. Chin. Formation of a sparse bus impedance matrix and its application to short circuit study. In Proceedings of the 8th Power Industry Computer Applications Conference, pages 63–69, Minneapolis, MN, June 1973. IEEE Power Engineering Society. Robert E. Tarjan and Mihalis Yannakakis. Addendum: Simple linear-time algorithms to test chordality of graphs, test acyclicity of hypergraphs, and selectively reduce acyclic hypergraphs. SIAM Journal on Computing, 14(1):254, 1985. Anh Tong and Jaesik Choi. Characterizing deep Gaussian processes via nonlinear recurrence systems. Proceedings of the AAAI Conference on Artificial Intelligence, 35(11):9915–9922, 2021. doi: 10.1609/aaai. v35i11.17191. URL https://ojs.aaai.org/index.php/AAAI/article/view/17191. Ivan Ustyuzhaninov, Ieva Kazlauskaite, Markus Kaiser, Erik Bodin, Neill Campbell, and Carl Henrik Ek. Compositional uncertainty in deep Gaussian processes. In Jonas Peters and David Sontag, editors, Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence, volume 124 of Proceedings of Machine Learning Research, pages 480–489. PMLR, 2020. URL https://proceedings.mlr.press/v124/ ustyuzhaninov20a.html. Martin J. Wainwright and Michael I. Jordan. Graphical models, exponential families, and variational inference. Foundations and Trends in Machine Learning, 1(1–2):1–305, 2008. doi: 10.1561/2200000001.

David Vernon Widder. The Laplace Transform, volume 6 of Princeton Mathematical Series. Princeton University Press, Princeton, NJ, 1941. Sam Witty, Kenta Takatsu, David Jensen, and Vikash Mansinghka. Causal inference using Gaussian processes with structured latent confounders. In Hal Daumé III and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 10313–10323. PMLR, 13–18 Jul 2020. URL https://proceedings.mlr.press/v119/witty20a. html. Mihalis Yannakakis. Computing the minimum fill-in is NP-complete. SIAM Journal on Algebraic Discrete Methods, 2(1):77–79, 1981. doi: 10.1137/0602010. Xun Zheng, Bryon Aragam, Pradeep K. Ravikumar, and Eric P. Xing. DAGs with no tears: Continuous optimization for structure learning. In Advances in Neural Information Processing Systems, volume 31, pages 9472–9483, 2018. URL https://proceedings.neurips.cc/paper/2018/hash/ e347c51419ffb23ca3fd5050202f9c3d-Abstract.html.

Appendix Contents A Graph-theoretic preliminaries and standing notation A.1 Directed acyclic graphs and reachability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Chains, antichains, and progressivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3 Standing probabilistic notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.4 Routes below an antichain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19 19 19 20 21

B Prior non-collapse from separating nodes B.1 Progressive antichain decomposition for DAGs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2 A conditional non-degeneracy bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.3 Proof of Theorem 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.4 Sparse separating families . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.5 Sufficient conditions for separation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.6 Recovery of the chain input-connection mechanism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.7 Conditional refresh from intermediate supervision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21 22 23 25 27 29 31 32

C Indegree and outdegree effects C.1 Layered radial blocks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.2 Indegree effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.3 Outdegree effects through branching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33 33 35 38

D Intermediate observations as stochastic skip connections D.1 Filtering on an observed antichain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . D.2 Retaining routes below the observed antichain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . D.3 Proof of Theorem 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . D.4 Source strengths for common conditional distributions . . . . . . . . . . . . . . . . . . . . . . . . . . .

43 43 45 47 48

E Variational inference E.1 ELBO derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E.2 Marginal ancestral sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E.3 Practical ELBO evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E.4 Chordal completion and elimination order . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E.5 Computational cost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E.6 Explaining away . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52 52 53 55 59 59 60

F Stochastic Deep Gaussian Processes over Graphs as a Special Case of DAG-DGP F.1 The DGPG model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . F.2 DGPG as a particular case of DAG-DGP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

63 63 64

G Experiments G.1 Branching-tree ELBO scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . G.2 Theory Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . G.3 Latent-Collider experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . G.4 Sachs Flow Cytometry Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . G.5 Heavy-ion collision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67 67 68 70 72 73

Graph-theoretic preliminaries and standing notation

This appendix provides the necessary background for the theoretical developments which follow. For clarity and completeness, we introduce graph-theoretic and order-theoretic notions, setting the general notation and fixing the conventions adopted. Whenever possible, we follow standard conventions; see, for example, Bang-Jensen and Gutin (2009, Ch. 1) for directed graphs, directed paths, and acyclic digraphs, and Stanley (2012, Ch. 3) for the basic language of finite partially ordered sets, including chains and antichains.

A.1

Directed acyclic graphs and reachability

A directed graph is a pair G = (V, E), where V is a finite set of vertices and E ⊆ V × V is a set of directed edges. We write u → v when (u, v) ∈ E. A directed path from v0 to vL is a sequence (v0 , v1 , . . . , vL ) such that vr−1 → vr for every r = 1, . . . , L. Its length is L. Paths of length zero are allowed. The graph is a directed acyclic graph (DAG) if it contains no directed cycle of positive length. For a node w ∈ V, its parent and child sets are Pa(w) := {v ∈ V : v → w},

Ch(w) := {v ∈ V : w → v}.

A root is a node with no parents. In the main text the root set is denoted by R, and U := V \ R denotes the set of non-root nodes. We write u ⊑ v if there exists a directed path from u to v, allowing the length-zero path. We write u ⊏ v when u ⊑ v and u = ̸ v. Because G is acyclic, ⊑ is a partial order on V. For any S ⊆ V, throughout the appendix Anc(S) := {u ∈ V : ∃v ∈ S such that u ⊏ v} denotes the set of strict ancestors of S. We also use the shorthand Anc(v) := Anc({v}), and occasionally write Anc(S) := Anc(S) ∪ S for the weak ancestral closure of S. Similarly, Desc(S) := {u ∈ V : ∃v ∈ S such that v ⊏ u} denotes strict descendants.

A.2

Chains, antichains, and progressivity

A chain in the reachability order is a set of vertices any two of which are comparable. An antichain is a set of vertices no two distinct elements of which are comparable. Equivalently, A ⊆ V is an antichain if there is no directed path from one element of A to another distinct element of A. A finite progressive antichain decomposition of a DAG is a partition V=

h−1 G

Aℓ

ℓ=0

into h ∈ N non-empty antichains such that Aℓ ⊆ Anc(Aℓ+1 ),

ℓ = 0, . . . , h − 2.

For asymptotic statements we use an infinite progressive antichain sequence (Aℓ )ℓ≥0 , or finite truncations thereof, satisfying the same successive ancestry condition. The antichains should be read as successive cross-sections of the DAG rather than as independent layers. Figure 8 gives a visualization of a simple DAG partitioned into three antichains.

Figure 8: Progressive antichain sequence on an example DAG. Each highlighted row is an antichain: there are no directed paths between distinct nodes in the same row. Solid arrows show direct parent–child relations, whereas dashed arrows indicate longer reachability relations. The progressive condition only requires every node in Aℓ to be a strict ancestor of some node in Aℓ+1 .

A.3

Standing probabilistic notation

The DAG-DGP prior is the one defined in Section 2. Throughout, let (Ω, A , P) be a probability space supporting the mutually independent nodewise GP modules {fw : w ∈ U}, and let P denote the joint law induced by the DAG-DGP prior recursion. For an antichain Aℓ , we use the module-generated filtration Fℓ := σ(fv : v ∈ Anc(Aℓ )). All GP modules {fw : w ∈ U } are mutually independent under the prior probability measure P, and E denotes expectation under P. Root states are deterministic (i.e. almost surely constant), and no GP module is associated with a root. Once two distinct cases a ̸= b ∈ [n] are fixed, we write ∆w := Fw(a) − Fw(b) ,

w ∈ V,

(i) (i) recalling that Fw = fw {Fp }p∈Pa(w) , and here i ∈ {a, b}. For an antichain Aℓ , the maximum contrast at depth ℓ is Mℓ := max ∥∆w ∥2 . w∈Aℓ

Let P0 denote the joint law of the DAG-DGP prior, including the independent nodewise GP modules and the latent states induced by the DAG recursion. For a non-root node w and two cases a, b, say Uw,a := (Fp(a) )p∈Pa(w) ,

Uw,b := (Fp(b) )p∈Pa(w)

for the parent-input configurations generated by the DAG-DGP prior. The two-point contrast covariance at w is defined as Γw (a, b) := CovP Fw(a) − Fw(b) Uw,a , Uw,b . Equivalently, by the finite-dimensional distributions of the fresh GP module fw , Γw (a, b) = Kw (Uw,a , Uw,a ) + Kw (Uw,b , Uw,b ) − Kw (Uw,a , Uw,b ) − Kw (Uw,b , Uw,a ). A non-root node w is v⋆ -separating for (a, b) if Γw (a, b) ⪰ v⋆ Idw

P-almost surely.

This condition means that, for P-almost every pair of parent-input configurations generated by the DAG-DGP prior, the fresh two-point GP contrast at node w, evaluated at those inputs, has covariance bounded below. This is a pathwise conditional non-degeneracy condition along the DAG-DGP prior. It is meant to abstract mechanisms that refresh the deep composition by ensuring that, after the parent inputs of w have been 20

realised, the fresh nodewise GP module still sees a non-degenerate two-case contrast. In the architectures that we study in this work, the refresh is provided by structural coordinates or kernel components whose separation is not washed out by the upstream stochastic composition, but other such mechanisms can be considered in specialised modelling settings. Developing similar results under weaker conditions, perhaps under specific kernel choices and fusion rules, would require a different analysis and is left beyond the scope of the present work. Here and below ⪰ denotes the Loewner order on symmetric matrices. Notice that v⋆ separation is not a purely local property but imposes some requirements on the ancestors of the node in question: two distinct root input case configurations must propagate through the earlier part of the graph and remain sufficiently distinguished that they can be separated at this node. In this sense, although v⋆ -separation is a useful concept, establishing it in particular graph topologies can be non-trivial. We provide some examples in Appendix B.5; these cover a number of important cases. For a progressive antichain sequence (Aℓ )ℓ≥0 , the prior non-collapse proofs use the module-generated filtration with elements Fℓ := σ(fv : v ∈ Anc(Aℓ )) , ℓ ≥ 0. (10) The posterior-refresh section uses a different, state-generated notation. Hℓ records strict ancestral state matrices, whereas Fℓ also includes the current antichain. These are motivated and defined precisely in Appendix D.1.

A.4

Routes below an antichain

The stochastic-skip arguments use directed routes whose off-route parents are already known at the source antichain. We define these route notions once here. Definition 3 (Admissible route). Fix ℓ1 ≥ ℓ0 . A directed path γ = (v0 , v1 , . . . , vL ) with v0 ∈ Aℓ0 and vL ∈ Aℓ1 is admissible below Aℓ0 if, for each r = 1, . . . , L, every parent of vr is either vr−1 , a root node, or belongs to Anc(Aℓ0 ). Its interior is int(γ) := {v1 , . . . , vL }. If L = 0, the route is degenerate and its interior is empty. Definition 4 (Disjoint routes below an antichain). A family of admissible routes γj = (vj,0 , . . . , vj,Lj ), j ∈ [s], is pairwise disjoint below Aℓ0 if int(γj ) ∩ int(γk ) = ∅, j ̸= k. Figure 9 provides an illustration of the two definitions.

Prior non-collapse from separating nodes

This appendix uses concepts introduced in Appendix A to establish Theorem 1. We first verify that finite DAGs admit progressive antichain decompositions, so depth can be represented by successive antichain cross-sections. This allows us to define a rigorous notion of depth across the DAG. We then prove the key probabilistic lemma, which gives a uniform lower bound on the conditional probability that a separating node preserves a non-trivial two-case contrast, and combine this lemma with a martingale averaging argument to prove Theorem 1. After the main proof, we relax the uniform per-antichain assumption to sparse families of separating nodes, obtaining both an almost-sure non-collapse statement and a quantitative lower-frequency bound. The second part of the appendix identifies concrete mechanisms that produce separating nodes. We give sufficient conditions based on root-retaining additive and ANOVA-type kernel decompositions, including the multi-fidelity kernels used in deep graphical multi-fidelity GPs. We then specialise the general result to single-node antichains and use it to formalise the input-connection mechanism for DGP chains. In particular, 21

(a) admissible route

(b) disjoint routes

Anc(Aℓ0 ) ∪ R Aℓ0

A ℓ0 γ

Aℓ0 γ1

γ2

int(γ)

int(γ1 ) Aℓ1

Aℓ1

int(γ2 )

int(γ) Aℓ1

Figure 9: Route conventions below an antichain. In panel (a), the thick path is an admissible route γ from Aℓ0 to Aℓ1 ; its off-route parents are roots or belong to Anc(Aℓ0 ). In panel (b), two admissible routes are disjoint below Aℓ0 , since their interiors are disjoint. In panel (c), the displayed path is inadmissible: the highlighted node u1 is an additional parent of a route node, but it lies below Aℓ0 and is therefore not available at the source antichain. Corollary 2 proves, for the squared-exponential input-connected chain considered in Dunlop et al. (2018, Remark 5(3)), that the two-input contrast does not collapse almost surely. Finally, we show that exactly observed internal nodes can play the same anchoring role under the conditional prior: once conditioned on, their states act as deterministic coordinates for downstream kernels and yield the same non-collapse conclusions.

B.1

Progressive antichain decomposition for DAGs

We begin with a structural fact on finite DAGs that will be used extensively: Endowing the vertex set with the reachability order turns a finite DAG into a finite poset. Following the level-decomposition logic underlying Mirsky’s theorem (Mirsky, 1971), yields a partition of the vertex set into progressive antichains. Proposition 1 (Finite DAGs admit progressive antichain decompositions). Let G = (V, E) be a finite DAG. Then there exist an integer h ≥ 1 and non-empty antichains A0 , A1 , . . . , Ah−1 ⊆ V such that V=

h−1 G

Aℓ

and

Aℓ ⊆ Anc(Aℓ+1 ),

ℓ = 0, . . . , h − 2.

ℓ=0

Proof. For each v ∈ V, define n o d(v) := max k ≥ 0 : ∃ v1 , . . . , vk ∈ V such that v ⊏ v1 ⊏ · · · ⊏ vk . Thus d(v) is the maximum number of strict comparability steps in a chain starting at v. Since V is finite, d(v) is well defined for every v. Let h := 1 + max d(v), v∈V

Aℓ := {v ∈ V : d(v) = h − 1 − ℓ},

Then (Aℓ )h−1 ℓ=0 is a partition of V. We next show that each Aℓ is non-empty. Let v0 ⊏ v1 ⊏ · · · ⊏ vh−1 22

ℓ = 0, . . . , h − 1.

be a chain of maximum length h − 1. For each j = 0, . . . , h − 1, there exists a tail chain vj ⊏ vj+1 ⊏ · · · ⊏ vh−1 of length h − 1 − j, so d(vj ) ≥ h − 1 − j. Conversely, if d(vj ) ≥ h − j, then there would exist a chain of length at least h − j starting at vj , and prepending v0 ⊏ · · · ⊏ vj−1 ⊏ vj would yield a chain starting at v0 of length at least h, contradicting the maximality of h − 1. Hence d(vj ) = h − 1 − j,

j = 0, . . . , h − 1.

Therefore every value in {0, . . . , h − 1} is attained by some d(v), so every Aℓ is non-empty. We claim that each Aℓ is an antichain. Indeed, suppose u, v ∈ Aℓ and u ⊏ v. Then any chain of length d(v) starting at v can be prepended by u, yielding a chain of length at least d(v) + 1 starting at u. Hence d(u) ≥ d(v) + 1, contradicting d(u) = d(v) = h − 1 − ℓ. Thus distinct vertices in Aℓ are incomparable, so Aℓ is an antichain. It remains to prove the progressive property. Fix ℓ ∈ {0, . . . , h − 2} and let v ∈ Aℓ . Then d(v) = h − 1 − ℓ ≥ 1. By definition of d(v), there exists a chain v = v0 ⊏ v1 ⊏ · · · ⊏ vd(v) of length d(v) starting at v. For w := v1 , the tail w = v1 ⊏ v2 ⊏ · · · ⊏ vd(v) shows that d(w) ≥ d(v) − 1. Conversely, if d(w) ≥ d(v), then prepending v would yield a chain of length at least d(v) + 1 starting at v, contradicting the definition of d(v). Therefore d(w) = d(v) − 1 = h − 1 − ℓ − 1 = h − 1 − (ℓ + 1), so w ∈ Aℓ+1 . Since v ⊏ w, it follows that v ∈ Anc(w) ⊆ Anc(Aℓ+1 ). As v ∈ Aℓ was arbitrary, we conclude that Aℓ ⊆ Anc(Aℓ+1 ),

ℓ = 0, . . . , h − 2.

This proves the result. Remark 1. Proposition 1 provides a progressive antichain decomposition for every finite DAG. In particular, for a fixed finite DAG, every non-empty progressive antichain family is necessarily finite. Accordingly, asymptotic statements in Section 4.1 should be interpreted over increasing sequences of graphs whose common components coincide.

B.2

A conditional non-degeneracy bound

The proofs of Theorem 1 and its sparse generalisation both rest on the same one-step conditional bound, which we state and prove here as a self-contained lemma. Throughout this subsection, P denotes the probability measure induced by the DAG-DGP prior of Section 2, and E the corresponding expectation. We use the strict-ancestor convention, the contrast notation ∆w , the maximum Mℓ , the contrast covariance Γw (a, b), and the module filtration Fℓ from Appendix A.3.

Lemma 1 (Conditional antichain bound). Under the DAG-DGP prior of Section 2, fix distinct cases a= ̸ b ∈ [n] and a progressive antichain sequence (Aℓ )ℓ≥0 . Then (Fℓ )ℓ≥0 is a filtration. Moreover, for every ℓ ≥ 0 and every ε > 0, if Gℓ ⊆ Aℓ is a deterministic set of nodes that are v⋆ -separating for (a, b) with common constant v⋆ > 0, then P Mℓ ≤ ε | Fℓ ≤ (1 − pε )|Gℓ | P-almost surely, (11) where

√ pε := 2 ΦN (−ε/ v⋆ ) ∈ (0, 1).

(12)

Furthermore, 1{Mℓ > ε} is Fℓ+1 -measurable. Proof. Let u ∈ Anc(Aℓ ). Then there exists v ∈ Aℓ such that u is a strict ancestor of v. Since v ∈ Aℓ ⊆ Anc(Aℓ+1 ), there exists t ∈ Aℓ+1 such that v is a strict ancestor of t. By transitivity of the strict ancestor relation, u is a strict ancestor of t, hence u ∈ Anc(Aℓ+1 ). Therefore Anc(Aℓ ) ⊆ Anc(Aℓ+1 ). Consequently Fℓ ⊆ Fℓ+1 , so (Fℓ )ℓ≥0 is increasing and forms a filtration. We next establish a measurability fact that will be used repeatedly: for every node v ∈ V and every case (i) i ∈ [n], the state Fv is measurable with respect to σ(fu : u ∈ Anc(v) ∪ {v}) . (i)

(13)

(i)

If v ∈ R, then Fv = xv is deterministic, so the claim is immediate. If v ∈ U, choose any topological ordering of the DAG and argue by induction along that ordering. For each parent p ∈ Pa(v), the induction hypothesis (i) gives that Fp is measurable with respect to σ(fu : u ∈ Anc(p) ∪ {p}) . Since every ancestor of p is also an ancestor of v, and p itself is a strict ancestor of v, one has Anc(p) ∪ {p} ⊆ Anc(v). (i)

Hence each parent state Fp

is measurable with respect to σ(fu : u ∈ Anc(v)). Because Fv(i) = fv {Fp(i) }p∈Pa(v) ,

(i)

it follows that Fv is measurable with respect to (13), as claimed. As a first consequence, we verify that Mℓ is Fℓ+1 -measurable. Let v ∈ Aℓ and i ∈ [n]. Since v ∈ Aℓ ⊆ Anc(Aℓ+1 ), the node v is a strict ancestor of Aℓ+1 . Also, every strict ancestor of v is a strict ancestor of Aℓ+1 , so Anc(v) ∪ {v} ⊆ Anc(Aℓ+1 ). (i)

By (13), Fv is therefore Fℓ+1 -measurable. Since this holds for every v ∈ Aℓ , the random variable Mℓ is Fℓ+1 -measurable, and so is 1{Mℓ > ε}. We now turn to the conditional bound. Fix ℓ ≥ 0 and w ∈ Gℓ ⊆ Aℓ . Recall Uw,a := (Fp(a) )p∈Pa(w) , Uw,b := (Fp(b) )p∈Pa(w) . If p ∈ Pa(w), then p is a strict ancestor of w, and since w ∈ Aℓ , it follows that p ∈ Anc(Aℓ ). Moreover, every ancestor of p is also a strict ancestor of w, hence of Aℓ , so Anc(p) ∪ {p} ⊆ Anc(Aℓ ). (a)

(b)

By (13), Fp and Fp are therefore Fℓ -measurable. Root coordinates are deterministic by construction. Hence the parent-input configurations Uw,a , Uw,b are Fℓ -measurable. Since the GP modules are mutually independent under P, and since w ∈ / Anc(Aℓ ) by the antichain property of Aℓ , the fresh module fw is independent of Fℓ . Hence, under the conditional law given Fℓ , the parent 24

inputs of fw are fixed at Uw,a , Uw,b , while the remaining randomness comes only from fw . Therefore ∆w is conditionally centred Gaussian with variance Vw := CovP (∆w | Fℓ ) = Γw (a, b). The last equality follows from the definition of Γw (a, b) and the Fℓ -measurability of Uw,a , Uw,b . Since w is v⋆ -separating for (a, b), we have Vw ⪰ v⋆ Idw (P-almost surely). Fix any deterministic unit vector ew ∈ Rdw . By the Cauchy–Schwarz inequality, |e⊤ w ∆w | ≤ ∥ew ∥2 ∥∆w ∥2 = ∥∆w ∥2 . Therefore Conditional on Fℓ ,

P ∥∆w ∥2 ≤ ε | Fℓ ≤ P |e⊤ w ∆w | ≤ ε | F ℓ .

(14)

2 e⊤ w ∆w ∼ N (0, σw ),

(15)

2 σw = e⊤ w Vw ew ≥ v⋆ .

Let Z ∼ N (0, 1). Then P |e⊤ w ∆w | > ε | Fℓ = P |Z| > ε/σw . Since the map σ 7→ P(|Z| > ε/σ) is increasing on (0, ∞), (15) implies √ P |e⊤ w ∆w | > ε | Fℓ ≥ P |Z| > ε/ v⋆ = pε .

(16)

Combining (14) and (16) gives P ∥∆w ∥2 ≤ ε | Fℓ ≤ 1 − pε .

(17)

We now pass from a single node to the whole antichain. The family {∆w }w∈Gℓ is conditionally independent given Fℓ . Indeed, for each w ∈ Gℓ , the variable ∆w is obtained by evaluating the module fw at the Fℓ measurable inputs ua , ub and then taking a difference, so conditional on Fℓ it is a measurable function of fw alone. Since the modules {fw }w∈Gℓ are mutually independent under the prior and each is independent of Fℓ , the conditional independence follows. Using conditional independence and (17), ! \ Y P Mℓ ≤ ε | F ℓ ≤ P {∥∆w ∥2 ≤ ε} Fℓ = P ∥∆w ∥2 ≤ ε | Fℓ ≤ (1 − pε )|Gℓ | . w∈Gℓ

w∈Gℓ

This establishes (11).

B.3

Proof of Theorem 1

Proof of Theorem 1. Fix distinct cases a = ̸ b ∈ [n] and an antichain sequence (Aℓ )ℓ≥0 with Aℓ ⊆ Anc(Aℓ+1 ). Let pε be as in (12). By hypothesis, every antichain Aℓ contains at least s ≥ 1 nodes that are v⋆ -separating for (a, b). Fix once and for all a deterministic ordering of the vertices, and let Gℓ be the first s separating nodes in Aℓ under this ordering. Applying Lemma 1 gives P Mℓ > ε | Fℓ ≥ 1 − (1 − pε )s a.s. (18) Define Dℓ := 1{Mℓ > ε} − E[1{Mℓ > ε}|Fℓ ], which is a bounded martingale difference sequence with respect to the filtration (Fl+1 )l≥0 . Indeed, Lemma 1 tells us that the first term is Fl+1 -measurable, and hence Dl . Also, E[Dℓ |Fℓ ] = 0 and |Dℓ | ≤ 1 follow by construction. The partial sums Sm :=

m−1 X

Dℓ ,

ℓ=0

m ≥ 1,

with S0 = 0, form a martingale with respect to (Fm )m≥1 . Indeed, Sm is Fm -measurable and "m # X E[Sm+1 | Fm ] = E Dℓ Fm ℓ=0 m−1 X

Dℓ + E[Dm | Fm ] = Sm .

ℓ=0

Furthermore, |Sm − Sm−1 | = |Dm−1 | ≤ 1 almost surely. By the Azuma–Hoeffding inequality (Azuma, 1967; Hoeffding, 1963), for every η > 0 and every m ≥ 1, 2 2 2 η m η m P |Sm − S0 | ≥ ηm ≤ 2 exp − = 2 exp − . (19) 2m 2 We now apply the first Borel–Cantelli lemma. Fix η > 0 and define Am := {|Sm | ≥ ηm},

m ≥ 1.

From (19), ∞ X

∞ X

2 η m 2 exp − P(Am ) ≤ < ∞. 2 m=1 m=1 Hence

P lim sup Am = 0. m→∞

Equivalently, for this fixed η, there exists an almost surely finite random integer Mη such that |Sm | < ηm

for all m ≥ Mη .

(20)

For each fixed integer k ≥ 1, apply (20) with η = 1/k. Thus there is an event Ek ⊆ Ω with P(Ek ) = 1 such that, for every outcome ω ∈ Ek , there exists an integer Mk (ω) < ∞ satisfying |Sm (ω)| <

m k

for all m ≥ Mk (ω).

Define E :=

∞ \

Ek .

(21)

(22)

k=1

Since the intersection in (22) is countable and each Ek has probability one, P(E) = 1.

(23)

We now prove convergence on this probability-one event. Fix an outcome ω ∈ E, and let δ > 0 be arbitrary. Choose an integer k > 1/δ. Since ω ∈ E ⊆ Ek , there exists an integer Mk (ω) < ∞ such that (21) holds. Hence, for every m ≥ Mk (ω), Sm (ω) 1 < < δ. m k Because δ > 0 was arbitrary, this proves that, for every ω ∈ E, Sm (ω) →0 m

as m → ∞.

(24)

P-almost surely.

(25)

Combining (23) and (24), we obtain Sm →0 m

Taking lower limits as m → ∞, using (25) and the lower bound (18), gives m−1

m−1

Sm 1 X P Mℓ > ε | F ℓ + m m

1 X lim inf 1{Mℓ > ε} = lim inf m→∞ m m→∞ ℓ=0

ℓ=0

1 m→∞ m

≥ lim inf

m−1 X

1 − (1 − pε )s

ℓ=0 = 1 − (1 − pε )s

B.4

almost surely.

Sparse separating families

Theorem 1 generalises to families of separating nodes that may be unevenly distributed across depths: Proposition 2 (Non-collapse from sparse separating families). Under the DAG-DGP prior of Section 2, fix distinct cases a ̸= b ∈ [n] and a progressive antichain sequence (Aℓ )ℓ≥0 . Let Gℓ ⊆ Aℓ be a deterministic set of v⋆ -separating nodes at each depth ℓ. P (i) If ℓ≥0 |Gℓ | = ∞, then P(Mℓ → 0 as ℓ → ∞) = 0, and for every fixed ε > 0, P Mℓ > ε for infinitely many ℓ = 1. (ii) For s ∈ N, let Ls := {ℓ ≥ 0 : |Gℓ | ≥ s}. Then, for every ε > 0, m−1

1 X |Ls ∩ {0, . . . , m − 1}| 1 {Mℓ > ε} ≥ 1 − (1 − pε )s lim inf m→∞ m m→∞ m

lim inf

a.s.,

ℓ=0

where pε is defined in (12). Proof. By Lemma 1, P(Mℓ ≤ ε | Fℓ ) ≤ (1 − pε )|Gℓ |

a.s.

(26)

For part (i), fix ε > 0, n ≥ 0, and N ≥ n, and define En,N :=

N \

{Mℓ ≤ ε}.

ℓ=n

Since Mℓ is Fℓ+1 -measurable by Lemma 1, the event En,N −1 is FN -measurable. Therefore, using the tower property and (26), P(En,N ) = E 1En,N −1 1{MN ≤ε} = E 1En,N −1 P(MN ≤ ε | FN ) ≤ (1 − pε )|GN | P(En,N −1 ). (27) Iterating (27) from N down to n yields P(En,N ) ≤

N Y

(1 − pε )|Gℓ | = (1 − pε )

PN ℓ=n

|Gℓ |

ℓ=n

Letting N → ∞ and using continuity from above gives PN P Mℓ ≤ ε for all ℓ ≥ n ≤ lim (1 − pε ) ℓ=n |Gℓ | . N →∞

(28)

P PN If ℓ≥0 |Gℓ | = ∞, then for every n ≥ 0 the tail sum ℓ=n |Gℓ | → ∞ as N → ∞, so the right-hand side of (28) is zero. Hence P Mℓ ≤ ε for all ℓ ≥ n = 0 for every n ≥ 0. Equivalently, P Mℓ > ε for infinitely many ℓ = 1.

(29)

Since (29) holds for any fixed ε > 0, it follows in particular that P(Mℓ → 0) = 0. For part (ii), fix s ∈ N and define Yℓ := 1{ℓ∈Ls } 1{Mℓ > ε},

Dℓ := Yℓ − E[Yℓ | Fℓ ].

Because Ls is deterministic and 1{Mℓ > ε} is Fℓ+1 -measurable, Yℓ is Fℓ+1 -measurable. Also, |Dℓ | ≤ 1 almost surely and E[Dℓ | Fℓ ] = 0. Thus the partial sums Sm :=

m−1 X

Dℓ

(30)

ℓ=0

form a martingale with respect to (Fm )m≥1 with bounded increments. Applying the pathwise conclusion proved in (23)–(25) to the martingale in (30), there exists an event Esparse ⊆ Ω with P(Esparse ) = 1 such that, for every ω ∈ Esparse , m−1 Sm (ω) 1 X Dℓ (ω) = −→ 0 as m → ∞. m m ℓ=0

Equivalently, m−1

1 X Dℓ → 0 m

P-almost surely as m → ∞.

(31)

ℓ=0

For every ℓ ∈ Ls , one has |Gℓ | ≥ s, so by (26), P(Mℓ > ε | Fℓ ) ≥ 1 − (1 − pε )s Therefore,

a.s.

E[Yℓ | Fℓ ] = 1{ℓ∈Ls } P(Mℓ > ε | Fℓ ) ≥ 1{ℓ∈Ls } 1 − (1 − pε )s .

Using 1{Mℓ > ε} ≥ Yℓ and the definition of Dℓ , we obtain m−1

m−1

1 X 1 X 1{Mℓ > ε} ≥ Yℓ m m ℓ=0

ℓ=0

1 m

m−1 X

m−1

E[Yℓ | Fℓ ] +

ℓ=0

1 X Dℓ m ℓ=0

m−1 |Ls ∩ {0, . . . , m − 1}| 1 X + Dℓ . ≥ 1 − (1 − pε )s m m ℓ=0

Taking lower limits as m → ∞ and using (31) proves the claim. Theorem 1 assumes that every antichain Aℓ contains at least s separating nodes, leading to a uniform positive lower bound on the asymptotic frequency of depths where the maximum contrast exceeds ε. In many DAGs, separating nodes may individually be less v⋆ - separating. Proposition 2 relaxes this uniformity. Part (i) shows that it suffices that the total number of separating nodes across all depths is infinite: then the prior does not collapse, in the sense that Mℓ does not converge to zero, and for every fixed ε > 0, the event Mℓ > ε occurs infinitely often almost surely. Part (ii) provides a quantitative refinement: the empirical frequency of large contrasts is bounded below by the product of the per-depth guarantee 1 − (1 − pε )s and the lower asymptotic density of depths that contain at least s separating nodes. Hence the non-collapse phenomenon persists even when separating nodes appear only sporadically, as long as the set of such depths has positive lower density. When every antichain contains at least s separating nodes, the lower density equals 1 and the bound reduces to that of the main proposition, so the latter is a special, more intuitive, case of the more general result. 28

B.5

Sufficient conditions for separation

A convenient class of fusion kernels that automatically produce separating nodes is given by those retaining a uniformly non-trivial root-only main effect. This provides the natural link between the separating condition required for prior non-collapse and the connection of intermediate latent nodes to the roots of the DAG: the natural generalization of input connection from DGPs to DAG-DGPs. Proposition 3 (ANOVA root retention implies separation). Fix a non-root node w and split its input as (x, z), where x collects one or more root coordinates. For the pair (a, b), write xa , xb for the corresponding values of these root-coordinate components. Assume that the kernel at w has the form   Jw Mw X X (0) Kw (x, z), (x′ , z ′ ) = Bw  βw,j κw,j (x, x′ ) + αw,m ρw,m (x, z), (x′ , z ′ )  , m=1

j=1

with Bw ⪰ σ 2 Idw , nonnegative weights satisfying Jw X

βw,j +

Mw X

αw,m = 1,

m=1

j=1 (0)

and scalar correlation kernels κw,j and ρw,m . If Jw X

βw,j ≥ γ0 > 0

(0)

sup κw,j (xa , xb ) ≤ r < 1,

and

j=1

then w is v⋆ -separating for (a, b) with v⋆ := 2 σ 2 1 − r̄ ,

r̄ := (1 − γ0 ) + γ0 r < 1.

Proof. For the pair (a, b), write ua = (xa , za ), and define ϑw :=

Jw X

ub = (xb , zb ), Mw X

(0)

βw,j κw,j (xa , xb ) +

αw,m ρw,m (ua , ub ).

m=1

j=1

(0)

Jw X

(0)

βw,j |κw,j (xa , xb )| +

Jw X

βw,j +

Jw X

Mw X

αw,m

m=1

j=1

αw,m |ρw,m (ua , ub )|

m=1

j=1

≤r

Mw X

βw,j + 1 −

j=1

= 1 − (1 − r)

Jw X

βw,j

j=1 Jw X

βw,j

j=1

≤ 1 − γ0 (1 − r) = (1 − γ0 ) + γ0 r = r̄. (0)

Moreover, since κw,j and ρw,m are scalar correlation kernels, evaluating on the diagonal gives Kw (ua , ua ) = Kw (ub , ub ) = Bw . 29

The cross-covariances are Kw (ua , ub ) = Bw ϑw ,

Kw (ub , ua ) = Bw ϑw .

Hence Γw (a, b) = 2(1 − ϑw )Bw . Since ϑw ≤ |ϑw | ≤ r̄ < 1, it follows that Γw (a, b) ⪰ 2(1 − r̄)Bw ⪰ 2 σ 2 (1 − r̄)Idw = v⋆ Idw . Thus w is v⋆ -separating. The simple root-retention decomposition used in the motivating discussion of Section 4.1 is recovered as the special case Jw = 1. The following remark shows that the same mechanism applies more broadly, including to multi-fidelity kernels that do not need to be written as a normalised convex combination. Remark 2 (Separation from an additive root-only kernel component). Suppose the kernel at a non-root node w decomposes as (int) (0) Kw (x, z), (x′ , z ′ ) = Kw (x, z), (x′ , z ′ ) + Kw (x, x′ ), (int)

(0)

where Kw is a valid matrix-valued positive-semidefinite kernel on the full input space and Kw is a valid matrix-valued positive-semidefinite kernel depending only on root coordinates. For the pair (a, b), write xa , xb for the corresponding values of these root-coordinate components. Since both summands are valid kernels, the contrast covariance decomposes as Γw (a, b) = Γ(int) (a, b) + Γ(0) w w (a, b), with each summand positive semidefinite. In particular, Γw (a, b) ⪰ Γ(0) w (a, b), where (0) (0) (0) (0) Γ(0) w (a, b) = Kw (xa , xa ) + Kw (xb , xb ) − Kw (xa , xb ) − Kw (xb , xa ).

(32)

If, in particular, (0) Kw (x, x′ ) = σ02 k̄(x, x′ )Idw ,

where k̄ is a scalar correlation kernel, then (32) becomes 2 Γ(0) w (a, b) = 2σ0 1 − k̄(xa , xb ) Idw . Therefore, whenever k̄(xa , xb ) < 1, Γw (a, b) ⪰ 2σ02 1 − k̄(xa , xb ) Idw =: v⋆ Idw , with v⋆ > 0, and w is separating. This condition holds, for example, for a squared-exponential root kernel whenever xa ̸= xb . This criterion applies to the multi-fidelity deep GP kernel of Cutajar et al. (2019) and its graphical generalisation in Ji et al. (2024). In their formulation, the kernel at each non-source node takes the form Kw [x, z], [x′ , z ′ ] = KSE,ρ (x, x′ ) KLIN (z, z ′ ) + KSE (z, z ′ ) + KSE,δ (x, x′ ), where z collects the parent-node outputs and x is the exogenous input, and the subscripts SE and LIN denote the squared-exponential and linear kernel, respectively. The discrepancy kernel KSE,δ (x, x′ ) depends only on root coordinates and is a squared-exponential kernel with variance σδ2 . Since KSE,δ (xa , xb ) < σδ2 whenever xa ̸= xb , the above argument gives #! " X v⋆ = 2σδ2 1 − exp − (xa,l − xb,l )2 /(2λ2l ) > 0, l

where the λl are the length-scale parameters of KSE,δ . Thus every non-source node in a graphical multi-fidelity DGP with this kernel is separating for any pair of distinct input cases, regardless of the parent latent states. 30

B.6

Recovery of the chain input-connection mechanism

The results of Section B.5 refer directly to the empirical evidence in the literature that connecting all latent nodes to the input space prevents the prior-collapse pathology. This observation goes back to Duvenaud et al. (2014), who provided the first theoretical result on prior collapse and investigated in detail the proposal, originally elaborated by Neal (1995, Chapter 2) as a general recommendation for arbitrarily deep Bayesian neural network, of connecting every latent layer of a DGP to the corresponding input. In their theoretical paper on this pathological prior behaviour, after formalising and proving that the pathology occurs for non-input-connected DGPs under RBF kernels, Dunlop et al. (2018, Remark 5(3)) conjecture that the same collapse result does not hold when every layer is input-connected, and give an intuition for this. With our general theory, we can formalise and prove that conjectured non-collapse mechanism as a special case. The following corollary shows that uniform v⋆ -separation is sufficient to avoid prior collapse in the case of chain DAGs (an unsurprising fact), whereas the subsequent corollary demonstrates that skip-connection is enough to ensure this property and hence non-collapse. Corollary 1 (Separating input-connected chains). Assume that Aℓ = {vℓ } for every ℓ ≥ 0, and that each vℓ is v⋆ -separating for the pair (a, b), with common constant v⋆ > 0. Then, for every ε > 0, m−1 o 1 X n (a) 1 ∥Fvℓ − Fv(b) ∥ > ε ≥ pε 2 ℓ m→∞ m

lim inf

a.s.,

ℓ=0

√ where pε = 2ΦN (−ε/ v⋆ ) ∈ (0, 1). Moreover, (b) P ∥Fv(a) − F ∥ > ε for infinitely many ℓ =1 vℓ 2 ℓ

for every ε > 0,

and in particular (b) P ∥Fv(a) − F ∥ → 0 as ℓ → ∞ = 0. 2 vℓ ℓ Proof. Set Gℓ := {vℓ }, so that |Gℓ | = 1 for every ℓ. Then

ℓ≥0 |Gℓ | = ∞, and Proposition 2 (i) gives

P ∥Fv(a) − Fv(b) ∥2 > ε for infinitely many ℓ = 1 ℓ ℓ

for every ε > 0.

Moreover, part (ii) with s = 1 and lower density equal to 1 yields m−1 o 1 X n (a) 1 ∥Fvℓ − Fv(b) ∥2 > ε ≥ pε ℓ m→∞ m

lim inf

a.s.

ℓ=0

(a)

(b)

(a)

(b)

Finally, if ∥Fvℓ − Fvℓ ∥2 → 0 as ℓ → ∞, then for every ε > 0 one must have ∥Fvℓ − Fvℓ ∥2 ≤ ε for all sufficiently large ℓ, which excludes the event that the norm exceeds ε for infinitely many ℓ. Since that event has probability one, the convergence-to-zero event has probability zero. Corollary 2 (Input-connected chain of Dunlop et al. (2018)). Consider the input-connected chain in Dunlop et al. (2018, Remark 5(3)), un+1 (x) = ξn+1 (un (x), x), 1 1 m m where ξn+1 = (ξn+1 , . . . , ξn+1 ) and the scalar fields ξn+1 , . . . , ξn+1 are independent copies of a centred Gaussian m d process on R × R with squared-exponential kernel ∥z − z ′ ∥22 + ∥x − x′ ∥22 h (z, x), (z ′ , x′ ) = σ 2 exp − . 2w2

Let x, x0 ∈ D with x ̸= x0 . Then, for every ε > 0, N −1

1 X 1{∥un (x) − un (x0 )∥2 > ε} ≥ pε N →∞ N n=0

lim inf

a.s.,

where



 ε

 pε = 2 ΦN − q ∈ (0, 1). 2σ 2 1 − exp −∥x − x0 ∥22 /(2w2 ) Moreover, P(∥un (x) − un (x0 )∥2 → 0 as n → ∞) = 0. Proof. This is an application of Corollary 1. The model is a chain, so the antichains are the singletons An = {vn }. Fix n ≥ 0 and condition on un . Set z := un (x) and z0 := un (x0 ). The matrix-valued kernel is K (z, x), (z ′ , x′ ) = h (z, x), (z ′ , x′ ) Im . Hence K (z, x), (z, x) = K (z0 , x0 ), (z0 , x0 ) = σ 2 Im , and

∥z − z0 ∥22 + ∥x − x0 ∥22 Im . K (z, x), (z0 , x0 ) = σ exp − 2w2

Therefore the conditional contrast covariance is Γn (x, x0 ) = K (z, x), (z, x) + K (z0 , x0 ), (z0 , x0 ) − K (z, x), (z0 , x0 ) − K (z0 , x0 ), (z, x) ∥z − z0 ∥22 + ∥x − x0 ∥22 2 = 2σ 1 − exp − Im . 2w2 Since ∥z − z0 ∥22 ≥ 0, the exponential is bounded above by exp(−∥x − x0 ∥22 /(2w2 )), so ∥x − x0 ∥22 2 Γn (x, x0 ) ⪰ 2σ 1 − exp − Im =: v⋆ Im . 2w2 Since x = ̸ x0 , one has v⋆ > 0, so every layer is v⋆ -separating for (x, x0 ). Applying Corollary 1 gives both the positive lower-frequency bound and the almost-sure non-collapse.

B.7

Conditional refresh from intermediate supervision

The following corollary formalises the claim from Section 4.1 that fully observed and noiseless internal nodes stabilise the conditional prior without any architectural modification. An observed internal node can be used by downstream kernels in the same way that input-connected kernels use root inputs. Accordingly, throughout this subsection let E ⊆ U be a set of non-root nodes that are observed without error, in the sense that Oe = [n] × [de ] and the realised observed value is Ye = Fe for every e ∈ E. Since exact observations of continuous latent variables should be interpreted through regular conditional laws, we write PE ( · ) := P( · |σ(Fe : e ∈ E)) evaluated at {Fe = Ye : e ∈ E} for a chosen regular conditional version of the prior given the exactly observed internal states. Under PE , (i) (i) each conditioned state Fe is almost surely equal to the observed value Ye . Hence any downstream kernel that depends on coordinates coming from nodes in E sees them as fixed anchor coordinates. Corollary 3 (Perfect intermediate supervision as conditional refresh). Fix distinct ̸ b ∈ [n], and let S cases a = (Aℓ )ℓ≥0 be antichains in U \ E with Aℓ ⊆ Anc(Aℓ+1 ). Assume that no node in ℓ≥0 Aℓ is an ancestor of any node in E. Let Gℓ ⊆ Aℓ be a deterministic family of nodes such that, under the conditional prior PE , each w ∈ Gℓ satisfies the assumptions of Proposition 3, with anchor coordinates taken from E, and with a common resulting lower bound v⋆ > 0. Then each w ∈ Gℓ is v⋆ -separating for (a, b) under PE , and the conclusions of Proposition 2 apply under the conditional prior PE . In particular, if every antichain contains at least s ≥ 1 such nodes, the conclusion of Theorem 1 applies under PE . 32

S Proof. For every node w ∈ ℓ≥0 Gℓ , any input split (c, z) in which c collects coordinates from nodes in E becomes, under PE , an input split with fixed anchor coordinates. These coordinates therefore play the same role that root inputs play in Proposition 3. By assumption, each w ∈ Gℓ satisfies the hypotheses of Proposition 3 with these conditioned anchor coordinates in place of the root coordinates. Proposition 3 therefore applies under PE , and the common lower-bound assumption yields that every such node is v⋆ -separating for (a, b) under the conditional prior. It remains to verify that the conditional independence structure required by Lemma 1 is preserved under PE . Under the prior P, the exactly observed variables are measurable with respect to the modules ( ) [ ME := fv : v ∈ Anc(e) ∪ {e} . (33) e∈E

By hypothesis, no node in and w ∈ / E, it follows that

ℓ≥0 Aℓ is an ancestor of any node in E. Since every w ∈ Aℓ is a non-root node

w∈ /

[

(Anc(e) ∪ {e}).

(34)

e∈E

Thus fw is one of the prior-independent GP modules outside ME . Now consider the filtration (Fℓ )ℓ≥0 from (10) and a separating node w ∈ Gℓ ⊆ Aℓ . Under P, the module fw is jointly independent of the sigma-field generated by {fv : v ∈ Anc(Aℓ )} ∪ ME . Indeed, w ∈ / Anc(Aℓ ) by the antichain argument in the proof of Lemma 1, and w is not contained in the module index set (33) by (34). Mutual independence of the GP modules therefore gives the joint independence. Consequently, after conditioning on the realised values of {Fe : e ∈ E}, the module fw remains independent of Fℓ and retains its prior GP law. Similarly, for distinct w, w′ ∈ Gℓ , the modules fw and fw′ remain conditionally independent given Fℓ under PE . Therefore every argument in the proof of Lemma 1 applies under PE without modification. Hence Proposition 2 holds under the conditional prior, and Theorem 1 holds under PE when its uniform per-antichain assumption is satisfied.

Indegree and outdegree effects

This section isolates how local graph degrees affect two-case contrasts when the fusion kernel is radial in the concatenated parent state. The indegree analysis shows that, in product-type radial blocks, aggregating many parents can attenuate expected squared contrasts and yields a local contraction criterion. The outdegree analysis gives the complementary mechanism: if a node has sufficiently many disjoint designated children, branching can sustain threshold-size contrasts with positive probability. These results are local to the analysed block and do not alter the separating-node non-collapse criterion above; additive root-retaining components remain separating regardless of the number of additional parents. (a) (b) We continue with the notation of Appendix A; in particular, ∆w = Fw − Fw for a fixed pair of distinct cases a ̸= b ∈ [n].

C.1

Layered radial blocks

To isolate the effect of indegree, we work on a layered block of the DAG in which all parents of a node in the current antichain lie in the immediately preceding antichain, and the nodewise fusion kernel is radial in the concatenated parent state. This is a local structural assumption on the block being analysed, not a global restriction on the whole DAG. Assumption 1 (Layered Laplace–radial block). Let (Aℓ )ℓ≥0Sbe non-root antichains such that Pa(w) ⊆ Aℓ−1 for every w ∈ Aℓ and every ℓ ≥ 1. Assume that all nodes in ℓ≥0 Aℓ have the same output dimension d, and that there exist τ > 0 and a probability measure ν on [0, ∞) such that Z ∞ Kw (u, u′ ) = τ 2 κ(∥u − u′ ∥22 )Id , κ(r) = e−sr ν(ds), 0

for every node w ∈

ℓ≥0 Aℓ .

A visualisation of such a layered block is given in Figure 10.

Aℓ−1

Aℓ

Aℓ+1

Figure 10: Layered block as in Assumption 1. Each row is an antichain, and all displayed edges in the analysed block run from one layer to the next. The highlighted node w ∈ Aℓ illustrates the local indegree mechanism: all its parents lie in Aℓ−1 , and its descendants lie in Aℓ+1 . Dashed arrows indicate omitted portions of the surrounding DAG above and below the displayed block. By Schoenberg’s theorem (Schoenberg, 1938), a continuous function κ : [0, ∞) → R with κ(0) = 1 yields a radial kernel (x, x′ ) 7→ κ(∥x − x′ ∥22 ) that is positive definite on Rm for every m if and only if κ is completely monotone. By the Bernstein–Widder theorem (Widder, 1941), this is equivalent to the Laplace-transform representation Z ∞

e−sr ν(ds)

κ(r) = 0

for a probability measure ν on [0, ∞). Therefore Assumption 1 covers isotropic radial kernels that are valid in every ambient dimension, including the squared-exponential and rational-quadratic families. In particular, κ is nonincreasing, takes values in [0, 1], with κ(0) = 1. The homogeneity across nodes is imposed only for notational simplicity; a layer-dependent version is stated in Remark 3 below. A relevant special case is product fusion of squared-exponential parent kernels. If ! ′ 2 ∥u − u ∥ p p 2 (p) Kw (up , u′p ) = exp − 2λ2 for every p ∈ Pa(w), and the fused kernel is normalised so that Kw (u, u) = τ 2 Id , then ! Y ∥up − u′p ∥22 ∥u − u′ ∥22 ′ 2 2 Kw (u, u ) = τ exp − Id = τ exp − Id , 2λ2 2λ2 p∈Pa(w)

where u = (up )p∈Pa(w) denotes the concatenated parent state. Hence Assumption 1 holds with ν = δ1/(2λ2 ) . For each non-root node w we define the squared parent contrast X ϱ2w = ∥∆p ∥22 , (35) p∈Pa(w)

which is the squared Euclidean distance between the concatenated parent states evaluated at cases a and b. Lemma 2 (Two-point law for a single radial module). Let w ∈ U be a non-root node whose kernel has the form Kw (u, u′ ) = τ 2 κ(∥u − u′ ∥22 )Id , and let ϱ2w be the squared parent contrast defined in (35). Conditional on the parent states of w, ∆w ∼ N 0, 2τ 2 1 − κ(ϱ2w ) Id . In particular, h i E ∥∆w ∥22 | {Fp(a) , Fp(b) }p∈Pa(w) = 2dτ 2 1 − κ(ϱ2w ) . 34

(36)

Proof. By the conditional Gaussianity of ∆w established in the proof of Lemma 1, the contrast ∆w is conditionally centred Gaussian with covariance Γw (a, b) given the parent states. Under the radial kernel assumption, the diagonal evaluations give Kw (ua , ua ) = Kw (ub , ub ) = τ 2 Id , and the cross-evaluations give Kw (ua , ub ) = Kw (ub , ua ) = τ 2 κ(ϱ2w )Id . Hence Γw (a, b) = 2τ 2 1 − κ(ϱ2w ) Id . For the second-moment identity, write Σw = 2τ 2 (1 − κ(ϱ2w ))Id for the conditional covariance. Since ∆w is conditionally centred, h i h i (a) (b) E ∥∆w ∥22 | {Fp(a) , Fp(b) }p∈Pa(w) = E ∆⊤ w ∆w | {Fp , Fp }p∈Pa(w) = tr(Σw ) = 2dτ 2 1 − κ(ϱ2w ) , via the standard identity E[Z ⊤ Z] = tr(Var(Z)) for any centred random vector Z.

C.2

Indegree effects

Throughout this subsection, recall from Section 4.2 that kℓ = maxw∈Aℓ |Pa(w)| denotes the maximal indegree on the ℓ-th antichain and Cℓ = maxw∈Aℓ E[∥∆w ∥22 ] the maximal expected squared contrast. The following proposition is the main result of this section and it shows how indegree relates to the passage of contrasts through the graph. Proposition 4 (Indegree controls the contrast recursion). Under Assumption 1, for every ℓ ≥ 1, we have with τ 2 the marginal variance and ν is the Schoenberg mixing measure: " # −dk/2 Z ∞ 2s Cℓ ≤ Ψkℓ (Cℓ−1 ), Ψk (u) = 2dτ 2 1 − 1+ u ν(ds) , d 0 Moreover, for every fixed u ≥ 0, the map k 7→ Ψk (u) is nondecreasing. Proof. Fix ℓ ≥ 1 and w ∈ Aℓ . Enumerate the parents of w as Pa(w) = {p1 , . . . , pm }, where m = |Pa(w)| ≤ kℓ . By Lemma 2, E[∥∆w ∥22 ] = 2dτ 2 1 − E[κ(ϱ2w )] . Using the Laplace-transform representation of κ and Tonelli’s theorem, since the integrand is nonnegative, gives Z ∞

E[e−sϱw ] ν(ds) .

E[∥∆w ∥22 ] = 2dτ 2 1 −

(37)

0 2

We therefore seek a lower bound on E[e−sϱw ]. Let Gw denote the sigma-field generated by the latent states at all strict ancestors of the parents of w,   [ Gw = σ Fu(i) : u ∈ Anc(p), i ∈ {a, b} . p∈Pa(w)

Conditional on Gw , the inputs of each parent module fpj are fixed. By the mutual independence of the GP modules, the random vectors ∆p1 , . . . , ∆pm are therefore conditionally independent given Gw . Applying Lemma 2 to each parent pj gives ∆pj | Gw ∼ N(0, ξj Id ), ξj = 2τ 2 1 − κ(ϱ2pj ) . Since ∥∆pj ∥22 is conditionally distributed as a scaled χ2d variable, for every s ≥ 0, h i 2 E e−s∥∆pj ∥2 Gw = (1 + 2sξj )−d/2 . 35

(38)

Combining (38) with conditional independence, m h i Y 2 E e−s∥∆pj ∥2 Gw

E[e−sϱw | Gw ] = =

j=1 m Y

(1 + 2sξj )−d/2 .

(39)

j=1

By the arithmetic–geometric mean inequality, m  m  m m m X X Y 2s 1 (1 + 2sξj ) = 1 + ξj  . (1 + 2sξj ) ≤  m m j=1 j=1 j=1

(40)

Since x 7→ x−d/2 is decreasing on (0, ∞), (39) and (40) give  2

E[e−sϱw | Gw ] ≥ 1 +

m 2s X

m j=1

−dm/2 ξj 

(41)

Define gs,m : R+ → R+ by: gs,m (x) :=

2s 1+ x m

−dm/2 .

Then it’s first two derivatives are: −dm/2−1 2s ′ gs,m (x) = − ds 1 + x m 2 −dm/2−2 dm dm 2s 2s ′′ gs,m (x) = +1 1+ x 2 2 m m

≤ 0,

(42)

≥ 0.

Thus gs,m is nonincreasing and convex. Equation (41) and Jensen’s inequality give      m m h i X X −sϱ2w −sϱ2w E[e ] = E E[e | Gw ] ≥ Egs,m  ξj  ≥ gs,m  E[ξj ] . j=1

j=1

We now bound E[ξj ]. By (36), E[∥∆pj ∥22 | {Fp(a) , Fp(b) }p∈Pa(pj ) ] = dξj . Taking expectations of both sides gives E[ξj ] =

Cℓ−1 1 E[∥∆pj ∥22 ] ≤ , d d

because pj ∈ Aℓ−1 . Thus m X

E[ξj ] ≤

j=1

mCℓ−1 . d

By the monotonicity of gs,m in (42), −sϱ2w

E[e

] ≥ gs,m

mCℓ−1 d

2s 1 + Cℓ−1 d

−dm/2 .

(43)

Substituting (43) into (37) gives " E[∥∆w ∥22 ] ≤ 2dτ 2

Z ∞

1− 0

2s 1 + Cℓ−1 d

−dm/2

# ν(ds) .

Since m ≤ kℓ and m 7→ (1 + 2sCℓ−1 /d)−dm/2 is nonincreasing, we further obtain " # −dkℓ /2 Z ∞ 2s 2 2 E[∥∆w ∥2 ] ≤ 2dτ 1 − 1 + Cℓ−1 ν(ds) = Ψkℓ (Cℓ−1 ). d 0 Taking the maximum over w ∈ Aℓ proves Cℓ ≤ Ψkℓ (Cℓ−1 ). Finally, for every fixed u, s ≥ 0, the map k 7−→

−dk/2 2s 1+ u d

is nonincreasing, so k 7→ Ψk (u) is nondecreasing. The next corollary gives a local contraction criterion from the recursion of Proposition 4, showing that the derivative of Ψk at the origin determines whether small contrasts decay geometrically. Corollary 4 (Local contraction criterion). Suppose Assumption 1 holds, and let k̄ = supℓ≥1 kℓ and µ1 = R∞ s ν(ds) < ∞. If 0 2dτ 2 k̄ µ1 < 1, then there exist ρ ∈ (0, 1) and u⋆ > 0 such that, whenever Cℓ0 ≤ u⋆ for some ℓ0 ≥ 0, Cℓ0 +m ≤ ρm Cℓ0 ,

m ≥ 0.

In particular, for the squared-exponential kernel with ν = δ1/(2λ2 ) , the condition becomes dτ 2 k̄/λ2 < 1. Proof. For every k ≥ 1, Ψk (0) = 0. Since µ1 < ∞, dominated convergence justifies differentiating under the integral sign and gives −dk/2−1 Z ∞ 2s ′ 2 Ψk (u) = 2dτ ks 1 + u ν(ds). d 0 In particular, Ψ′k (0) = 2dτ 2 kµ1 , and by assumption, Ψ′k̄ (0) = 2dτ 2 k̄µ1 < 1. Choose any ρ ∈ (Ψ′k̄ (0), 1). Since Ψk̄ is differentiable at 0 with Ψk̄ (0) = 0, there exists u⋆ > 0 such that Ψk̄ (u) ≤ ρu

for all u ∈ [0, u⋆ ].

Whenever Cℓ−1 ≤ u⋆ , Proposition 4 and the monotonicity of Ψk in k give Cℓ ≤ Ψkℓ (Cℓ−1 ) ≤ Ψk̄ (Cℓ−1 ) ≤ ρCℓ−1 . In particular, Cℓ ≤ ρCℓ−1 ≤ u⋆ , so the bound propagates. Iterating from ℓ0 onward gives Cℓ0 +m ≤ ρm Cℓ0 ,

m ≥ 0.

For the squared-exponential case, µ1 = 1/(2λ2 ), so the condition 2dτ 2 k̄µ1 < 1 reduces to dτ 2 k̄/λ2 < 1. The corollary shows that, once the expected squared contrast enters a sufficiently small neighbourhood of zero, it decays geometrically fast, with the rate governed by Ψ′k̄ (0). The criterion 2dτ 2 k̄µ1 < 1 makes the role of indegree clear: higher indegree k̄ tightens the condition, because the concatenated parent input lives in a space of dimension k̄ · d and the kernel sees a larger total squared distance. For a concrete illustration, consider a layered block of squared-exponential modules with output dimension d = 2, marginal variance τ 2 = 1, and length-scale λ = 2. The local contraction criterion then reads k̄ < λ2 /d = 2, so a chain (k̄ = 1) contracts locally, while a block with maximal indegree k̄ = 2 or higher may fail to satisfy this sufficient contraction condition.

Remark 3 (Layer-dependent dimensions and kernel parameters). The homogeneity assumptions in Proposition 4 are only used to keep the notation compact. Suppose instead that every node in Aℓ−1 has common output dimension dℓ−1 , every node in Aℓ has common output dimension dℓ , and every node in Aℓ uses a common radial kernel Z ∞ Kw (u, u′ ) = τℓ2 κℓ (∥u − u′ ∥22 )Idℓ , κℓ (r) = e−sr νℓ (ds). 0

The proof of Propostion 4, mutatis mutandis, gives " # −dℓ−1 kℓ /2 Z ∞ 2s 2 Cℓ−1 Cℓ ≤ 2dℓ τℓ 1 − 1+ νℓ (ds) . dℓ−1 0 Thus the argument extends verbatim to layer-dependent dimensions and kernel hyperparameters. The contraction result of Proposition 4 is specific to product-type fusion. Under fusion mechanisms that have an additive root-only summand, the separating property is robust to indegree: Proposition 5 (Additive kernel decompositions preserve separating contributions). Let w ∈ U be a non-root P (S) (S) node with kernel Kw = S∈Sw Kw , where Sw ⊆ 2Pa(w) \ {∅} and each Kw is a valid positive-semidefinite kernel on the coordinates indexed by S. Then, for any two distinct cases a ̸= b, X Γ(S) Γw (a, b) = Γ(S) w (a, b), w (a, b) ⪰ 0 for every S ∈ Sw . S∈Sw ⋆ In particular, if some subfamily Sw ⊆ Sw satisfies X Γ(S) w (a, b) ⪰ v⋆ Idw ⋆ S∈Sw

almost surely, then w is v⋆ -separating for (a, b), irrespective of the remaining summands. (a)

(b)

Proof. Write ua = (Fp )p∈Pa(w) and ub = (Fp )p∈Pa(w) for the concatenated parent states, and let ua,S , ub,S denote their restrictions to the coordinates indexed by S. For each S ∈ Sw , define (S) (S) (S) (S) Γ(S) w (a, b) := Kw (ua,S , ua,S ) + Kw (ub,S , ub,S ) − Kw (ua,S , ub,S ) − Kw (ub,S , ua,S ).

Then, by the additive structure of Kw , Γw (a, b) =

Γ(S) w (a, b).

S∈Sw (S)

(S)

For each S ∈ Sw , validity of the kernel Kw implies Γw (a, b) ⪰ 0. Therefore X X Γw (a, b) = Γ(S) Γ(S) w (a, b) ⪰ w (a, b) ⪰ v⋆ Idw , ⋆ S∈Sw

S∈Sw

⋆ where the first inequality drops the positive-semidefinite summands outside Sw . Thus w is v⋆ -separating.

C.3

Outdegree effects through branching

We now isolate a branching mechanism inside the layered block of Assumption 1. Let Bℓ ⊆ Aℓ be non-empty subsets and, for each w ∈ Bℓ , choose a set Ch⋆ (w) ⊆ Bℓ+1 of designated children such that w ∈ Pa(v) for every v ∈ Ch⋆ (w). Assume that |Ch⋆ (w)| = b for all w and that the families {Ch⋆ (w) : w ∈ Bℓ } are pairwise disjoint for each ℓ. The designated edges then form a rooted b-ary branching subgraph. Nodes in this branching subgraph are allowed to have additional parents outside the designated edges. This can only strengthen the lower bound used below. To see this, consider a designated edge w → v. Since w ∈ Pa(v), the term corresponding to w appears in the sum defining ϱ2v , so X ϱ2v = ∥∆p ∥22 ≥ ∥∆w ∥22 . p∈Pa(v)

Therefore, on the event {∥∆w ∥2 ≥ t}, ϱ2v ≥ t2 .

(44)

Because κ is nonincreasing, (44) implies 2τ 2 1 − κ(ϱ2v ) Id ⪰ 2τ 2 1 − κ(t2 ) Id =: σ 2t Id , where σ 2t is the smallest conditional contrast variance parameter compatible with the event that one designated parent already has contrast at least t. For each t > 0, define  P χ2d ≥ t2 /σ 2t , σ 2t > 0, pt := 0, σ 2t = 0. Figure 11 shows a binary instance (b = 2) of this designated branching pattern inside the layered block.

Aℓ

Bℓ

Aℓ+1

Bℓ+1

Aℓ+2

Bℓ+2

Figure 11: Binary designated branching pattern (b = 2) inside the layered block. The horizontal contours indicate the antichains Aℓ , Aℓ+1 , and Aℓ+2 . The dashed inner contours indicate the selected node subsets Bℓ ⊆ Aℓ , Bℓ+1 ⊆ Aℓ+1 , and Bℓ+2 ⊆ Aℓ+2 . Solid arrows denote the designated branching edges between these subsets; dashed arrows indicate omitted incoming or outgoing edges in the surrounding DAG. Proposition 6 (Outdegree can sustain contrasts through branching). Assume Assumption 1. Let (Bℓ )ℓ≥0 and {Ch⋆ (w)} be as above, with common branching factor b ≥ 1, and fix a seed w0 ∈ B0 . If P(∥∆w0 ∥2 ≥ t) > 0 and bpt > 1, then P(Mℓ ≥ t for all ℓ ≥ 0) > 0. Proof. Fix t > 0. By the definition above, the assumption bpt > 1 implies pt > 0, hence σ 2t > 0. Therefore pt = P χ2d ≥ t2 /σ 2t . We prove the stronger statement that the threshold persists already on the designated branching subgraph, namely P max ∥∆w ∥2 ≥ t for all ℓ ≥ 0 w∈Bℓ

> 0.

(45)

Since Bℓ ⊆ Aℓ for every ℓ, (45) immediately implies the proposition. For each ℓ ≥ 0, define the sigma-field Hℓ := σ(fu : u ∈ Anc(Bℓ+1 )) . Because |Ch⋆ (w)| = b ≥ 1 for every w ∈ Bℓ , each w ∈ Bℓ has at least one designated child in Bℓ+1 , and therefore Bℓ ⊆ Anc(Bℓ+1 ). By the same measurability argument used in Lemma 1, every parent state of a node in Bℓ+1 is Hℓ -measurable, and ∆w is Hℓ -measurable for every w ∈ Bℓ . 39

We first derive a one-step lower bound. Fix ℓ ≥ 0, w ∈ Bℓ , and v ∈ Ch⋆ (w). The squared parent contrast at v is X ϱ2v = ∥∆p ∥22 . p∈Pa(v)

The parent states of v are Hℓ -measurable. Therefore, by Lemma 2, under the conditional law given Hℓ , ∆v | Hℓ ∼ N 0, 2τ 2 1 − κ(ϱ2v ) Id . On the event Aw := {∥∆w ∥2 ≥ t}, equation (44) gives ϱ2v ≥ t2 . Therefore, on Aw , 2τ 2 1 − κ(ϱ2v ) ≥ σ 2t . Since

we obtain

∥∆v ∥22 2 2τ (1 − κ(ϱ2v )) P(∥∆v ∥2 ≥ t | Hℓ ) = P χ2d ≥

Hℓ ∼ χ2d ,

t2 2τ 2 (1 − κ(ϱ2v ))

t2 2 ≥ P χd ≥ 2 = pt σt

on Aw .

(46)

S We next record the relevant conditional independence. Let v = ̸ v ′ be distinct nodes in u∈Bℓ Ch⋆ (u) ⊆ Bℓ+1 . Conditional on Hℓ , the contrasts ∆v and ∆v′ are measurable functions of the independent GP modules fv and fv′ , respectively. Hence ∆v and ∆v′ are conditionally independent given Hℓ . In particular, the indicators 1{∥∆v ∥2 ≥ t} are conditionally independent across distinct designated children. We now build an embedded Galton–Watson process by thinning these threshold exceedances. On a product extension of the original probability space, let {Uv }v∈S Bℓ be an i.i.d. family of uniform random variables in ℓ≥1

(0, 1), independent of the DAG-DGP prior. This auxiliary extension does not change the marginal DAG-DGP probabilities. For each ℓ ≥ 0, define ! ℓ [ b Hℓ := σ Hℓ , {Uv : v ∈ Bm } . m=1

bℓ contains the past uniforms up to level ℓ, but not the current uniforms on level ℓ + 1. Thus H We recursively construct random sets Dℓ ⊆ Bℓ , to be interpreted as a distinguished surviving population. Set ( {w0 }, ∥∆w0 ∥2 ≥ t, D0 := ∅, ∥∆w0 ∥2 < t. b0 -measurable and Then D0 is H

|D0 | = 1{∥∆w0 ∥2 ≥ t}.

bℓ -measurable. Let Suppose Dℓ has been constructed and is H [ Iℓ := Ch⋆ (w). w∈Dℓ

Because the designated child families are pairwise disjoint, each v ∈ Iℓ belongs to the designated child set of a unique w ∈ Dℓ . For v ∈ Iℓ , define rv := P(∥∆v ∥2 ≥ t | Hℓ ) . bℓ -measurable. If v ∈ Iℓ , then v ∈ Ch⋆ (w) for some w ∈ Dℓ , and by The variable rv is Hℓ -measurable, hence H construction ∥∆w ∥2 ≥ t. Therefore (46) gives rv ≥ pt > 0

for every v ∈ Iℓ .

Hence pt /rv ∈ [0, 1]. Define

Yv := 1{∥∆v ∥2 ≥ t}1

pt Uv ≤ rv

v ∈ Iℓ ,

and set Dℓ+1 := {v ∈ Iℓ : Yv = 1}. bℓ+1 -measurable and every retained child is above the threshold. By construction, Dℓ+1 is H bℓ , the family {Yv : v ∈ Iℓ } is independent with each Yv ∼ Bernoulli(pt ). We claim that, conditional on H Let v1 , . . . , vm ∈ Iℓ be distinct. For each j, the variable ∆vj is measurable with respect to σ(Hℓ , fvj ), because the parent inputs of vj are Hℓ -measurable. The GP modules fv1 , . . . , fvm are mutually independent and jointly bℓ . The current uniforms Uv , . . . , Uv are i.i.d. and independent independent of the past uniforms entering H 1 m of both the DAG-DGP prior and the past uniforms. Hence the pairs (∆v1 , Uv1 ), . . . , (∆vm , Uvm ) bℓ , and so are the variables Yv , . . . , Yv . are conditionally independent given H 1 m bℓ ), we obtain For the conditional success probability, using that Uv is independent of σ(∆v , H pt b b Hℓ P(Yv = 1 | Hℓ ) = E 1{∥∆v ∥2 ≥ t}1 Uv ≤ rv pt b = E 1{∥∆v ∥2 ≥ t} H ℓ rv pt bℓ . = P ∥∆v ∥2 ≥ t | H rv

(47)

bℓ beyond Hℓ consists only of past uniforms, which are independent of the The additional information in H current GP variables. Therefore bℓ = P(∥∆v ∥2 ≥ t | Hℓ ) = rv . P ∥∆v ∥2 ≥ t | H Substituting into (47) gives bℓ ) = pt . P(Yv = 1 | H For each w ∈ Dℓ , define Nw :=

Yv .

v∈Ch⋆ (w)

bℓ , each Nw has law Because |Ch⋆ (w)| = b and the designated child sets are pairwise disjoint, conditional on H Bin(b, pt ), and the family {Nw : w ∈ Dℓ } is conditionally independent. Define Zℓ := |Dℓ |. Then X Zℓ+1 = Nw . w∈Dℓ

Thus, on the event {Zℓ = n}, the variable Zℓ+1 is the sum of n independent Bin(b, pt ) variables. Its conditional law depends only on n, and not on the earlier history. Hence (Zℓ )ℓ≥0 is a Galton–Watson process with offspring distribution Bin(b, pt ) and random initial state Z0 = 1{∥∆w0 ∥2 ≥ t}. Its generating function is g(s) = (1 − pt + pt s)b ,

s ∈ [0, 1],

and its mean offspring number is g ′ (1) = bpt > 1. Hence the process is supercritical. By the classical extinction-probability theorem for Galton–Watson processes, the extinction probability ξt is the smallest nonnegative solution of g(s) = s. Since g ′ (1) = bpt > 1, one has ξt < 1, and the survival probability πt := 1−ξt is strictly positive (Athreya and Ney, 1972, Ch. I, Sec. 5, Thm. 1). Conditional on Z0 = 1, P(Zℓ > 0 for all ℓ ≥ 0|Z0 = 1) = πt > 0.

Since Z0 = 1{∥∆w0 ∥2 ≥ t}, P(Zℓ > 0 for all ℓ ≥ 0) = P(Zℓ > 0 for all ℓ ≥ 0|Z0 = 1) P(Z0 = 1) = πt P(∥∆w0 ∥2 ≥ t) > 0. On the event {Zℓ > 0 for all ℓ ≥ 0}, every set Dℓ is non-empty, and every w ∈ Dℓ satisfies ∥∆w ∥2 ≥ t. Since Dℓ ⊆ Bℓ , max ∥∆w ∥2 ≥ t for all ℓ ≥ 0. w∈Bℓ

This proves (45). Since Bℓ ⊆ Aℓ , on the same event Mℓ ≥ max ∥∆w ∥2 ≥ t w∈Bℓ

for all ℓ ≥ 0.

Therefore P(Mℓ ≥ t for all ℓ ≥ 0) > 0, which proves the claim. Corollary 5 (Scalar threshold probability). Under the additional assumption d = 1, the threshold probability pt in Proposition 6 takes the form  2 Φ − t , σ 2 > 0, N t σt pt :=  0, σ 2t = 0. Proof. If d = 1 and σ 2t > 0, then t2 t 2 pt = P χ1 ≥ 2 = P |Z| ≥ , σt σt Hence

t pt = 2ΦN − . σt

If σ 2t = 0, the conclusion is immediate from the definition of pt .

Z ∼ N(0, 1).

Intermediate observations as stochastic skip connections

This appendix proves Theorem 2 and the following corollary thereto, which characterises a simple setting in a readily interpretable way. Corollary 6 (Gaussian refresh at an observed source). Let u ∈ Aℓ0 be scalar, assume Ou = {a, b} × {1}, and suppose that Ku (x, x) = τu2 for all parent inputs x. If Yu(i) = Fu(i) + ξu(i) ,

i.i.d.

ξu(i) ∼ N (0, σu2 ), (a)

i ∈ {a, b},

(b)

then, conditional on the parent states and on (Yu , Yu ), on the event Γu (a, b) > 0, 2σu2 Γu (a, b) Γu (a, b) (a) (b) . Fu(a) − Fu(b) ∼ N Y − Y , u Γu (a, b) + 2σu2 u Γu (a, b) + 2σu2 Consequently, if mu and vu denote the mean and variance in the display above, then the local filtered law √ (a) (b) assigns probability at least 1 − ΦN (ε − |mu |)/ vu to the event |Fu − Fu | > ε. If Γu (a, b) = 0, the conditional contrast is degenerate and this source strength is zero. (a)

Thus the Gaussian source strength is governed by the observed discrepancy Yu

(b)

− Yu , the noise variance

σu2 , and the local two-point prior variance Γu (a, b). The argument is organised so that the main theorem follows from three ingredients. First, filtering on an observed antichain yields a conditional product structure across the nodes of that antichain. Second, a refreshed contrast at one such node can be transported downstream by additive kernel components. Third, disjoint retaining routes remain conditionally independent, so their effects combine multiplicatively. Throughout this section the dataset is fixed and only the latent states remain random. We work under the standing well-posedness condition that every posterior or local conditional distribution displayed below is well defined, i.e. that the corresponding normalising constant is finite and strictly positive. As in the previous sections, root states are deterministic and are omitted from latent sigma-fields.

D.1

Filtering on an observed antichain

Consider a progressive antichain sequence (Aℓ )ℓ≥0 of non-root nodes. For each ℓ ≥ 0, let U≤ℓ := Anc(Aℓ ) ∪ Aℓ . The corresponding filtering posterior is Y Πℓ d{Fw }w∈U≤ℓ ∝ p0 dFw | {Fp }p∈Pa(w) pw (Yw | Fw , Ow ) , w∈U≤ℓ

where pw (Yw | Fw , Ow ) ≡ 1 when Ow = ∅. We also write Hℓ := σ(Fv : v ∈ Anc(Aℓ )) ,

Fℓ := σ(Fv : v ∈ U≤ℓ ) .

Thus Hℓ records the strict latent ancestors of the current antichain, whereas Fℓ also contains the current antichain variables. If w ∈ Aℓ , every latent parent of w is Hℓ -measurable. We define the local prior kernel Pwℓ (dFw | Hℓ ) := p0 dFw | {Fp }p∈Pa(w) . The one-node filtered kernel is the probability kernel Qℓw (dFw | Hℓ ) ∝ pw (Yw | Fw , Ow )Pwℓ (dFw | Hℓ ), For ℓ1 ≥ ℓ0 , we write Πℓ0 →ℓ1 for the predictive law obtained by first drawing from Πℓ0 and then propagating the DAG-DGP prior forward from Aℓ0 to Aℓ1 , without assimilating observations beyond level ℓ0 . 43

Lemma 3 (Conditional factorisation of the filtered antichain). For every ℓ ≥ 0, the conditional law of the current antichain under Πℓ factorises as O Πℓ (d{Fw }w∈Aℓ |Hℓ ) = Qℓw (dFw | Hℓ ). w∈Aℓ

Proof. The proof makes a repetitive use of important conditional independent properties. Along the proof, we mainly refer to Section 3 and 4 of Dawid (1979). Fix ℓ ≥ 0. By definition, Hℓ = σ(Fv : v ∈ Anc(Aℓ )) , so each state matrix Fv , v ∈ Anc(Aℓ ), is Hℓ -measurable. Equivalently, under any regular conditional law given Hℓ , these strict-ancestor states are degenerate at their realised values. Since Aℓ is an antichain, no node in Aℓ is a strict ancestor of any other node in Aℓ . Therefore, for each w ∈ Aℓ , the parent states of w are either deterministic roots or are Hℓ -measurable. We first justify the conditional-independence step. Temporarily regard the local observations as random e w , w ∈ Aℓ , generated from the nodewise conditional distributions pw (· | Fw , Ow ), and define variables Y e w ), Bw := (Fw , Y

w ∈ Aℓ .

Given Hℓ , the blocks {Bw : w ∈ Aℓ } are jointly conditionally independent, as the latent states are generated by distinct independent GP modules evaluated at Hℓ -measurable inputs, and the observation variables are then generated nodewise from their corresponding latent states. Consider first two disjoint subcollections I, J ⊂ Aℓ , and write BI = (Bw )w∈I , BJ = (Bw )w∈J , with e I and Y e J . We have analogous notation for Y BI ⊥⊥ BJ | Hℓ . e I is a measurable function of BI , the conditioning-stability property in (Dawid, 1979, Sec. 4, Since Y Lemma 4.2(ii)) gives e I. BI ⊥⊥ BJ | Hℓ , Y Using the symmetry of conditional independence (Dawid, 1979, Sec. 3.1, Theorem 3.1), applying again Dawid e J of BJ , and then using symmetry once more, gives (1979, Sec. 4, Lemma 4.2(ii)) to the measurable function Y e I, Y e J. BI ⊥⊥ BJ | Hℓ , Y By the closure of conditional independence under measurable transformations (Dawid, 1979, Sec. 4, Lemma 4.2(i)), applied to the maps BI 7→ FI and BJ 7→ FJ , this implies e I, Y e J. FI ⊥⊥ FJ | Hℓ , Y The same argument applied inductively, using the joint-independence convention following (Dawid, 1979, Sec. 4, Lemma 4.3), gives joint conditional independence of {Fw : w ∈ Aℓ } after conditioning on Hℓ and on the realised local observations at the current antichain. The preceding conditional-independence argument shows that the conditional law of {Fw : w ∈ Aℓ }, given Hℓ and the realised current-antichain observations, factorises over nodes. It remains to identify the corresponding nodewise conditional factor, and to verify that it is exactly Qℓw (· | Hℓ ). The conditional prior law of {Fw : w ∈ Aℓ } given Hℓ is O Pwℓ (dFw | Hℓ ), (48) w∈Aℓ

because distinct current-antichain states are obtained by evaluating distinct independent GP modules at Hℓ -measurable inputs. The conditional distribution factors for the current-antichain observations also factorise nodewise: Y pw (Yw | Fw , Ow ), (49) w∈Aℓ

with factors equal to one at unobserved nodes. All conditional distribution factors and prior factors associated with strict ancestors are Hℓ -measurable and hence do not affect the conditional distribution of {Fw : w ∈ Aℓ } beyond the realised value of Hℓ . ⋆ Let Bw be a measurable set in the state space of Fw , one for each w ∈ Aℓ . By Bayes’ rule, (48), and (49), ⋆ Πℓ (Fw ∈ Bw for all w ∈ Aℓ |Hℓ ) RQ Q pw (Yw | Fw , Ow )P ℓ (dFw | Hℓ ) ⋆ w∈Aℓ

w∈Aℓ

ℓ w∈Aℓ pw (Yw | Fw , Ow )Pw (dFw | Hℓ )

(50)

Since the integrands are non-negative products of nodewise terms, Tonelli’s theorem gives Z Y pw (Yw | Fw , Ow )Pwℓ (dFw | Hℓ ) Q w∈Aℓ

⋆ Bw w∈A

ℓ

Y Z

w∈Aℓ

⋆ Bw

pw (Yw | Fw , Ow )Pwℓ (dFw | Hℓ ),

(51)

and Z

pw (Yw | Fw , Ow )Pwℓ (dFw | Hℓ )

w∈Aℓ

Y Z

pw (Yw | Fw , Ow )Pwℓ (dFw | Hℓ ).

(52)

(53)

w∈Aℓ

Substituting (51) and (52) into (50) yields

w∈Aℓ

Since the rectangles

D.2

⋆ w∈Aℓ Bw generate the product sigma-field, (53) proves the stated factorisation.

Retaining routes below the observed antichain

We recall that disjointness conventions are given in Appendix A.4 and now formalise a downstream transport mechanism. Definition 5 (ε-retaining route). An admissible route γ = (v0 , v1 , . . . , vL ) below Aℓ0 is ε-retaining with factor ρ ∈ [0, 1] if, under the predictive law propagated from Πℓ0 , n o (b) (b) (a) − F ∥ > ε F ≥ ρ 1 ∥F − F ∥ > ε Πℓ0 -a.s. P ∥Fv(a) ℓ 2 2 v v v 0 0 0 L L Definition 6 (Additive route-retention coefficients). Let γ = (v0 , . . . , vL ) be an admissible route below Aℓ0 , and fix ε > 0. We say that γ satisfies the additive ε-retention condition if, for each r = 1, . . . , L, the kernel at vr , as a function of the distinguished parent coordinate vr−1 , admits the decomposition Kvr (u, z), (u′ , z ′ ) = Kv(rest) (u, z), (u′ , z ′ ) + ςr2 k̄r (u, u′ )Idvr , r (rest)

where Kvr

is positive semidefinite, ςr > 0, and k̄r is a scalar correlation kernel satisfying sup

k̄r (u, u′ ) ≤ r̄r (ε) < 1.

∥u−u′ ∥2 >ε

For such a route, define βr (ε) := 2 ΦN − p

ε 2ςr2 (1 − r̄r (ε))

and ργ (ε) :=

L Y

r = 1, . . . , L,

βr (ε),

r=1

with the empty product interpreted as one. Proposition 7 (Additive routes are retaining). If an admissible route γ below Aℓ0 satisfies the additive ε-retention condition, then γ is ε-retaining with factor ργ (ε). Proof. Start considering γ = (v0 , . . . , vL ). If L = 0, then vL = v0 , so n o P ∥Fv(a) − Fv(b) ∥2 > ε Fℓ0 = 1 ∥Fv(a) − Fv(b) ∥2 > ε , 0 0 L L and the claim holds with the empty product equal to one. Assume now that L ≥ 1. For r = 0, . . . , L, define Gr := Fℓ0 ∨ σ(Fv1 , . . . , Fvr ),

Er := {∥Fv(a) − Fv(b) ∥2 > ε}, r r where G0 = Fℓ0 . We first prove the one-step inequality

P(Er | Gr−1 ) ≥ βr (ε)1Er−1 ,

r = 1, . . . , L.

(54)

Consider r ∈ {1, . . . , L}. Write ua := Fv(a) , r−1

ub := Fv(b) r−1

for the two values of the distinguished parent coordinate, and write za , zb for the remaining parent coordinates of vr at cases a, b. By admissibility, every parent of vr other than vr−1 is either a root node or belongs to U≤ℓ0 . Hence za , zb , ua , ub are Gr−1 -measurable. The regular conditional law of ∆vr := Fv(a) − Fv(b) r r given Gr−1 is centred Gaussian with covariance Γvr (a, b) = Kvr ((ua , za ), (ua , za )) + Kvr ((ub , zb ), (ub , zb )) − Kvr ((ua , za ), (ub , zb )) − Kvr ((ub , zb ), (ua , za )). The additive decomposition of the kernel gives Γvr (a, b) = Γ(rest) (a, b) + 2ςr2 1 − k̄r (ua , ub ) Idvr , vr (rest)

where Γvr

(a, b) ⪰ 0. On the event Er−1 , one has ∥ua − ub ∥2 > ε, and therefore k̄r (ua , ub ) ≤ r̄r (ε).

Consequently, on Er−1 , Γvr (a, b) ⪰ 2ςr2 (1 − r̄r (ε))Idvr . Let h ∈ Rdvr be any deterministic unit vector. On Er−1 , the scalar conditional distribution of h⊤ ∆vr given Gr−1 is Gaussian with mean zero and variance at least 2ςr2 (1 − r̄r (ε)). Since Er = {∥∆vr ∥2 > ε} ⊇ {|h⊤ ∆vr | > ε}, we obtain (54).

(55)

P(EL−1 | GL−2 ) ≥ βL−1 (ε)1EL−2 .

Combining (55) and (56),

(56)

P(EL | GL−2 ) ≥ βL (ε)βL−1 (ε)1EL−2 .

Repeating this backward induction along the route yields L Y

P(EL | Fℓ0 ) = P(EL | G0 ) ≥

! βr (ε)

1 E0 ,

r=1

which is the retaining-route condition with factor ργ (ε).

D.3

Proof of Theorem 2

Proof of Theorem 2. For each target node vj , define Gj := {∥Fv(a) − Fv(b) ∥2 > ε}. j j Since Mℓ1 is the maximum contrast over the whole antichain Aℓ1 , s [

{Mℓ1 > ε} ⊇

Gj .

(57)

j=1

Ss It is therefore enough to lower-bound the predictive probability of j=1 Gj . We first work conditionally on Fℓ0 . By the ε-retaining property of γj , P(Gj | Fℓ0 ) ≥ ρj 1Ej ,

j = 1, . . . , s.

(58)

We next verify conditional independence of the target events under the same conditioning. For a route γj , all non-route parents of interior nodes are either deterministic roots or belong to U≤ℓ0 , hence are Fℓ0 measurable. Thus, under the predictive law, the variables generated along the interior of γj are measurable with respect to Fℓ0 together with the GP modules attached to int(γj ). The route interiors are pairwise disjoint, so these collections of GP modules are disjoint. Since distinct GP modules are mutually independent under the DAG-DGP prior, the route-generated random elements are conditionally independent given Fℓ0 . The events G1 , . . . , Gs , being measurable functions of these route-generated random elements and of Fℓ0 , are therefore conditionally independent given Fℓ0 . Hence   s s \ Y P Gcj Fℓ0  = 1 − P(Gj | Fℓ0 ) j=1

j=1

≤

s Y

1 − ρ j 1 Ej ,

(59)

j=1

where the inequality uses (58) and ρj 1Ej ∈ [0, 1]. We now pass from conditioning on Fℓ0 to conditioning on Hℓ0 . By the tower property,       s s \ \ P Gcj Hℓ0  = EP Gcj Fℓ0  Hℓ0  j=1

j=1

 ≤ EΠℓ0 

s Y

 (1 − ρj 1Ej ) Hℓ0  .

j=1

(60)

By Lemma 3, the source states Fu1 , . . . , Fus are conditionally independent given Hℓ0 under Πℓ0 . Since each Ej depends only on Fuj , the source events E1 , . . . , Es are conditionally independent given Hℓ0 . Thus   s s Y Y EΠℓ0  (1 − ρj 1Ej ) Hℓ0  = EΠℓ0 1 − ρj 1Ej Hℓ0 j=1

j=1

s Y

(1 − ρj Πℓ0 (Ej | Hℓ0 ))

j=1 s Y

(1 − ρj qj ),

(61)

j=1

because the conditional law of Fuj given Hℓ0 is Qℓu0j (· | Hℓ0 ), and hence Πℓ0 (Ej | Hℓ0 ) = qj . Combining (60) and (61) gives   s s \ Y c   P Gj Hℓ0 ≤ (1 − ρj qj ). j=1

j=1

Taking complements gives  P

s [

 Gj Hℓ0  ≥ 1 −

j=1

s Y

(1 − ρj qj ).

j=1

Finally, taking expectation under Πℓ0 and using (57) yields  Πℓ0 →ℓ1 (Mℓ1 > ε) ≥ EΠℓ0 1 −

s Y

 (1 − ρj qj ) .

j=1

D.4

Source strengths for common conditional distributions

The theorem is agnostic about how the source probabilities qj are obtained. We now give source-strength calculations for common nodewise conditional distributions. The main result is a bounded-curvature source bound for one-parameter exponential-family conditional distributions in canonical form (e.g., (Wainwright and Jordan, 2008, Sec. 3.2)). For a scalar source node u ∈ Aℓ0 and two cases a ̸= b, we use the local notation for simplicity Fa := Fu(a) ,

Fb := Fu(b) ,

Ya := Yu(a) ,

Yb := Yu(b) .

When the parent inputs of u at cases a and b are denoted by xa and xb , and Ku (x, x) = τu2 for all parent inputs x, we write Γu (a, b) := 2 τu2 − Ku (xa , xb ) , Ωu (a, b) := 2 τu2 + Ku (xa , xb ) . For a one-parameter exponential-family conditional distribution in canonical form, p(y | θ) = h(y) exp{T (y)θ − A(θ)}, define dy := T (Ya ) − T (Yb ). For Γ > 0 and curvature constants 0 ≤ mA ≤ LA < ∞, set αL (Γ) := Γ−1 + LA /2,

αm (Γ) := Γ−1 + mA /2.

Proposition 8 (Bounded-curvature source bound). Consider a scalar observed node u ∈ Aℓ0 with Ou = {a, b} × {1}. Assume the scalar-source notation above, and suppose that Γu (a, b) > 0 and Ωu (a, b) > 0. If the nodewise conditional distribution belongs to a one-parameter exponential family in canonical form and satisfies 0 ≤ mA ≤ A′′ (θ) ≤ LA < ∞, 48

θ ∈ R,

then, for every ε > 0, Qℓu0 |Fu(a) − Fu(b) | > ε Hℓ0 ≥ qEF (ε, dy , Γu (a, b), mA , LA ), where, for any Γ > 0, s

( ) d2y d2y αm (Γ) qEF (ε, dy , Γ, mA , LA ) := exp − αL (Γ) 8αL (Γ) 8αm (Γ) " p |dy | × 1 − ΦN αL (Γ) ε − 2αL (Γ) # p |dy | + ΦN − αL (Γ) ε + . 2αL (Γ)

(62)

Proof. Set ∆ := Fa − Fb ,

S := Fa + Fb ,

ry := T (Ya ) + T (Yb ).

For readability, also write Γ := Γu (a, b),

Ω := Ωu (a, b),

ᾱL := αL (Γ),

ᾱm := αm (Γ).

Under the two-point GP law conditional on Hℓ0 , ∆ and S are centred Gaussian random variables with variances Γ and Ω, respectively. The constant-diagonal assumption gives Cov(∆, S | Hℓ0 ) = Var(Fa | Hℓ0 ) − Var(Fb | Hℓ0 ) = 0. Since the pair is jointly Gaussian, ∆ and S are conditionally independent given Hℓ0 . The change of variables S+∆ S−∆ Fa = , Fb = 2 2 has constant Jacobian. Ignoring the factor h(Ya )h(Yb ), which is constant in (S, ∆), the joint posterior density of (S, ∆) is proportional to 2 δ s2 dy ry s+δ s−δ exp − − + δ+ s−A −A . 2Γ 2Ω 2 2 2 2 Thus, after integrating out s, the marginal posterior density of ∆ is proportional to 2 dy δ exp − + δ Ry (δ), 2Γ 2 where

s2 ry s+δ s−δ exp − + s−A −A ds. 2Ω 2 2 2 R

Z Ry (δ) :=

The bounded-curvature assumption controls the symmetric second difference of A. By Taylor’s theorem with Lagrange remainder (Rudin, 1976, Theorem 5.15), for all x, c ∈ R, mA c2 ≤ A(x + c) + A(x − c) − 2A(x) ≤ LA c2 . Taking x = s/2 and c = δ/2 gives LA δ 2 mA δ 2 exp − Ry (0) ≤ Ry (δ) ≤ exp − Ry (0). 4 4 Let Bε := {|∆| > ε}. The lower bound in (63) gives the lower bound Z ᾱL dy Ry (0) exp − δ 2 + δ dδ 2 2 |δ|>ε 49

(63)

(64)

for the unnormalised posterior mass of Bε . The upper bound in (63) gives the upper bound Z ᾱm 2 dy Ry (0) exp − δ + δ dδ 2 2 R for the full normalising constant. Taking the ratio of (64) and (65), and cancelling Ry (0), yields o n R d exp − ᾱ2L δ 2 + 2y δ dδ |δ|>ε n o Qℓu0 (Bε | Hℓ0 ) ≥ R . dy ᾱm 2 exp − δ + δ dδ 2 2 R

(65)

(66)

Completing the square, α dy α − δ2 + δ = − 2 2 2

dy δ− 2α

2 +

d2y . 8α

(67)

Using (67) in the denominator of (66) gives ( )r d2y ᾱm 2 dy 2π δ + δ dδ = exp . exp − 2 2 8ᾱm ᾱm R

(68)

Using (67) in the numerator gives ᾱL 2 dy exp − δ + δ dδ 2 2 |δ|>ε ( )r d2y 2π = exp P(|ZEF | > ε), 8ᾱL ᾱL

where

ZEF ∼ N

dy 1 , 2ᾱL ᾱL

(69)

Substituting (68) and (69) into (66) gives Qℓu0 (Bε | Hℓ0 ) ≥

( ) d2y d2y ᾱm exp − P(|ZEF | > ε). ᾱL 8ᾱL 8ᾱm

(70)

Writing the two-sided Gaussian tail in (70) explicitly, and using |dy | because the event is symmetric, yields (62). Proof of Corollary 6. For Gaussian observations, ξ ∼ N (0, σu2 ),

Y = F + ξ, the conditional density has canonical form with T (y) = y/σu2 ,

A(θ) = θ2 /(2σu2 ).

Hence mA = LA = σu−2 , so the upper and lower bounded-curvature inequalities in Proposition 8 are equalities. Equivalently, and more directly, the prior contrast ∆u := Fu(a) − Fu(b) has conditional distribution ∆u | Hℓ0 ∼ N (0, Γu (a, b)). The observation difference satisfies Yu(a) − Yu(b) = ∆u + ηu , 50

ηu ∼ N (0, 2σu2 ),

with ηu independent of ∆u . On {Γu (a, b) > 0}, the one-dimensional Gaussian conditioning formula gives ∆u | Yu(a) , Yu(b) , Hℓ0 ∼ N (mu , vu ), where mu =

Γu (a, b) Y (a) − Yu(b) , Γu (a, b) + 2σu2 u

vu =

2σu2 Γu (a, b) . Γu (a, b) + 2σu2

If Γu (a, b) = 0, the conditional prior of the contrast is degenerate at zero, and so is the filtered contrast. Finally, if Z ∼ N (mu , vu ), then {|Z| > ε} ⊇ {sign(mu )Z > ε}, with either sign used when mu = 0. This gives P(|Z| > ε) ≥ 1 − ΦN

ε − |mu | √ vu

which is the claimed lower bound. Binomial refresh. We say that a scalar observed source node has a binomial conditional distribution with N trials and canonical parameter θ if N p(y | θ) = exp{yθ − N log(1 + eθ )}, y ∈ {0, . . . , N }. y The Bernoulli conditional distribution corresponds to N = 1. This is the canonical form of the binomial one-parameter exponential family; the Bernoulli case is one of the standard one-parameter exponential-family examples( e.g., Wainwright and Jordan (2008, Table 3.1)). Corollary 7 (Bernoulli and binomial refresh). In the setting of Proposition 8, suppose that the scalar observed source node has a binomial conditional distribution with N trials and canonical parameter θ. Then N . Qℓu0 |Fu(a) − Fu(b) | > ε Hℓ0 ≥ qEF ε, Ya − Yb , Γu (a, b), 0, 4 Proof. Here T (y) = y and A(θ) = N log(1 + eθ ). Therefore A′′ (θ) = N

eθ . (1 + eθ )2

The function eθ /(1 + eθ )2 is nonnegative and bounded above by 1/4, with maximum at θ = 0. Hence 0 ≤ A′′ (θ) ≤

N , 4

The result follows from Proposition 8.

θ ∈ R.

Variational inference

E.1

ELBO derivation

The ELBO derivation follows the Doubly Stochastic VI line of research (Salimbeni and Deisenroth, 2017; Ustyuzhaninov et al., 2020; Lindinger et al., 2020). Here we recall that, since the GP modules are independent a priori, the augmented joint distribution factorises as Y p(D, F, U) = p0 (Uw ) p0 Fw | {Fp }p∈Pa(w) , Uw pw (Yw | Fw , Ow ), (71) w∈U

where, as in the main text, pw (Yw | Fw , Ow ) = 1 whenever Ow = ∅. The variational family is Y q(F, U) = q(U) p0 Fw | {Fp }p∈Pa(w) , Uw ,

(72)

w∈U

where we have not yet imposed any resitrctions upon q(U). Starting from the marginal likelihood, Jensen’s inequality gives Z log p(D) = log p(D, F, U) dF dU Z p(D, F, U) dF dU = log q(F, U) q(F, U) p(D, F, U) ≥ Eq(F,U) log =: L(q). (73) q(F, U) Substituting Eqs. (71) and (72) into Eq. (73), we obtain " # Q w∈U p0 (Uw ) p0 Fw | {Fp }p∈Pa(w) , Uw pw (Yw | Fw , Ow ) Q L(q) = Eq(F,U) log q(U) w∈U p0 Fw | {Fp }p∈Pa(w) , Uw " # X X = Eq(F,U) log pw (Yw | Fw , Ow ) + log p0 (Uw ) − log q(U) w∈U

w∈U

! =

Eq(Fw ) [log pw (Yw | Fw , Ow )] − KL q(U)

p0 (Uw )

w∈U

! =

Eq(Fw ) [log pw (Yw | Fw , Ow )] − KL q(U)

p0 (Uw ) .

(74)

w∈U

w: Ow ̸=∅

The last equality removes the nodes without observations, since their observational distribution factors are identically one. Recall that for DAG-SVI we can write our distribution in Eq. (6) equivalently as: q(U) = qH (U) = N (m, Λ−1 ),

Λvw = 0

{v, w} ∈ / H.

(75)

We denote by qH (Fw ) the marginal distribution induced by propagating qH (U) through the GP conditionals in the DAG: Z Y qH (F) = qH (U) p0 Fw | {Fp }p∈Pa(w) , Uw dU. (76) w∈U

This induced distribution is not available in closed form, so the expectation terms in Eq. (74) are estimated by Monte Carlo, drawing samples of Fw by following the topological order of the DAG: starting from the roots and propagating samples through each non-root node’s GP conditional given the (already sampled) parent values. The mean-field DAG-VI objective is obtained by restricting the inducing posterior to factorise across nodes, Y q(U) = qw (Uw ), (77) w∈U

which yields LVI =

EqVI (Fw ) [log pw (Yw | Fw , Ow )] −

KL (qw (Uw ) ∥ p0 (Uw )) ,

(78)

w∈U

w: Ow ̸=∅

Q where qVI (Fw ) is induced by v∈U qv (Uv ) through the same GP conditionals. Finally, since each observational distribution pw (Yw | Fw , Ow ) factorises over cases, the term in Eq. (74) decomposes over observations for any choice of q(U) (whether the structured qH , the mean-field qVI , or any other family of the form Eq. (72)) as h i X X (i) (i) (i) Eq(F (i) ) log pw Yw | Fw(i) , Ow , Ow := {j : (i, j) ∈ Ow }. (79) w

w: Ow ̸=∅ i: O (i) ̸=∅ w

(i)

Two practical consequences follow. First, evaluating the ELBO requires only the per-case marginals q(Fw ); correlations between latent values at different cases i = ̸ i′ never enter, so they need not be tracked during ancestral sampling. Second, the outer sum over i admits unbiased mini-batch estimation, recovering the standard doubly stochastic scheme of Salimbeni and Deisenroth (2017) at the case level while the structured q(U) preserves cross-node coupling at the inducing level.

E.2

Marginal ancestral sampling

The ELBO in Eq. (74) requires expectations under the latent law induced by the variational family, Z Y q(F) = q(U) p0 Fw | {Fp }p∈Pa(w) , Uw dU. w∈U

For chain DGPs, structured Gaussian posteriors over inducing outputs can be marginalised recursively while retaining dependencies between latent processes (Lindinger et al., 2020). We use the same Gaussianconditioning principle along a topological ordering of the DAG. Fix a topological ordering of the non-root nodes U. We write v < w whenever v appears before w, and define F<w := {Fv : v < w}, U<w := {Uv : v < w}. When these collections appear in matrix expressions, they are understood as the corresponding vectorised concatenations in the chosen topological order. We consider the global Gaussian inducing posterior q(U) = N (U; m, Σ), and for any subset S ⊆ U we write mS and ΣS,S for the sub-vector and sub-block of m and Σ indexed by the inducing entries in {Uv : v ∈ S}; cross-blocks ΣS,T are defined analogously. In particular, mw is the marginal mean of Uw , Σww its marginal covariance, and Σw,<w the cross-covariance between Uw and U<w . b <w ) denotes a realised value of the corresponding random variable. We denote Throughout, a hat (e.g. F by Av and Rv the finite-dimensional GP conditional mean map and residual covariance at node v, evaluated at these realised parent-state inputs, so that b p }p∈Pa(v) , Uv = N (Fv ; Av Uv , Rv ). p0 Fv | {F Equivalently, in unwhitened inducing coordinates, −1 Av = Kv,F Z Kv,ZZ ,

−1 Rv = Kv,F F − Kv,F Z Kv,ZZ Kv,ZF ,

where the kernel blocks are computed from the nodewise kernel Kv , the current parent-state inputs, and the inducing locations Zv . For a fixed node w, define the stacked prefix matrices A<w := blockdiag(Av : v < w),

R<w := blockdiag(Rv : v < w).

b <w . All Av and Rv in these blocks are evaluated along the realised F 53

Proposition 9 (Marginal ancestral factorisation). Consider the DAG-DGP variational family in Eq. (72) with global Gaussian inducing posterior q(U) = N (U; m, Σ). Fix a topological ordering of U, and define C<w := A<w Σ<w,<w A⊤ <w + R<w . Then the induced latent law factorises along the topological order as Y q(F) = q(Fw | F<w ), w∈U

Eq. (9) immediately yields an exact sampler from q(F): traverse U in topological order and, at each node w, b <w ) using the realised F b <w of previously sampled latents. draw Fw from q(Fw | F Proof. Start by considering a fixed node w. By the variational family, Y q(F, U) = q(U) p0 Fv | {Fp }p∈Pa(v) , Uv . v∈U

If we integrate out all future nodes v > w in reverse topological order, their conditional densities integrate to one. Hence Z b b p }p∈Pa(w) , Uw q(U | F b <w ) dU. q(Fw | F<w ) = p0 Fw | {F b <w ). It remains to characterise the marginal conditional law of Uw under q(U | F Again by the variational family, Y b <w ) ∝ q(U) b v | {F b p }p∈Pa(v) , Uv . q(U | F p0 F v<w

Using the definition of Av and Rv , the product of the already-visited local conditionals can be written as Y b v | {F b p }p∈Pa(v) , Uv = N F b <w ; A<w U<w , R<w . p0 F v<w

Thus, conditional on the realised prefix, the previously sampled states act as a linear-Gaussian observation of U<w . b <w ) induced by this linear-Gaussian observation is jointly Under q(U) = N (U; m, Σ), the pair (Uw , F Gaussian with moments b <w ] = A<w m<w , Eq [Uw ] = mw , Eq [F b <w ) = Σw,<w A⊤ , Covq (Uw , F <w and b <w ) = A<w Σ<w,<w A⊤ + R<w = C<w . Covq (F <w Gaussian conditioning (Bishop, 2006, Sec. 2.3, Eqs. (2.81)–(2.82)) therefore gives b <w ] = mw + Σw,<w A⊤ C −1 F b <w − A<w m<w , Eq [Uw | F <w <w and b <w ) = Σww − Σw,<w A⊤ C −1 A<w Σ<w,w . Covq (Uw | F <w <w 54

Finally, at node w, the local GP conditional is b p }p∈Pa(w) , Uw = N (Fw ; Aw Uw , Rw ). p0 Fw | {F The remaining integration is the linear-Gaussian marginalisation (Bishop, 2006, Sec. 2.3.3, Eqs. (2.113)–(2.115)) Z N (Fw ; Aw Uw , Rw )N (Uw ; mw|<w , Σww|<w ) dUw = N Fw ; Aw mw|<w , Rw + Aw Σww|<w A⊤ w . Therefore, b <w ) = N Fw ; Aw mw|<w , Rw + Aw Σww|<w A⊤ . q(Fw | F w Applying this identity at every node in the chosen topological order gives the chain-rule factorisation Y q(F) = q(Fw | F<w ). w∈U

Therefore, sampling each node from the displayed conditional distribution in topological order yields a sample from the induced marginal law q(F).

E.3

Practical ELBO evaluation

The marginal ancestral factorisation in Prop. 9 gives the distribution that must be sampled in order to estimate the part of the ELBO that involves the observational distributions in logarithmic form. In practice, this observational distribution factorises over data indexes, so the estimator only requires the marginal law of the latent variables across DAG nodes for each individual data index and Monte Carlo sample. Cross-index latent correlations do not enter the likelihood estimator and are therefore not materialised. We describe the scalar-output case. Vector-valued nodes are obtained by replacing the interpolation row vectors below by block interpolation matrices and the residual variances by residual covariance blocks. We use lowercase letters for pointwise quantities: aw,b is a single interpolation vector, rw,b a single residual variance, 0 µw,b a scalar mean, and vuv,b a scalar base covariance. Their uppercase counterparts in Prop. 9 denote stacked or matrix quantities. Whitened pointwise GP conditionals. As often in sparse DGP inference, we work in whitened inducing coordinates. We keep the notation Uw for the whitened inducing vector at node w, so that the prior is Uw ∼ N (0, I). For a Monte Carlo sample s ∈ {1, . . . , S} and a minibatch data index i ∈ B, write b = (s, i). Given the already-sampled parent values for the same sample-index pair b, the sparse GP conditional at node w has the pointwise linear-Gaussian form p0 Fw,b | {Fp,b }p∈Pa(w) , Uw = N Fw,b ; a⊤ (80) w,b Uw , rw,b , where aw,b ∈ RMw is the whitened interpolation vector and rw,b > 0 is the corresponding diagonal residual variance. Both aw,b and rw,b depend on the current parent-state input to node w, and hence on the ancestral samples already drawn for the same sample-index pair. E.3.1

Dense implementation

The dense implementation first materialises the full covariance Σ = Λ−1 of the structured inducing posterior qH (U) = N (m, Λ−1 ). During the ancestral pass, for each sample-index pair b, it maintains the already-sampled latent values Fb<w,b , their base means, and their base covariance matrix under the Gaussian model induced by Σ. For any two nodes u and v whose pointwise GP conditionals have already been constructed for sample-index pair b, define the base pointwise covariance 0 vuv,b = a⊤ u,b Σuv av,b + 1{u = v} ru,b ,

(81)

Algorithm 1 Dense ancestral sampling for DAG-SVI Require: Minibatch of data indexes B, number of Monte Carlo samples S, structured posterior parameters (m, Λ), local GP modules 1: Restrict all root design matrices, observations, and masks to indexes i ∈ B 2: Form the dense covariance Σ = Λ−1 b←∅ 3: Initialise sampled latent values F 4: Initialise the running pointwise Gaussian state for each sample-index pair b = (s, i) 5: for each non-root node w ∈ U in topological order do 6: Build the current node inputs from root values and sampled parent values 7: Compute interpolation vectors and residual variances (aw,b , rw,b ) from Eq. (80) 8: Compute base means and covariances using Eqs. (81)–(82) 2 9: Compute µw,b and σw,b by Gaussian conditioning, Eqs. (83)–(84) q 10: Draw Fbw,b = µw,b + σ 2 εw,b , with εw,b ∼ N (0, 1) w,b

11: Store Fbw,b , aw,b , rw,b , and update the running Gaussian state 12: end for b 13: return sampled latent values F

and the base mean µ0u,b = a⊤ u,b mu .

(82)

0 0 Here Σuv is the inducing covariance block between nodes u and v. By vw,<w,b , v<w,<w,b , and µ0<w,b , we denote the row vector, covariance matrix, and mean vector obtained by stacking these scalar quantities over the nodes preceding w in the chosen topological order. At node w, Gaussian conditioning gives −1 0 0 µw,b = µ0w,b + vw,<w,b v<w,<w,b Fb<w,b − µ0<w,b , (83) −1 0 2 0 0 0 σw,b = vww,b − vw,<w,b v<w,<w,b v<w,w,b . (84)

The latent value is then sampled using the pathwise reparametrisation q 2 ε Fbw,b = µw,b + σw,b εw,b ∼ N (0, 1). w,b , This is the pointwise implementation of Prop. 9: for each data index, it samples the latent variables across DAG nodes in topological order, while avoiding cross-index covariance terms that do not enter the likelihood estimator. E.3.2

Sparse implementation

The sparse implementation represents the same Gaussian posterior in canonical form, 1 qH (U) ∝ exp − U⊤ ΛU + h⊤ U , h = Λm, 2

(85)

with Λ sparse on the chordal graph H. The chordal completion ensures that sparse Cholesky elimination can be carried out without introducing fill-in outside H under a perfect elimination order (Lauritzen, 1996; Rue and Held, 2005). The key observation is that conditioning on an already-sampled latent value adds a Gaussian site involving only the corresponding inducing block. If Fbv,b has been sampled and Fbv,b | Uv ∼ N (a⊤ v,b Uv , rv,b ),

then the canonical parameters are updated by av,b a⊤ v,b , rv,b av,b Fbv,b ∆hv,b ← ∆hv,b + . rv,b

∆Λvv,b ← ∆Λvv,b +

(86) (87)

Thus the off-diagonal sparsity pattern is unchanged. For sample-index pair b, let e b = Λ + ∆Λb . Λ Rather than solving with the full information vector Λm + ∆hb , the implementation uses the centred identity e −1 (Λm + ∆hb ) = m + Λ e −1 (∆hb − ∆Λb m) . Λ b b

(88)

The right-hand side in the second term is nonzero only at blocks that have already contributed sites. At node w, define a block vector gw,b by ( aw,b , v = w, [gw,b ]v = 0, v ̸= w. Let e −1 (∆hb − ∆Λb m) , yb = Λ b

e −1 gw,b . zb = Λ b

(89)

The conditional moments needed to sample Fw,b are then µw,b = a⊤ w,b (mw + [yb ]w ) ,

(90)

2 σw,b = rw,b + a⊤ w,b [zb ]w .

(91)

e b . After Both yb and zb are obtained by sparse triangular solves using the current sparse Cholesky factor of Λ sampling Fbw,b , the site update in Eqs. (86)–(87) is added. The site is local in the precision matrix; numerically, the Cholesky factor is updated over the affected part of the elimination tree. KL term.

Both implementations use the same inducing KL. In whitened coordinates, with P = ! Y 1 tr(Σ) + m⊤ m − P + log |Λ| , Σ = Λ−1 . KL qH (U) N (0, I) = 2

w∈U Mw ,

(92)

w∈U

e b used The KL is computed from the variational precision Λ, not from the temporary site-updated precisions Λ −1 inside the ancestral sampler. The dense implementation obtains tr(Σ) after explicitly forming Σ = Λ . The sparse implementation obtains logP |Λ| from the sparse Cholesky factor of Λ, and obtains the diagonal covariance blocks Σww needed for tr(Σ) = w tr(Σww ) by selected-inverse, or Takahashi, recursions (Takahashi et al., 1973; Erisman and Tinney, 1975). Thus the sparse implementation computes the KL without materialising the full covariance matrix. E.3.3

DAG-VI

DAG-VI uses the same pointwise GP conditionals and the same topological ancestral pass, but restricts the inducing posterior to factorise across DAG nodes, Y qVI (U) = N (Uw ; mw , Sw ). w∈U

Consequently, no conditioning on previously sampled node values is performed at the inducing level. At node w, for sample-index pair b, µw,b = a⊤ w,b mw ,

2 σw,b = rw,b + a⊤ w,b Sw aw,b .

(93)

The sampled parent values still enter downstream GP inputs, so DAG-VI propagates marginal uncertainty through the DAG, but it removes posterior coupling between distinct node mechanisms. 57

Algorithm 2 Sparse ancestral sampling for DAG-SVI Require: Minibatch of data indexes B, number of Monte Carlo samples S, structured posterior parameters (m, Λ), chordal graph H, local GP modules 1: Restrict all root design matrices, observations, and masks to indexes i ∈ B 2: Compute the sparse Cholesky factor of Λ on H b←∅ 3: Initialise sampled latent values F 4: For each sample-index pair b = (s, i), initialise site terms ∆Λb = 0, ∆hb = 0, and the corresponding sparse Cholesky state 5: for each non-root node w ∈ U in topological order do 6: Build the current node inputs from root values and sampled parent values 7: Compute interpolation vectors and residual variances (aw,b , rw,b ) from Eq. (80) 8: Using sparse triangular q solves, compute the block readouts in Eqs. (89)–(91) 9: Draw Fbw,b = µw,b + σ 2 εw,b , with εw,b ∼ N (0, 1) w,b

10:

Add the site updates ∆Λww,b ← ∆Λww,b + aw,b a⊤ w,b /rw,b ,

∆hw,b ← ∆hw,b + aw,b Fbw,b /rw,b

11: Update the sparse Cholesky state over the affected elimination-tree region 12: end for b 13: return sampled latent values F

Algorithm 3 Ancestral sampling for DAG-VI Require: Minibatch of data indexes B, number of Monte Carlo samples S, mean-field posterior Q N (mw , Sw ), local GP modules w∈U 1: Restrict all root design matrices, observations, and masks to indexes i ∈ B b←∅ 2: Initialise sampled latent values F 3: for each non-root node w ∈ U in topological order do 4: Build the current node inputs from root values and sampled parent values 5: Compute interpolation vectors and residual variances (aw,b , rw,b ) 2 6: Compute µw,b and σw,b from Eq. (93) q 7: Draw Fbw,b = µw,b + σ 2 εw,b , with εw,b ∼ N (0, 1) w,b

8: end for b 9: return sampled latent values F

E.3.4

Stochastic ELBO estimator

Let

(i) (i) ℓ(i) , w (f ) := log pw Yw | f, Ow (i)

(i)

with ℓw (f ) = 0 whenever Ow = ∅. For a minibatch of data indexes B ⊂ [n] of size B, sampled uniformly, the estimator is S n 1 X X X (i) b(s,i) ℓw Fw . (94) Lbobs = B S s=1 i∈B w∈U

The stochastic ELBO estimate is Lb = Lbobs − K,

(95)

where K is the inducing KL in Eq. (92) for DAG-SVI, or the sum of nodewise KL terms for DAG-VI. If additional per-index latent variables are used, their KL terms are added to K, with the corresponding (i) minibatch scaling. For Gaussian distributions, the expectation of ℓw under the final Gaussian conditional can be evaluated analytically; this Rao–Blackwellised variant reduces Monte Carlo variance but leaves the objective unchanged. 58

Algorithm 4 One stochastic training step Require: Dataset D, minibatch size B, Monte Carlo samples S, variational family method ∈ {DAG-VI, DAG-SVI dense, DAG-SVI sparse} 1: Draw a minibatch of data indexes B ⊂ [n] uniformly, typically without replacement within an epoch 2: if method = DAG-SVI dense then 3: Draw ancestral samples using Algorithm 1 4: else if method = DAG-SVI sparse then 5: Draw ancestral samples using Algorithm 2 6: else 7: Draw ancestral samples using Algorithm 3 8: end if 9: Estimate the observation term using Eq. (94) 10: Compute the analytic KL term K b = Lbobs − K 11: Form L 12: Update variational parameters, inducing locations, and kernel hyperparameters by a stochastic gradient step

E.4

Chordal completion and elimination order

We construct the chordal graph H over the inducing-variable blocks associated with the non-root ancestors of observed nodes. Specifically, we consider An ({w ∈ U : Ow ̸= ∅}) ∩ U, where ancestors include the observed nodes themselves. We then moralise this induced DAG, adding undirected edges between each retained parent–child pair and between all retained co-parents. Root nodes are excluded because they do not carry inducing points. Non-root nodes outside this ancestral set are retained as isolated vertices. When the resulting moral graph is not chordal, we chordally complete it using the MCS-M minimaltriangulation heuristic (Berry et al., 2004). The fill edges returned by MCS-M give an inclusion-minimal triangulation, meaning that no added edge can be removed while preserving chordality. This is a local minimality guarantee and does not imply minimum treewidth, minimum maximum-clique size, or minimum computational cost. Finding an optimal chordal completion under such criteria is computationally intractable in general, with minimum fill-in and bounded-treewidth formulations being classical NP-complete problems (Yannakakis, 1981; Arnborg et al., 1987). If the moral graph is already chordal, no fill edges are added. Clique sizes are therefore induced by moralisation and chordal completion. Alternative completion or ordering heuristics designed to control clique size could be incorporated within the same general construction. After constructing H, we compute a deterministic perfect elimination order (PEO) π = (v1 , . . . , vJ ) using maximum cardinality search on H (Tarjan and Yannakakis, 1985), with ties broken by the fixed declaration order of the non-root nodes. This PEO is computed on the completed graph. For each vi , its later neighbours in π form a clique, and the corresponding elimination front consists of vi together with these later neighbours. DAG-SVI-sparse eliminates blocks in ascending PEO and performs selected inversion and triangular sampling in reverse PEO.

E.5

Computational cost

We report the leading cost of one stochastic ELBO evaluation. Let J = |U| be the number of non-root inducing blocks, let M = maxw Mw , and let K = SB, where S is the number of Monte Carlo samples and B is the minibatch size. For DAG-SVI-sparse, let c be the maximum number of inducing blocks in a clique of the chordal graph H. We assume comparable block sizes and scalar node outputs. DAG-VI factorises across inducing blocks, giving O(JM 3 + KJM 2 ) time and O(JM 2 ) memory. DAG-SVI-dense materialises the full covariance Σ = Λ−1 . Its cost is O((JM )3 + KJ 2 M 2 + KJ 3 ) time and O((JM )2 ) memory. The first term is the global dense Gaussian computation, while the remaining terms come from dense cross-node covariance contractions and Gaussian conditioning during the ancestral pass. 59

DAG-SVI-sparse instead exploits the chordal precision structure and never materialises Λ−1 . Sparse Cholesky, log-determinants, and the selected-inverse/Takahashi recursions used for the KL trace term cost O(Jc2 M 3 ) time and O(JcM 2 ) memory. The pointwise GP propagation through the DAG contributes O(KJM 2 ). The collapsed sampler also performs local Gaussian conditioning updates along the elimination tree. IfP ρt denotes the number of Cholesky block columns, or supernodes, affected at ancestral step t, and J RH := t=1 ρt , then these updates add O(KcM 2 RH ) to the sparse cost. Hence the sparse implementation costs O(Jc2 M 3 + KJM 2 + KcM 2 RH ), with persistent memory O(JcM 2 ), up to minibatch-specific temporary storage. This expression shows how DAG-SVI-sparse benefits from local graph structure. In particular, the base Gaussian computation scales with the maximal clique size c, while the collapsed-sampling overhead is governed by the elimination-tree update profile RH . In locally sparse DAGs, c is small and RH grows slowly with J. For example, in the balanced branching-tree setting of Fig. 2, each node has at most one parent, so the moralised graph remains a tree and c = 2. With a balanced elimination profile, the affected update regions grow with the tree depth, giving RH = O(J log J), and hence O(JM 3 + KJ log(J)M 2 ), which is substantially below the dense scaling as J grows. Finally, recall that for a chain DGP, H reduces to a block-tridiagonal structure and our structural family specialises to Ustyuzhaninov et al. (2020, Sec. 4.1). The corresponding marginalised conditionals can then be computed analytically, saving computation.

E.6

Explaining away

Compared with a standard chain DGP, a DAG-DGP can exhibit posterior coupling between independent mechanisms that share an observed child. This behaviour is the classical explaining-away effect (Pearl, 1988; Lauritzen, 1996). The linear Gaussian case gives a simple illustration. Linear Gaussian case.

Consider the scalar collider

2 A ∼ N (0, σA ),

2 B ∼ N (0, σB ),

2 C | A, B ∼ N (αA + βB, σC ),

with A and B independent a priori. Once a value c of the child variable is observed, the posterior over (A, B) is Gaussian with precision −2 −2 −2 σA + α2 σC αβσC QAB|c = (96) −2 −2 −2 . αβσC σB + β 2 σC Inverting (96) gives Cov(A, B | c) = −

2 2 αβ σA σB 2 2 2 . σC + α 2 σA + β 2 σB

(97)

Thus, whenever αβ > 0, the two parents become negatively correlated a posteriori. Intuitively, once one branch explains a substantial part of the observed value c, less support is needed from the other. The observation at the child therefore induces posterior dependence between parents that are independent a priori. Explaining away in DAG-DGPs.

The same mechanism appears in the DAG-DGP collider w1 → w3 ← w2 ,

where w1 and w2 are two parent mechanisms and observations are attached only to the child node w3 . For i ∈ [n], write the latent recursion as (i) (i) (i) (i) (i) (i) F(i) = f F , F = f F , F = f F , F (98) w w w w1 w2 w3 w1 w2 . 1 2 3 Pa(w1 ) Pa(w2 ) For notational simplicity, suppressing the root inputs and any other possible upstream variables outside the collider, the prior latent law factorises as p0 (Fw1 , Fw2 , Fw3 ) = p0 (Fw1 ) p0 (Fw2 ) p0 (Fw3 | Fw1 , Fw2 ).

Suppose that observations are available only at the child, through a nodewise conditional distribution pw3 (Yw3 | Fw3 , Ow3 ). The posterior marginal over the two parent branches is then p(Fw1 , Fw2 | Yw3 ) ∝ p0 (Fw1 ) p0 (Fw2 ) Z × pw3 (Yw3 | Fw3 , Ow3 ) p0 (Fw3 | Fw1 , Fw2 ) dFw3 .

(99)

The integral in (99) depends jointly on Fw1 and Fw2 , and therefore does not factorise in general. This is the DAG-DGP analogue of explaining away: once the observed child is partly accounted for by one branch, the posterior mass over the other branch shifts accordingly. The effect is especially transparent when the child kernel contains separate contributions from the two parents, for instance under additive fusion. In this case the child receives two distinct nonlinear contributions whose combined effect is constrained by the observations at w3 . The posterior therefore induces dependence between the two parent branches even though the GP modules fw1 and fw2 are independent under the prior. DAG-VI cannot retain explaining away. We now show that DAG-VI cannot represent the posterior coupling required by explaining away. Consider again the collider discussed in the main paper w1 → w3 ← w2 , and assume that observations are available only at the child node w3 . Under DAG-VI, the inducing posterior factorises across nodes, Y qDAG-VI (U) = qw (Uw ). w∈U

Together with the DAG-DGP conditionals, this gives qDAG-VI (Fw1 , Uw1 , Fw2 , Uw2 , Fw3 , Uw3 ) = p(Fw1 | Uw1 )qw1 (Uw1 ) p(Fw2 | Uw2 )qw2 (Uw2 )

(100)

× p(Fw3 | Uw3 , Fw1 , Fw2 )qw3 (Uw3 ). Marginalising the child variables yields qDAG-VI (Fw1 , Uw1 , Fw2 , Uw2 ) = p(Fw1 | Uw1 )qw1 (Uw1 ) p(Fw2 | Uw2 )qw2 (Uw2 ) Z Z × qw3 (Uw3 ) dUw3 p(Fw3 | Uw3 , Fw1 , Fw2 ) dFw3 .

(101)

Both integrals are equal to one. Hence qDAG-VI (Fw1 , Uw1 , Fw2 , Uw2 ) = qDAG-VI (Fw1 , Uw1 ) qDAG-VI (Fw2 , Uw2 ).

(102)

In particular, qDAG-VI (Fw1 , Fw2 ) = qDAG-VI (Fw1 ) qDAG-VI (Fw2 ).

(103)

Thus DAG-VI preserves marginal uncertainty along each branch, but assigns no posterior dependence between the two parents. It therefore cannot represent the explaining-away dependence induced by observations at w3 . A topological directed factorisation is also insufficient. One might instead consider a directed variational posterior that follows the topological order of the DAG, Y qtop (U) = qw Uw | {Up }p∈Pa(w) , (104) w∈U

which is the direct DAG analogue of autoregressive or layer-wise posterior factorisations used in chain settings (Ustyuzhaninov et al., 2020; Ober and Aitchison, 2021). This approach is also unable to capture explaining-away. For the same collider w1 → w3 ← w2 , (104) gives qtop (Uw1 , Uw2 , Uw3 ) = qw1 (Uw1 ) qw2 (Uw2 ) qw3 (Uw3 | Uw1 , Uw2 ). 61

Marginalising the child gives Z qtop (Uw1 , Uw2 ) = qw1 (Uw1 )qw2 (Uw2 )

qw3 (Uw3 | Uw1 , Uw2 ) dUw3

(105)

= qw1 (Uw1 )qw2 (Uw2 ). The child conditional can model how Uw3 depends on its parents, but it does not create marginal posterior dependence between the co-parents once the child is integrated out. Explaining away requires exactly such dependence: conditioning on observations at w3 couples the plausible contributions of w1 and w2 . This is why the structured approximation is based instead on the moralised ancestral graph. For the collider w1 → w3 ← w2 , moralisation adds the co-parent edge {w1 , w2 }, allowing the approximate posterior over inducing variables to retain the posterior dependence induced by the observed child.

Stochastic Deep Gaussian Processes over Graphs as a Special Case of DAG-DGP

Stochastic Deep Gaussian Processes over Graphs (DGPG) (Li et al., 2020) were introduced for a modelling task different from ours, namely learning maps between input and output signals defined on the vertices of a fixed graph. In that setting, the graph indexes the components of each signal and specifies which neighbouring components are used by each graph-indexed GP module. We show that this construction is nevertheless contained in the DAG-DGP framework. After unrolling the base graph across depth, any DGPG model can be represented as a DAG-DGP on a layered DAG, with deterministic roots F 0 = X, concatenation fusion rule at each non-root node, and observations restricted to the terminal layer. Under this identification, the DGPG ELBO corresponds to the DAG-VI objective. This observation positions DGPG as one particular graph-structured architecture within the broader DAG-DGP class. Thus, the DAG-DGP framework is strictly more general at the modelling level, since it allows arbitrary DAG structures, heterogeneous fusion rules, and observations at arbitrary nodes. It is also more general at the inferential level, since the same embedded DGPG model can be equipped with our DAG-SVI objective, which applies directly to this model.

F.1

The DGPG model

Setup. Following Li et al. (2020), the dataset is D = {G, Ψ, Φ}, where G = ⟨V, E⟩ is a graph with vertices V = {v1 , . . . , vK } and edges E ⊆ V × V . For vk ∈ V , we write PaG (k) := {j : (vj , vk ) ∈ E} for the parent indices of vk in the base graph, with self-loops allowed. Each input ψ ∈ Ψ is a graph signal ψ : V → Rdin , and each output ϕ ∈ Φ is a graph signal ϕ : V → Rdout . The learning task is to infer a map h : Ψ → Φ, taking an input graph signal to an output graph signal. With N training signals, the inputs and outputs are stacked row-wise as X = (x1 , . . . , xN )⊤ and Y = (y1 , . . . , yN )⊤ , where xi ∈ RKdin and yi ∈ RKdout concatenate the vertex-wise features of the i-th input and output signals. Generative model. DGPG stacks L layers of graph-indexed GP mappings. For layer l the latent matrix is F l ∈ RN ×Kdl (with per-node dimension dl , d0 = din , dL = dout ) and F 0 = X; each layer is augmented with M inducing inputs Z l and inducing outputs U l . For a matrix M ∈ {X, Y, F, Z, U }, M l,k denotes the sub-block of layer l associated with vertex vk , M l,PaG (k) the concatenated sub-block at the base-graph parents of vk , and Mil,k the i-th row of M l,k . Assuming Gaussian GP priors and inducing outputs U l,k that are independent across layers and vertices, the joint density factorises over observations, layers, and vertices as p Y, {F

l,k

}l,k =

N Y K Y

p ynk | FnL,k

n=1 k=1

L Y K Y

(106) p F

l,k

l−1,Pa(k)

p U

l,k

l−1,Pa(k)

l=1 k=1

with F 0 = X. Crucially, the GP module at node k in layer l acts on the concatenation F l−1,Pa(k) of its graph-parent signals from the previous layer. The notation of Li et al. (2020) follows the standing convention of Salimbeni and Deisenroth (2017), according to which the semicolon separates fixed inputs and kernel-design quantities from random quantities being conditioned on. In our DAG-DGP notation, this fixed dependence is kept implicit. Variational family and recursive sampling. DGPG uses the doubly-stochastic family of Salimbeni and Deisenroth (2017), retaining the GP conditionals and a factorised Gaussian inducing posterior, L Y K Y q {F l,k , U l,k }l,k = p F l,k | U l,k ; F l−1,Pa(k) , Z l−1,Pa(k) q U l,k ,

(107)

l=1 k=1

where q(U l,k ) := N U l,k | ml,k , S l,k . Marginalising each U l,k yields per-node Gaussian marginals q(F l,k ) = N (F l,k | µ̃l,k , Σ̃l,k ) whose moments depend only on the parent states. Consequently the terminal-layer 63

marginal q(FiL,k ) depends only on the ancestors of ⟨L, k⟩ and can be drawn recursively across depth by the reparameterisation, q l−1,Pa(k) l−1,Pa(k) b l−1,Pa(k) l,k l,k b b Fi = µml,k ,Z l−1,Pa(k) Fi + ϵi ⊙ ΣS l,k ,Z l−1,Pa(k) Fbi , Fi , (108) where ϵl,k i ∼ N (0, Idl ) and µ and Σ are the sparse variational predictive mean and covariance. Evidence lower bound. LDGPG =

N X K X

With (106)–(107), the DGPG ELBO is L X K X Eq(FnL,k ) log p(ynk | FnL,k ) − KL q(U l,k ) p(U l,k ; Z l−1,Pa(k) ) .

n=1 k=1

F.2

(109)

l=1 k=1

DGPG as a particular case of DAG-DGP

e is a layered DAG and that DGPG coincides with the Proposition 10 shows that the depth-unrolled graph G DAG-DGP supported on it, under a concatenation fusion rule and observations available only at the terminal layer. We first define the depth-unrolling of the base graph and the corresponding notion of layered DAG. We then prove that the unrolling of any base graph is layered. Finally, we show that, under this specialization, the DGPG augmented joint and ELBO are recovered exactly as the corresponding DAG-DGP joint and DAG-VI objective. Figure 12 shows an example of this. Throughout this subsection, parent sets in the original DGPG e are denoted by Pa . base graph G are denoted by PaG , whereas parent sets in the depth-unrolled DAG G e G Definition 7 (Depth-unrolling of the base graph). Given a DGPG model with base graph G = ⟨V, E⟩ and e = (Ve , E) e with depth L, the depth-unrolling of G is the directed graph G e = (l − 1, j) → (l, k) : 1 ≤ l ≤ L, j ∈ PaG (k) , Ve = (l, k) : 0 ≤ l ≤ L, k ∈ V , E so that PaG e (l, k) = {(l − 1, j) : j ∈ PaG (k)},

l ≥ 1.

e = {(0, k) : k ∈ V } the roots, Ue = {(l, k) : 1 ≤ l ≤ L, k ∈ V } the non-roots, and Al = {(l, k) : k ∈ V } We call R the l-th depth slice. Definition 8 (Layered DAG). A DAG G = (V, E) is layered if there exists a progressive antichain decomposition V = A0 ⊔ A1 ⊔ · · · ⊔ AL such that every edge connects consecutive antichains: E⊆

L−1 [

Aℓ × Aℓ+1 .

ℓ=0

Lemma 4 (The depth-unrolling of any DGPG base graph is a layered DAG). For any DGPG base graph G e of Definition 7 is a layered DAG. and depth L, the depth-unrolling G e takes the form Proof. By construction, every edge of G (l − 1, j) → (l, k),

1 ≤ l ≤ L,

j ∈ Pa(k),

and therefore increases the depth index by one. Hence every directed path strictly increases the depth index, e is acyclic and no two vertices in the same slice Al are comparable. Thus each Al is an antichain. so G Moreover, the antichains satisfy Ve = A0 ⊔ A1 ⊔ · · · ⊔ AL , and, again by construction, e⊆ E

L [

Al−1 × Al .

l=1

Therefore A0 , . . . , AL form a progressive antichain decomposition satisfying the consecutive-layer condition. e is a layered DAG. Hence G 64

Proposition 10 (DGPG is a DAG-DGP and DAG-VI recovers its ELBO). For any DGPG base graph G, the e of Definition 7, with deterministic DGPG model of (Li et al., 2020) is the DAG-DGP on the layered DAG G 0 roots F = X, observations supported only on the terminal antichain-layer AL , and concatenation fusion at each non-root node w = (l, k), meaning that the parent tuple is treated as a single concatenated input and the local kernel is obtained by applying a standard kernel to this concatenated parent state. Under this identification, the DGPG joint (106) is the augmented DAG-DGP joint, and its objective equals the corresponding DAG-VI bound, that is, LDGPG = LVI . e is a layered DAG, so the DAG-DGP construction on it is well defined. We identify Proof. By Lemma 4, G e with the DGPG GP module at layer l and graph vertex vk . Under this each non-root node w = (l, k) of G identification, Uw = U l,k ,

Zw = Z l−1,PaG (k) ,

l−1,PaG (k) {Fp : p ∈ PaG . e(w)} = F

The root nodes are deterministic and given by F(0,k) = X 0,k , so that F 0 = X. For every non-root node (l, k), the concatenation fusion rule makes the DAG-DGP parent input the DGPG input F l−1,PaG (k) . e is The augmented DAG-DGP joint on G Y p(D, F, U ) = p0 (Uw ) p0 Fw | {Fp }p∈Pa (w) , Uw pw (Yw | Fw , Ow ), G e e w∈U where pw ≡ 1 whenever Ow = ∅. We set O(l,k) = ∅ for l < L and p(L,k) (Y(L,k) | F(L,k) ) =

N Y

p(ynk | FnL,k )

n=1

at the terminal layer. Hence, substituting w = (l, k) in the preceding augmented DAG-DGP joint factorisation gives p0 (U(l,k) ) = p U l,k ; Z l−1,PaG (k) , and p0 F(l,k) | {Fp }p∈Pa ((l,k)) , U(l,k) = p F l,k | U l,k ; F l−1,PaG (k) , Z l−1,PaG (k) . G e Therefore p(D, F, U ) =

L Y K Y

K Y N Y p(ynk | FnL,k ), p U l,k ; Z l−1,PaG (k) p F l,k | U l,k ; F l−1,PaG (k) , Z l−1,PaG (k) k=1 n=1

l=1 k=1

which is the DGPG joint in (106), up to reordering of factors. e gives For the variational family, the DAG-VI construction (here on G) Y q(F, U ) = q(U ) p0 Fw | {Fp }p∈Pa (w) , Uw . G e e w∈U Choosing the mean-field inducing posterior (DAG-VI) Y q(U ) = qw (Uw ), qw (Uw ) = N (U l,k | ml,k , S l,k ), e w∈U recovers the DGPG variational family (107). The ancestral sampling recursion is also the same, since the DAG-DGP parent input at (l, k) is precisely the concatenated DGPG state F l−1,PaG (k) .

v10

v20

v30

v11

v21

v31

v12

v22

v32

(a) G

(b) DGPG

e (c) DAG-DGP on G

Figure 12: Three views built from one base graph. (a) The base graph G on K = 3 vertices, with a single self-loop on v1 (highlighted) and a feedback pair v2 ↔ v3 , hence cyclic. (b) DGPG: a standard chain DGP whose state at every layer is a signal on G—shown as an identical small copy of G inside each layer node—and e each node v ℓ whose layer map f ℓ is wired by G. (c) The DAG-DGP defined on the depth-unrolled graph G: k denotes the copy of the original node vk at unrolled layer ℓ. In both (b) and (c), diamonds denote root/input nodes or layers, circles denote latent nodes or layers, and squares denote observed nodes or layers. It remains to match the objectives. Starting from the DGPG ELBO, we have LDGPG =

K X N X

Eq(FnL,k ) log p(ynk | FnL,k )

−

L X K X

k=1 n=1

l=1 k=1

K X N X

Eq(FnL,k ) log p(ynk | FnL,k ) −

k=1 n=1

KL q(U l,k ) ∥ p(U l,k ; Z l−1,PaG (k) )

KL qw (Uw ) ∥ p0 (Uw )

e w∈U

K X

X EqVI (F(L,k) ) log p(L,k) (Y(L,k) | F(L,k) , O(L,k) ) − KL qw (Uw ) ∥ p0 (Uw ) k=1 e w∈U X X = EqVI (Fw ) log pw (Yw | Fw , Ow ) − KL qw (Uw ) ∥ p0 (Uw ) w: Ow ̸=∅ e w∈U =

= LVI .

Remark 4 (DGPG is a strict subclass of DAG-DGP). Let MDGPG denote the DGPG model class. Let MDAG-DGP denote the general DAG-DGP model class as presented in this paper. Proposition 10 shows that every DGPG specification is a DAG-DGP specification, i.e. MDGPG ⊆ MDAG-DGP . The inclusion is strict because the embedding fixes several choices that are free in the general DAG-DGP formulation. In particular, DGPG fixes the fusion rule to concatenation, whereas DAG-DGPs allow heterogeneous additive, product, or domain-specific fusion rules; and it fixes the observation pattern to the terminal antichain, whereas DAG-DGPs allow observations at internal nodes. Both degrees of freedom are used in the multi-fidelity and protein-signalling models of the main text. The simplest separation, however, relates to the graphical structure. By Lemma 4, every depth-unrolling produces a layered DAG. Hence any non-layered DAG is outside the DGPG domain. For example, consider the three-nodes DAG H with edges a → b, b → c, and the skip edge a → c. This DAG is not layered. Indeed, recall that in a layered DAG every edge must connect two consecutive antichains. The path a → b → c would require a, b, c to lie in three consecutive antichains, whereas the edge a → c would require a and c to lie in consecutive antichains. These two requirements are incompatible. Nevertheless, a DAG-DGP is well defined directly on H for any admissible kernels and fusion rule.

Experiments

Experiments were run in float64 on NVIDIA TITAN RTX and GeForce RTX 2060 SUPER GPUs. To indicate the computational scale of the real-data experiments, we report representative wall-clock training times. For Sachs, on the extrapolation task and using the current 10000-step joint-training protocol, a single split required 29.2 ± 0.3 minutes for DAG-VI and 33.0 ± 0.6 minutes for DAG-SVI on an NVIDIA TITAN RTX, averaged over the completed projection splits. For the published HeavyIon split, averaged over the five reported seeds, end-to-end runs required 6.3 ± 0.7 minutes for DAG-VI and 8.8 ± 1.0 minutes for DAG-SVI. We now provide more details on our experiments.

G.1

Branching-tree ELBO scaling

Experimental design. Figure 2 benchmarks the cost of ELBO evaluation on a synthetic binary branchingtree DAG, used to compare DAG-VI with the dense and sparse structured backends of DAG-SVI. The DAG has one two-dimensional root input X and depth d ∈ {2, . . . , 9}. Every non-root node has exactly one parent and splits into two children, so the graph contains 2d+1 − 2 latent non-root nodes in total, with the 2d leaves observed. For each depth, we generate 20 training cases. Root inputs are sampled uniformly from [−2, 2]2 . Each latent node is then generated from a one-parent nonlinear GP draw using random Fourier features with variance 1, lengthscale 1.2 at the first latent layer and 0.9 thereafter, followed by a mild Gaussian perturbation and a tanh squashing nonlinearity. The observed leaves are obtained by adding independent Gaussian noise with standard deviation 0.05. Models and timing protocol. All methods use the same nodewise GP architecture, with an ARD RBF kernel at each latent node and 20 inducing points per latent node, initialized from training parent-state inputs. We compare three ELBO evaluators: DAG-VI, DAG-SVI with dense marginalization, and DAG-SVI with sparse marginalization. To isolate the marginalization cost only, the dense and sparse structured models are created from the same initialized parameter state; they differ only in the backend used to evaluate the ELBO. For each depth, we evaluate the full-batch ELBO on all 20 training cases using a single Monte Carlo sample, and record wall-clock time on GPU. Relation to the main implementation claims. This benchmark is intended to validate that the sparse backend can substantially reduce ELBO-evaluation cost relative to dense structured inference as the inducing dimension grows, while DAG-VI is the most efficient algorithm that we proposed, being its variational structure limited. The branching-tree topology gives a sparsity structure that the algorithm uses.

G.2

Theory Validation

Simulation design. Figure 3 reports finite-depth Monte Carlo checks of the prior and posterior non-collapse statements in Section 4. Panels (a)–(c) are generated under the DAG-DGP prior of Eq. (2), by ancestral sampling on the simulated DAG in topological order. Rather than sampling whole GP paths, we sample the finite two-case Gaussian conditionals induced at each node by its realised parent states, which is the finite-dimensional prior construction studied in Section 4 and formalised in Appendix A.3. Throughout, we refer to the quantities defined there. Panel (d) is generated from the filtering and forward-predictive laws used in Section 4.3 and Appendix D. For simplicity, all experiments use scalar node outputs and squared-exponential kernels, with the hyperparameters stated below. Panel (a): repeated separating nodes. Panel (a) validates Theorem 1 on a layered DAG with depth L = 60 and antichain width W = 16. We induce separation through root connections. The two cases have a single deterministic root coordinate, fixed at r(a) = 0 and r(b) = 0.3. For each Monte Carlo realisation, every non-root node selects one parent uniformly from the preceding antichain. Non-separating nodes use a squared-exponential kernel on this parent coordinate, with variance 1 and lengthscale 1. Separating nodes use the same parent component, additively fused with a squared-exponential root component with the same hyperparameters. This means that the simulated separating nodes are precisely instances of the additive root-retention mechanism in Remark 2, with the associated separation constant obtained from Eq. (32). The one-step antichain probability used in the theorem curve is therefore the pε of Eq. (12). For each s ∈ {0, . . . , 10}, exactly s separating nodes are chosen in each antichain. The case s = 0 is included only as a no-separation baseline. For each s, we simulate 1500 independent prior realisations and compute the finite-window version of the frequency appearing in Theorem 1, L−1

ρ̂ε =

X 1 1{Mℓ > ε}, L − L0

L0 = 30,

ε = 0.30.

(110)

ℓ=L0

The plotted points and error bars are the empirical mean and standard deviation of ρ̂ε across realisations. The curve is the lower bound in Eq. (1), evaluated with the above root-retaining separation constant. Panels (b) and (c): topology effects. Panel (b) validates the indegree recursion of Proposition 4 in the layered radial block of Assumption 1. We use depth 50, width 128, kernel variance 1, parent lengthscale 2.0, root lengthscale 0.8, and root gap 1.0. The first antichain is initialised by applying the same two-case Gaussian prior rule to the two fixed root inputs; hence the initial contrast variance is determined by the root gap and the root lengthscale. All later node values are propagated using the two-point Gaussian law in Lemma 2. For a fixed indegree k ∈ {1, 2, 4, 8}, every node receives k parents from the previous antichain. To keep all nodes statistically symmetric and avoid boundary effects, we assign these parents by taking k consecutive nodes and wrapping around at the edge of the layer. This is only a convenient regular layered DAG used to isolate the effect of indegree. Under product fusion, the product of squared-exponential parent correlations is the radial squared-exponential kernel on the concatenated parent state, so the simulation is the squared-exponential special case of Assumption 1. We estimate Cℓ by averaging the simulated squared contrasts over nodes and over 2400 Monte Carlo draws. The curves illustrate the monotonicity in k in Proposition 4 and the corresponding contraction discussion in Corollary 4: larger indegree broadens the recurrence relative to the chain case. Panel (c) validates the outdegree mechanism in Proposition 6. We simulate the designated b-ary branching subgraph of Appendix C.3, with b ∈ {1, 2, 3, 5}, maximum depth 9, threshold t = 0.45, variance 1, lengthscale 1.0, and root gap 1.0. Children are conditionally independent given their parent contrast and are sampled using Lemma 2. For each branching factor and each L ∈ {3, 6, 9}, the plotted value estimates the finite-depth survival event corresponding to Proposition 6, namely that the antichain maximum remains above threshold at every depth up to L. The estimates are based on 700 Monte Carlo draws. This isolates the claim that larger outdegree creates more parallel opportunities for a threshold-size contrast to survive. Panel (d): intermediate observations as stochastic skip connections. Panel (d) validates Theorem 2 using the Gaussian refresh setting of Eq. (9) and Corollary 6. At each observed source, the two-case source 68

contrast is sampled from the Gaussian filtering law determined by source contrast variance 1, observation-noise variance 0.20, observed gap 1.0, and threshold ε = 0.35. The corresponding source strength is the qj appearing in Theorem 2. Downstream routes use squared-exponential route kernels with variance 0.60 and lengthscale 0.55; their retention factors are computed from Definition 6 and Proposition 7. We compare two geometries. The first is a chain, giving the s = 1 case of Theorem 2. The second is a disjoint-route DAG, as in the geometry of Fig. 1, with two observed sources and two pairwise interior-disjoint admissible routes, in the sense of Definitions 3 and 4. In this case, the theorem bound is computed as one minus the probability that both routes fail. Thus, each refreshed route contributes its own success probability, and the two contributions combine through the “at least one route succeeds” term in Eq. (8). The empirical curves estimate the forward-predictive probability Πℓ0 →ℓ0 +h (Mℓ0 +h > ε) for h = 0, . . . , 6, using 50000 posterior-predictive draws. In the disjoint-route simulation, each target receives its main parent along the designated route and a weak additive parent from the opposite refreshed source, with variance 0.02 and lengthscale 0.55. This makes the simulated target a simple multi-parent DAG node, while the theorem bound is evaluated using only the designated disjoint routes. Since the extra parent enters additively, it only contributes additional non-negative contrast variance and is not needed to certify route retention. The plotted theorem curve is therefore the lower bound for the designated routes. The panel validates both parts of the stochastic-skip statement: noisy intermediate observations refresh source contrasts, and multiple disjoint refreshed routes increase the downstream probability that at least one contrast survives.

Figure 13: Illustration of the prior-theory experiments.

G.3

Latent-Collider experiment

Data-generating process.

We consider the collider DAG x1 → w1 → w3 ← w2 ← x2 ,

where only the child node C is observed. The root inputs x1 and x2 are placed on fixed asymmetric onedimensional grids with n1 = 15 and n2 = 11 points, respectively, yielding 15 × 11 = 165 observed locations for C. The latent parent functions w1 (x1 ) and w2 (x2 ) are sampled independently from zero-mean Gaussian processes with RBF kernels of lengthscale 0.4 and variance 1. The child latent surface is then sampled from a GP defined on the parent pair (w1 , w2 ), with fusion kernel KC (a, b), (a′ , b′ ) = kRBF (a + b, a′ + b′ ), using child lengthscale 0.60 and variance 1. Observations are generated as y(x1 , x2 ) = w3 (x1 , x2 ) + ε,

ε ∼ N (0, 0.0152 ).

We deliberately use a small observation-noise level because the experiment is designed to test explaining-away behaviour. Indeed, as suggested by the linear-Gaussian covariance of the child in Eq. (97), lower observation noise is expected to induce stronger posterior dependence between the latent parents. This construction induces a genuine explaining-away geometry because multiple parent configurations (w1 , w2 ) can produce nearly the same child value through the fused direction w1 + w2 . Model training and evaluation. We compare DAG-VI and DAG-SVI on the same DAG, kernels, inducing locations, and optimisation setup; only the variational family differs. Both methods use 15 inducing points for each parent node and a 7 × 7 inducing grid for the child node in the latent (A, B) plane. Kernel hyperparameters, the observation-noise variance, and inducing locations are fixed throughout training, such that the comparison isolates the effect of the variational family at a fixed training budget. Training is full-batch for 3000 gradient steps with learning rate 0.01 and 4 Monte Carlo samples per ELBO estimate. Summary statistics are aggregated over 3 random seeds. Visualisation of explaining-away geometry. The left panel of Fig. 14 is designed to isolate posterior geometry from posterior scale, already represented in Fig. 4 of the main paper. We therefore select a small set of representative interior observation locations on the child surface, chosen to span the input domain and to cover a range of local w3 -level-set orientations while avoiding boundary cases where the geometry is less informative. For each selected case (x1,i , x2,j ), we extract posterior draws of the corresponding parent states (w1 , w2 ). To display these local posterior clouds in the original input space, we map latent perturbations back to input space using the local sensitivities of the true parent functions, (s)

(s)

(s) ∆x1 =

w1,ij − w̄1,ij , ∂w1 (x1,i )/∂x1

(s) ∆x2 =

w2,ij − w̄2,ij , ∂w2 (x2,j )/∂x2

where w̄1,ij and w̄2,ij denote the posterior means at the selected case. Each cloud is then recentered at its corresponding observation location and rescaled isotropically for display. This removes absolute scale differences between DAG-VI and DAG-SVI, so the panel emphasises the orientation of posterior uncertainty relative to the true w3 -level set passing through that observation. The goal is to visualise whether the posterior captures the explaining-away direction: DAG-SVI aligns local uncertainty with the child level set, whereas DAG-VI remains close to isotropic. Visualisation of compositional uncertainty. The right panel of Fig. 14 complements the geometry plot by restoring posterior scale on the natural domains of the parent functions. For each method, we draw posterior samples of the latent parent functions w1 (x1 ) and w2 (x2 ), and summarise them by pointwise credible bands and posterior means. Unlike the explaining-away panel, no recentering or display rescaling is applied here, so the figure reflects the actual posterior spread learned by each variational family. This panel is therefore intended to visualise compositional uncertainty: not having observations, the latent functions for w1 and w2 70

Latent w1

Explaining Away in Input Space avg |posterior corr(w1, w2)| DAG-VI: 0.02 DAG-SVI: 0.93

DAG-VI: posterior var(w1) = 1.01e-04 DAG-SVI: posterior var(w1) = 1.60e-01

1.0

1 0

0.5

DAG-VI DAG-SVI True w3 level set

0.0

1.00

0.75

0.50

0.25

0.00

0.25

0.50

0.75

1.00

Latent w2

4 3 2

0.5

1 0

1.0 1.00

1 0.75

0.50

0.25

0.00

0.25

0.50

0.75

1.00

DAG-VI: posterior var(w2) = 6.67e-05 DAG-SVI: posterior var(w2) = 1.65e-01 1.0

DAG-VI

0.5

0.0

DAG-SVI x2

0.5

1.0

Figure 14: Explaining-away and compositional uncertainty in the synthetic V-structure (Section G.3). Left: projected and rescaled local posterior samples of (A, B) in input space. Under DAG-SVI, samples concentrate along the true C-level set, recovering the explaining-away geometry; under DAG-VI, they remain nearly isotropic. Rescaling is used only to expose the geometric behaviour of the two posteriors on a comparable scale—the raw variability of DAG-VI is in fact very small, as shown on the right. Right: DAG-SVI retains uncertainty over the latent factorisation of C, whereas DAG-VI collapses to a single point-like representation. should preserve posterior uncertainty, as explained by Ustyuzhaninov et al. (2020) for the standard DGPs. DAG-SVI preserves substantially more posterior variance over w1 and w2 , whereas DAG-VI collapses toward an almost deterministic decomposition.

G.4

Sachs Flow Cytometry Experiment

Interpolation and extrapolation tasks on Sachs. For the interpolation task, we use the standard random 80/20 train–test split described in the main text. For the extrapolation task, we adapt the projection-based protocol of Lindinger et al. (2020) to the Sachs DAG. This induces a train–test split in which the test set lies farther from the training support, requiring generalization beyond the region covered by the training data. We use 32 inducing points for each latent node. Consider xi = (PKCi , Plcgi ) ∈ R2 as the observed inputs with index i. We standardize these two coordinates over the full dataset, sample a random unit direction w ∈ R2 , and project each case onto the scalar coordinate zi = x̃⊤ i w. We then order the observations by zi , train on the central 80%, and test on the lower and upper 10% tails, yielding an 80/20 split with train and test separated along a one-dimensional projection of the input space. Note that, as in the interpolation task, PKA and Raf remain observed in both training and test sets. We report results averaged over five random projection seeds. Model training and evaluation. All nodes use an RBF kernel initialised with an ARD lengthscale of 0.1, and observed nodes have a Gaussian likelihood with noise initialised to 0.01. We use Adam (Kingma and Ba, 2015) with a fixed learning rate of 1e − 3. For both DAG-VI and DAG-SVI we first train node-wise with a meanfield approximate posterior per node for 1000 steps. And then jointly train for 10000 steps. In the joint training phase we hold the inducing locations and likelihood noise for 40% of the steps. For the extrapolation task, we keep the same model and optimisation settings, but train both methods jointly for 10000 steps without the node-wise pretraining stage. Table 3: Full test results for the Sachs extrapolation task, averaged over five random projection splits. Lower is better for RMSE, CRPS, and NLPD, while higher is better for PICP. Model DAG-VI DAG-SVI

RMSE ↓

CRPS ↓

NLPD ↓

PICP ↑

0.7371 ± 0.0655 0.6423 ± 0.0307

0.3252 ± 0.0232 0.2995 ± 0.0065

0.9104 ± 0.0985 0.8274 ± 0.0496

0.7980 ± 0.0139 0.8395 ± 0.0115

Here RMSE and CRPS are the standard metrics already defined in the main text, while NLPD denotes the negative log predictive density of the posterior predictive distribution on the test target and PICP the empirical prediction interval coverage probability.

G.5

Heavy-ion collision

Model training. We use the elicited multi-fidelity DAG of Figure 6 on a shared nine-dimensional input space, with two lower-fidelity latent nodes L1 and L2 feeding the high-fidelity node H. In all protocols, L1 and L2 use 200 inducing points each, while H uses all available high-fidelity training locations in the current split or fold: 25 on the published split, 20 in each fold of the repeated 5-fold protocol, and 90 in each fold of the pooled 10-fold protocol. Inducing locations are fixed after initialization. For H, we consider two initialization schemes. In the predmean variant, inducing inputs are initialized at sampled high-fidelity input locations augmented with the predictive means of pre-trained lower-fidelity models. In the free variant, the high-fidelity inducing inputs are initialized at sampled high-fidelity input locations, with the lower-fidelity coordinates initialized freely rather than from lower-fidelity predictive means. We follow the multifidelity inducing-input construction described by (Cutajar et al., 2019) (Section 4.4), which is motivated by the practical difficulty of freely optimizing such augmented inducing representations. All runs use full-batch Adam in double precision with learning rate 0.01, likelihood learning-rate multiplier 0.1, and 10 Monte Carlo samples per ELBO estimate. DAG-VI first pre-trains each node independently as an SVGP and then optimizes the joint DAG objective. DAG-SVI follows the same nodewise warmup, then trains an auxiliary mean-field model for 2000 joint steps to initialize the structured posterior, and finally optimizes the structured objective. On the published split used in Table 2, both methods use 500 pretraining steps per node; DAG-VI then runs 12000 joint steps, whereas DAG-SVI runs 10000 structured steps after the 2000-step mean-field warm start. Results are reported over five random seeds. In the high-fidelity-scarce protocol of Table 4, we perform 10 repeated 5-fold cross-validation runs over the 25 high-fidelity observations only, using 2000 pretraining steps for L1 and L2 , and 10000 joint steps for both methods; fold-wise metrics are averaged within each repetition and then summarized across the 10 split seeds. In the pooled joint 10-fold protocol of Tables 5 and 6, we use 200, 200, and 100 observations at L1 , L2 , and H, respectively, with held-out observations across all fidelities, 3000 pretraining steps for L1 and L2 , 1000 pretraining steps for H, and the same 12000-step DAG-VI / 2000 + 10000-step DAG-SVI budgets as above. The main paper reports the free initialization for this protocol, while Table 6 gives the corresponding validation run with predmean initialization. We use the free initialization only in the pooled 10-fold protocol, where the larger training set supports this more flexible procedure. Predictive metrics. RMSE and CRPS are reported with their standard definitions, see Ji et al. (2024)); it is worth noting that, following Ji et al. (2024), the normalized RMSE (N-RMSE) is defined as 1 − RMSE/RMSEbase , where RMSEbase is the RMSE of a constant sample-mean baseline predictor, so that larger values indicate better performance. In the published-split experiment, the baseline mean is computed from the available training targets in that protocol; in the cross-validation protocols, it is recomputed from the corresponding training fold. Joint predictive metrics. Let F denote the joint posterior predictive distribution of the concatenated held-out vector yval = {L1 (xi )}i∈IL1 , {L2 (xi )}i∈IL2 , {H(xi )}i∈IH . We evaluate the multivariate weighted energy score 1 ESW (F, yval ) = E ∥X − yval ∥W − E ∥X − X ′ ∥W , 2

iid

X, X ′ ∼ F,

with ∥z∥W =

1 1 1 ∥zL1 ∥22 + ∥zL2 ∥22 + ∥zH ∥22 2 2 2 nL 1 σL n σ n σ L H 2 L2 ,tr H,tr 1 ,tr

!1/2 ,

where nL1 , nL2 , nH are the held-out block sizes in the current fold, and σL1 ,tr , σL2 ,tr , σH,tr are the empirical standard deviations computed from the observed training data of that fold only. Results on the other two experiments.

Model

RMSE (×10−2 ) ↓

N-RMSE ↑

CRPS (×10−2 ) ↓

2.44 ± 0.30 2.50 ± 0.34 1.74 ± 0.26

0.754 ± 0.032 0.749 ± 0.032 0.826 ± 0.024

1.68 ± 0.29 1.48 ± 0.20 1.06 ± 0.19

d-GMGP DAG-VI DAG-SVI

Table 4: Predictive performance on the heavy-ion collision emulation task using 200 observations at each lower fidelity and 25 high-fidelity observations, with validation folds formed on the high-fidelity set only. Results are reported over 10 repeated 5-fold cross-validation runs as mean ± std. RMSE and CRPS are reported in units of ×10−2 ; lower is better for both, while higher is better for N-RMSE.

Model

Weighted ES ↓

H RMSE (×10−2 ) ↓

H CRPS (×10−2 ) ↓

H N-RMSE ↑

DAG-VI DAG-SVI

0.1624 ± 0.0352 0.1583 ± 0.0317

1.15 ± 0.58 1.06 ± 0.55

0.62 ± 0.27 0.58 ± 0.26

0.883 ± 0.051 0.891 ± 0.051

Table 5: Joint 10-fold cross-validation results on the multi-fidelity heavy-ion task using the free initialization for the high-fidelity inducing inputs in both DAG-VI and DAG-SVI. Results are reported as mean ± std across folds. Lower is better for weighted energy score, RMSE, and CRPS; higher is better for N-RMSE.

Model

Weighted ES ↓

H RMSE (×10−2 ) ↓

H CRPS (×10−2 ) ↓

H N-RMSE ↑

DAG-VI DAG-SVI

0.1635 ± 0.0363 0.1621 ± 0.0353

1.15 ± 0.56 1.10 ± 0.55

0.62 ± 0.26 0.59 ± 0.26

0.883 ± 0.050 0.888 ± 0.050

Table 6: Joint 10-fold cross-validation results on the multi-fidelity heavy-ion task using the predmean initialization for the high-fidelity inducing inputs in both DAG-VI and DAG-SVI. Results are reported as mean ± std across folds. Lower is better for weighted energy score, RMSE, and CRPS; higher is better for N-RMSE.

Pretraining effect.

H-Pretraining vs Joint H pretrain

H pretrain

Joint

H pretrain Joint training

20 Validation CRPS on H (×10 2)

Validation RMSE on H (×10 2)

0 500

25000 1k

4k Training progress

12k

500

25000 1k

4k Training progress

12k

Figure 15: Effect of high-fidelity pretraining in the pooled 10-fold heavy-ion protocol. Validation RMSE and CRPS on H are shown during local H-pretraining and subsequent joint training, aggregated over folds. The same pretraining strategy is used for both DAG-VI and DAG-SVI: nodewise pretraining phase stabilizes the high-fidelity branch before full DAG optimisation. The joint training that targets the DAG-DGP variational objective is fundamental to achieve better and stable performance in the task.

Record · ID 361434 · SHA-256 02484f92ddc8c8ce
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.