A Bayesian Take on Gaussian Process Networks
Authors: Enrico Giudice
Jack Kuipers
Giusi Moffa
Venue: NeurIPS 2023 poster (2023)
Keywords: Bayesian networks
structure learning
graphical models
gaussian processes
Bayesian inference
MCMC sampling
importance sampling

Abstract: Gaussian Process Networks (GPNs) are a class of directed graphical models which employ Gaussian processes as priors for the conditional expectation of each variable given its parents in the network. The model allows the description of continuous joint distributions in a compact but flexible manner with minimal parametric assumptions on the dependencies between variables. Bayesian structure learning of GPNs requires computing the posterior over graphs of the network and is computationally infeasible even in low dimensions. This work implements Monte Carlo and Markov Chain Monte Carlo methods to sample from the posterior distribution of network structures. As such, the approach follows the Bayesian paradigm, comparing models via their marginal likelihood and computing the posterior probability of the GPN features. Simulation studies show that our method outperforms state-of-the-art algorithms in recovering the graphical structure of the network and provides an accurate approximation of its posterior distribution.

A Bayesian Take on Gaussian Process Networks
Enrico Giudice
Jack Kuipers
Dep. Mathematics and Computer Science University of Basel, Basel, Switzerland enrico. waaicc.gguiiin giudiceuunibas
Dep. of Biosystems Science and Engineering
ETH Zurich, Basel, Switzerland jack. kuipers@bsse. ethz. ch
Giusi Moffa
Dep. of Mathematics and Computer Science, University of Basel, Basel, Switzerland and Division of Psychiatry, University College London, London, UK giusi.moff www.mmffffmmbaa.cccmmmmm
Abstract
Gaussian Process Networks (GPNs) are a class of directed graphical models which employ Gaussian processes as priors for the conditional expectation of each variable given its parents the network. The model allows the description of continuous joint distributions in compact but flexible manner with minimal parametric assumptions on the dependencies between variables. Bayesian structure learning of GPNs requires computing the posterior over graphs of the network and is computationally infeasible even in low dimensions. This work implements Monte Carlo and Markov Chain Monte Carlo methods to sample from the posterior distribution network structures. As such, the approach follows the Bayesian paradigm, comparing models via their marginal likelihood and computing the posterior probability of the GPN features. Simulation studies show that our method outperforms state-of-tt the-rrtertt algorithms in recovering the graphical structure of the network and provides an accurate approximation of its posterior distribution.
1 Introduction
Bayesian networks (BNs) are a powerful tool for compactly representing joint distributions and the underlying relationships among large set of variables 33 These relationships are described via directed acyclic graph (DAG), with each node in the graph representing a random variable. The joint distribution random variables X (X1,... X,) factorizes into conditional distributions for each variable given its parents in the DAG:
p(X) II p(X, Pax.). 1=1
The DAG provides visual description of the dependency structure among the variables, where missing edges encode conditional independence relations among pairs of variables. Thanks their inherent flexibility and their ability to combine expert knowledge with data, Bayesian networks have been applied to large range of domains (2).
In the absence of full prior knowledge of the underlying graph, inference on the structure is necessary. For this purpose, a plethora of structure learning algorithms have been proposed, which involve either learning the structure from the conditional independence relations in the data (i.e. constraint-hased) or assigning score each DAG and searching the graph space for high-scoring networks. recent overview of current algorithms appears in Kitson et al. X3 while reproducible workflow for
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

benchmark studies is available in Rios et al. [37. Hybrid methods combining both constraint- and score-based aspects generally offer current state-of-the thetthhhtttttt performance.
Typical scoring functions rely on measures of goodness of fit of the graphical model for the given data, such the Bayesian Information Criterion (BIC) and the Mutual Information Test (MIT) Another common choice is to use the graph marginal likelihood that we can derive by integrating over the parameter prior on a graphical model. The vast majority of score-based methods implement searches aiming to maximise the scoring function and return one single optimal structure. By complementing the marginal likelihood with graph prior, one can score graphs according to their posterior distribution and sample from it. By fully characterizing the posterior distribution over graphs possible to go beyond optimal point estimates and enable Bayesian model averaging. In this context, approaches have emerged to capture both network and parameter uncertainties across a variety of research areas 27
24 B
Likewise, we focus on Bayesian methods which sample DAGs according to their posterior distribution. Unlike approaches based on maximizing the score function, which return single graph estimate, Bayesian approaches provide a full posterior distribution over DAGs, which we can use to perform inference on the network's features of interest. However, Bayesian structure inference typically requires a parametric model for the conditional distributions of each variable given its parents in the network compute the posterior. Variations of the BDe score 2 are the usual choice for discrete-variable discrete- variable networks. Other information theory-based score functions All B] are not Bayesian
since they do not estimate marginal likelihood but rather provide a measure of fitness of a DAG to the data based on the minimum description length principle.
In the case of continuous-variablee variable networks, most of the current research has focused the linear- Gaussian case, due to the availability of the closed-form BGe score 12 26 When relaxing the assumption, current approaches fall outside of the Bayesian framework since they do not target a posterior probability for DAGs but rather search the DAG space for a high-scoring network by minimizing a generic loss function. For example, Elidan 66 employs rank correlation as a goodness-ot of fit measure between a graph and the observed data; Sharma and van Beek 4Q propose modeling the relations among variables via regression splines and scoring the networks via cross- validated score. Using the formulation by Zheng et al. [52] of the structure learning problem as a constrained continuous optimization, Yu al. 55 model the conditional distributions via graph neural network.
In this work, we develop sampling scheme perform fully Bayesian structure inference on generic continuous-variable networks with potentially non-linear relations among variables. We follow the strategy of Friedman and Nachman [9] and model the functional dependencies between each variable and its parents via Gaussian process priors. We extend their approach to the Bayesian framework by making structural inference based on the posterior distribution, which also involves treating the hyperparameters of the model as random. After a short review in Section 2 of Bayesian structure inference, Gaussian processes and how they can be used to parameterize BNs, Section proposes proposes a sampling scheme to perform fully Bayesian inference on nonlinear, continuous networks. In Sections we evaluate and compare our algorithm existing approaches on simulated and real data.
and 5
2
Background
2.1
Bayesian Structure Inference
When modeling set of variables X with Bayesian network, a typical objective is to estimate the probability of a generic feature of interest V. Examples of such a feature are the presence and direction of certain edges, the topological ordering of nodes conditional independence relations among variables. Given a complete set D of observations of X, we can obtain the posterior probability of vp by integrating the probability or presence of the feature over the posterior distribution of graphs:
p D) = Σ (()
D( x p(D 9) PG))
(2)
wherein key component is the likelihood of the data integrated over the prior distribution of the parameters 0 of the conditional distributions:
p(D ( = -f p(DIg, p00 d8.
(3)
2

---

The marginal likelihood p(D () for given DAG & is available in closed form for specific combina- tions of priors on parameters and likelihood functions
The number of DAGs grows super--xponentiall in the number of nodes 38) hence exact posterior inference is still exponentially complex, making it effectively intractable for large networks [47]. Approximate inference for p(V D) is possible if one can obtain samples GI,..., 9M from the posterior distribution over graphs p(g D).
Markov Chain Monte Carlo (MCMC) methods have successfully tackled the problem of posterior inference in the large and complex space of DAGs. The rationale consists of defining Markov chain whose stationary distribution the posterior distribution of interest p(g D). Madigan et al. [0) suggested Metropolis-Hastiee accepts or rejects it according to probability a. The value of a depends on the relative posterior D) of the proposal graph compared to that of the current DAG:
sampler: at each step, the algorithm proposes
DAG G' and it
probability p(g')
min (1, p(g'l D) Q() P(g D) Q(G'
(4)
where Q(:-)) are the transition probabilities from one DAG to another. In its original formulation, the algorithm proposes new graph by adding or deleting one edge from the current DAG, with Q uniformly distributed in the set of neighbouring DAGs (including the current one).
Several refinements to the original method have been made over the years, providing for example The order MCMC approach by Friedman and Koller S introduced an important variation by sampling in the space of node orderings, which is smaller and smoother than the space of DAGs. Because of these attractive properties of the sample space, the order MCMC algorithm achieves superior mixing and convergence results compared to regular structure MCMC. Sampling orders however introduces a bias in the resulting posterior since node orders do not induce uniform coverage in the DAG space 77. The partition MCMC algorithm (25] corrects this bias in order MCMC by sampling in the larger space of ordered partitions of nodes to achieve unbiased samples from the posterior.
more efficient proposal distributions 14 20
Current fate-of-theeart methods rely on conditional independence tests to obtain a first coarse estimate of the structure to use as a starting point for score-based method 4S More sophisticated approaches involve iteratively expanding the search space to correct for errors its initial estimate 2X To achieve efficient posterior inference
[4 28 in the selected search space it helps to precompute all combinations of scores potentially needed in the chain so it can run at a very low computational cost.
Key to enabling efficient DAG sampling is decomposability, property ensuring that we can express the posterior probability of a DAG G as a product of local scores that depend only on each variable and its parents in G:
p D) II S(x,, D).
(5)
Decomposability guarantees that we need recompute only those scores whose parent sets have changed, or we can precompute all needed combinations of scores efficiently (28).
2.2 Gaussian Processes
Gaussian processes are flexible tool used in machine learning for regression and classification tasks. Formally, a Gaussian Process (GP) is a distribution over functions such that every finite collection of its function values {f(x1) f(22), .... f (2+)} has multivariate Gaussian distribution 36. A GP therefore fully specified by its mean function (xx) = [[f(xx) and covariance function (=), f(='). ]. Due to their ability to model large range of functional behaviours, GPs find common use as priors over regression functions f(X) = E(Y X). A common GP regression model assumes independent Gaussian additive noise:
k(x, x²)
Y = fX) + f(X) ~ GP(O,K(z,*')) k(x,,')). x²)), E ~~ Niu,
(6)
Notably, GP models admit closed-form expression of the marginal likelihood, in this case, the likelihood of the N observations y marginalised over the prior distribution of
(u (K + - (2-11) m1
p(y) (2T)
K + o²
exp
where K is the N x N Gram matrix Ki3 k(I1,1z).
3

---

2.3 Gaussian Process Networks
Gaussian process networks (GPNs) refer to Bayesian networks whose conditional distributions are modeled via Gaussian process priors 19]. The structural equation model defining the distribution of each variable X; given its parents in GPN
X; = (Pax,) E
(8)
where has normal distribution (u,, (2) independent of the data and a Gaussian process prior placed on the function f.. Thanks the nonparametric nature of GPs, the model { can capture a wide range of functional dependencies between variables while maintaining the closed-form expression for the marginal likelihood.
The covariance function & of the Gaussian process is typically parameterized by set of hyperpa- rameters 0 which determine the behaviour of samples of Sa such smoothness, shape periodicity, When these hyperparameters (together with the noise mean and variance σ²) are unknown they are typically learned by maximizing the marginal likelihood. Since the marginal likelihood and its derivatives are available in closed form, the optimization can be performed efficiently via gradient ascent 36 The GPN model enables structure learning of networks on continuous variables without the need to make strict parametric assumptions on the distributions of the variables. In scenarios even with low sample sizes, Friedman and Nachman [ have shown that searching for the highest- scoring structure can accurately reconstruct the underlying DAG under different functional relationships.
When estimating the score of GPN, the common approach learning the hyperparameter' by simple maximization can however lead to problematic estimates, since many local optima may exist 4. In addition, the resulting plug-in estimate of the marginal likelihood would not correctly reflect the uncertainty in the data of the full posterior. Bayesian model averaging provides a natural solution by integrating over the prior distribution of the hyperparameters obtain a true marginal likelihood p( D 9), which we can then use to perform posterior inference on the structure.
3 Bayesian Structure Inference for GPNs
In this section, we describe a method to sample from the posterior distribution v( D) of generic feature y of a GPN for continuous, non-linear data. To implement a fully Bayesian approach, we place priors over the hyperparameters 0 of the kernel function and the Gaussian noise.
= f(Pax) tG, E~ ~~((,,,0)
f GP(O, ko(...) 0 (0).
(9)
1.0 ~ T(H,o).
The priors ensure that the uncertainty in the functional relationships between each variable and its parents fully accounted for. On the other hand, a maximum likelihood approach to estimating the hyperparameters could yield overly confident score functions and in turn misrepresent the posterior p(g) D). Under the GPN model, the score remains decomposable, necessary condition for its efficient evaluation. As in Section 2. making inference about network features T of interest is possible by sampling graphs from the posterior:
1 M
p V
D)
Σ (V))))
9j p(9, D).
(10)
M
Sampling graphs hinges on computing the scores of all variables given different combinations of their possible parent sets (see equation 5
Let θ = {1, 0} be the d-dimensional set of hyperparameters for a given node X and its parents Pax. Unless stated otherwise, throughout the rest of the text we assume a uniform prior over all structures p(G) x 1. The score function is then the likelihood 7 of the observations x marginalized with respect to the hyperparameter priors:
5(X, Pax) = -f |
Pax, 0) () T( Pax) dO.
(11)
variable X has no parents then the Gram matrix of the kernel is zero and the score function reduces to Gaussian marginal likelihood. Since the above score function is generally intractable,

---

Algorithm 1 GP network sampling scheme
Input Data D of variables, feature of interest y Output Posterior probability of the feature p( D)
1: for j E {l.... ,...... ) do
2: Sample DAG Gj according its Laplace approximate posterior ((9) D).
Eq.
3: for i E [2, do
4:
Compute SX,,, P));) via MC estimation.
Eq. 12
5: Compute posterior p(9) D).
> Eq. 5
6: Compute posterior probability of via importance sampling.
Eq.
(4
one option to use Monte Carlo (MC) approaches such as bridge sampling G to approximate it. Bridge sampling employs Gaussian proposal distribution and bridge function h chosen to
minimize the MSE of the resulting estimator. The bridge sampling estimator of is then defined as
ΣN p(a Pax, O2) x (O; Pax) ((O; Pax)
S(X, 2
M
(12)
N2 1 Pax) h()) Pax)
θ ~ ~ Pax),
0 ~ p Pax). i= 1, N1, 1=1,..,,N2
The estimator is also function of samples from the posterior of the hyperparameters p(O) x, Pax). which can easily be obtained via MCMC sampling. One can show that other approaches such as importance sampling or harmonic mean are special cases of bridge sampling [[8]. MC methods based on bridge sampling provide consistent estimators for the marginal likelihood but may be biased for finite samples [50]. Nonetheless, such methods can become computationally expensive in high dimensions (i.c. for large parent sets).
Since sampling DAGs via MCMC requires computing large number of scores of potential parent sets, which may not all be represented in the final sample, we avoid computing these expendable scores by first running the MCMC algorithm using Laplace approximation of the score [I1 around the MAP value of O
(2m)
Sr(X, p(x Pax, θ) ( Pax) H/2
(13)
with é argmax p(a: Pax. θ) T(O) Pax). and Hij 8pp(1 Pax, +) ( Pax)
00,00,
@=
We denote the resulting posterior probability of DAG G from this Laplace approximated score as q(G D) to distinguish it from the true posterior p(g D).
The Laplace approximate score provides an approximation of the posterior at lower computational cost, speeding up considerably the running time of the MCMC algorithm used to sample graphs. After sampling M graphs in the first step with the Laplace approximate score, we can make inference with respect to the true posterior by re-computing the scores and performing importance sampling. To estimate the posterior probability of a feature of interest V via importance sampling we evaluate
P(V D) ~ ΣØ 1 wj
E3=1"} , wy = p(9, D)
(14)
<(9, D)
where p(G) D) and (( D) are the posterior probabilities of DAG Gj computed respectively with the bridge sampling MC estimate (2) and the Laplace approximation (33) for the score. The two posteriors do not need to be normalized for valid posterior inference V. This is because the normalizing constants in the importance sampling weights (4) cancel out. The procedure is
summarized as pseudo-code in Algorithm 1.
Re-computing the score functions in a second step and implementing importance sampling compu- tationally advantageous compared to running the MCMC algorithm directly with the MC estimates of the scores. The advantage is simply due to the number unique scores in the final chain being much lower than those evaluated or needed during the chain itself (see figure A.3 in the supplementary material). Regardless the MCMC algorithm used, we expect a substantial improvement run-time compared to using the MC scores directly.
5

---

3.1 Implementation Details
Bayesian inference and optimization of the hyperparameters were performed via the Stan interface RStan 43.. The library offers highly efficient C++ implementation of the No U-turn sampler [22]. providing ofat--of-the-rr posterior inference on the hyperparameters. We performed MC estimation of the marginal likelihood via the bridgesampling package [19], which can easily be combined with fitted Stan models.
In the implementation of Algorithm 1 we use order or partition MCMC to generate the network samples ((; D); the procedure can however accommodate a variety of sampling methods, long as they result in samples from the posterior. Given its good performance in benchmarking studies 45 hybrid implementation for MCMC inference on the graph. The
[37, we use the recent BiDAG
hybrid sampler requires an initial search space which is then improved upon to correct for possible estimation errors 28 As initial search space we take the output of the dual PC algorithm [3] For
the bridge sampling estimator (12) of the marginal likelihood we used N1 N2 = 300 particles
from the proposal and posterior distribution over the hyperparameters. The proposal function g was set to a normal distribution, with its first two moments chosen to match those of the posterior distribution. R code to implement Algorithm and reproduce the results in Sections and is available at htttps://githhn hwwpttttu...omeeeiieeeeeee http///rrroogo.dddc////aaaaaaaaa
3.2 Score Equivalence
Without any assumptions the functional form of structural models, observational data cannot generally distinguish between any two DAGs within the same Markov equivalence class [34], since the joint distribution always factorizes according to either DAG. Scoring functions like the BGe that assign the same score to all DAGs in the same class satisfy score equivalence. Imposing parametric models on the local probability distributions of Bayesian network may however break score equivalence. Specifically for GPNs, an alternative factorization may not admit a representation that follows the structural equation model in B Consequently, DAGs belonging the same Markov equivalence class may display different scores and become identifiable beyond their equivalence class. For this reason, the GP score generally differentiates models according the direction of any of their edges.
Importantly, the identifiability of otherwise score equivalent DAGs is direct consequence of the assumptions underpinning the GPN model. Because the functions fo are not generally invertible, the likelihood in 7 will assign higher values to functional dependencies that admit the GPN structural equation model. Further, Peters et al. 35 demonstrated that the asymmetry distinguishing between Markov equivalent factorizations holds beyond the case of non-invertible functions. Indeed, with the exception of the linear case, all structural models that take the form in ( with additive Gaussian noise violate score equivalence and may identify DAG beyond its Markov equivalence class. GPNs may encompass cases where equivalent DAGs are indistinguishable, such as when the functions Si are linear. In this case, the computed scores of equivalent DAGs will be similar, as long as the GP prior allows learning sufficiently linear relations. The numerical experiments in supplementary Section show that the Ga-dd score function displays near score equivalent behaviour when the data come from a joint Gaussian distribution, coherently with theoretical considerations (35
Experimental Results
We evaluate the Bayesian GP network inference scheme on data generated from known random networks with n = 10 nodes. The DAG structures are constructed from an Erdós- Kénis model, where each node has an independent probability of 0.2 to be connected with another with higher topological ordering. For every node X, in each randomly generated network, we then sample 100 observations non-linear function of its parents. The nonlinear data generated by transforming the parents' instances using a weighted combination of six Fourier components
X, Σ
B
Wi.j. X3
Σ
D.j. sin (k:X)) + W.g.. cos (1xx)) eg.
(15)
JJX, Eaax
The weights Ui.jA and Wi.j. are sampled from Dirichlet distribution with concentration parameters k/2
equal to seven positive exponentially decreasing values 7/k = E.e 67x, for k = {0, (0,...66}. 6}. The

---

parameter À controls the rate of exponential decay, with values of À close to zero providing mostly linear effects between variables and higher values resulting increasingly non-linear relationships.
The edge coefficients Bi, determine the strength of the dependencies between X; and its parents and are sampled from a uniform distribution on (-2, - U (22 2); the noise variable e; has a standard normal distribution. Instances of root nodes (nodes without parents) are also sampled from standard normal. The linear-Gaussian case corresponds asymptotically to À = 0: in this case, we set the weights jij,k and Wi,j. zero except for wi.j.--
We compare all structure learning algorithms in terms of structural Hamming distance (SHD)
48
which compares estimated graphs G with the ground truth graph G*. Following Lorch et al. 29 compute the average SHD of the samples weighted by the posterior:
E- E-SHD(p. 9*)
Σ SHD(G, P(G)D)
(16)
The E-SHD summarizes how close the estimated posterior is to the true DAG; lower values therefore desirable. is however not direct measure of the similarity the estimated posterior to the true posterior.
4.1 Choice of Priors
Different choices kernels for the GP prior result in different behaviours of the conditional expecta- tion in equation 99) of variable X given its parents Pax. A simple model employs additive kernel
k(...) Σ Ko. (...)
(17)
which corresponds to modeling each variable X as sum of the individual contributions of each of its parents. The additive model serves to reduce the computational burden calculating the scores by keeping the number of parameters small possible while preserving non-linearity in the model. The approach can however easily accommodate more complex relationships at a higher computational cost, such an additive kernel with all first-order interactions
[5]:
Paxl k(.,.) = n Σ ko, (. + 72
Paax Paxl
Σ
Σ
ko, (..)) (. Roy (...)
(18)
For each parent Z = (Pax); we used a squared exponential kernel function Kg, (z, 2') =
||===111") with each O; measuring the degree of non-linearity along the Z-th dimension. The kernel function has unit variance since the data are always normalized in the structure learn- ing process. We assign the following independent prior distributions to the hyperparameter set θ (n.0, A :i e. l,... [axl]
exp
20
| of the GPN model P
0, 1G(,,2),
NNO,, 1),
IG.
(19)
The overse-gamma priors for each lengthscale 0; and the noise variance suppress values near zero, which in either case would result in overfitting and degenerate behaviour of samples of [13]. The additional parameters and 72 of the kernel function (88 were assigned independent IG(1, 1) priors.
4.2 Results
We compare different versions of the GP-based structure learning scheme: the first employs the order MCMC algorithm for sampling DAGs, and the second uses partition MCMC. We then compare both versions of the GPN sampling scheme with order and partition MCMC with the BGe score 28 To reduce the computational burden of the simulations, we employ the simpler additive kernel for the GP score.
As an additional benchmark, we include the DiBS+ algorithm 22 which models the adjacency matrix probabilistically. using particle variational inference to approximate posterior over structures. In the simulations, we parameterize DiBS+ by neural network with one hidden layer with nodes. We also consider methods which rely on the constraint-basee PC algorithm 42], which
7

---

2=1
=.5
==
15
15
15
E S-SHD
10
E-SHD
10
E-SHD
10
GP, partiion
GP: ordar
B0e, partiton
XPC-HSIC
KPC-DC
DIBS.
GP, parmon
BGc, partiion
B0e onder
IPC-HSIC
APC-DC
DIBS-
GP partiton
GP onder
80e, partiion
BGe order
IPC-HBIC
MPC-DO
DIBS.
GP. partitio
GP, onder
BGo, parition
BGe, order
XPO-HBIG
APC-00
DBS+
Figure 1: Distribution of E SHD values for all the different algorithms. 0 corresponds to
linear-Gaue Gaussian data while higher values increase the degree of non non-linearity of the relations among variables.
learns network by testing for conditional independence among the variables. Since these not Bayesian methods, we bootstrap the data to obtain different graph estimates via the PC algorithm
10]. To account for the non-linear dependencies between the variables apply the kernel PC (kPC) algorithm 17 to the resampled data using two different independence tests: HSIC [66] and distance correlation 20
Figure D shows, for three values A. the distribution of E-SHD values for the different algorithms over the 100 generated structures. For non-linear data (A 1), both versions of the GP samplers outperform existing methods: the median E-SHD for either our algorithms is equal or lower than the bottom quartile of the BGe samplers. For linear (A 0) and slightly non-linear data (A 0.5), the BGe score-based sampling algorithms. All EHD values computed with respect to DAGs; the results comparing completed partially directed acyclic graphs (CPDAGs) are available in figure A.3 in the supplementary material.
GPN sampler performs competitively with the ofate-of-the
Figure displays the run-times for the same methods as in figure I The run-times are stable across the different levels of non-linearity parameterized by X, with the k-PC algorithms being the most computationally expensive, followed by the GP score-based approach, DiBS+ and finally the BGe score-based MCMC methods which are the fastest.
1
0
1000
1000
1000
Time (8)
Time (s)
Time
GP pariiion
GP. order
BDa, pertiion
IPC HBI-
IPC-DC
DBS+
GP. order
partiion
BGa order
APC-HSIC
DES
BP partiion
GP. order
BCs perliion
8Ge order
IPc. HSIC
APC-DC
DBS
IPC.
BGe
partition
APC-HBIC KPC-0O
Figure 2: Distribution of run-times for all the different algorithms. À = O corresponds to linear- Gaussian data while higher values increase the degree of non-linearity of the relations among variables.

---

Proteetor protablltty
82-
Reverse K-L divergence
0.2
Ra-08
Postertor probability
0.1
0e.00
400
0.0
1000
3000
10000
100
Samples
DAG
BGe
Laplace
weighted
Figure 3: Left: reverse K-L divergence between the true posterior and the BGe posterior (green), the Laplace approximate posterior (blue) and the posterior obtained via importance sampling (red) as function of the number of sampled DAGs. Right: the true posterior (gray) together with the BGe posterior (green), the Laplace approximate posterior (blue) and the posterior obtained via importance sampling (red). The majority DAGs have a very low true posterior probability and are therefore never sampled by the MCMC algorithms (see inset).
Besides estimating accurately the structure, our approach can quantify the uncertainty in its estimates via the sampled graphs from the posterior distribution over GPNs. To evaluate the accuracy of the sampling approach in estimating the general posterior over structures, we compare its estimated posterior distribution over DAGs with the "true" posterior, obtained by enumerating every possible structure and computing its score directly with equation (12). Due to the exceedingly large number of DAGs, this approach is only feasible with small structures. The left panel of figure 3 shows the reverse Kullback-Leible (K-L) divergence between the estimated and true posteriors for given network with n 4 nodes, as function of the number samples M in the partition MCMC algorithm. The estimated ("GP, weighted")) posterior probability p(g D) for generic DAG G is obtained by setting p( 92) [g) in equation (4). Reverse K-L divergence was chosen as a metric since the algorithms assign probability of 7ero to DAGs that were not sampled.
The plot includes the divergence between the Laplace approximate posterior q( ()) in equation (14) and the true posterior, as well as between the posterior obtained with the BGe score and the true posterior. The divergence of the Laplace approximation is reduced by roughly one order of magnitude by weighting the samples via importance sampling. For reference, the DiBS+ algorithm yields reverse K-L divergence of 27.8 with 1000 samples, i.e. two orders of magnitude higher than our approach, despite also being allocated longer run-time unntimmes (see figure A..4in the supplementary material), Sampling a higher number of graphs with DiBS+ quickly becomes infeasible since its run-time scales quadratically with the number of samples 2 while MCMC sampling scales linearly with M. The right panel of figure shows the posterior distributions over the 543 possible DAGs, for M 10 obtained after complete enumeration and by sampling with different methods. The plots confirm that as the number of sampled DAGs increases, our approach can accurately reflect the full posterior uncertainty. The results also underline the importance of sampling from the hyperparameter'' priors to obtain an accurate representation of the posterior, as the Laplace approximation of the marginal likelihood results in highly biased approximation even when sampling large number of DAGs.
5 Application on Protein Signaling Networks
We also applied the GP score-based structure inference to the flow cytometry dataset of Sachs et al. (39 to learn protein signaling pathways. The authors provide an accepted consensus network, which we used as reference. We considered the first experimental condition, consisting of 853 single-cell observations of 11 phosphoproteins and phospholipids in human T cells. The first three columns of table show the performances of all the different algorithms in reconstructing the consensus network. The results include the GP model using the additive kernel (18) with all first-order interactions, denoted as GP2. We measure the different algorithms in terms of the E SHD,

---

Table I: Performance of the different algorithms in reconstructing the consensus network from flow cytometry data. The last two columns show the posterior probabilities of the two features experimentally validated by Sachs et al. 391 where the edge on the left should be present and the one on the right absent; up/down arrows indicate higher/lower values are better.
E SHD 4
E
TP
E-FP 4
Erk - Akt t
Erk A PKA 1
GP, partition
14.5
7.2
4.7
0.75
GP, order
14.6
6.8
4.4
1
0.42
GPå partition
16.7
6.6
6.3
I
0
GP order
16.7
6.7
6.4
0.34
BGe, partition
15.9
4.3
3.2
0
BGe, order
15.3
4.4
2.7
0.17
0.98
kPC-HSIC
17
5.5
5.5
0.69
0.26
kPC-DC
16.9
6
5.8
0.72
0.28
DiBS+
16.1
4.2
0.45
0.7
as well as the EPP and E-FP, the absolute number of TP and FP edges the DAGs weighted by the posterior, obtained by replacing the SHD in equation T6 by TP or FP, respectively.
One of the benefits Bayesian structure inference is the possibility of deriving posterior probabilities of specific edges of interest in the network. In their work, Sachs et al. (39 experimentally tested two relationships between the proteins by intervening directly on the cells. By means of small interfering RNA inhibition, they concluded that the inhibition of Ekt has direct effect on Akt, while there was a weaker and non-significant effect on PKA. We therefore expect an edge Erk - Akt but no directed edge (or path) from Erk to PKA in the learned networks. The last two columns of table display the posterior probabilities of these two features according to the different algorithms. GP- and BGe-based methods perform best, with the GP learner correctly assigning the highest probability to the edge Erk - Akt, while the lack of the edge Erk to PKA predicted with a lower probability. The sparser BGe score-based methods assign lower probability the first edge, but correctly predict the absence of the second edge.
Conclusions
In this work, we have proposed procedure that efficiently performs Bayesian structural inference on Gaussian process networks and estimates the posterior of generic features of interest. Building on the original GPN idea by Friedman and Nachman we now embed it in a fully Bayesian framework. In particular, our approach involves placing priors on the hyperparameters of the score function and sampling DAGs via MCMC. Although more computationally expensive than greedy search over the DAG space for high-scoring network, the Bayesian approach allows one to accurately quantify uncertainty features of the network. This is made feasible by minimizing the number of scores to compute in MCMC chain via importance sampling. as well as harnessing the advances in MC that have been made since the introduction of GP networks.
and MCMC methods 19
The flexible nature of GPs allows the corresponding BNs to model a large range of functional dependencies between continuous variables. In the case of linear dependencies, our method remains competitive with ifate-of-the-art methods and exhibits desirable properties such as (approximate) score equivalence. The versatility of the model is particularly useful domains where the nature of relations among variables in unknown and strict parametric assumptions on the distributions of variables want to be avoided. Based on promising simulation results and the convenient properties of the proposed method, we believe that it holds potential for making accurate inference on the underlying structure BNs in complex domains.
Acknowledgements
The authors are grateful to acknowledge partial funding support for this work from the two Cantons of Basel through project grant PMB-02-- 02-18 granted by the ETH Zurich.
10

---

References
[1] Federico Castelletti and Guido Consonni. Discovering causal structures in Bayesian Gaussian directed acyclic graph models. Journal of the Royal Statistical Society: Series A (Statistics in Society), 183:1727-1745, 2020.
[2] Rónán Daly, Qiang Shen, and Stuart Aitken. Learning Bayesian networks: approaches and issues. The Knowledge Engineering Review, 26:99-157, 2011.
[3] Luis M. de Campos. A scoring function for learning Bayesian networks based mutual information and conditional independence tests. Journal of Machine Learning Research, 2 2149-2187, 2006.
[4] David Duvenaud, James Lloyd, Roger Grosse, Joshua Tenenbaum, and Ghahramani Zoubin. Structure discovery in nonparametric regression through compositional kernel search. In Proceedings of the 30th International Conference on Machine Learning, volume 28, pages 1166-115 1174, 2013.
[5] David Duvenaud, Hannes Nickisch, and Carl Rasmussen. Additive Gaussian processes, In Advances in Neural Information Processing Systems, volume 24, pages 226-234, 2011.
[6] Gal Elidan. Lightning structure learning of nonlinear continuous networks. In Proceedings of the Fifteenth International Conference on Artificial Intelligence and Statistics, volume 22, pages 355-363 2012.
[7] Byron Ellis and Wing Hung Wong. Learning causal Bayesian network structures from experi- mental data. Journal of the American Statistical Association, 103:778-789, 2008.
[8] Nir Friedman and Daphne Koller, Being Bayesian about network structure: a Bayesian approach to structure discovery in Bayesian networks. Machine Learning, 50:95-125, 2003.
[9] Nir Friedman and Iftach Nachman. Gaussian process networks. In Proceedings of the Sixteenth Conference on Uncertainty in Artificial Intelligence, pages 211-219, 2000.
[10] Nir Friedman, Moises Goldszmidt, and Abraham Wyner. Data analysis with Bayesian networks: a bootstrap approach. In Proceedings of the Fifteenth Conference on Uncertainty in Artificial Intelligence, pages 196-205, 2013.
[11] Dan Geiger and David Heckerman. Learning Gaussian networks. In Proceedings of the Tenth Conference Uncertainty in Artificial Intelligence, pages 235-243, 1994.
[12] Dan Geiger and David Heckerman. Parameter priors directed acyclic graphical models and the characterization of several probability distributions. The Annals of Statistics, 30:1412-1440, 2002.
[13] Enrico Giudice, Jack Kuipers. and Giusi Moffa. The dual PC algorithm and the role of Gaussianity for structure learning of Bayesian networks. International Journal of Approximate Reasoning, 161:108975, 2023.
[14] Paolo Giudici and Robert Castelo. Improving Markov chain Monte Carlo model search for data mining. Machine Learning, 50:127-158, 2003.
[15] Robert B. Gramacy and Herbert K. H. Lee. Bayesian treed Gaussian process models with an application to computer modeling. Journal of the American Statistical Association, 103: 11 1111-130. 2008.
[16] Arthur Gretton, Kenji Fukumizu, Choon Teo, Le Song, Bernhard Schölkopf, and Alex Smola. A kernel statistical test of independence. In Advances in Neural Information Processing Systems, volume 20, pages 585-592, 2007.
[17] Arthur Gretton, Peter Spirtes, and Robert Tillman. Nonlinear directed acyelic structure learning with weakly additive noise models. Advances Neural Information Processing Systems, volume 22, pages 1847-1855, 2009.

---

[18] Quentin F. Gronau, Alexandra Sarafoglou, Dora Matzke, Alexander Ly, Udo Boehm, Maarten Marsman, David S. Leslie, Jonathan J. Forster, Eric-Jan Wagenmakers, and Helen Steingroever. A tutorial on bridge sampling. Journal of Mathematical Psychology, 81:80-97, 2017.
[19] Quentin F. Gronau, Henrik Singmann, and Eric-Jan Wagenmakers. Bridgesampling: An package for estimating normalizing constants. Journal of Statistical Software, 92:1-29, 2020.
[20] Marco Grzegorezyk and Dirk Husmeier. Improving the structure MCMC sampler for Bayesian networks by introducing
edge reversal move. Machine Learning, 71:265-305, 2008
[21] David Heckerman and Dan Geiger. Learning Bayesian networks: a unification for discrete and Gaussian domains. Proceedings of the Eleventh Conference Uncertainty in Artificial Intelligence, pages 274-284, 1995.
[22] Matthew D. Homan and Andrew Gelman. The No-U-tum sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15:1593-1623, 2014.
[23] Neville Kitson, Anthony Constantinou, Guo Zhigao, Yang Liu, and Kiattikun Chobtham. A survey of Bayesian network structure learning. Artificial Intelligence Review, 56:1 2023.
[24] Daphne Koller and Nir Friedman. Probabilistic graphical models: principles and techniques adaptive computation and machine learning. The MIT Press, 2009.
[25] Jack Kuipers and Giusi Moffa. Partition MCMC for inference acyclic digraphs. Journal of the American Statistical Association, 112:282-299, 2017.
[26] Jack Kuipers, Giusi Moffa, and David Heckerman. Addendum the scoring of Gaussian directed acyclic graphical models. The Annals of Statistics, 42:1689-1691 2014.
[27] Jack Kuipers, Thomas Thurneer, Giusi Moffa, Polina Suter, Jonas Behr, Ryan Goosen, Gerhard Christofori, and Niko Beerenwinkel. Mutational interactions define novel cancer subgroups. Nature Communications, 9:4353, 2018.
[28] Jack Kuipers, Polina Suter, and Giusi Moffa. Efficient sampling and structure learning of Bayesian networks, Journal of Computational and Graphical Statistics, 31:639-650, 2022.
[29] Lars Lorch, Jonas Rothfuss, Bernhard Schölkopf, and Andreas Krause. DiBS: differentiable Bayesian structure learning. In Advances in Neural Information Processing Systems, volume 34, pages 241 24111-200 11-222,,, 2021.
[30] David Madigan, Jeremy York, and Denis Allard. Bayesian graphical models for discrete data. International Statistical Review, 63:215-232. 1995.
[31] Xiao-Li Meng and Wing Hung Wong. Simulating ratios of normalizing constants via a simple identity: theoretical exploration. Statistica Sinica, 6:831-860, 1996.
[32] Giusi Moffa, Jack Kuipers, Giuseppe Carrà, Cristina Crocamo, Elizabeth Kuipers, Matthias Angermeyer, Traolach Brugha, Mondher Toumi, and Paul Bebbington. Longitudinal symp- tomatic interactions in long-standing schizophrenia: a novel five-point is based on directed acyelic graphs. Psychological Medicine, pages 1-8, 2021.
[33] Judea Pearl. Probabilistic reasoning in intelligent systems: networks of plausible inference. Morgan Kaufmann Publishers Inc., 1988.
[34] Judea Pearl. Causality: Models, reasoning, and inference. Cambridge University Press, 2000.
[35] Jonas Peters, Joris Mooij, Dominik Janzing, and Bernhard Schölkopf. Causal discovery with continuous additive noise models. Journal of Machine Learning Research, 15:2009-2053 2013.
[36] Carl Edward Rasmussen, Gaussian processes machine learning. Springer Berlin Heidelberg, 2004.
12

---

[37] Felix L. Rios, Giusi Moffa, and Jack Kuipers. Benchpress: a scalable and platform- independent workflow for benchmarking structure learning algorithms for graphical models. UrXiv:2107.03866 2021.
[38] Robert W. Robinson. Counting labeled acyclic digraphs. In New Directions in Graph Theory, pages 239-273. New York: Academic Press, 1973.
[39] Karen Sachs, Omar Perez, Dana Pe'er er, Douglas A. Lauffenburger, and Garry P. Nolan. Causal protein-sinnaling networks derived from multiparameter single-cell data. Science, 308:523-529, 2005.
[40] Charupriya Sharma and Peter van Beek. Scalable Bayesian network structure learning with splines. In Proceedings of The 11th International Conference on Probabilistic Graphical Models, volume 186, pages 181-192, 2022.
[41] Tomi Silander, Teemu Roos, and Petri Myllymäki. Locally minimax optimal predictive modeling with Bayesian networks. In Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics, volume 5, pages 504-511, 2009.
[42] Peter Spirtes, Clark Glymour, and Richard Scheines. Causation, prediction, and search. Springer New York, 1993.
[43] Stan Development Team. RStan: the R interface to Stan, R package version 2.21.8, 2023. URL https::/mcc-stan.org
[44] Chengwei Su and Mark Borsuk. Improving structure MCMC for Bayesian networks through Markov blanket resampling. Journal of Machine Learning Research, 17:1-20 2016.
[45] Polina Suter, Jack Kuipers, Giusi Moffa, and Niko Beerenwinkel. Bayesian structure learning and sampling of Bayesian networks with the R package BiDAG. Journal of Statistical Software, 105:1-31, 2023.
[46] Gábor J. Székely, Maria L. Rizzo, and Nail K. Bakirov. Measuring and testing dependence by correlation of distances. The Annals of Statistics, 35:2769-2794, 2007.
[47] Topi Talvitie, Aleksis Vuoksenmaa, and Mikko Koivisto. Exact sampling of directed acyclic graphs from modular distributions. In Proceedings of The 35th Uncertainty in Artificial Intelligence Conference, pages 965-974, 2020.
[48] Ioannis Tsamardinos. Laura Brown, and Constantin Aliferis. The max-min hill-climbing Bayesian network structure learning algorithm. Machine Learning, 65:31-78, 2006.
[49] Jussi Viinikka, Antti Hyttinen, Johan Pensar, and Mikko Koivisto. Towards scalable Bayesian learning of causal DAGs. In Advances in Neural Information Processing Systems, volume 33, pages 1588465544 2020.
[50] Jackie S.T. Wong, Jonathan J. Forster, and Peter W.F. Smith. Properties of the bridge sampler with focus splitting the MCMC sample. Statistics and Computing. 30:799 816, 2020.
[51] Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG structure learning with graph neural networks. In Proceedings of the 36th International Conference on Machine Learning, volume 97, pages 7154-7163, 2019.
[52] 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 9492-9503, 2018.
13


Reviews:
[{"review_id": "Qkrau1TPMM", "reviewer": "Reviewer_JYiJ", "summary": "Within the broader task of causal network structure identification from observational data, i.e. finding the probabilistic graphical model which best explains (has highest likelihood) observations of a set of variables, this paper focus on a type of network known as GPN. GPNs are adapted to scenarios where all variables are continuous, and are built from a simple building block in which every child random variable's distribution, conditional on its parents, is a GP regression with independent Gaussian noise. For this task and this type of network, the paper proposes a fully Bayesian inference procedure based on MC and MCMC to sample PGM structure (i.e. edge presence and direction), and compute graph-related statistics of interest (such as the probability of presence of an edge between two given variables). This procedure holds the promise of more accurately representing the graph structure distribution than score-based methods.", "questions": "- To make the point lines 168 ff., would it be worth counting the number of score evaluations?\n- Do some real-world scenarios require a different choice from the uniform prior over graph structures (line 141) ?\n- Table 1, rightmost column: can the poor performance of your method be explained by the fact that there is no reason to predict edge absence, when absence of correlation can be captured with a fairly wide conditional distribution, i.e. low covariance in GP?", "limitations": "The limitations tied to algorithm choices are well discussed. It is worth motivating better why it is reasonable to focus on GPNs with respect to real-world needs, where all variable might not be continuous, where GP conditional distributions might be too loose, or where neural-network-parameterized conditional distributions might provide a more adequate model.\n\nI see no interesting limitations to discuss related to applications or societal impact.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "Originality: The paper's approach consists of \"pushing through\" with a strong MCMC sampling procedure for the stated problem. The problem statement is not original per se, but there is moderate inventivity in putting together the required technical steps to overcome each hurdle: speeding up eq 5 (graph posterior) thanks to eq 12 (Laplace approximation of per-variable score), importance sampling for feature posteriors in eq 13. Originality is not the paper's strong suit however.\nQuality: Derivations and implementations are all correct. Experimental validation on synthetic data is a necessity, especially where full enumeration of network structures is needed; the validation on the Sachs dataset is interesting, though one would wish for further real-data experiments. I appreciate that score equivalence/graph identifiability is addressed thoroughly. Evaluation methods and prior choices seem well-motivated and comparable to other literature.\nClarity: The paper is exceptionally clear. Citations are almost all accurately chosen (with the exception of a few in the introduction, cf below). \nSignificance: The paper tackles its stated problem heads-on. I believe we will see further Bayesian takes on network structure inference, maybe with smarter graph posterior approximations or even variational approaches for certain types of networks. The choice of GPN (and the particular further choices made for experimental evaluation) is rather flexible as it supports varied functional dependencies of child upon parent variables. The experimental evaluation places variants (order/partition) of the paper's methods at better or on par with BGe, and better than other methods. It is clearly superior in its sampling behaviour (fig 2). Its running time, unfortunately hidden in the supplementary material fig A8, is 30x that of BGe however, which might restrict applicability. ", "weaknesses": "Significance, as discusse above. \nMore real-world experiments would help assess the impact, in particular if they were analysed with the level of detail of synthetic experiments sec4.\n\nMinor errors and typos I found\n- line 162: perform, not make\n- line 208: data is generated\n- line 266: define CPDAG\n- line 289: increase, our approach\n- lines 24 and 29, citations seem quite cherry-picked and do not provide a variety of viewpoints as I would expect; I would hope there are no biased self-citations here.\n- supplementary line 479: do you mean TP instead of TN? eq 19, what does the lowercase p stand for in FPRp?\n- table 1: I would add the usual arrow-up, arrow-down to indicate \"higher/lower is better\" for easier reading", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their constructive comments and suggestions, and will incorporate them into the manuscript. \n* *Minor errors and typos I found [...] lines 24 and 29, citations seem quite cherry-picked and do not provide a variety of viewpoints as I would expect; I would hope there are no biased self-citations here.\nsupplementary line 479: do you mean TP instead of TN? eq 19, what does the lowercase p stand for in FPRp?*\n\nThanks for spotting the typos! The citations are based on papers performing Bayesian analyses with Bayesian networks, and while Bayesian networks are wide-spread, Bayesian inference for Bayesian networks is less so. For example, citation 25 is the only Bayesian network paper cited in a recent review of Bayesian statistics (doi.org/10.1038/s43586-020-00001-2). The lowercase p means we normalise by P rather than N so that FPs and FNs count the same in the comparisons.\n\n* *To make the point lines 168 ff., would it be worth counting the number of score evaluations?*\n\nThis would indeed be an interesting additional result. Figure 1 in the additional pdf in the global response shows the number of scores that were evaluated for the same simulations of Figure 1 in the main text. In red is the number of scores computed via the Laplace approximated posterior (“Laplace approximate”), in light green the number of additional scores evaluated via bridgesampling for the partition MCMC algorithm (“Bridgesampling, partition”) and in dark green the same quantity for order MCMC. \nNaively running the MCMC algorithm without our importance sampling approach would require computing the number of scores in red with the computationally expensive bridge sampling method. On the other hand, our importance sampling approach requires computing via bridge sampling only the number of scores shown in light/dark green, leading to much more efficient procedure.\n***\n* *Do some real-world scenarios require a different choice from the uniform prior over graph structures (line 141) ?*\n\nYes, and knowledge of the underlying graph can be easily incorporated into the graph prior. For example, in the real world data example in Section 5, prior knowledge of the edges that were experimentally validated by Sachs et al. (2015) would normally be included by up- or down-weighting the edges in question. Of course, to avoid circular reasoning we didn’t use them in the experiments on the Sachs data, and in fact in all of our experiments (unless stated otherwise) we employ a uniform prior when examining the performance of our approach. For real-world analyses typically one chooses a non-uniform prior informed by the setting, and indeed the cited application references [1,25,30] all have this (using either a sparsity cut-off, information from a separate biological database, or excluding edges on the basis of temporal restrictions, respectively).\n***\n* *Table 1, rightmost column: can the poor performance of your method be explained by the fact that there is no reason to predict edge absence, when absence of correlation can be captured with a fairly wide conditional distribution, i.e. low covariance in GP?* \n\nEdge absence implies a conditional (including marginal) independence relation between the two variables in question. However, edge absence implies neither zero correlation nor zero partial correlation in general (see for example A. Vargha, 1996). Predicting edge absence can therefore be more challenging, and we agree that the absence of an edge in the “ground truth” network could have different interpretations. In particular we find stronger evidence for an edge when we allow non-linear relationships than in the linear cases, which could be an indication of more complex dependencies. Of course we didn’t want to make too strong statements, but this aspect may be worth discussing in more detail and we could adapt accordingly. \n"}, {"review_id": "R8sJHCrxnA", "reviewer": "Reviewer_tnR2", "summary": "This paper proposes an MCMC algorithm for Bayesian inference for Gaussian Process Networks where you model the distribution of a node as a function of its parents plus noise, with the function a Gaussian process.\n\nThe sampler uses a Laplace approximation to make informed moves between Networks.\n\nThe sampler is demonstrated on a simulated and real example.", "questions": "It would be helpful to see the robustness of the results of your method as you vary the prior.\n\nIt would be good to see more thorough evaluation of the methods -- in particularly looking at understanding when the method works well/when it does not.\n\nIt would be helpful to have (easy to use) code available and details of how to reproduce the results. ", "limitations": "N/A", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "This is a fine paper -- which presents a sensible method for a general learning problem. It compares with some other methods and empirically gives some evidence of better performance.\n\nThe presentation is generally good.", "weaknesses": "The ideas come across as rather incremental -- taking standard ideas and putting them together.\n\nThe question then becomes how important and useful is the resulting method. Here is where I think the paper is a bit lacking. The empirical study is somewhat limited. Also, with their method there a numerous tuning parameters (i.e. the specification of the prior) -- and it is unclear how robust their results are to different choices of the priors. It would also be good to see some investigation of when and why the method works well relative to other methods. I can see that the Bayesian paradigm potentially gives advantages in terms of averaging over model uncertainty and being able to quantify uncertainty -- but these come with certain disadvantages (potential dependence on prior; extra computational cost; how reliable are measures of uncertainty when the model is incorrect).\n\nThe results could be a bit clearer -- e.g. more information captions; make it clearer the E-SHD is smaller for better methods; perhaps clearer about which method is which.\n\nIs there any code available for the method?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their comments, and acknowledge the time they took to assess the paper. \n* *It would be helpful to see the robustness of the results of your method as you vary the prior.*\n\nThe priors we use for the hyperparameter set are standard, non-informative priors (see for instance Gramacy & Lee, 2009). Nevertheless, we agree that it would be interesting to see to what extent the results are dependent on the choice of prior. In Figure 2 of the new additional pdf in the global response we compare the results in terms of E-SHD for our GP score-based approach for different prior distributions on the length-scale hyperparameter $\\theta$. In addition to the default IG(2,2) prior, we explore a IG(1,3) prior that favors higher lengthscales, as well as a IG(3,1) prior favoring lower values of $\\theta$. The results show that varying the prior hyperparameter has a limited effect on the results. \nAs for the prior over the graphs, we explore different choices when generating the ROC-like plots in Figure A.2 in the supplementary material.\n***\n* *It would be good to see more thorough evaluation of the methods -- in particularly looking at understanding when the method works well/when it does not.*\n\nIn our simulations, we look at different levels of non-linearity for specifically this purpose - comparing the performance of our approach under situations when linear models such as BGe are expected to work better. \n***\n* *[...] make it clearer the E-SHD is smaller for better methods*\n\nWe agree this could provide additional clarity, and will add a short comment under equation (15) mentioning this, though, as we previously noted in the manuscript, the E-SHD is not a direct measure of distance; the true posterior has a non-zero distance, so overly small values are not best.\n***\n* *Is there any code available for the method?*\n\nAs mentioned in lines 187-188, code to implement the method and reproduce the results of Sections 4 and 5 is publicly available.\n***\nWe hope the above comments have addressed the reviewer’s concerns and will be happy to hear back from them.\n"}, {"review_id": "76BrplZp4L", "reviewer": "Reviewer_uGD4", "summary": "The paper proposes methodology to perform Bayesian inference on the (hyper) parameters of a so-called Gaussian Process Networks (GPNs), which are sets of functional equations with Gaussian Process (GP) priors on the functions relating a variable to its parents, as well as inference on the graph structure and graph posterior expectations.\nPrevious works in the GPN literature that consider a posterior over graphs did not put priors over hyperparameters of the GP kernels and perform associated posterior inference. \nTo address the challenges, the authors crucially exploit the decomposability property of p(G | D): \" scores \" (conditionals of X_i given PA_i ) only need to be recomputed if the parents have changed, so that scores can be reused across different graphs. \nSince scores are intractable integrals over (hyper) parameters of the GPs, the authors use bridge sampling to approximate the integrals, and make use of Laplace approximations  in a two-step procedure to reduce computational cost compared to running MCMC directly with the estimated scores. \nThe methodology is evaluated on a good set of synthetic and real-data examples. \n", "questions": "- How is the Gaussian proposal for the numerator of Eq. (11) chosen ? I would have expected a Laplace - IS here too ? Perhaps clarify what is the choice made here and possibly why.    \n - Could you comment theoretically/intuition-wise/experimentally on which of these approximations are the ones that affects accuracy most and how ? For example, should we focus on getting better proposals $q(G|D)$ or in approximating the scores with MC (bridge sampling, in your case) more efficiently to approximate the target (numerator of the IS weights, in IS terminology) better ? Further, could you provide a couple of example $\\Phi$ of interest already in the methodology (Section 3) for concreteness ? \n    - The comment in lines 290-293 is important in my view, but should be expanded upon. Why/how does marginalizing the hyperparameters mitigate the bias of the Laplace approximation? Why can we not \"just\" not marginalize them and simply use more MCMC samples ? \n- Regarding score equivalency, could you clarify why it sounds like it's not good if DAGs from the same MEC obtain different scores ? Does the property make it harder to approximate the true $p(G|D)$ somehow ?  \n- Does any Bernstein von Mises - type result hold for $p(G|D)$ ? \n", "limitations": "The authors do not have an explicit limitations subsection, although limitations are not hard to infer from the text, including mainly computational cost of marginalizing the hyperparameters, and the unclear effects on the performance of the approximations performed by bridge sampling; it is also not clear whether the authors can use kernels appropriate for with non-stationary functions. \n", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "Overall a technically solid paper on an interesting topic for the NeurIPS community. Main strengths are as follows: \n- GPNs are gaining more attention in particular due to their application in causal structure learning (von Kügelgen et al 2019), so the authors address a relevant problem of performing \"fully\" Bayesian inference with these models \n- The scope and objective of the paper and proposed methods are clear and targeted \n- The paper demonstrates an excellent knowledge of the surrounding relevant literature on MCMC for structure learning, and describes the relationship with the current work clearly \n- Well explained and solid experimental section  \n\n", "weaknesses": "- [**Presentation of Section 3**] The main weakness of the paper is the presentation of the approximations in Section 3. Even as a reader well familiar with bridge sampling, Laplace approximations and MCMC, the part from lines 143 to 172 is unclear at first read and too condensed. I am guessing the final goal is to approximate p( \\Phi | D) of Eq (13). Rather than start from the various approximations needed to get there, I would suggest starting from there and write down the integrals that arise, then explain how to approximate each step. Otherwise, while I'm reading I keep thinking \"Why am I approximating the scores ? What is the final goal ?\". Some related questions/comments about clarity of this part below.\n    - Is the reason that you approximate the scores both **(a)** via Laplace  **and**  **(b)** via bridge sampling because  **(a)** is used to build the proposal while  **(b)** is needed to compute the (self-normalized) estimator in Eq. (13), specifically, the unnormalized target ? If so you should note that the estimator in Eq (13) has an **additional bias** ( in addition to the standard bias coming from self-normalized IS)  due to the fact that the bridge sampling estimator is not unbiased. You should also state that because you are using self-normalized IS, you do not need to know the normalizing constants of either q(G|D) or p(G|D), as the ratio of these two constants cancels out. \n- [**Difficulty introduced by the additional marginalization**]: while usually people assume there is inherent interest in a \"fully\" Bayesian treatment and marginalize hyperparameters as well when feasible (in this case, kernel hyperparameters), the authors should make it clearer why it is beneficial to do so specifically in the setting of GPNs with interest in p(G|D). There is a comment about this in lines 290-293 but it is quite vague.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their constructive comments and interesting questions. \n* *How is the Gaussian proposal for the numerator of Eq. (11) chosen ? I would have expected a Laplace - IS here too ? Perhaps clarify what is the choice made here and possibly why.*\n\nThe proposal function g() was set to a normal distribution, with its first two moments chosen to match those of the posterior distribution. This is a standard choice in the literature and the default in the bridge sampling implementation that we used (Gronau et al., 2017). We will add a clarification of this when discussing implementation details in Section 3.1.\n***\n* *Could you comment theoretically/intuition-wise/experimentally on which of these approximations are the ones that affects accuracy most and how ? For example, should we focus on getting better proposals $q(G|D)$ or in approximating the scores with MC (bridge sampling, in your case) more efficiently to approximate the target (numerator of the IS weights, in IS terminology) better ?*\n\nThis is an interesting and relevant question which is however challenging to address in all generality. Improving the proposal distribution can provide more uniform weights, leading to a more efficient exploration of the posterior and reducing variance. More accurate proposals would however require more computational power and could undermine the point of having a proposal distribution in the first place. Improving the MC estimates of the scores would on the other hand reduce the bias deriving from the bridge sampling approach (as you noted in your comments). The approximations and trade-offs used in the experiments seemed to perform quite well, though determining the best approach to perform posterior inference would be an interesting area of experimentation in this domain.\n***\n* *Further, could you provide a couple of example $\\Phi$ of interest already in the methodology (Section 3) for concreteness ?*\n\nThree examples of $\\Phi$ are already provided in Section 2.1, lines 61-62.\n***\n* *Why/how does marginalizing the hyperparameters mitigate the bias of the Laplace approximation?*\n\nThe Laplace approximation is meant as an approximation of the score function, since it does not fully take into account the uncertainty regarding the hyperparameters, and only provides a good approximation if the posterior over the hyperparameters is approximately Gaussian. \nMarginalizing the hyperparameters on the other hand provides an accurate estimation of the score because it directly follows the definition of the score in equation (10). This process of marginalizing is however computationally costly, leading us to use the Laplace approximation in the first step of our approach. Finally we correct the approximation with importance sampling.\n***\n* *Does any Bernstein von Mises - type result hold for $p(G|D)$ ?*\n\nWhile we are not aware of any such result, we agree it would be interesting if any such results could be established and an interesting line of research.\n***\n* *Regarding score equivalency, could you clarify why it sounds like it's not good if DAGs from the same MEC obtain different scores ?*\n\nTwo DAGs from the same MEC are theoretically indistinguishable from observational data. The lack of score equivalence in our case is caused by the GPN assumption, which introduces asymmetry into the factorizations of the joint distribution via the structural equation model (8). Because score equivalence is a common, and in some cases desirable, property to avoid strong distributional assumptions, we simply wished to emphasize that the identifiability of otherwise score equivalent DAGs is a consequence of the assumptions underpinning the GPN model. \n***\n* *[...] it is also not clear whether the authors can use kernels appropriate for with non-stationary functions.*\n\nWhile we didn’t explicitly mention non-stationary kernels, the method is applicable to any kernel function as long as a marginal likelihood can be computed.\n"}, {"review_id": "cNcAItslZi", "reviewer": "Reviewer_sMvJ", "summary": "This paper proposes a Bayesian structure learning of the GPNs framework that is claimed to be less computationally. To address this, the approach presented in this work utilizes Monte Carlo and importance sampling to sample from the posterior distribution of network structures. This approach compares models using their marginal likelihood and computes the posterior probability of GPN features. Some simulations have been discussed. ", "questions": "A few questions that come to mind: \n- Laplace approxiamtion is sensitive to the MAP estiamte. as the score s are being approximated this may cause poor approxiamtion and bias and thus it may impact the perfomrance of the subsequenct imporatnce sampling and causing high variance of imporance weights and thus inefficient sampling and inaccurate estmation. could the paper elaborate on this issue. \n- Imporatnce sampling can introduce bias if the proposal does not adequetly cover the suppot of the target posterior. Also, high variance of the weights can lead to unstable estiamtes. \n- Althogh use of Lapace approxiamtion and imporatance sampling may in general reduce the computaional intensity  but overall performance dependes on model and the size of dataset (specifically when approxiamting the score and so on). It is well-known that as dimension of the netwrok increases or dataset grows larger the computational cost of Laplace approxiamtion and imporatnece samling can be significant. How is the relationship between effective sample size, variabaility of the weights and growth of the netwrok, dataset, and dimentionality justified? ", "limitations": "- Laplace approximation is only valid if the distribution is Gaussian or close to Gaussian. This assumption can simply be violated. \n- Lack of discussion on computational intensity and how this framework overcomes those known computational issues. \n- Lack of validation, sensitivity analysis and comparison to alternative models. ", "rating": 3, "confidence": 5, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- This paper studies a rather important problem.\n- It is easy to follow the thoughts. ", "weaknesses": "- The contribution of the paper is not clearly justified and seems moderate at best. It is a combination of a few approaches (Gaussian processes, Laplace approximation, and Importance sampling).\n- The claims of the paper are not completely justified. It provides an approach but does not discuss how the main concerns are addressed. It is not clear how efficient this method is relative to others. \n- Some references are missing, for instance: \n     * Efficient and scalable structure learning for Bayesian networks: algorithms and applications by Zhu 2020\n     * Efficient structure learning of Bayesian networks using constraints by de Campos \n     * MCMC algorithms for Bayesian analysis of phylogenetic trees\n- Lack of comparison/validation to other state-of-the-art algorithms and also lack of discussion on how efficient this framework is by providing computational intensity analysis. \n- Lack of theoretical discussion regarding the convergence and issues rising from Laplace approximation and Importance sampling. \n- Lack of validation, sensitivity analysis and comparison to alternative models. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their comments, and acknowledge the time they took to assess the paper. \n* *The contribution of the paper is not clearly justified and seems moderate at best. It is a combination of a few approaches (Gaussian processes, Laplace approximation, and Importance sampling).*\n\nThe objective of our work is addressing an important problem, as recognised also in the review, indeed drawing on technical tools originating in different domains of the statistical inference literature and we combine them effectively in a novel whole. The Bayesian approach accounts for uncertainty both in the parameters and structure, allowing us to compute posterior probabilities of network features in a principled manner. \n***\n* *The claims of the paper are not completely justified. [...] It is not clear how efficient this method is relative to others.*\n\nThe supplementary material already included analyses of the run-times illustrating the computational costs of the procedure, which are expected to be higher compared to simpler methods. The added value is that this is the only method correctly targeting the posterior distribution (as illustrated for example in Figure 2 in the main text). \n***\n* *Some references are missing, for instance [...]*\n\nThe references suggested are not entirely pertinent to the topic covered in the paper. The first two references do not qualify as Bayesian methods, since they do not provide a posterior distribution over DAGs and/or do not compute the marginal likelihood of the data, opting for a generic score function instead. The last reference concerns phylogenetic trees and is not directly related to our problem of Bayesian network structure learning. However, we agree that to better embed the work in the literature, adding some text on non-Bayesian point estimate approaches, and then modifying the current sentence to better emphasize the Bayesian side would work well, and we will modify accordingly.\n***\n* *Lack of comparison/validation to other state-of-the-art algorithms and also lack of discussion on how efficient this framework is by providing computational intensity analysis.*\n\nThere are few methods that currently are able to perform Bayesian structure learning for BNs with non-linear, continuous data. Again, as for the references, we have to distinguish between Bayesian and non-Bayesian approaches. We of course welcome concrete suggestions for other algorithms to include. As for the computational intensity analysis, we already examined this in Figures A.3, A.5 and A.8 in the supplementary material.\n***\n* *Lack of validation, sensitivity analysis and comparison to alternative models.*\n\nWe compared the performance of our method to five different state-of-the-art algorithms, and provided a comparison of run-times in Figures A.3, A.5 and A.8 in the supplementary material. Furthermore, we have now added new figures in the additional pdf in the global response, which include an analysis of our method’s performance for different choices in the prior hyperparameters.\n***\n* *Laplace approximation is only valid if the distribution is Gaussian or close to Gaussian. This assumption can simply be violated.*\n\nThe Laplace approximation is indeed only exact if the posterior distribution is Gaussian; we are however not relying directly on this approximation to compute the marginal likelihoods. We are using the Laplace approximation as a proposal distribution for subsequent importance sampling. The exact posterior is then computed by marginalizing the likelihood with respect to the prior using bridge sampling. The results so obtained finally provide valid samples.\n***\n* *Laplace approxiamtion is sensitive to the MAP estiamte. as the score s are being approximated this may cause poor approxiamtion and bias and thus it may impact the perfomrance of the subsequenct imporatnce sampling and causing high variance of imporance weights and thus inefficient sampling and inaccurate estmation.*\n\n* *Imporatnce sampling can introduce bias if the proposal does not adequetly cover the suppot of the target posterior. Also, high variance of the weights can lead to unstable estiamtes.*\n\nThe same Laplace approximation is used both in the target for graph sampling and in the importance sampling. Even if the approximation is poor, these terms “cancel” out and do not induce bias, though a poor approximation may increase variance. A high variance in the weights will only be observed when there is a large discrepancy between the proposal and true posterior. Using a Laplace approximation as a proposal distribution for importance sampling is commonplace in the literature (see e.g. Kuk, 1998 or Bek et al., 2018) since it provides a reasonable approximation of the target distribution. There are also results showing convergence in Hellinger distance of the Laplace approximation to the posterior (Schillings et al., 2020). \n***\n* *It is well-known that as dimension of the netwrok increases or dataset grows larger the computational cost of Laplace approxiamtion and imporatnece samling can be significant. How is the relationship between effective sample size, variabaility of the weights and growth of the netwrok, dataset, and dimentionality justified?*\n\nWe are unfortunately not entirely sure what is meant with this comment, but we will be happy to discuss during the discussion period.\n***\nGiven the response comments above, and that many of the concerns were already addressed in the original submission, we would appreciate a reconsideration of the perceived weaknesses of the paper and hopefully an adjustment of the score accordingly."}, {"review_id": "yXBMQCoFuD", "reviewer": "Reviewer_9qjJ", "summary": "Gaussian Process (GP) networks are directed graphical models for continuous data, where the function mapping from parent node values to parameters of child node is a Gaussian process. Given the graphical network structure and a dataset of observations, learning the GP model for each node simply reduces to a standard supervised GP regression using parents nodes as input $x$ values and child node as output $y$ values. Likewise marginal likelihood evaluation and hyperparameter optimization/marginalization can use the standard methods for GP regression.\n\nIn this work, as I understand, the authors goal is to marginalise over network structures, each network structure yields a different set of parents-child links and each one requires marginalizing GP hyperparameters again, if done naively, this is very very expensive. The authors thus propose methods to speed up this process hence making network structure marginalisation much more computationally efficient. These methods include\n- using Metropolis Hasting like algorithms to make small steps through the space of graphs thereby removing the need for excessive recomputation.\n- for each graph in the Metropolis Hastings proposal, the authors use a Laplace approximation to marginalize GP hyperparameters instead of full MCMC marginalisation.\n- once a dataset of sampled graphs is collected, the expensive MCMC marginalisation is used to weight each sampled graph with importance sampling.\n\n\n\nExperiments on both synthetic toy example and a protein interaction network example are presented.", "questions": "- as mentioned above regarding $\\lambda=0$, is my understanding correct? Can the authors comment on this?\n- Figure 1: I assume the graph hamming distance shows high \"tight\" the posterior over graphs is clustered around the true graph, tighter is better, whereas with 100 observations in 10D space, this seems very sparse dataset to learn any complex model and I would expect the posterior to be rather broad. What is the hamming distance of uniformly randomly generated graphs? What is the average hamming distance of the \"true\" posterior from the true generating graph?", "limitations": "- the authors provide an explanation of the Score equivalence issue (which hadn't occurred to me until it was mentioned)", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- I feel that the proposed methods are very standard approaches and I feel are well justified in this use case and the overall model and algorithm is is well designed.\n- personally coming from a traditional GP background, I felt this paper was accessible and easy to digest.", "weaknesses": "- Everything from line 101 up until line 159 is focused on model fitting and hyperparameter marginalization for a single GP and, as far as I can tell, has nothing to do with networks, hence the section title felt a bit misleading and rearranging material might make this easier. Equations 8 - 12 could hypothetically be written using $x, y$ instead of $Pa_{X}, X$ and moved to Section 2. Then section 3 could contain a description about how graph structure determines each nodes \"local GP training set\" and how the model is a summation over GPs over individual training sets. and the fact that there are so many justifies a hyper-parameter approximation.\n- Equation (10) appears to be a standard marginal likelihood with an extra layer of marginalisation over kernel hyperparameters (with or without hyperpriors), this is standard practice in GP regression [Code here](https://github.com/rmgarnett/gpml_extensions/) and [this paper](https://proceedings.neurips.cc/paper/2010/hash/4d5b995358e7798bc7e9d9db83c612a5-Abstract.html). The authors cite bridge sampling which I am not personally familiar with, but there exist many, many algorithms to solve this problem (this is not a weakness of the paper as any MCMC algorithm will suffice, but at first glance seems a rather strange choice)\n\n- In Section 4, the authors repeatedly state that $\\lambda =0$ removes higher frequencies and is a linear case and so $w_{i,0}$ is the only non-zero coefficient. If I take equation 14 and make the suggested substitution, I have a constant with noise\n$$\nX_i = \\sum_{Z\\in Pa_{X_i}} \\beta_i w_{i,0} + \\epsilon_i\n$$\nWhich implies no dependency between $X$ and $Z$. (The only way I could see a linear relationship between $X$ and $Z$ is if $\\sin(jZ)$ for very small $j$ where $\\sin()$ is approximately linear near the origin, but given normalised $Z$ and $j=1$ this isn't likely). Hence the $\\lambda=0$ case is generating data that ignores graph structure, and the hamming distance is meaningless as the models will learn the empty graph from the data. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their constructive comments and suggestions. \n\n* *as mentioned above regarding $\\lambda=0$, is my understanding correct? Can the authors comment on this?*\n\nThere is indeed a typo in equation (14), the linear terms should have been included and the index $j$ of the second sum should start at 1 instead:\n$$ X_i \\\\, = \\sum_{Z \\\\, \\in \\\\, \\textrm{Pa}_{X_i}} \\beta_0 \\\\, w\\_{i,0} \\\\,Z \\\\,+\\\\, \\\\left\\\\{\\sum\\_{j=1}^6 \\beta_i \\\\,v\\_{i,j} \\sin{(j Z)} \\\\,+\\\\, \\beta_i\\\\, w\\_{i,j} \\cos{(j Z)} \\\\right\\\\} \\\\,+\\\\, \\epsilon_i. $$\nIt was however coded correctly in the experiments. We thank the reviewer for this careful observation, and will update equation (14) in the manuscript.\n***\n* *Figure 1: I assume the graph hamming distance shows high \"tight\" the posterior over graphs is clustered around the true graph, tighter is better, whereas with 100 observations in 10D space, this seems very sparse dataset to learn any complex model and I would expect the posterior to be rather broad. What is the hamming distance of uniformly randomly generated graphs? What is the average hamming distance of the \"true\" posterior from the true generating graph?*\n\nTighter is indeed better, but only up to a point, because as noted the posterior is spread and its average distance is non-zero. To check convergence to the true posterior we need measures like the K-L divergence in Figure 2 of the main text, but this requires enumerating all DAGs and computing their scores. The current best complexity for exact computations scales like 3^(# of nodes) times the score computation time, and is only feasible for small graphs (as in Figure 2). The number of DAGs with 10 nodes is over 4 quintillion, making it infeasible to compute the exact E-SHD of the true posterior distribution. However, since we have samples from the posterior, the “GP partition” results in Figure 1 of the main text are Monte Carlo estimates of the average Hamming distance between the true posterior and the true generating graph (and the “BGe partition” for a linear-Gaussian model). Exact computations for the uniform distribution over DAGs are cheaper (there is no score component) and MC estimates are simple, so this could be added, but as an average DAG with 10 nodes has more than 20 edges, the distances would, as we might expect, be high compared to the scale of the figure.\n***\n* *the authors provide an explanation of the Score equivalence issue (which hadn't occurred to me until it was mentioned)*\n\nWe’re glad the reviewer appreciated this aspect, as it is often overlooked. It is unclear to us why this is a limitation though.\n"}]

Author rebuttal:
Following the reviewers' suggestions, we attach a pdf containing additional simulation results.

Meta-review: Four knowledgeable referees support to accept. A fifth reviewer recommends rejection. Having looked at this reviewer's comments and the rebuttal, I side with the authors on two important points: 1) the paper, indeed, is related to a different problem than the one alluded to in the references provided by the reviewer; 2) the reviewer does not provide a concrete list of methods that could be used as baselines in the experiments, rather refers to a lack of comparisons against state-of-the-art approaches. Several other concerns by the reviewer were already addressed in the original submission. I therefore recommend accepting. Authors must follow the reviewer's suggestions, in particular, addressing their comments on the presentation of the approximations in Section 3.
Decision: Accept (poster)