## 4 Application to Gaussian Graphical Models
### Hyper Inverse-Wishart Induced Cholesky Factor Density  
We first consider the case where \(G\) is decomposable. The hyper-inverse Wishart (HIW) distribution [34] for \(\Sigma=\Omega^{-1}\) and the corresponding induced class of distributions [14] for \(\Omega\) are attractive choices for priors. Given \(G\), we place a hyper-inverse Wishart prior \(\operatorname{HIW}_{G}(\delta,\Lambda)\) on \(\Omega=\Sigma^{-1}\), where \(\delta>2\) is the degrees of freedom and \(\Lambda\in\mathbb{S}^{p}_{>0}\) is fixed. The HIW distribution is defined over the cone of \(p\times p\) positive definite matrices, with the corresponding density:  
\[f(\Omega\mid G)\propto\,|\Omega|^{(\delta-2)/2}e^{-\operatorname{tr}(\Omega \Lambda)/2}. \tag{10}\]  
Despite being able to sample from the posterior distribution, \(\operatorname{HIW}_{G}(\delta+n,\Lambda+B)\), we cannot directly carry out the EPSOM-Hyb algorithm. Since samples are drawn from a sub-manifold of \(\mathbb{R}^{p\times p}\), proceeding to obtain a partition over \(\mathbb{R}^{p\times p}\) does not guarantee that a given point in the partition can be reconstructed to form a valid covariance matrix. As such, we circumvent this issue by working with an alternative representation of \(\Sigma\) given by the Cholesky decomposition \(\Sigma^{-1}=\Omega=\phi^{\prime}\phi\). Provided that the vertices of \(G\) are enumerated according to a perfect vertex elimination scheme, the upper triangular matrix \(\phi\) observes the same sparsity as \(\Omega\). Consequently, we need only compute the induced prior on the nonzero elements of \(\phi\), which, together with the likelihood in Equation (9), gives us anexplicit expression for the negative log posterior \(\Psi\). The determinant of the Jacobian matrix \(J\) of the transformation \(\Omega\to\phi\) is given by \(|J|=2^{p}\prod_{i=1}^{p}\phi_{ii}^{\
u_{i}+1}\), where the \(i\)-th row of \(\phi\) has exactly \(v_{i}+1\) nonzero elements [14]. Further, the distribution of \(\Sigma\) has the strong hyper-Markov property [35], so we can ascertain mutual independence of the rows of \(\phi\). Then, we can specify the joint density of the free elements of \(\phi\) as follows,  
\[\pi(\phi)=\left[\prod_{i=1}^{p}\frac{2^{-(\delta+v_{i})/2}}{\Gamma((\delta+v_{ i})/2)}\phi_{ii}^{\delta+v_{i}-2}e^{-\frac{1}{2}\phi_{ii}^{2}}(2\phi_{ii}) \right]\times\Bigg{[}\prod_{(r_{s})>r,(v_{s},r_{t})\in E}\frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2}\phi_{rs}^{2}}\Bigg{]}. \tag{11}\]  
Recall that we also require expressions for the gradient and Hessian of \(\Psi\), both of which are provided in Appendix D. Putting these ideas together, we can employ the EPSOM-Hyb estimation framework by drawing samples from the posterior distribution, taking the Cholesky factor of each sample, and computing \(\Psi(\phi)\) using the likelihood in Equation (9) and the prior in Equation (11). While this procedure appears to be quite simple, the implications are significant; if we have a different prior on \(\Sigma\) for which we can do the posterior computation, all other aspects of the algorithm would remain the same. All that is required is a way to sample from the posterior of \(\Omega\) and an expression for the Jacobian of the corresponding transformation. In our simulations, we emulate a high-dimensional setting by stacking the adjacency matrix of the decomposable graph \(G_{9}\) in Figure 1 multiple times along a block diagonal to construct larger graphs, e.g., \(A_{G_{p}}=I_{r}\otimes A_{G_{9}}\), where \(A_{G_{p}}\) is the adjacency matrix of \(G_{p}\) and \(r\) is the number of times we stack the graph. For \(r=8,10\), we form the graphs \(G_{72}\) and \(G_{90}\), which have \(d=200,250\) free elements, respectively. Drawing data that satisfy their corresponding GGMs and taking \(\delta=3,n=100\), and \(\Lambda=I_{p}\), we can sample from the posterior distribution and compute the marginal likelihood estimates.