## 1 Introduction  
In Bayesian inference, many tasks, such as model selection, averaging, and comparison, rely on being able to compute the marginal likelihood (or model evidence) in order to form the Bayes factor. The marginal likelihood is the normalizing constant of a posterior distribution, which plays a particularly important role in probabilistic graphical models. Suppose \(\gamma\) is a general Gibbs distribution defined on \(\mathcal{U}\subseteq\mathbb{R}^{d}\), \(\gamma(u)=\mathcal{Z}^{-1}e^{-\Psi(u)}\), with  
\[\mathcal{Z}=\int_{\mathcal{U}}e^{-\Psi(u)}du<+\infty,\quad\mathcal{U}\subseteq \mathbb{R}^{d}. \tag{1}\]  
In a Bayesian setup, \(\Psi:\mathbb{R}^{d}\rightarrow\mathbb{R}\) is the negative log posterior. Since \(\Psi\) is typically complicated, and the space over which we are integrating tends to be high-dimensional, the resulting integral is usually intractable. This calculation is exacerbated further for Gaussian graphical models (GGM) [1, 2] because the integral is taken over a specialized subset of positive definite matrices. Generic algorithms for approximating the marginal likelihood include [PERSON]'s method [3], the harmonic mean estimator [4, 5], the corrected arithmetic mean estimator [6], annealed importance sampling [7], [PERSON]'s method [8], (warp) bridge sampling [9, 10], and nested sampling [11]. More recent developments include the integrated nested Laplace approximation [12] and the expectation propagation-approximate Bayesian computation algorithm [13]. While these methods can be applied to many marginal likelihood estimation problems, they do not provide a straightforward way to deal with the specificity of GGMs and the positive definite restriction on the precision matrices. Additionally, many of these aforementioned methods tend to rely on high-quality MCMC samples and have numerous problem-specific model settings that limit the overall practicality. While there have been developments that specifically target the marginal likelihood of GGMs ([PERSON], 2016), the stringent prior restrictions (inverse-Wishart and hyper-inverse Wishart) prevent many of these approaches from being used in broader contexts. Furthermore, there exists substantial literature on inference for both decomposable ([PERSON], 2016) and non-decomposable ([PERSON], 2016; [PERSON], 2016; [PERSON], 2016) graphs; however, the need for more general and scalable algorithms is ever present. This has consequently led to dedicated methods for sampling from the G-Wishart distribution and conducting model comparison ([PERSON], 2016; [PERSON], 2016; [PERSON], 2016). Also relevant is the availability of the G-Wishart normalizing constant in closed form ([PERSON], 2016), but the viability of these results is limited when considering most high-dimensional graphs of practical interest. [PERSON] et al. (2017) propose an application of [PERSON]'s method and a telescoping decomposition of the precision matrix that simplifies the ensuing marginal likelihood calculations, but the main advantages of this approach are seen with element-wise priors, and in most cases, the time complexity is on par with other GGM-specific methods. A recent approach that combines some of the probabilistic ideas from MCMC-based methods while drawing inspiration from quadrature is the hybrid estimator ([PERSON], 2016). This method moves away from an over-reliance on MCMC samples, as these can be time-consuming to obtain if the likelihood is expensive to evaluate. Rather, the MCMC samples are utilized during a preliminary step to learn a partition of the parameter space in order to identify regions of posterior concentration. By making local approximations to the negative log posterior over these partition sets, the the hybrid algorithm simplifies the integral and forms a piecewise estimate of the marginal likelihood. The high-probability partition of the parameter space yields multiple benefits. First, the areas of the parameter space that admit finer partitions encourage precise approximations to the log posterior.