Continuous Bayesian model selection for multivariate causal discovery

In brief

• Observational causal discovery typically requires rigid functional assumptions (like additive noise models) to ensure structure identifiability.
• Continuous Bayesian model selection relaxes these rigid constraints by scoring graph hypotheses via the marginal likelihood.
• Using Causal Gaussian Process Conditional Density Estimators (CGP-CDE), the method continuously optimises adjacency weights while enforcing acyclicity.

This retrospective explainer is dated to the November 2024 preprint release. It discusses Continuous Bayesian model selection for multivariate causal discovery by Anish Dhir, Ruby Sedgwick, Avinash Kori, Ben Glocker and Mark van der Wilk.

Imagine being handed observational data for temperature, humidity, and barometric pressure, and being asked to draw the causal arrows between them without running controlled experiments.

Because multiple directed acyclic graphs (DAGs) can produce the identical conditional independencies—a phenomenon known as Markov equivalence—traditional causal discovery algorithms rely on strict mathematical assumptions to pick a unique graph. Standard methods might assume that the causal mechanisms are linear with non-Gaussian noise (LiNGAM) or that the error terms enter strictly additively ($y = f(x) + \epsilon$). If real-world physics or biology violates those rigid equations, the mathematical guarantees evaporate, and the predicted graph can be completely wrong.

This paper proposes an alternative route for the multivariate setting: Continuous Bayesian Model Selection. Instead of enforcing brittle functional forms, it treats causal discovery as a Bayesian hypothesis comparison problem over continuous graph relaxations.

Escaping rigid additive noise assumptions

Consider two causal variables $X$ and $Y$. In an Additive Noise Model (ANM), we assume:

\[Y = f(X) + \epsilon_Y, \quad \epsilon_Y \perp\!\!\perp X.\]

If the true relationship contains heteroscedastic noise (the noise variance depends on $X$) or non-linear interactions ($Y = f(X, \epsilon_Y)$ where $\epsilon_Y$ cannot be neatly separated), ANMs fail.

Bayesian model selection offers a principled escape. Rather than declaring that an inverted causal direction is mathematically impossible, it asks: under which causal ordering is the observed dataset more probable after integrating over all possible non-linear mechanisms? By evaluating the marginal likelihood:

\[p(\mathbf{X} \mid G) = \int p(\mathbf{X} \mid \boldsymbol{\theta}, G) \, p(\boldsymbol{\theta} \mid G) \, d\boldsymbol{\theta},\]

the Bayesian Occam’s razor penalises overly complex models automatically. In the bivariate case ($X \to Y$ vs $Y \to X$), the marginal likelihood overwhelmingly favours the true generative direction with only a negligible error probability.

Scaling to multivariate graphs continuously

Scaling Bayesian model selection to $d$ variables is notoriously difficult because the space of DAGs grows super-exponentially. Evaluating $p(\mathbf{X} \mid G)$ for every candidate discrete graph quickly becomes intractable.

The paper overcomes this by continuously relaxing the discrete search space using Causal Gaussian Process Conditional Density Estimators (CGP-CDE). Let $W \in \mathbb{R}^{d \times d}$ denote a continuous weighted adjacency matrix where entry $W_{ij} \ge 0$ governs the influence of variable $j$ on variable $i$. Each conditional density is modeled as a non-parametric Gaussian Process:

\[p(X_i \mid \mathbf{X}_{\text{pa}(i)}) = \mathcal{N}\left(X_i; \mu_i(\mathbf{X}_{\text{pa}(i)}), \sigma_i^2(\mathbf{X}_{\text{pa}(i)})\right).\]

To ensure that the learned matrix $W$ corresponds to a valid Directed Acyclic Graph rather than a cyclic network with feedback loops, the algorithm adopts the continuous algebraic characterisation of acyclicity:

\[h(W) = \mathrm{Tr}\left(e^{W \odot W}\right) - d = 0,\]

where $\odot$ is the Hadamard (element-wise) product and $e^M$ is the matrix exponential. When $h(W) = 0$, the graph contains zero directed cycles.

The overall learning objective becomes:

\[\max_{W \ge 0} \quad \log p(\mathbf{X} \mid W) - \lambda h(W) - \rho \|W\|_1,\]

where $\log p(\mathbf{X} \mid W)$ is the continuous marginal likelihood, $\lambda$ penalises cycles, and $\rho$ encourages structural sparsity.

Figure from the paper. Structural Hamming Distance (SHD ↓), Structural Intervention Distance (SID ↓), and F1 score (↑) for 3-variable causal discovery benchmarks comparing NOTEARS, SDCD, CAM, NOGAM, SCORE, DiBS, and the proposed CGP-CDE and DGP-CDE frameworks. Source: Anish Dhir, Ruby Sedgwick, Avinash Kori, Ben Glocker and Mark van der Wilk. Extracted from the source manuscript; click to enlarge.

Empirical performance on synthetic and real graphs

The framework is evaluated on synthetic and benchmark DAGs against continuous optimization methods (such as NOTEARS) and non-linear causal discovery methods (such as CAM, SCORE, and DiBS):

1. Lower Structural Error: As shown in the figure, the proposed Gaussian Process models (CGP-CDE and Deep Gaussian Process DGP-CDE) achieve near-zero Structural Hamming Distance (SHD) and Structural Intervention Distance (SID) on 3-variable systems, significantly outperforming NOTEARS and heuristic baselines.
2. Robustness to Heteroscedastic Noise: Because the conditional density estimator models input-dependent variance, the framework correctly recovers causal directionality in noise regimes where standard additive-noise assumptions fail completely.
3. Continuous Tractability: The smooth relaxation allows standard gradient-based optimisers (such as Adam) to navigate the combinatorial landscape of causal topologies efficiently.

Paper and sources

• Anish Dhir, Ruby Sedgwick, Avinash Kori, Ben Glocker and Mark van der Wilk. Continuous Bayesian model selection for multivariate causal discovery. arXiv:2411.10154 (2024).