Decomposable graphs are used in probabilistic graphical models for mathematical convenience to describe the dependence structure of a multivariate probability distribution.
We demonstrate two applications of the \smcabbr\, sampler from \cite{olsson2019}, related to structure learning in graphical models, where the expander and collapser kernels  $\{\bk{\jtind}\}_{\jtind=1}^{\p-1}$ and $\{\prop{\jtind}\}_{\jtind=1}^{\p-1}$ respectively,  and the factorisation formula for \(\mu\) are employed as presented in Section~\ref{sec:smc}. 

The first example shows an application of the particle Gibbs approach presented in \cite{olsson2019}  to Bayesian graph learning. 
In the second example, we provide unbiased estimates of the cardinality of the space of decomposable graphs of a given order.

The Python code used for generating the examples in this paper is contained in the \emph{trilearn} library available at \url{https://github.com/felixleopoldo/trilearn}.


\subsection{Bayesian structure learning in Gaussian graphical models}
It is well known that the Gaussian distribution encodes conditional independence statements by the non-zero elements in its precision matrix, which can be represented as a graph.
The model is sometimes termed a Gaussian graphical graphical model; see \emph{e.g.} \citep{10.2307/2241271}.

In this example, we approximate the graph posterior distribution of a zero mean Gaussian decomposable graphical model with graph \(\graph=(\graphnodeset,\graphedgeset)\) and adjacency matrix displayed in Figure~\ref{fig:AR1-5}.
To generate a  precision matrix \(\precmat\), we first sampled a \(\p\times\p\) matrix \(\Omega\) as
\begin{align*}
 \Omega_{ij}=\Omega_{ji} \sim \begin{cases}
        \unif[{(-\upsilon,-\epsilon)\cup (\epsilon,\upsilon)}], &\text{ if } (i, j)\in \graphedgeset\\
        \unif[{(\epsilon,\upsilon)}], &\text{ if } i=j,
    \end{cases}
\end{align*}
where \((\epsilon,\upsilon) \in [0,1]^2\).
Then, to ensure positive definiteness, we let \(\precmat = \Omega + \delta I_\p \), where \(\delta = \argmin_{\delta'\ge 0} (\Omega + \delta' I_\p \in \posdefsp_{\graph})\), \(I_\p\) is the \(\p\times\p\) identity matrix and \(\posdefsp_{\graph}\) denotes the space of positive definite matrices with zero elements corresponding to non-edges in \(\graph\).
We have sampled $n=100$  data vectors of dimension \(\p=30\) from this model with \(\epsilon=1\) and \(\upsilon=0.5\). 
The observations are stored in an \(n \times \p \) data matrix \(\y\).
In the following, for any matrix $p \times p$ matrix $Q$ and
    $\genseta \subseteq \nn[\p]$, we denote by $Q_\genseta$ the $|\genseta|\times |\genseta|$ matrix
    obtained by extracting the elements $(Q_{ij})_{(i, j) \in \genseta\times\genseta}$ from $Q$.

A fully Bayesian formulation of the structure learning problem requires two types of prior distributions, one defined on \(\graphsp[\p]\) and the other on \(\posdefsp_{\graph}\) for any \(\graph \in \graphsp[\p]\). 
For simplicity we assigned the uniform prior on \(\graphsp[\p]\) and for \(\posdefsp_{\graph}\) we use a hyper Wishart distribution having \(\df=\p\) degrees of freedom and scale matrix \(\scalmat \) being set to \(I_\p\). 
Since we are merely interested in the graph posterior, we consider the marginal graph likelihood \(f(\y\,|\,\graph)\), obtained after properly integrating out \(\precmat \in \posdefsp_{\graph}\) with respect to the latter prior.
Thus, using the notation from Section~\ref{sec:smc}, we are interested in approximating  \(\graphtarg(\graph) = f(\graph\, |\, \y) \propto f(\y\,|\,\graph)\).
One may show that for any \(\comb{m}\in \combsp{m}, m\in \nn[\p] \) and \( \graph \in \graphsp[\comb{m}],\) we obtain the  clique-separator factorisation form
 \begin{align}\label{eq:graphlikelihood}
    \ungraphtarg[\comb{m}](\graph) = 
    \frac{\prod_{\jtnodeind=1}^k
    \frac{
    J_{\jtnode_{\jtnodeind}}(\scalmat_{\jtnode_{\jtnodeind}} + (\y^\trans \y)_{\jtnode_{\jtnodeind}},\, \df + n)}{
    J_{\jtnode_{\jtnodeind}}(\scalmat_{\jtnode_{\jtnodeind}},\, \df)}
    }
    {\prod_{\jtnodeind=2}^k
    \frac{
    J_{\sep_{\jtnodeind}}(\scalmat_{\sep_{\jtnodeind}}  + (\y^\trans \y)_{\sep_{\jtnodeind}},\, \df+n)
    }{
     J_{\sep_{\jtnodeind}}(\scalmat_{\sep_{\jtnodeind}},\, \df)
    }
    },
 \end{align}
where \( J_{\jtnode_{\jtnodeind}}(\scalmat_{\jtnode_{\jtnodeind}},\, \df)\) is the clique-specific likelihood corresponding to the normalising constant in a hyper Wishart distributions with scale matrix \(\scalmat_{\jtnode_{\jtnodeind}}\) and \(\df_{\jtnode_{\jtnodeind}}\) degrees of freedom and so on; see \emph{e.g.} \cite{dawid1993, Roverato01032000, Atay-Kayis01062005}.


Instead of directly approximating \(f(\graph\,|\,\y)\) using  \smcabbr, we employed the particle Gibbs  sampler derived in full detail in \cite{olsson2019}, where a conditional version of the \smcabbr\, update in Algorithm~\ref{alg:SMC:update} is used as a subroutine.
From the particle Gibbs sampler we may extract a Markov chain of graphs \((\graph[\p]^\ell)_{\ell=1}^\Nmcmc, \graph[\p]^\ell\in \graphsp[\p]\),  having \(\graphtarg\) as its stationary distribution as \(\Nmcmc\rightarrow \infty\). 
An estimator of \(\graphtarg\) is thus obtained after a burn-in period,
as
\begin{align*}
    \graphtargMCMC{\Nmcmc}(\graph) \eqdef \frac{1}{\Nmcmc}
    \sum_{\ell = 1}^\Nmcmc  \1_{\graph}
    (\graph[p]^\ell)), 
    
    
\end{align*}
for all \(\graph \in \graphsp[\p]\).


In this example, 10 particle Gibbs trajectories of length $\Nmcmc=1000$ were sampled, of which the first 300 samples were removed as burn-in period.
To optimise the choice of $\subtreecont$ and $\subtreenonempty$, (where optimality we mean obtaining the best mixing properties) we have evaluated every combination of the tuning parameters \((\subtreecont,\subtreenonempty) \in [0.3,0.5,0.8]^2\), and found that the choice $\subtreecont=0.8$ and $\subtreenonempty=0.5$ we were able to obtain a heatmap pattern of the marginal posterior edge probabilities that best mimics the true one.
The heatmaps together with the maximum aposteriori graph estimates most similar to the true graph are displayed in Figure~\ref{fig:AR1-5_est}.

\begin{figure}
    \centering
    \begin{minipage}[b]{0.4\textwidth}
        \includegraphics[width=\textwidth]{figures/randomgraph/agraph.eps}
    \end{minipage}
    ~
    \begin{minipage}[b]{0.38\textwidth}
        \includegraphics[width=\textwidth]{figures/randomgraph/adjmatcrop.eps}
    \end{minipage}
    \caption{A decomposable graph on $\p = 30$ nodes along with its adjacency matrix.}
    \label{fig:AR1-5}
\end{figure}
\begin{figure}
    \centering
    \begin{minipage}[t]{0.38\textwidth}
        
        \includegraphics[width=\textwidth]{{figures/randomgraph/pgibbs/pgibbs_graph_trajectory_ggm_posterior_n_50_p_30_prior_scale_1.0_shape_x_length_1000_N_50_alpha_0.5_beta_0.5_radius_30_map_0_crop}.eps}
    \end{minipage}
    ~
    \begin{minipage}[b]{0.45\textwidth}
        
        \includegraphics[width=\textwidth,clip]{{figures/randomgraph/pgibbs/pgibbs_graph_trajectory_ggm_posterior_n_50_p_30_prior_scale_1.0_shape_x_length_1000_N_50_alpha_0.5_beta_0.5_radius_30_heatmap_cbar_0_burnin_0_crop}.eps}
    \end{minipage}

    \caption{MAP graph estimate (left) and heatmap (right) for the particle Gibbs sampler.}
    \label{fig:AR1-5_est}
\end{figure}

\subsection{Estimating  $|\graphsp[\p]|$ }
The main motivating factor for using decomposable graphical models is that the clique-specific marginal likelihoods in \eqref{eq:graphlikelihood} have closed form expressions obtained from known distributions.
However, this restriction also ignores the vast majority of graphs, some of which might serve better for expressing the true underlying independence structure, even though recent results by \citet{Fitch20122430} shows that the assumption are anyway reasonable, showing that, asymptotically the posterior will converge to graphical structures that are so-called decomposable covers of the true graph.

Even though the restriction to decomposable graphs is a serious reduction of the graph space, the space is still immense. 
For $\p\le 12$, the cardinality $|\graphsp[\p]|$ can be calculated exactly using the expression provided in \citep{countchordal}.
For larger $\p$, the exact number is intractable and can be regarded as unknown. 
For very large values of $\p $, an expression for the asymptotic behavior is derived in \citep{countchordal} as \(|\graphsp[\p]| \sim \sum_{r=1}^{\p}{\p \choose r} 2^{r(\p-r)}\).

In this example we provide unbiased estimates of $|\graphsp[\p]|$ and compare them to the space of all undirected graphs with the same number of vertices using the \smcabbr\, construction provided in Section~\ref{sec:smc}. 
Specifically, we consider Algorithm~\ref{alg:SMC:update} for approximating $\graphtarg =\unif[{\graphsp[\p]}]$, with for any \(\comb{\jtind}\in \combsp{\jtind}, \jtind\in \nn[\p] \) and \(\graph \in \graphsp[\comb{\jtind}]\),
\begin{align*}
    \ungraphtarg[\comb{\jtind}](\graph) = 1.  
\end{align*}
With this formulation, an unbiased estimate \(\widehat{|\graphsp[\p]|}\) 
of $|\graphsp[\p]|$, 
is obtained directly using the estimator \eqref{eq:constest}.

Table~\ref{tab:numpropdec} shows means and standard errors from 10 
 \smcabbr\, estimates  for $\p \in \nn[15]$.
The upper panel of the table shows  $\widehat{|\graphsp[\p]|}$ while the lower panel shows $\widehat{|\graphsp[\p]|}/2^{{\p(\p-1)/ 2}}$, i.e. estimates of the fraction of undirected graphs which are decomposable.
For $\p\le 12$ the exact enumerations are given in the second column. 
To obtain the estimates with minimal empirical variance, we ran the  \smcabbr\, sampler on the grid of tuning parameters \((\subtreecont,\subtreenonempty) \in [0.3,0.5,0.8]^2\) and found the optimal combination as $\subtreecont=0.5$ and $\subtreenonempty=0.5$. 
The number of particles was set to $\N=10000$.
Figure~\ref{fig:asymtpotics} displays the asymptotic behavior of $|\graphsp[\p]|$ and \( \widehat{|\graphsp[\p]|}\) for $\p\le50$, along with the exact values for $\p \leq 12$, justifying a concordance with the exact results.


\begin{table}[ht!]
\centering
\begin{tabular}{lllll}
$\p $ &\(|\graphsp[\p]|\)  &\(\widehat{|\graphsp[\p]|}\) & std. err. \\
\hline

4       &   61    &   60.6    &   0.72 \\

5       &   822    &   814    &   18.3 \\

6       &   1.82    &   1.78    &   0.48 &$\times 10^{4}$\\

7       &   6.18    &   5.94    &   0.18 &$\times 10^{5}$\\

8       &   3.09    &   2.96    &   0.10 &$\times 10^{7}$\\

9       &   2.19    &   2.09    &   0.12 &$\times 10^{9}$\\

10      &   2.15    &   2.02    &   0.13 &$\times 10^{11}$\\

11      &   2.88    &   2.66    &   0.21 &$\times 10^{13}$\\

12      &   5.17    &   4.77    &   0.42 &$\times 10^{15}$\\

13      &   -    &   1.13    &   0.12 &$\times 10^{18}$\\

14      &   -    &   3.56    &   0.40 &$\times 10^{20}$\\

15      &   -    &   1.46    &   0.17 &$\times 10^{23}$\\
&&&\\
$\p $ & $|\graphsp[\p]|/2^{{\p(\p-1)/ 2}}$  &$\widehat{|\graphsp[\p]|}/2^{{\p(\p-1)/ 2}}$ & std. err. & \\
\hline

4   &   9.53    &   9.47    &   0.11 & $\times 10^{-1}$\\

5   &   8.03    &   7.95    &   0.18 &$\times 10^{-1}$\\

6   &   5.54    &   5.43    &   0.15 &$\times 10^{-1}$\\

7   &   2.95    &   2.83    &   0.08 &$\times 10^{-1}$\\

8   &   1.15    &   1.10    &   0.04 &$\times 10^{-1}$\\

9   &   3.19    &   3.05    &   0.18 &$\times 10^{-2}$\\

10  &   6.12    &   5.73    &   0.36 &$\times 10^{-3}$\\

11  &   7.99    &   7.38    &   0.57 &$\times 10^{-4}$\\

12  &   7.00    &   6.46    &   0.57 &$\times 10^{-5}$\\

13  &   -    &   3.75    &   0.41 & $\times 10^{-6}$\\

14  &   -    &   1.44    &   0.16 & $\times 10^{-7}$\\

15  &   -    &   3.61    &   0.43 &$\times 10^{-9}$\\


 \end{tabular}
\vspace{0.6cm}


\caption{Sequential Monte Carlo estimation of the number decomposable graphs and the fraction of graphs which are decomposable.}
\label{tab:numpropdec}
\end{table}
\begin{figure}[!ht]
    \centering
    \includegraphics[width=0.6\textwidth]{figures/n_chordal_p50_N_10000_T_10}
    \caption{The number of decomposable graphs as a function of the number of vertices.}
    \label{fig:asymtpotics}
\end{figure}
{"config": "arxiv", "file": "1806.00584/simulations.tex"}