# Parallel sampling of decomposable graphs using Markov chain on junction
trees

Mohamad Elmasri Department of Statistical Sciences, University of Toronto, 100
St. George Street, Toronto, ON M5S 3G3, Canada; E-mail:
<EMAIL_ADDRESS>
###### Abstract

Bayesian inference for undirected graphical models is mostly restricted to the
class of decomposable graphs, as they enjoy a rich set of properties making
them amenable to high-dimensional problems. While parameter inference is
straightforward in this setup, inferring the underlying graph is a challenge
driven by the computational difficultly in exploring the space of decomposable
graphs. This work makes two contributions to address this problem. First, we
provide sufficient and necessary conditions for when multi-edge perturbations
maintain decomposability of the graph. Using these, we characterize a simple
class of partitions that efficiently classify all edge perturbations by
whether they maintain decomposability. Second, we propose a new parallel non-
reversible Markov chain Monte Carlo sampler for distributions over junction
tree representations of the graph, where at every step, all edge perturbations
within a partition are executed simultaneously. Through simulations, we
demonstrate the efficiency of our new edge perturbation conditions and class
of partitions. We find that our parallel sampler yields improved mixing
properties in comparison to the single-move variate, and outperforms current
methods. The implementation of our work is available in a Python package.

Keywords — Conditional independence graph; Bayesian structure learning; model
determination; distributed learning.

##  1 Introduction

A graphical model is a collection of joint probability distributions of a
vector of random variables $Y=(Y_{1},\dots,Y_{p})$ that satisfy conditional
independence constraints specified by a graph having a set of vertices
$\left\\{1,\dots,p\right\\}$. A prominent class of graphical models considers
only the case when the underlying graph $G$ is undirected and decomposable
(chordal). Bayesian structure learning (or model determination) is the process
of simultaneously inferring the underlying conditional independence graph and
model parameters, from observed data, through prior specification. In this
work, we develop Markov chain Monte Carlo methods for computational inference
in this setting.

The decomposability assumption is a server restriction on the space of
possible graphs. Less than $8\times 10^{-5}$ of graphs with 12 vertices are
decomposable, and exact enumeration exists only for graphs with up to 15
vertices (Olsson et al., 2022; Wormald, 1985). Bayesian methods requiring
quantification of the prior normalization constant are therefore impractical.
While general Bayesian structure learning methods can take minutes to converge
for medium-sized problems Mohammadi et al. (2021), assuming decomposability
shrinks this convergence time to seconds. In sparse high-dimensional problems,
this inferential efficiency is vital.

The line of work in Frydenberg and Lauritzen (1989); Giudici and Green (1999);
Thomas and Green (2009); Green and Thomas (2013) has culminated in development
of an efficient Metropolis–Hastings (Hastings, 1970) algorithm for
decomposable graphical models. They exploit the  _junction tree_
representation of decomposable graphs to construct efficient proposals that
maintain decomposability across the chain states. The junction tree sampler
proposed in Green and Thomas (2013) is commendably fast, even in large
dimensions. However, it exhibits a high degree of within-sample correlation,
shown in Section 7.4 and in Olsson et al. (2019). The class of multi-edge
perturbations in Green and Thomas (2013) only sufficiently maintains
decomposability. This leads to sub-optimal proposals that unnecessarily
increase the sampler’s rejection rate.

Our first contribution provides the sufficient and necessary conditions under
which multi-edge updates maintain the decomposability of the underlying graph.
We do so by building intuition from a rather simple algorithm that generates
decomposable graphs from random walks on trees (Sec. 3). This algorithm
exhibits properties similar to that of junction trees. We exploit such
properties to characterize a partitioning of junction trees (Sec. 4). That is,
for a given graph vertex, all elements of the tree are partitioned into three
mutually disjoint sets: those that the vertex can connect to, disconnect from,
and everything else. These sets enumerate all possible edge perturbations
pertaining to the vertex that maintain decomposability.

Our second contribution (Sec. 5) proposes a hierarchical latent sampler of
junction trees, in which we project junction trees into higher dimensions
while preserving their desirable factorization properties (Sec. 5.2). Updates
in this latent space do not necessarily translate to graph updates. However,
coupled with our multi-edge partitioning sets, the hierarchical sampler
carries simultaneous parallel updates on the graph. To illustrate the
advantage of parallelism, we first implement a single-move reversible
Metropolis–Hastings Markov chain (Sec. 5.2) and contrast it with a parallel
non-reversible Metropolis–Hastings Markov chain (Sec. 6).

It is now well-understood that non-reversible chains offer significant
advantages over their reversible counterparts (Neal, 1998; Rey-Bellet and
Spiliopoulos, 2015; Bierkens, 2016; Duncan et al., 2016). We illustrate this
advantage in a simulated study (Sec. 7.3) under a Gaussian graphical model. In
this study, the parallel sampler achieves significantly faster convergence and
better mixing properties than the single-move counterpart. To assess the
improvements attributable to our newly developed class of decomposability-
maintaining perturbations, we compare our parallel sampler to Green and Thomas
(2013)’s junction tree sampler (Sec. 7.4). In so doing, we achieve a
significant reduction in within-sample correlation, faster and more stable
convergence, and up to a five-fold increase in the sampler’s acceptance rate
when compared to Green and Thomas (2013)’s. A Python package to accompany this
work is available111https://github.com/melmasri/parallelDG.

Our characterization of the multi-edge partitioning sets and the parallel
sampling approach are likely applicable to different statistical problems,
such as the sampling of directed acyclical graphs and possibly also to
computational graph theory. As with Giudici and Green (1999); Green and Thomas
(2013); Olsson et al. (2019), our latent junction tree sampler can be used
with proper parameter-inference methods in different statistical models. An
important example is the multinomial models for discrete data (Tarantola,
2004).

##  2 Decomposable graphs and junction trees

This work is framed in the classic graphical modelling literature, where
Lauritzen (1996) and Cowell et al. (2006) are excellent references on the
topic.

An _undirected graph_ $G=(\mathcal{V},\mathcal{E})$ is composed of a set of
_vertices_ $\mathcal{V}$ connected by a set of undirected _edges_
$\mathcal{E}$. For any subset $V\in\mathcal{V}$, $G_{V}$ denotes the induced
subgraph with the vertex set $V$. A path $v_{0}\sim v_{k}$ in $G$ consists of
a sequence of vertices $\left(v_{0},v_{1},\dots,v_{k}\right)$, such that
$(v_{i},v_{i+1})\in\mathcal{E}$. A graph is _connected_ when there is a path
between every pair of vertices, and it is _complete_ (a clique) if an edge
exists between every pair of vertices. A complete subgraph is called a
_maximal clique_ if it is not a subgraph of any other clique. A graph is
called a _tree_ if there is a unique path between any pair of nodes.

A graph $G$ is _decomposable_ if and only if the set of maximal cliques of $G$
can be ordered as $(C_{1},\dots,C_{c})$, such that, for each $i=1,\dots,c$, if

$S_{i}=C_{i}\cap\bigcup_{j=1}^{{i-1}}C_{j}\quad\text{then}\quad S_{i}\subset
C_{k},\text{ for some }k<i;$ (2.1)

$S_{i}$ may be empty. The relation in (2.1) is called the running intersection
property, and the sequence $(C_{1},\dots,C_{c})$ is called the _perfect
ordering sequence_. The set
$\texttt{sep}(G)=\left\\{S_{1},\dots,S_{c}\right\\}$, formed by (2.1), is
known as the _minimal separators_ of $G$, and define
$\texttt{cl}(G)=\left\\{C_{1},\dots,C_{c}\right\\}$. Maximal cliques are
unique to a decomposable graph, while separators can repeat in (2.1). The
edge-relation $(C_{i},C_{k})$ that $S_{i}$ forms between $C_{i}$ and $C_{k}$
in (2.1) leads to a tree representation of $G$. A (reduced) _junction tree_ of
$G$ is a tree with a vertex set as the maximal cliques and an edge set as the
minimal separators of $G$ written as $J=(\texttt{cl}(G),\texttt{sep}(G))$. It
follows that $G$ is decomposable if and only if it admits a junction tree
representation. A decomposable graph admits multiple (reduced) junction tree
representations. We denote $\mathcal{J}^{G}$ to be the set of (reduced)
junction trees of $G$. The maximal cardinality search algorithm of Tarjan and
Yannakakis (1984) allows a junction tree representation to be found in time of
order $|\mathcal{V}|+|\mathcal{E}|$.

As in Green and Thomas (2013), we adopt the convention that we allow
separators to be empty, with the effect that every junction tree is connected.
We reserve the terms vertex and edge to refer to the elements of $G$, and use
cliques and separators for the elements of junction trees. To explain this
difference: a vertex represents a single element while a clique, or a
separator, can represent multiple vertices.

Junction trees are general graphical objects defined independently of a
decomposable graph. A tree $T$ with vertices as subsets of $\mathcal{V}$ is
called a _junction tree_ , if it holds that for any pair of vertices $C_{1}$,
$C_{2}$ of $T$, and any tree vertex $C$ in the unique path $C_{1}\sim C_{2}$
in $T$,

$C_{1}\cap C_{2}\subseteq C.$ (2.2)

The relation in (2.2) is called the _junction property_. It is equivalently
stated that for any subset $V\subseteq\mathcal{V}$, the set of vertices in $T$
containing $V$ form a connected subtree. Under this definition, any tree
$T=(\mathcal{C},{\mathcal{S}})$ formed over subsets of cliques $\mathcal{C}$
of $G$ satisfying (2.2) is a junction tree. Given this, it is possible to form
junction trees of a decomposable graph $G$, where the set $\mathcal{C}$
includes, but is not restricted to, the maximal cliques $\texttt{cl}(G)$.

123233234$T$$J=\texttt{reduce}(T)$1243Decomposable graph $G$123234

Figure 1: A decomposable graph $G$, its unique junction tree
$J=(\texttt{cl}(G),\texttt{sep}(G))$ and an expanded tree
$T=(\mathcal{C},{\mathcal{S}})$, where $\texttt{cl}(G)\subset\mathcal{C}$.

The decomposable graph $G$ in Fig. 1 admits the junction tree
$J=(\texttt{cl}(G),\texttt{sep}(G))$, as well the expanded junction tree $T$
that follows from the junction property. The junction property makes it
possible to encode extra, albeit probabilistically superfluous, relations in
$T$. Here, $T$ has four cliques: two maximal $\left\\{1,2,3\right\\}$ and
$\left\\{2,3,4\right\\}$, and two non-maximal, $\left\\{3\right\\}$ and
$\left\\{2,3\right\\}$.

Moving forward, we always refer to junction trees through their junction
property, i.e. having $\texttt{cl}(G)\subseteq\mathcal{C}$, and denote them as
$T$, unless otherwise specified, in which case we refer to them as reduced
junction trees and denote them as $J$. We will refer to elements of junction
trees as cliques and separators, rather than as vertices and edges, as they
represent subsets of vertices. It is always possible to compress a junction
tree $T$ into a reduced junction tree
$J=\texttt{reduce}(T)=(\texttt{cl}(G),\texttt{sep}(G))$, by iteratively
removing cliques in $T$ that are subsets to an adjacent clique, while rewiring
their associated edges to the adjacent superset clique and omitting self-
loops. We denote this operation as $\texttt{reduce}(T)$.

We simplify notations by letting $C\in T$ denote a clique in $T$, with
$(C,C^{\prime})\in T$ as a separator in $T$. We reserve the notation
$\mathcal{C}$ and ${\mathcal{S}}$ to sets of cliques and separators of $T$,
i.e. $T=(\mathcal{C},{\mathcal{S}})$. For a graph vertex $v\in\mathcal{V}$,
$T_{v}$ is the induced subtree of $T$, where every $C\in T_{v}$ includes $v$,
noted as $v\in C$. For $C\in T$, $C\setminus T_{v}$ is the clique $C$
excluding all graph vertices in $T_{v}$. For a specific junction tree $T$, we
let $g(T)$ be the unique decomposable graph represented by $T$, which is
formed by connecting vertices in cliques of $T$. That is to say, by adding an
edge $(v,u)$ to $G$ for every $v,u\in C$, $C\in T$. We let ${\sf
deg}(C,T)=\left|{\sf nei}(C,T)\right|$ as the number of adjacent cliques to
$C$ in $T$. All such notations apply to $J$ as well.

Decomposability allows for the specification of the joint probability
distributions of a random vector through products of marginals over the set
$\texttt{cl}(G)$, known as clique-separator factorization. For a random vector
$Y$ that has a decomposable conditional independence graph $G$, its
distribution factorizes as

$p(Y)=\frac{\prod_{C\in\texttt{cl}(G)}p(Y_{C})}{\prod_{S\in\texttt{sep}(G)}p(Y_{S})},$
(2.3)

where $Y_{A}$ is a subvector of $Y$ indexed by the set $A$.

##  3 Decomposable graphs from random-walks on trees

We start by a simple hierarchical model. Generate an arbitrary tree skeleton
$t$, and conditional on $t$, sample a junction tree $T$ in the hierarchical
schemes

$\displaystyle t$ $\displaystyle\sim\pi(t),$ (3.1) $\displaystyle T{\,|\,}t$
$\displaystyle\sim\pi(T{\,|\,}t),$

while ensuring that the sampling of $T$ adheres to the junction property
(2.2). Here $\pi$ is an underlying probability measure on the space of the
specified random variable. One such sampling method is, given $t$ of $n$
vertices, to initiate $n$ random-walks on $t$. Each walk is initiated at an
arbitrary vertex of $t$, while continuously attempting to visit adjacent
vertices in $t$ of already visited vertices, with each adjacent vertex
attempted only once. The visit-trace of the $i$th random walk is a connected
subtree $t_{i}$ of $t$. The union of all subtrees $\left\\{t_{i}\right\\}$ is
a junction tree of an $n$-vertex decomposable graph. We refer to this
algorithm as an attempt-once random-walk on trees, and detail it in Algorithm
1. We refer to the elements of $t$ as vertices and edges, as they represent
single units.

Input: An arbitrary tree $t$ with $\left\\{1,2,\dots,n\right\\}$ vertices.

1 for  _ $i\leftarrow 1$ to $n$_ do

2 draw $j\sim\mathsf{Uniform}(\left\\{1,2,\dots,n\right\\})$;

$w=\left\\{j\right\\}$; // visited vertices

$a=\left\\{\right\\}$; // attempted vertices

3 while  _${\sf nei}(t_{w},t)\setminus\left\\{w\cup a\right\\}\not=\emptyset$_
do

4 $k\sim\textsf{Uniform}({\sf nei}(t_{w},t)\setminus\left\\{w\cup
a\right\\})$;

5 if  _$\mathsf{Uniform}([0,1])\geq 1/2$_ then add $k$ to $w$ else add $k$ to
$a$;

6

7 $t_{i}$ $\leftarrow$ a copy of $t$, where the label of each vertex $k\in w$
is set to $i$.

return $\bigcup_{i\leq n}t_{i}$ // union over labels.

Algorithm 1 Attempt-once random-walk on trees.

Figure 2 showcases an example from Algorithm 1 with four attempt-once random
walks initiated (independently) on skeleton $t$. The resulting junction tree
$T$, $J=\texttt{reduce}(T)$ and the decomposable graph are in Fig. 1.

123233234$T$$J=\texttt{reduce}(T)$1243Decomposable graph
$G$123234122333324$t$$t_{1}$$t_{2}$$t_{3}$$t_{4}$

Figure 2: Example of a 4-vertex skeleton tree $t$, and 4 attempt-once random-
walks $\left\\{t_{i}\right\\}$, for $i=1,\dots,4$, generated via Algorithm 1;
visited vertices are labelled by the walk’s count (1-4) with solid-edge
subtrees. The resulting junction tree $T$, $\texttt{reduce}(T)$, and
decomposable graph $G$ are in Fig. 1.

From the junction property (2.2), we know that for a fixed junction tree $T$
and graph vertex $v\in\mathcal{V}$, $T_{v}$ is a connected subtree. In fact,
if we remove all vertices in $T_{v}$ except $v$, the result is the subtree
skeleton $t_{v}$, as
$t_{v}=T_{v}\setminus\left\\{{\mathcal{V}\setminus\left\\{v\right\\}}\right\\}$.
Thus, updating $t_{v}$ as done in Algorithm 1 ensures that $T$ is a junction
tree, as confirmed in the following theorem.

######  Theorem 3.1.

A tree $T$ generated by Algorithm 1 is a junction tree.

###### Proof.

With a direct proof, every $C,C^{\prime}\in T$ are connected via a unique path
$C\sim C^{\prime}$. If $v\in C\cap C^{\prime}$, then $C,C^{\prime}\in T_{v}$,
and so is every clique on the path $C\sim C^{\prime}$. ∎

Algorithm 1 is both simple and quite revealing of how decomposable graphs can
be constructed. It illustrates three essential properties that motivate our
work: $(i)$ updating a junction tree can be done solely through vertex-induced
subtrees $\left\\{T_{v}\right\\}$ for $v\in\mathcal{V}$; when conditional
independence is not the focus, $(ii)$ a random decomposable graph can be
constructed from the union of independently generated subtrees
$\left\\{t_{i}\right\\}$ (shown in Fig. 2). When conditional independence is
relevant, $(iii)$ each subtree $T_{v}$ can be expanded and contracted in
multiple directions in parallel by expanding $t_{v}$; that is to say, $(iiia)$
expanding $T_{v}$ can be carried out by expanding the $v$th random-walk,
$t_{v}$, (in parallel) to adjacent vertices in $t$ (lines 5-7 in Alg. 1); and
$(iiib)$ contracting $T_{v}$ is a rollback of expansion, by contracting
$t_{v}$ from its leaves.

Observations $(i)-(iii)$ are the bedrock of our analysis, as they enable a
parallelized method to work with junction trees through their junction
property. The next section connects the expansion and contraction methods
described in $(iiia)-(iiib)$ to decomposable graph updates.

##  4 Decomposable graph updates

For a given junction tree $T=(\mathcal{C},{\mathcal{S}})$, and graph vertex
$v\in\mathcal{V}$, we partition $\mathcal{C}$ into three clique sets: cliques
that $v$ can be added (_connect_) to, referred to as $\mathcal{N}_{{v}}$,
cliques that $v$ can be removed from (_disconnect_), referred to as
$\mathcal{L}_{{v}}$, and everything else. Graphically, $\mathcal{N}_{{v}}$ is
composed of cliques in $T$ that are adjacent to the subtree $T_{v}$, and
$\mathcal{L}_{{v}}$ are the leaf cliques of $T_{v}$, as

$\displaystyle\mathcal{N}_{{v}}(T):=\left\\{C\in T:(C,C^{\prime})\in
T\,,C^{\prime}\in T_{v}\,,C\not\in T_{v}\right\\},$ (4.1)
$\displaystyle\mathcal{L}_{{v}}(T):=\left\\{C\in T_{v}:{\sf
deg}(C,T_{v})=1\right\\}.$

The terms add and remove are used when referring to updates on $T$, and
(dis)connect when referring to updates on the underlying graph $G$. In $T$,
adding vertex $v$ to an adjacent clique $C\in\mathcal{N}_{{v}}$ implies
updating $C$ to $C\cup\left\\{v\right\\}$. Removing $v$ from a leaf clique
$C\in\mathcal{L}_{{v}}$ implies updating $C$ to
$C\setminus\left\\{v\right\\}$. Although we work solely on $T$, in $G=g(T)$,
these operations translate to the following. Adding $v$ to $C$ entails
connecting edges $(v,u)$ in $G$, where $u\in C\setminus C_{{\sf adj}}$; here,
$C_{{\sf adj}}$ is the unique clique in $T_{v}$ that is adjacent to $C$ (i.e.
$(C,C_{{\sf adj}})\in T$). Removing $v$ from a clique $C\in\mathcal{L}_{{v}}$
entails disconnecting all graph edges $(v,u)$ from $u\in C\setminus C_{{\sf
adj}}$, where $(C,C_{{\sf adj}})\in T_{v}$. In both cases, there exists a
unique clique $C_{{\sf adj}}\in T_{v}$ that is adjacent to the clique of
interest. To gain better understanding about those sets, consider the
following graphical example.

123233234$T$$J=\texttt{reduce}(T)$1243Decomposable graph
$G$123234122333324$t$$t_{1}$$t_{2}$$t_{3}$$t_{4}$9876123451110$C_{1}$$C_{2}$$C_{3}$$C_{4}$$C_{5}$$C_{6}$$C_{7}$12$C_{8}$

Figure 3: 12-vertex decomposable graph of 8 maximal cliques
$\texttt{cl}(G)=\left\\{C_{1},\dots,C_{8}\right\\}$.

######  Example 4.1.

Figure 3 illustrates a 12-vertex, 8-maximal-clique
$\left\\{C_{1},\dots,C_{8}\right\\}$, decomposable graph having a junction
tree $T$ in Fig. 4 (left). The tree $T$ includes one clique that have no graph
vertices, hence shown as $\emptyset$, and two non-empty non-maximal cliques
labeled as $(8,7)$ and $3$ to indicate the vertices included in them. The
induction of $T$ to vertex $6$ ($T_{6}$), shown in Fig. 4 (right), is composed
of 3 cliques $\left\\{C_{2},C_{3},C_{4}\right\\}$ connected with solid edges.
The set $\mathcal{N}_{{6}}(T)$ is coloured in blue and $\mathcal{L}_{{6}}(T)$
in green.

123233234$T$$J=\texttt{reduce}(T)$1243Decomposable graph
$G$123234122333324$t$$t_{1}$$t_{2}$$t_{3}$$t_{4}$9876123451110$C_{1}$$C_{2}$$C_{3}$$C_{4}$$C_{5}$$C_{6}$$C_{7}$12$C_{8}$$C_{1}$$C_{2}$$C_{3}$$C_{4}$$C_{5}$$C_{6}$$C_{7}$$C_{8}$$\emptyset$$8,7$$T$3

123233234$T$$J=\texttt{reduce}(T)$1243Decomposable graph
$G$123234122333324$t$$t_{1}$$t_{2}$$t_{3}$$t_{4}$9876123451110$C_{1}$$C_{2}$$C_{3}$$C_{4}$$C_{5}$$C_{6}$$C_{7}$12$C_{8}$$C_{1}$$C_{2}$$C_{3}$$C_{4}$$C_{5}$$C_{6}$$C_{7}$$C_{8}$$\emptyset$$8,7$$T$3$C_{1}$$C_{2}$$C_{3}$$C_{4}$$C_{5}$$C_{6}$$C_{7}$$C_{8}$$\emptyset$$8,7$$T_{6}$3

Figure 4: A junction tree $T$ (left) of decomposable graph in Fig. 3, its
induction by vertex $6$ (right) in solid lines, with neighbouring
$\mathcal{N}_{{6}}$ (blue) and leaf $\mathcal{L}_{{6}}$ (green) cliques. None-
maximal cliques are labeled with vertices contained in them, as
$\emptyset,\left\\{8,7\right\\}$ and $\left\\{3\right\\}$.

Cliques in $\mathcal{N}_{{6}}$ and $\mathcal{L}_{{6}}$ are nonadjacent in $T$,
except $C_{4}\in\mathcal{L}_{{6}}$, which is adjacent to
$\emptyset,C_{5},C_{8}\in\mathcal{N}_{{6}}$. Removing $6$ from $C_{4}$ would
hence remove $\emptyset,C_{5}$ and $C_{8}$ from $\mathcal{N}_{{6}}$.
Similarly, adding $6$ to any of $\emptyset,C_{5}$ or $C_{8}$, would remove
$C_{4}$ from $\mathcal{L}_{{6}}$. Therefore, we term a clique
$C\in\mathcal{N}_{{v}}\cup\mathcal{L}_{{v}}$ a partition-divisor if it is
adjacent to other cliques in those sets. A partition-divisor occurs only when
a leaf clique $C\in\mathcal{L}_{{v}}$ is adjacent to a clique in
$\mathcal{N}_{{v}}$.

As implicitly mentioned above, the sets $\mathcal{N}_{{v}}$ and
$\mathcal{L}_{{v}}$ exhibit a congruence relation across updates, where leaf
cliques become neighbouring cliques, and vice versa, if updated. In Example
4.1, adding vertex $6$ to the neighbouring clique
$C_{5}\in\mathcal{N}_{{6}}(T)$ to form the clique
$C^{\prime}_{5}=C_{5}\cup\left\\{6\right\\}$ in a new junction $T^{\prime}$
renders $C^{\prime}_{5}\in\mathcal{L}_{{6}}(T^{\prime})$. Similarly, if 6 is
removed from $C_{4}$, the new $C^{\prime}_{4}$ becomes a neighbouring clique
to $T^{\prime}_{6}$. This congruence relation, demonstrated in (4.2), will
come in handy in carrying graph updates across a Markov chain, as shown in
subsequent sections.

$C\in\mathcal{N}_{{v}}\Longleftrightarrow
C\cup\left\\{v\right\\}\in\mathcal{L}_{{v}}.$ (4.2)

It is clear, from Fig. 4, that $\mathcal{L}_{{v}}$ and $\mathcal{N}_{{v}}$ are
both vertex-wise symmetric with respect to graph updates. That is to say, for
an unconnected vertex pair $v,u\in\mathcal{V}$, $u\in C$ for some
$C\in\mathcal{N}_{{v}}$ if and only if $v\in C^{\prime}$ for some
$C^{\prime}\in\mathcal{N}_{{u}}$. If $v$ and $u$ are in a single clique of
$G$, then $u\in C$ for some $C\in\mathcal{L}_{{v}}$ if and only if $v\in
C^{\prime}$ for some $C^{\prime}\in\mathcal{L}_{{u}}$. These findings match
the single-edge decomposability-preserving conditions illustrated in
Frydenberg and Lauritzen (1989) and Giudici and Green (1999).

Extension to a multi-vertex set $V\in\mathcal{V}$, where $V$ is a clique in
$G$, follows directly by using $T_{V}$ in (4.1). For improved readability, the
partition sets (4.1) are stated for singleton vertices.

Now that we have analytically and visually demonstrated the class of
partitions generated by (4.1), the following theorems illustrate the necessary
and sufficient in order for a perturbation of the edge set of a decomposable
graph to maintain decomposability.

######  Theorem 4.2.

Let $G=(\mathcal{V},\mathcal{E})$ be a decomposable graph, and let $V,U$ be
two disjoint subsets of $\mathcal{V}$, each forming complete graphs in $G$.
Adding edges from every vertex in $V$ to every vertex in $U$ results in a
decomposable graph $G^{\prime}=(\mathcal{V},\mathcal{E}^{\prime})$ if and only
if $U$ can be partitioned into mutually exclusive subsets of vertices ordered
as $\left(U^{(1)},\dots U^{(K)}\right)$, where $U^{(i)}\cap U^{(j)}=\emptyset$
if $i\neq j$. Such that, for some $J\in\mathcal{J}^{G}$, there exists path of
maximal cliques $\left(C_{1},\dots,C_{K}\right)$ in $J$ satisfying the
recursive relation that $U^{(1)}\subseteq C_{1}$ for
$C_{1}\in\mathcal{N}_{{V}}(J)$, and $U^{{(k)}}\subseteq C_{k}$ for
$C_{k}\in\mathcal{N}_{{U^{(k-1)}}}(J)$ where $k=2,\dots,K$.

######  Theorem 4.3.

Let $G=(\mathcal{V},\mathcal{E})$ be a decomposable graph, and let $V,U$ be
two disjoint subsets of $\mathcal{V}$ that are completely connected, i.e.
$V\cup U$ is complete in $G$. Disconnecting all edges $(v,u)$, where $v\in V$
and $u\in U$, results in a decomposable graph
$G^{\prime}=(\mathcal{V},\mathcal{E}^{\prime})$ if and only if $U$ can be
partitioned into mutually exclusive subsets of vertices ordered as
$\left(U^{(1)},\dots U^{(K)}\right)$, where $U^{(i)}\cap U^{(j)}=\emptyset$ if
$i\neq j$. Such that, for some $J\in\mathcal{J}^{G}$, there exists a path of
maximal cliques $\left(C_{1},\dots,C_{K}\right)$ in $J$ satisfying the
recursive relation that $U^{(1)}\subseteq C_{1}$ for
$C_{1}\in\mathcal{N}_{{V}}(J)$, and $U^{{(k)}}\subseteq C_{k}$ for
$C_{k}\in\mathcal{N}_{{U^{(k-1)}}}(J)$ where $k=2,\dots,K$.

A special case in Theorem 4.2 and 4.3 is when the whole vertex set $U$ is a
subset of $C_{1}$, meaning $U=U^{(1)}$ and $\left\\{U^{(k)}\right\\}$ are
empty, for $k=2,\dots,K$. This implies that $C_{1}$ is a neighbouring clique
($C_{1}\in\mathcal{N}_{{V}}(J)$) in Theorem 4.2, or a leaf clique
($C_{1}\in\mathcal{L}_{{V}}(J)$) in Theorem 4.3. When the sets
$\left\\{U^{(k)}\right\\}$, for $k=1,\dots,K$, are non-empty, this implies
that the set $U$ is spread out over multiple adjacent cliques in $J$ that form
a path. Metaphorically, a disconnect move can be viewed as an act of peeling a
bandage, where $V$ is disconnected sequentially from $U^{(1)}$, then
$U^{(2)}$, and so on. Here, $U^{(k)}$ is always contained in a leaf clique of
the junction tree formed by disconnecting $V$ from $\cup_{j<k}U^{(j)}$. A
connect move is the reverse operation, where $U^{(k)}$ is always contained in
a neighbouring clique to the induced junction tree $J_{V}$ formed by
connecting $V$ to $\cup_{j<k}U^{(j)}$. These operations follow directly from
the junction property (2.2).

The main consequence of Theorem 4.2 and 4.3 is that they enumerate all
decomposability-preserving updates pertaining to a vertex in terms of
partition sets (4.1). Theorem 4.2 and 4.3 are stated in terms of graph
updates, since it is the relevant structure in graphical models. Similar
statements that govern updates on general junction trees follow directly from
the junction properly (2.2). The next sections begin by proposing a single-
move sampler, and then generalizing it to a parallel sampler over junction
trees.

##  5 Single-move sampler

###  5.1 Proposal probability

This section presents a single-move sampler based on partition sets (4.1). For
a state of the chain $T$, updating clique $C$, to result in state
$T^{\prime}$, occurs with probability

$q(T,T^{\prime})=\frac{1}{2}\frac{1}{|\mathcal{V}|}\frac{1}{\left|\mathcal{P}_{{v}}(T)\right|},$
(5.1)

where the ratio $1/\left|\mathcal{P}_{{v}}(T)\right|$ accounts for the
probability of selecting $C$ out of the set $\mathcal{P}_{{v}}$, where
$\mathcal{P}_{{v}}=\mathcal{N}_{{v}}$ for the addition, or
$\mathcal{P}_{{v}}=\mathcal{L}_{{v}}$ for the removal move. The ratio $1/2$
accounts for the update-type choice, and $1/\left|\mathcal{V}\right|$ accounts
for the vertex choice.

For the addition move, the reverse proposal requires enumerating
$\mathcal{L}_{{v}}(T^{\prime})$ of the new junction tree $T^{\prime}$, which
can be computed from $\mathcal{L}_{{v}}(T)$, as a result of the congruence
relation between the sets (4.1) defined in (4.2). If $v$ is added to
$C\in\mathcal{N}_{{v}}(T)$ to form $C^{\prime}$, it is clear that
$C^{\prime}\in\mathcal{L}_{{v}}(T^{\prime})$ as a leaf clique of $T^{\prime}$.
No other clique in $\mathcal{L}_{{v}}(T)$ changes status, except if the clique
adjacent to $C$ in $T_{v}$ is a partition-divisor. In Example 4.1, if $6$ is
added to $C_{1}$, the number of leaves of $T^{\prime}_{6}$ is one greater than
$T_{6}$, since $C_{2}$ is not a partition-divisor. However, if $6$ is added to
$C_{5}$, or any neighbour of the partition-divisor $C_{4}$, $C_{4}$ would lose
its status as a leaf clique in $T^{\prime}_{6}$, resulting in
$\mathcal{L}_{{6}}(T^{\prime})=\mathcal{L}_{{6}}(T)$. Let $C_{{\sf adj}}$ be
the unique clique in $T_{v}$ that is adjacent to $C$, $(C_{{\sf adj}},C)\in
T$, then

$\left|\mathcal{L}_{{v}}(T^{\prime})\right|=\left|\mathcal{L}_{{v}}(T)\right|+{\bm{1}}\\{C_{{\sf
adj}}\not\in\mathcal{L}_{{v}}(T)\\}.$ (5.2)

For the removal move, the reverse proposal requires enumerating
$\mathcal{N}_{{v}}(T^{\prime})$, which, similarly, can be computed from
$\mathcal{N}_{{v}}(T)$. If $v$ is removed from $C$, and $C$ is a partition-
divisor, all neighbours of $C$ in $\mathcal{N}_{{v}}(T)$ lose their status as
neighbouring cliques to $T_{v}$. For example, if $6$ is removed from $C_{4}$
(in Fig. 4), all neighbours (in blue) of $T_{6}$
($\left\\{\emptyset,C_{5},C_{8}\right\\}$) are not in
$\mathcal{L}_{{6}}(T^{\prime})$, and $C_{4}$ would become a leaf node of
$T^{\prime}_{6}$. Otherwise,
$\left|\mathcal{N}_{{v}}(T^{\prime})\right|=\left|\mathcal{N}_{{v}}(T)\right|+1$.
This results in

$\left|\mathcal{N}_{{v}}(T^{\prime})\right|=\left|\mathcal{N}_{{v}}(T)\right|+2-{\sf
deg}(C,T).$ (5.3)

In all cases, quantifying the reverse proposal ratio $q(T^{\prime},T)$ only
requires the enumeration of the sets $\mathcal{L}_{{v}}(T)$ and
$\mathcal{N}_{{v}}(T)$. This can be carried out efficiently with a single pass
over $T_{v}$, and its neighbours, without an actual modification to $T$.

###  5.2 Prior and posterior probabilities

Following (3.1), we propose a hierarchical sampler that samples a junction
tree $T$, and the underlying skeleton $t$, iteratively from the posterior
$\pi(T{\,|\,}t)$ and $\pi(t{\,|\,}T)$, respectively. Assuming a skeleton tree
of $n$ vertices (since a decomposable graph of $n$ vertices has at most $n$
maximal cliques), from Cayley (1889) we know that there is $n^{n-2}$ possible
trees of size $n$. A uniform prior on $t$ would be $\pi(t)\propto n^{-(n-2)}$.

Turning to prior specification over the space of junction trees. Byrne et al.
(2015); Green and Thomas (2018) have illustrated an important family of graph
laws through algebraic characterization, known as the  _clique-separator_
factorization laws. Laws in this family have densities $\pi(T)$ that factorize
as

$\pi(T)\propto\frac{\prod_{C\in\mathcal{C}}\phi(C)}{\prod_{S\in{\mathcal{S}}}\psi(S)},$
(5.4)

over the class of decomposable graphs, and consequently over the class of
junction trees, for positive functions $\phi,\psi$. As noted by Green and
Thomas (2018, Thm. 1), a graph law over the class of decomposable graphs,
whose support is all such class, is weakly structural Markov if and only if it
is a clique-separator factorization law. Specifying decomposability through
junction trees does not alter the clique-separator factorization of joint
distribution defined in (2.3), as stated in the following proposition.

######  Proposition 5.1.

For a random vector $Y$ that has a decomposable conditional independence graph
$G$, where $T=(\mathcal{C},{\mathcal{S}})$ is a junction tree of $G$ having
$\texttt{cl}(G)\subseteq\mathcal{C}$ and
$\texttt{sep}(G)\subseteq{\mathcal{S}}$, the distribution of $Y$ factorizes as

$p(Y)=\frac{\prod_{C\in\mathcal{C}}p(Y_{C})}{\prod_{S\in{\mathcal{S}}}p(Y_{S})}=\frac{\prod_{C\in\texttt{cl}(G)}p(Y_{C})}{\prod_{S\in\texttt{sep}(G)}p(Y_{S})}.$
(5.5)

###### Proof.

Let $C\in\mathcal{C}$ be any non-maximal clique in $T$, i.e. there exists a
$(C^{\prime},C)\in{\mathcal{S}}$ such that $C\subseteq C^{\prime}$. The terms
in (5.5) that are specific to $C$ and $C^{\prime}$ are

$\frac{p(Y_{C})p(Y_{C^{\prime}})}{p(Y_{C\cap
C^{\prime}})}=\frac{p(Y_{C})p(Y_{C^{\prime}})}{p(Y_{C})}=p(Y_{C^{\prime}}).$

Carrying such cancellations over non-maximal cliques results in the second
equality in (5.5). ∎

Following a similar argument to Proposition 5.1, the density $\pi(T)$ in (5.4)
factorizes in terms of $(\texttt{cl}(G),\texttt{sep}(G))$, where $G=g(T)$, if
$\psi=\phi$, as

$\pi(T)\propto\frac{\prod_{C\in\mathcal{C}}\phi(C)}{\prod_{S\in{\mathcal{S}}}\psi(S)}=\frac{\prod_{C\in\texttt{cl}(G)}\phi(C)}{\prod_{S\in\texttt{sep}(G)}\psi(S)}.$
(5.6)

Otherwise, a multiplicative factor of
$\prod_{C\in\left\\{\mathcal{C}\setminus\texttt{cl}(G)\right\\}}\phi(C)/\prod_{S\in\left\\{{\mathcal{S}}\setminus\texttt{sep}(G)\right\\}}\psi(S)$
exists on the right-hand side of (5.6).

Turning to posterior specification, we have $\pi(t{\,|\,}T)\propto 1$, as the
likelihood and prior terms cancel out, leaving a symmetric proposal, which
follows exactly as the tree randomization step performed in Green and Thomas
(2013), detailed in (Thomas and Green, 2009, Sec. 5).

For a posterior sample from $\pi(T{\,|\,}t)$, we adopt a standard
Metropolis–Hastings step (Hastings, 1970). Given a junction tree $T$, we
uniformly sample a vertex $v\in\mathcal{V}$, compute the partitions
$\mathcal{L}_{{v}}(T)$ and $\mathcal{N}_{{v}}(T)$, and uniformly sample an
update-type and a partition $C$ to update. Let $C_{{\sf adj}}$ be the unique
clique in $T_{v}$ that is adjacent to $C$ (i.e. $(C_{{\sf adj}},C)\in T$). The
acceptance probability of the proposal of updating $C$ to form $C^{\prime}$
and junction tree $T^{\prime}$ is

$\displaystyle\alpha(T,T^{\prime}{\,|\,}t)$
$\displaystyle=\min\left\\{1,\frac{\pi(T^{\prime}{\,|\,}t)q(T,T^{\prime})}{\pi(T{\,|\,}t)q(T^{\prime},T)}\right\\}$
(5.7) $\displaystyle=\min\left\\{1,\frac{p(Y_{C^{\prime}})p(Y_{C\cap C_{{\sf
adj}}})}{p(Y_{C^{\prime}\cap C_{{\sf
adj}}})p(Y_{C})}\frac{\phi(C^{\prime})\psi(C\cap C_{{\sf
adj}})}{\psi(C^{\prime}\cap C_{{\sf
adj}})\phi(C)}\frac{\left|\mathcal{P}_{{v}}(T^{\prime})\right|}{\left|\mathcal{P}_{{v}}(T)\right|}\right\\}.$

If the update-type is an addition, $\mathcal{P}_{{v}}(T)=\mathcal{N}_{{v}}(T)$
and $\mathcal{P}_{{v}}(T^{\prime})=\mathcal{L}_{{v}}(T^{\prime})$. Otherwise,
$\mathcal{P}_{{v}}(T)=\mathcal{L}_{{v}}(T)$ and
$\mathcal{P}_{{v}}(T^{\prime})=\mathcal{N}_{{v}}(T^{\prime})$. Both quantities
can be computed by (5.2) and (5.3).

The step in (5.7) is reversible by the congruence relation (4.2), across the
chain states. Moreover, $\pi(T{\,|\,}t)>0$. Therefore, we can always find a
path of junction trees from any state of the chain to the single-clique
junction tree. Since the state space of the chain is finite, it follows from
irreducibility that this chain is ergodic.

Computing the acceptance probability in (5.7), while not requiring a
modification of $T$ in practice, still requires enumeration of the two state-
dependent sets (4.1). Section 6 proposes an alternative junction tree sampler
with a state-invariant proposal that only requires enumeration of a single
partition set for every proposal.

##  6 Parallel sampler

###  6.1 Proposal probability

Partition sets (4.1) are mutually independent, i.e. , for a specific vertex
$v\in\mathcal{V}$, a clique in $T$ belongs either to the leaf set, the
neighbouring set, or neither. This allows us to perform all proposals within a
partition set in parallel. For example, all proposals generated by the blue
cliques in Fig. 4 can be carried out simultaneously, without interference,
i.e., in a single Markov chain step, $T$ can be updated into one of $2^{5}$
possible states ($\left|\mathcal{N}_{{6}}\right|=5$). The choice of which
depends on the acceptance probability of each proposal, independently. Such
localization of updates renders the clique choice in the proposal probability
of (5.1) unnecessary, since all updates are carried out. This leads to a
proposal probability that depends only on the vertex choice and update type,
as

$\frac{1}{2}\frac{1}{|\mathcal{V}|}.$ (6.1)

By eliminating the need to account for which clique to update, we only need to
enumerate a single partition set, $\mathcal{N}_{{v}}$ or $\mathcal{L}_{{v}}$,
depending on the chosen update type. Contrarily, the single-move sampler (Sec.
5) required the enumeration of both sets.

###  6.2 Prior and posterior probabilities

We follow the iterative sampling method and prior setup of the single-move
sampler (Sec. 5), leading to $\pi(t{\,|\,}T)\propto 1$. We adopt a
parallelized Metropolis–Hastings (Hastings, 1970) step to sample from the
posterior of $\pi(T{\,|\,}t)$. Given a junction tree $T$, we uniformly sample
a vertex $v\in\mathcal{V}$, and an update-type. Let $\mathcal{P}_{{v}}(T)$ be
the partition set corresponding to the update-type. For every
$C_{i}\in\mathcal{P}_{{v}}(T)$, let $C_{{\sf adj}(i)}\in T_{v}$ be the unique
clique adjacent to $C_{i}$, such that $(C_{{\sf adj}(i)},C_{i})\in T$. Let
$T^{\prime}_{i}$ be the junction tree formed by updating clique $C_{i}$ to
$C^{\prime}_{i}$. Accept every $C^{\prime}_{i}$, with a standard
Metropolis–Hastings acceptance probability $\alpha_{i}$, as

$\displaystyle\alpha_{i}(C_{i},C^{\prime}_{i}{\,|\,}t)$
$\displaystyle=\min\left\\{1,\frac{\pi(T^{\prime}_{i}{\,|\,}t)}{\pi(T{\,|\,}t)}\right\\}$
(6.2) $\displaystyle=\min\left\\{1,\frac{p(Y_{C^{\prime}_{i}})p(Y_{C_{{\sf
adj}(i)}\cap C_{i}})}{p(Y_{C_{{\sf adj}(i)}\cap
C^{\prime}_{i}})p(Y_{C_{i}})}\frac{\phi(C^{\prime}_{i})\psi(C_{{\sf
adj}(i)}\cap C_{i})}{\psi(C_{{\sf adj}(i)}\cap
C^{\prime}_{i})\phi(C_{i})}\right\\},$

where $\phi$ and $\psi$ as in (5.4). Although the first equality in (6.2) is
defined in terms of $T$, the final ratio only involves terms associated with
$C_{i}$, its update $C^{\prime}_{i}$, and its neighbour $C_{{\sf adj}(i)}$.
Since we do not update the skeleton $t$ after every update, all updates in
$\mathcal{P}_{{v}}(T)$ are carried out simultaneously. The proposed parallel
sampling scheme is straightforward to implement in practice, with some
computational considerations discussed in the next section.

The chain constructed in (6.2) is irreducible and aperiodic over a finite
state space, since it is possible to reach the single-clique junction tree
with a finite number of steps, from any state, and backward, and return to the
single-clique junction tree in any number of steps. Moreover,
$\pi(T_{i}{\,|\,}t)>0$. Hence, a unique stationary distribution exists.
However, the chain is not reversible, as it only satisfies the partial detail
balance equations (Whittle, 1985).

Conditional on the skeleton $t$, the marginal state space of clique $C_{i}$ is
reversible by the congruence relation in (4.2) across the chain states, where
leaf cliques become neighbouring cliques, and vice versa, if updated. Hence,
the congruence relation tells us that all updated cliques in
$\mathcal{P}_{{v}}(T)$ are now a subset of $\mathcal{P}_{{v}}(T^{\prime})$.
However, the latter also contains other cliques. Updating all cliques in
$\mathcal{P}_{{v}}(T^{\prime})$ in a reversal move would result in a state
$T^{\prime\prime}$ that might not be the starting state $T$. Hence,
reversibility is partially guaranteed under the set $\mathcal{P}_{{v}}(T)$, or
when
$\left\\{C\cup\left\\{v\right\\}:C\in\mathcal{P}_{{v}}(T)\right\\}=\mathcal{P}_{{v}}(T^{\prime})$
in an addition move and
$\left\\{C\setminus\left\\{v\right\\}:C\in\mathcal{P}_{{v}}(T)\right\\}=\mathcal{P}_{{v}}(T^{\prime})$
in a removal move.

###  6.3 Computational considerations

Our proposed single-move and parallel samplers rely on the junction property
of trees (Sec. 2), where in practice, $T=(\mathcal{C},{\mathcal{S}})$ can be
much larger than $\texttt{reduce}(T)=(\texttt{cl}(G),\texttt{sep}(G))$,
$G=g(T)$. Two aspects influence the number of non-empty cliques of
$\mathcal{C}$, the number of vertices of the underlying skeleton $t$, and the
prior over $T$, in (5.4).

The proposed samplers can be initialized with a skeleton $t$ of an arbitrary
number of vertices $n$. This can lead $\mathcal{C}$ to have as much as $n-1$
non-empty cliques, since the complete graph is a single clique. Very large $n$
can lead to a slowdown in convergence, since many updates can be redundant in
the graph space. Similarly, conservative $n$ can lead to very few accepted
updates. We find it intuitive to initialize $t$ with $n=|\mathcal{V}|$
vertices, since the no-edge graph has the greatest number of cliques
($\left|\mathcal{V}\right|$) over the vertex set.

Prior specification can control the number of non-empty cliques in
$\mathcal{C}$. Very diffused priors saturate $\mathcal{C}$ with many cliques
as non-empty and non-maximal. Tight priors are more conservative in proposing
moves, which reduces the effective size of $\mathcal{C}$ along the chain
iterations, even when $t$ is large. We formalize this trade-off in the
following example. Consider the $i$th acceptance ratio in (6.2). Proposition
5.1 enables us to factor out likelihood ratios corresponding to non-maximal
cliques, as if the likelihood ratio is computed over
$\texttt{reduce}(T)=(\texttt{cl}(G),\texttt{sep}(G))$. The prior ratio cancels
all terms associated with non-maximal cliques if $\phi=\psi$. Therefore, we
will assume that $\phi\neq\psi$.

If $C\in\mathcal{C}$ is a non-maximal clique, then $C$ is contained in one of
its neighbours in $T$, say $C_{{\sf adj}}$. Therefore, the contribution of $C$
to the prior ratio in the acceptance probability simplifies to the
multiplicative factor $\phi(C)/\psi(C_{{\sf adj}}\cap C)=\phi(C)/\psi(C)$.
Updating $C$ to the proposal $C^{\prime}$ can cause it to change its maximal
status, for example, $C^{\prime}$ being maximal while $C$ is non-maximal, and
vice versa. This leads to four cases, (a)-(d), indicated in Table 1. Each case
in this table indicates the contribution of the non-maximal clique to the
acceptance ratio, in a multiplicative factor. In case (a) for example, when
both $C$ and the proposal $C^{\prime}$ are both maximal cliques, the non-
maximal clique contribution to the acceptance probability ratio is simply 1.
If $C$ is a maximal clique while the proposal $C^{\prime}$ is non-maximal,
then other than the contribution of $C$ to the acceptance probability,
$C^{\prime}$ contributes with the prior ratio
$\psi(C^{\prime})/\phi(C^{\prime})$, as indicated in case (b). In the special
case (d), the acceptance probability is fully driven by the prior ratio, since
no other clique in $T$ would change status.

Table 1: Multiplicative contribution of non-maximal cliques to the acceptance ratio |  $C^{\prime}$ maximal |  $C^{\prime}$ non-maximal
---|---|---
$C$ maximal | (a) : 1 | (b) : $\frac{\phi(C^{\prime})}{\psi(C^{\prime})}$
$C$ non-maximal | (c) : ${\frac{\psi(C)}{{\phi(C)}}}$ | (d) : ${\frac{\phi(C^{\prime})\psi(C)}{\psi(C^{\prime})\phi(C)}}$

Adjusting the contribution factors in Table 1 can be effective in controlling
the number of non-maximal cliques in $\mathcal{C}$. For example, a tight prior
that penalizes separators by strictly specifying $\psi(C)>\phi(C)$ can
encourage proposals of maximal cliques. In this setup, the likelihood of
accepting the non-maximal proposal in case (b) is curtailed, while accepting
the maximal proposal in case (a) is encouraged. Case (d) depends on the prior
specification. For example, under the clique exponential family priors (Bornn
and Caron, 2011; Green and Thomas, 2018), where $\phi(C)=\exp(\alpha|C|)$ and
$\psi(S)=\exp(\beta|S|)$ for constants $\alpha,\beta>0$, the ratio
${\phi(C)}/{{\psi(C)}}=\exp\left\\{|C|(\alpha-\beta)\right\\}$. Hence, case
(d) translates to $\exp\left\\{\xi(\alpha-\beta)\right\\}$, where $\xi=1$ if
$|C^{\prime}|>|C|$ otherwise $-1$. If $\beta>\alpha$, case (d) favours smaller
non-maximal cliques.

##  7 Numerical performance of the new samplers

###  7.1 Introduction

To compare our proposed samplers and alternative models, we focus on the
Gaussian decomposable graphical model setup (Giudici and Green, 1999).
Alternative setups, for example, under a log-linear model, are detailed in
Tarantola (2004) and many others. In Supplementary Material (SM) Section A, we
derive the conjugate posterior of a Gaussian graphical model under a Wishart
prior for covariance matrix. This is a classic derivation, where the
covariance matrix is integrated out, resulting in a posterior formed in terms
of the graph priors and Wishart hyperparameters (Dawid and Lauritzen, 1993;
Giudici and Green, 1999).

Under the Gaussian decomposable graphical model, in Section 7.3, we compare
the single-move sampler (Sec. 5) to the parallel sampler (Sec. 6). A step in
the parallel sampler can involve multiple proposals carried out in parallel.
For example, updating the chain from $T$ to $T^{\prime}$ can involve $k>0$
(order-invariant) parallel updates, as $T_{i},T_{i1},\dots,T_{ik}=T^{\prime}$.
To clarify the effect of parallelism, every chain is processed post-sampling
to extract a serial chain that is indexed by every proposal. A step of $k$
parallel updates is $k$ steps in the serial chain.

Updates of parallel and single-move samplers occur on the junction tree $T$,
though not all updates on $T$ translate to updates in the underlying
decomposable graph $G=g(T)$. To assess sampler performance on the graph space,
every chain is processed post-sampling to create one with the graph as the
state variable. Metrics measured on this processed chain are referred to as
graph-updates metrics, in lieu of junction-updates metrics.

Lastly, we compare our parallel sampler to the Metropolis–Hastings sampler
over junction trees of Green and Thomas (2013) in Section 7.4.

###  7.2 Study setup for Gaussian graphical model

We simulate 10 datasets from a Gaussian graphical model, each dataset being
generated as follows. We randomly sample a $p=50$ vertex decomposable graph
$G$ using an auto-regressive process with lag varying between 1 and 5. For a
real scalar $\sigma^{2}>0$ and $\rho$, we define the nonnegative definite
matrix $\Theta^{-1}\in\mathbb{R}^{p\times p}$ as

$\left(\Theta^{-1}\right)_{ij}=\begin{cases}\sigma^{2}&\text{if }i=j\\\
\sigma^{2}\rho&\text{if }\left\\{i,j\right\\}\in G.\end{cases}$ (7.1)

and ${\Theta}_{ij}=0$ if $\left\\{i,j\right\\}\not\in G$. We set $\Theta^{-1}$
as the covariance matrix of a $p$-variate zero-mean (for simplicity) Gaussian
distribution, for which $G$ encodes the conditional independence graph. We set
$\sigma^{2}=1$, $\rho=0.9$ and sample 100 data vectors.

This process results in 10 datasets, each having its own structure dependence
graph. Figure 5 illustrates an example of a sampled decomposable graph $G$
alongside its adjacency matrix. This model is known as the graphical Gaussian
intraclass model, used in Green and Thomas (2013); Olsson et al. (2019). From
Grone et al. (1984), since $G$ is decomposable and $\Theta$ restricted to each
clique is positive definite, $\Theta$ exists and is unique.

Figure 5: Simulated 50-vertex auto-regressive decomposable graph and its
adjacency matrix.

To avoid using our knowledge of the generative model (mimicking real-world
conditions), in the posterior sampling, we assign a hyper Wishart prior to
$\Theta$, where for each clique $C$ the degree of freedom is set $\delta=5$,
with the $\left|C\right|$-dimensional identity matrix used as the scale
matrix. This prior leads to a conjugate posterior (as in Appendix A). Unless
otherwise specified, we assign a uniform prior over $T$ by setting
$\pi(T)\propto 1$.

###  7.3 Single-move and parallel samplers

This section compares the single-move sampler (Sec. 5) to the parallel sampler
(Sec. 6). We run a Metropolis–Hastings Markov chain sampler for one million
steps for the single-move sampler, and 500,000 for the parallel sampler. For
both cases, we initiate the sampler with a random $50$-vertex skeleton $t$ and
set it as the junction $T$, this corresponds to junction tree of a no-edge
50-vertex graph. We sample a posterior of the skeleton $\pi(t{\,|\,}T)$ every
100 Metropolis–Hastings samples of $\pi(T{\,|\,}t)$, as recommended by Green
and Thomas (2013). Sampling of $\pi(t{\,|\,}T)$ follows the randomization
approach in Thomas and Green (2009).

Using the dataset associated with Fig. 5, each sampler is initiated
arbitrarily ten times. Although sample performance varies between runs, all
parallel chains have converged within the first 50,000 iterations, and 8 out
of the 10 single-move chains converged by iteration 500,000, with the slowest
converging by iteration 800,000. We observed consistently faster convergence
to the chain’s stationary state, less within-sample correlation and improved
computation speed with the parallel sampler versus the single-move sampler.
For more details, refer to SM Section B, which also includes other graphical
diagnostics of the chains.

First, the computational speed improvement is due to parallelization, where
the parallel sampler performed approximately $\approx$ 710,000 total updates
when run for 500,000 samples (42% improved efficient). Second, the acceptance
probabilities of the parallel sampler do not require the calculation of the
reverse proposal, as is the case with the single-move sampler, which, at a
minimum, requires the enumeration of both partition sets in (4.1) and carrying
out (5.2) and (5.3).

Figure 6 illustrates the adjacency matrix of the maximum a posteriori graph
and the heatmap of one of the parallel sampler chains, all computed over the
last 300,000 samples. There is a strong correspondence between the estimated
and the true adjacency matrix (Fig. 5).

Figure 6: Maximum a posteriori (left) and heatmap (right) of the last 300,000
(out of 500,000) samples of the parallel junction tree sampler.

Figure 7 illustrates the cumulative acceptance rate, as the average number of
accepted proposals, for one of the chains, for various setups. Lines (a) and
(b) consider a proposal as accepted if the underlying decomposable graph is
altered (graph-updates). Line (a) is the parallel chain and Line (b) is the
serial chain. The gap between (a) and (b) illustrates the improvement
resulting from parallelism. Line (c) considers junction-updates of the serial
chain. The gap between (c) and (b) showcases the proportion of junction-only
updates, those that involve non-maximal cliques. As shown, the gap is large
when the chain starts, but diminishes as the chain converges. This suggests
that the advantage of using the expanded junction tree, as we have done in
this work, outweighs the added computational cost associated with junction-
only updates. Line (d) illustrates the acceptance rate of the single-move
sampler, which appears to be inferior when compared to the parallel and serial
samplers. This is a result of the extra proposal ratio used in the single-move
sampler, see (5.1), and to the difference in sampling methods between the two
samplers. The parallel sampler mimics that of a breadth-first sampling
structure, where for every $v\in\mathcal{V}$, all possible updates are
attempted. The single-move sampler follows a depth-first sampling structure,
where for every $v\in\mathcal{V}$ a single clique is attempted. Refer to SM
Section B for more diagnostic plots of the samplers.

Figure 7: Cumulative acceptance rate of the parallel sampler (a)-(c) and the
single-move sampler (d). Measured by graph-updates in (a)-(b) and junction-
updates in (c). Computed over the parallel sampler in (a), and over the serial
sampler in (b) and (c).

Bayesian inference for decomposable graphical models is commendably fast when
compared to their nondecomposable counterparts. For example, under a Gaussian
graphical model setup with $p=150$ vertices, when the underlying graph is
assumed nondecomposable, the state-of-the-art Bayesian structure learning
method executes 1000 iterations in about 12 minutes (Mohammadi et al., 2021,
Fig. 1). When we assume decomposability, the herein proposed parallel sampler
executes 1000 iterations in 0.3 seconds, while the single-move sampler does so
in 0.5 seconds. This 3 to 4 magnitudes improvement in execution time is main
advantage in assuming decomposability of the conditional independence graph
$G$. Refer to SM Fig. 14 for more details on the execution time of the
proposed methods.

Lastly, the average acceptance rate of graph-updates for the uniform prior is
approximately 4.3%, and 5.1% for junction-updates, with a gap of 0.8%. Under
the same setup and initial conditions, we ran the parallel sampler with a
modified version of the exponential family prior, as
$\phi(C)=\exp\left(\alpha\left\\{\left|C\right|-1\right\\}\right)$ and
$\psi(S)=\exp\left(\beta\left|S\right|\right)$ for parameters $\alpha=2$ and
$\beta=4$. The average acceptance rate for the exponential prior is 1.4% over
graph-updates and 1.43% over junction-updates, with a gap of 0.04%. This
confirms that conservative priors reduce the number of non-maximal cliques of
$T$ (discussed in Sec. 6.3). All averages are computed over the serial chain.
Refer to SM Fig. 13 for a visual comparison of the commutative acceptance
rates of those priors.

###  7.4 Comparison to alternative junction tree samplers

We compare our parallel sampler to the junction tree sampler in Green and
Thomas (2013), who proposed a single-move multi-edge junction tree sampler of
decomposable graphs. Unlike our junction-based sampler, their sampler is
graph-based, i.e. at every step they reduce $T$ to $\texttt{reduce}(T)$. Their
multi-edge proposal updates an arbitrary number of edges from their
characterization of decomposability-preserving updates, which is unlike our
all-or-none sampling method, where $v$ is added or removed from a clique.

Figure 8: Likelihood traceplots (top) and within-sample correlation of the
number of graph edges (bottom) of 10 (colour-matched) trajectories for the
junction tree sampler of Green and Thomas (2013) (left) and the proposed
parallel sampler (bottom).

Using the 10 simulated datasets in Section 7.3 for both samplers, we use a
uniform graph prior and assign a Wishart prior to $\Theta$ with degrees of
freedom $\delta=5$. Similarly, both samplers share the same junction tree
randomization step. The difference is conceptual: while Green and Thomas
(2013) refer to this step as junction tree randomization, we designed it as a
hierarchical sampling scheme as in (3.1). This step is carried once every 100
graph samples, for both. We ran the parallel sampler for 500,000 steps
resulting in $\approx$ 700,000 total updates. To have a similar number of
iterations, the Green and Thomas (2013) sampler is run for a million steps.

First, we find that all 10 parallel chains have converged within the first
100,000 steps, while 3 of Green and Thomas (2013) chains did not converge, as
illustrated by the likelihood traceplots in Fig. 8 (top panel). Second, the
within-sample correlation of the number of edges of the underlying
decomposable graph of each chain is almost negligible for the parallel
sampler. For example, the correlation drops below 0.2 for all chains before
lag 2500, and remains close to zero thereafter (Fig. 8 bottom-right). The
within-sample correlation of Green and Thomas (2013) chains remains high for
samples that are more than 80,000 steps away, for most convergent chains (Fig.
8 bottom-left). The within-sample correlation is calculated for the last
300,000 samples of the parallel sampler and 500,000 samples for Green and
Thomas (2013) sampler. Refer to SM Fig. 15 for traceplots of the number of
graph edges for both samplers.

Figure 9: Boxplot of acceptance rates for the Green and Thomas (2013) and the
proposed parallel and serial samplers, over the 10 simulated datasets of
Section 7.2, on graph-updates.

To conclude, parallelism is not the only advantage of our proposed sampler
over that of Green and Thomas (2013). In fact, all proposals by our samplers
are decomposability-preserving. Rejections are only the result of the data
likelihood, under the uniform graph prior. This is a consequence of the
partition sets in (4.1), as they illustrate the sufficient and necessary
conditions for a proposal to preserve decomposability (Thm. 4.2 and 4.3). In
Green and Thomas (2013), the disconnect move (Section 3.2 conditions (a)-(d))
can lead to proposals being rejected as they break the decomposability of the
graph, resulting in sub-optimal mixing. This phenomenon is illustrated in
Boxplot 9, which shows the acceptance rates of each sampler over the 10
datasets, over graph-updates. For example, the average acceptance rate of the
Green and Thomas (2013) sampler over the 10 datasets is $\approx$ 3%, while
the serial sampler has an average of $\approx$ 9% (three times as large). The
parallel sampler achieves an average of $\approx$ 12%.

##  8 Discussion

Computational bottlenecks still remain in the quest to design Markov chain
samplers for high-dimensional decomposable graphs. The junction tree data
structure, although proved more efficient computationally than the graph
structure, is not unique. Counting the number of junction tree representations
for a fixed decomposable graph is known to be computationally expensive
(Thomas and Green, 2009), which led Green and Thomas (2013) to propose a
modified Metropolis–Hastings chain that is inferior to the classical one, to
avoid such computation. Our proposed method trades the need to quantify the
space of junction trees with a form of latent dimensional expansion coupled
with hierarchical sampling of junction trees. Here, the latent expansion
represents non-maximal cliques of underlying decomposable graph.

While some updates in this latent space are probabilistically superfluous, as
they involve non-maximal cliques, this work has demonstrated that such latent
expansion is promising for two reasons: (a) controlling the size of this
expansion

the possible number of non-maximal cliques can be done efficiently and
intuitively by restricting the size of the latent junction trees to the number
of graph vertices; it is intuitive since a no-edge graph has the maximum
number of maximal cliques over the set of vertices; and (b) parallel sampling,
which significant improves mixing properties, easily adapts to our proposed
framework, as we do not require any junction tree modification after an
update. Such modification is required in other junction tree based sampler.

Exploiting the junction tree representation led to the largest gain in
sampling efficiency, as done in the work of Green and Thomas (2013); Giudici
and Green (1999); Olsson et al. (2022, 2019), and herein. The junction tree
data structure helped in the handling of decomposable graphs for computational
purposes. A greater part of this efficiency gain is driven mainly by an
improvement in graph proposals that the junction property explicitly exposed,
which this work further improves.

Olsson et al. (2019) have proposed a sequential Monte Carlo sampler for
decomposable graphical models. Their method is analytically more complex than
a Metropolis–Hastings approach, as proposed here and in Green and Thomas
(2013), and more computationally intensive. As demonstrated in Olsson et al.
(2019), they establish mixing properties that compare to those of this work,
and mixing properties that are better than that of Green and Thomas (2013). We
believe that the superiority gap between the sequential Monte Carlo methods
and the approach of Green and Thomas (2013) can be diminished with a more
efficient disconnect proposal for the latter.

Green and Thomas (2013) have also opted to use one multi-edge update per
proposal, where the set of edges are selected uniformly from the superset of
valid edges. This occurs with probability proportional to $2^{m}-1$, where $m$
is the cardinality of the superset. It remains unclear whether such a multi-
edge update is more efficient than the alternatives, e.g. attempting all
possible subsets of the superset. The average acceptance rate of their multi-
edge proposal is inferior to that of their single-edge proposal (Green and
Thomas, 2013, Supp. Fig. 3). This work opted for an all-or-none multi-edge
update, which is found to be more intuitive with respect to the junction
property, and which outperforms a single-edge proposal.

Exact formulas for the normalizing constants of Wishart priors for
nondecomposable graphical models have recently been illustrated by Uhler et
al. (2018). This unlocks a major analytical bottleneck faced by graphical
models when discarding the decomposability assumption (Jones et al., 2005;
Wang and Carvalho, 2010; Mitsakakis et al., 2011; Roverato, 2002; Atay-Kayis
and Massam, 2005; Dellaportas et al., 2003). However, computational barriers
still remain. As shown in Mohammadi et al. (2021), by using this formula, the
state-of-the-art Bayesian structure learning for graphical models, in
moderate-dimensional problems, can take up to 15 minutes to execute a thousand
iterations. Such execution takes less than a second for decomposable models,
as done here and by others (Green and Thomas, 2013; Giudici and Green, 1999).

Consequently, it still remains unclear whether Bayesian inference for
graphical models can be achieved for large-dimensional problems, e.g. beyond
1000 vertices, with no reliance on the decomposability framework.

## Acknowledgements

The author wishes to acknowledge Felix Rios (Royal Institute of Technology in
Stockholm) for providing necessary code through the trilearn Python package,
and the structure learning benchmarking workflow Benchpress, and for his
comments, which improved the quality of this work. The author was supported by
an NSERC postdoctoral fellowship.

## References

  * Atay-Kayis and Massam (2005) Atay-Kayis, A. and Massam, H. (2005).  A Monte Carlo method for computing the marginal likelihood in nondecomposable Gaussian graphical models.  Biometrika 92 317–335.
  * Bierkens (2016) Bierkens, J. (2016).  Non-reversible Metropolis-Hastings.  Statistics and Computing 26 1213–1228.
  * Bornn and Caron (2011) Bornn, L. and Caron, F. (2011).  Bayesian clustering in decomposable graphs.  Bayesian Analysis 6 829–846.
  * Byrne et al. (2015) Byrne, S., Dawid, A. P. et al. (2015).  Structural Markov graph laws for Bayesian model uncertainty.  The Annals of Statistics 43 1647–1681.
  * Cayley (1889) Cayley, A. (1889).  A theorem on trees.  Quart. J. Math. 23 376–378.
  * Cowell et al. (2006) Cowell, R. G., Dawid, P., Lauritzen, S. L. and Spiegelhalter, D. J. (2006).  Probabilistic Networks and Expert Systems: Exact Computational Methods for Bayesian Networks.  Springer Science & Business Media.
  * Dawid and Lauritzen (1993) Dawid, A. P. and Lauritzen, S. L. (1993).  Hyper Markov laws in the statistical analysis of decomposable graphical models.  The Annals of Statistics 21 1272–1317.
  * Dellaportas et al. (2003) Dellaportas, P., Giudici, P. and Roberts, G. (2003).  Bayesian inference for nondecomposable graphical Gaussian models.  Sankhyā: The Indian Journal of Statistics 43–55.
  * Duncan et al. (2016) Duncan, A. B., Lelievre, T. and Pavliotis, G. A. (2016).  Variance reduction using nonreversible Langevin samplers.  Journal of statistical physics 163 457–491.
  * Frydenberg and Lauritzen (1989) Frydenberg, M. and Lauritzen, S. (1989).  Decomposition of maximum likelihood in mixed graphical interaction models.  Biometrika 76 539–555.
  * Giudici and Green (1999) Giudici, P. and Green, P. (1999).  Decomposable graphical Gaussian model determination.  Biometrika 86 785–801.
  * Green and Thomas (2013) Green, P. J. and Thomas, A. (2013).  Sampling decomposable graphs using a Markov chain on junction trees.  Biometrika 100 91–110.
  * Green and Thomas (2018) Green, P. J. and Thomas, A. (2018).  A structural Markov property for decomposable graph laws that allows control of clique intersections.  Biometrika 105 19–29.
  * Grone et al. (1984) Grone, R., Johnson, C. R., Sá, E. M. and Wolkowicz, H. (1984).  Positive definite completions of partial Hermitian matrices.  Linear algebra and its applications 58 109–124.
  * Hastings (1970) Hastings, W. K. (1970).  Monte Carlo sampling methods using Markov chains and their applications.  Biometrika 57 97–109.
  * Jones et al. (2005) Jones, B., Carvalho, C., Dobra, A., Hans, C., Carter, C. and West, M. (2005).  Experiments in stochastic computation for high-dimensional graphical models.  Statistical Science 20 388–400.
  * Lauritzen (1996) Lauritzen, S. L. (1996).  Graphical Models.  Oxford University Press.
  * Mitsakakis et al. (2011) Mitsakakis, N., Massam, H. and Escobar, M. D. (2011).  A Metropolis-Hastings based method for sampling from the G-Wishart distribution in Gaussian graphical models.  Electronic Journal of Statistics 5 18–30.
  * Mohammadi et al. (2021) Mohammadi, R., Massam, H. and Letac, G. (2021).  Accelerating Bayesian structure learning in sparse Gaussian graphical models.  Journal of the American Statistical Association 1–14.
  * Neal (1998) Neal, R. M. (1998).  Suppressing random walks in Markov chain Monte Carlo using ordered overrelaxation.  In Learning in graphical models. Springer, 205–228.
  * Olsson et al. (2019) Olsson, J., Pavlenko, T. and Rios, F. L. (2019).  Bayesian learning of weakly structural Markov graph laws using sequential Monte Carlo methods.  Electronic Journal of Statistics 13 2865–2897.
  * Olsson et al. (2022) Olsson, J., Pavlenko, T. and Rios, F. L. (2022).  Sequential sampling of junction trees for decomposable graphs.  Statistics and Computing .
  * Rey-Bellet and Spiliopoulos (2015) Rey-Bellet, L. and Spiliopoulos, K. (2015).  Irreversible Langevin samplers and variance reduction: a large deviations approach.  Nonlinearity 28 2081\.
  * Roverato (2002) Roverato, A. (2002).  Hyper inverse wishart distribution for non-decomposable graphs and its application to bayesian inference for gaussian graphical models.  Scandinavian Journal of Statistics 29 391–411.
  * Speed and Kiiveri (1986) Speed, T. P. and Kiiveri, H. T. (1986).  Gaussian Markov distributions over finite graphs.  The Annals of Statistics 138–150.
  * Tarantola (2004) Tarantola, C. (2004).  MCMC model determination for discrete graphical models.  Statistical Modelling 4 39–61.
  * Tarjan and Yannakakis (1984) Tarjan, R. E. and Yannakakis, M. (1984).  Simple linear-time algorithms to test chordality of graphs, test acyclicity of hypergraphs, and selectively reduce acyclic hypergraphs.  SIAM J. Comput. 13 566–579.
  * Thomas and Green (2009) Thomas, A. and Green, P. J. (2009).  Enumerating the junction trees of a decomposable graph.  Journal of Computational and Graphical Statistics 18 930–940.
  * Uhler et al. (2018) Uhler, C., Lenkoski, A. and Richards, D. (2018).  Exact formulas for the normalizing constants of Wishart distributions for graphical models.  The Annals of Statistics 46 90–118.
  * Wang and Carvalho (2010) Wang, H. and Carvalho, C. M. (2010).  Simulation of hyper-inverse Wishart distributions for non-decomposable graphs.  Electronic Journal of Statistics 4 1470–1475.
  * Whittle (1985) Whittle, P. (1985).  Partial balance and insensitivity.  Journal of Applied Probability 22 168–176.
  * Wormald (1985) Wormald, N. C. (1985).  Counting labelled chordal graphs.  Graphs and combinatorics 1 193–200.

##  Appendix A Proof of Theorem 4.2 and 4.3

###  A.1 Theorem 4.2

We consider the simplest case, where $V=\left\\{v\right\\}$ and
$U\subset\mathcal{V}$ is a subset of vertices that form a complete graph in
$G=(\mathcal{V},\mathcal{E})$. There is no edge between $V$ and any node in
$U$, and $G$ is connected. Let $J=(\texttt{cl}(G),\texttt{sep}(G))$, be a
junction tree of $G$. Let $G^{\prime}=(\mathcal{V},\mathcal{E}^{\prime})$ be a
graph such that $\mathcal{E}\subset\mathcal{E}^{\prime}$, and $G^{\prime}$ is
formed by connecting edges $(v,u)$, for every $u\in U$. Let
$C_{V}\in\texttt{cl}(G)$ be a maximal clique containing $V$. Similarly, let
$C_{U}\in\texttt{cl}(G)$ such that $U\subseteq C_{U}$. when $U$ is partitioned
into mutually exclusive sets $\left\\{U^{(1)},\dots,U^{(K)}\right\\}$, then
none of those sets are empty, as it contradicts the fact that $U$ is
completely connected. It is sufficient to prove Theorem 4.2 for the simplest
case when $U=U^{(1)}$, since the same argument can be applied to connect to
$U^{(2)}$, once connected to $U^{(1)}$, and so forth.

123233234$T$$J=\texttt{reduce}(T)$1243Decomposable graph
$G$123234122333324$t$$t_{1}$$t_{2}$$t_{3}$$t_{4}$9876123451110$C_{1}$$C_{2}$$C_{3}$$C_{4}$$C_{5}$$C_{6}$$C_{7}$12$C_{8}$$C_{1}$$C_{2}$$C_{3}$$C_{4}$$C_{5}$$C_{6}$$C_{7}$$C_{8}$$\emptyset$$8,7$$T$3$C_{1}$$C_{2}$$C_{3}$$C_{4}$$C_{5}$$C_{6}$$C_{7}$$C_{8}$$\emptyset$$8,7$$T_{6}$3$\left\\{U\right\\}_{\perp}$$J$:$S\cup
V\cup
U$$\left\\{V\right\\}_{\perp}$$\left\\{V,U\right\\}_{\perp}$$\left\\{U\right\\}_{\perp}$updated:$S\cup
V$$S\cup V$$\left\\{V,U\right\\}_{\perp}$$S\cup V\cup U$$C^{\prime}$$S\cup
U$$S\cup
U$$\left\\{V\right\\}_{\perp}$$\left\\{U\right\\}_{\perp}$$J^{\prime}$:$S\cup
V$$C_{V}$$S$$\left\\{V,U\right\\}_{\perp}$$S$$S$$S\cup
U$$C_{U}$$\left\\{V\right\\}_{\perp}$disconnectconnect

Figure 10: Updating junction tree $J$ to $J^{\prime}$ (and backward) by
(dis)connecting vertices in $V,U\in\mathcal{V}$. The notation
$\left\\{A\right\\}_{\perp}$ indicates a subtree branch that does not contain
$A\subset\mathcal{V}$. For a disconnect, modify $J$ to include two dummy
cliques $S\cup\left\\{V\right\\}$ and $S\cup\left\\{U\right\\}$ in the
junction between $S\cup V\cup U$ and both branches
$\left\\{U\right\\}_{\perp}$ and $\left\\{V\right\\}_{\perp}$. The final tree
$J^{\prime}$ is formed by removing $V\cup U$ from the $S\cup V\cup U$.
$J^{\prime}$ adheres to the running intersection property and forms a
decomposable graph that does not include edges between $V$ and $U$. A connect-
move is the opposite. Circles are cliques, boxes are separators.

  * •

If $C_{U}\in\mathcal{N}_{{V}}(J)$, then $G^{\prime}$ is decomposable. This
case follows directly from (Green and Thomas, 2013, Prop. 1). Here we present
a different proof. The steps are illustrated graphically in Fig. 10. We know
that $(C_{V},C_{U})\in J$. Let $C^{\prime}=V\cup S\cup U$, where $S=C_{V}\cap
C_{U}$. Create $J^{\prime}$ as follows, add $C^{\prime}$ between $C_{V}$ and
$C_{U}$ by replacing the edge $(C_{V},C_{U})$ with $(C_{V},C^{\prime})$ and
$(C^{\prime},C_{U})$. $J^{\prime}$ satisfies the junction property (2.2).
Since for any clique $C_{i},C_{j},\in J$, if their path passes the edge
$(C_{V},C_{U})$, then $C_{i}\cap C_{j}\subseteq S\subset C^{\prime}\in
J^{\prime}$. Also, it is easily verified that
$\left(C_{1},\dots,C_{V},C^{\prime},C_{U},\dots,C_{c}\right)$ is a perfect
ordering sequence that generates $J^{\prime}$, when
$\left(C_{1},\dots,C_{V},C_{U},\dots,C_{c}\right)$ is the perfect ordering
sequence of $G$, as defined in (2.1). Now $C^{\prime}$ or $C_{U}$ might not be
maximal in $J^{\prime}$, but then we can reduce $J^{\prime}$ to
$\texttt{reduce}(J^{\prime})=(\texttt{cl}(G^{\prime}),\texttt{sep}(G^{\prime}))$,
where $G^{\prime}=g(J^{\prime})$.

  * •

We will show the necessary part of the condition by showing that $G^{\prime}$
contains a non-chordal cycle of 4 vertices. We will pick those vertices from
the recursive simplicial subsets over a path in $J$. Assume $G^{\prime}$ is
decomposable, and there does not exist a junction tree $J$ of $G$ such that
$(C_{V},C_{U})\in J$. Assume that there exists a
$J=(\texttt{cl}(G),\texttt{sep}(G))$ such that the path $C_{V}\sim C_{U}$
contains a single maximal clique $C$, i.e. $C_{V}\sim C_{U}=(C_{V},C,C_{U})$.
This path, by definition is part of the perfect ordering sequence of $G$.
Define the following: $S_{V}=C_{V}\cap C$ and $S_{U}=C\cap C_{U}$. Now, we
know that: (i) $S_{U}\setminus C_{V}\neq\emptyset$, otherwise $S_{U}\subseteq
C_{V}$ and a junction tree where the edge $(C_{V},C_{U})$ exists, leading to a
contradiction; and similarly (ii) $S_{V}\setminus C_{U}\neq\emptyset$. Then,
let $b\in S_{U}\setminus C_{V}$ and $a\in S_{V}\setminus C_{U}$. Then all the
following edges exits in $G^{\prime}$, $(v,a),(a,b)(b,u),(v,u)$, where where
$v\in V,u\in U$. However, neither edges $(v,b)$, nor $(a,u)$ exists in
$G^{\prime}$, hence $G^{\prime}$ has a non-chordal cycle of the 4 vertices
$(v,a,b,u,v)$, leading to a contradiction.

Now if the path $C_{V}\sim C_{U}=(C_{V},C,C^{\prime},C_{U})$, one can find a
vertex $d$ in $C\cap C^{\prime}\setminus C_{V}$ such that the non-chordal
cycle $(v,a,d,b,u)$ exist in $G^{\prime}$. Here, $b\in C_{U}\cap
C^{\prime}\setminus C$. It is possible that $d$ equals $a$ or $b$. In all
cases, $G^{\prime}$ is not chordal.

By induction the results of Theorem 4.2 follows.

###  A.2 Theorem 4.3

Under similar setup in proof of Theorem 4.2, in Section A.1. Here $U$ is fully
connected to $V$ in $G$ and $G^{\prime}$ is a graph formed from $G$ by
removing edges $(v,u),$ where $v\in V,u\in U$. Again, it is sufficient to
prove Theorem 4.3 for the simplest case when $U=U^{(1)}$, since the same
argument can be applied to disconnect from $U^{(2)}$, once disconnected from
$U^{(1)}$, and so forth.

  * •

If $U\in C$ for some $C\in\mathcal{L}_{{V}}(T)$, then $V\cup U$ are contained
in only one maximal clique $C$, and $G^{\prime}$ is decomposable. This follows
from (Green and Thomas, 2013, Prop. 2). We present a proof based on the
running intersection property (2.1). This process is demonstrated visually in
Fig. 10. Define $C^{\prime}=V\cup S\cup U$, where $S\subset\mathcal{V}$. Let
$C^{\prime}_{V},C^{\prime}_{U}\in\texttt{cl}(G)$ be such that
$(C_{1},\dots,C^{\prime}_{V},C^{\prime},C^{\prime}_{U},\dots,C_{c})$ is the
perfect ordering sequence of $T$ formed by (2.1). Here
$(C^{\prime}_{V},C^{\prime})\in T$ and $(C^{\prime},C^{\prime}_{U})\in T$, by
definition of junction trees edges with relation to (2.1). Now remove
$C^{\prime}$ from the ordering and add $S$ instead, alongside the cliques
$C_{V}=S\cup V$ and $C_{U}=S\cup U$, to form the new ordering
$(C_{1},\dots,C^{\prime}_{V},C_{V},S,C_{U},C^{\prime}_{U},\dots,C_{c})$ which
respects (2.1) and generates $G^{\prime}$ and
$T^{\prime}=(\texttt{cl}(G^{\prime}),\texttt{sep}(G^{\prime}))$. To see this,
let $C_{i},C_{j}\in\texttt{cl}(G)$ if the path $C_{i}\sim C_{j}$ in $T$ passes
through $(C^{\prime}_{V},C^{\prime})$ or $(C^{\prime},C^{\prime}_{U})$, then
$C_{i}\cap C_{j}\subseteq S$, where $S$ is contained in $C_{V}$ and $C_{U}$ as
well. From (2.1), we have
$C_{V}\cap\left\\{C_{1},\dots,C^{\prime}_{V}\right\\}\subseteq
C^{\prime}_{V}$,
$S\cap\left\\{C_{1},\dots,C^{\prime}_{V},C_{V}\right\\}\subseteq C_{V}$, and
$C_{U}\cap\left\\{C_{1},\dots,C^{\prime}_{V},S\right\\}\subseteq S$, while
everything else remains the same. If $T^{\prime}$ contains non-maximal
cliques, reduce it as $\texttt{reduce}(T^{\prime})$.

  * •

To show the necessary part of the conditions. Assume that $G^{\prime}$ is
decomposable but there does not exist a $C\in\texttt{cl}(G)$ such that
$C\in\mathcal{L}_{{V}}(T)$, for any junction tree $T$ of $G$. We will show by
contradiction, that $G^{\prime}$ has a non-chordal cycle of 4 vertices, by
picking vertices of recursively simplical subsets of a path in $T$ that
contains $V\cup U$. For simplicity, let $V\cup U$ be contained in only two
maximal cliques $C_{1},C_{2}\in\texttt{cl}(G)$, then by the run intersection
property (2.1), there exists a junction tree $T$, of $G$, such that the edge
$(C_{1},C_{2})\in T$. Let $C_{1}=A\cup V\cup U\cup S$ and $C_{2}=S\cup V\cup
U\cup B$, such that $C_{1}\cap C_{2}=V\cup U\cup S$, for
$A,B,S\subset\mathcal{V}$. Here $A$ and $B$ are no connected, and not empty,
i.e. $A\neq\emptyset$ and $B\neq\emptyset$. Otherwise, either $C_{1}$ or
$C_{2}$ would be non-maximal. Disconnecting edges $(v,u),v\in V,u\in U$, in
$G$ would result in the cycle $(a,v,b,u,a)$ for $a\in A$ and $b\in B$. Since
the edges in this cycle associated with $a$ and $b$ still exist. This
contradicts the fact that $G^{\prime}$ is decomposable.

This completes the proof of Theorem 4.3.

## Supplementary Materials

##  Appendix A Inference setup for Gaussian graphical model

Consider a $p$-dimensional zero-mean (for simplicity) Gaussian random vector
$Y\in\mathbb{R}^{n}$ that is globally Markov with respect to a decomposable
graph $G$. Its precision matrix (inverse covariance) $\Theta$ belongs to the
set

$A_{G}=\left\\{\Theta\in\mathcal{M}_{p}^{+}:\Theta_{{ij}}=0\text{ for
all}\left\\{i,j\right\\}\not\in G\right\\},$

where $\mathcal{M}_{p}^{{+}}$ is the space $p\times p$ positive definite
matrices. Conditional independence between $i$th and $j$th variable in
Gaussian models is equivalent to having $\Theta_{ij}=0$ (Speed and Kiiveri,
1986). Let $y=\left\\{y_{i}\right\\}$, for $i=1,\dots,n$, be $n$ observations
from this model. Following the factorization law in (2.3), the full likelihood
of $y$ can be specified by the clique marginals $\texttt{cl}(G)$ of $G$, and
their separators $\texttt{sep}(G)$. For clique $C\in\mathcal{C}$, the
likelihood marginal is

$f(y_{C}\mid\Theta_{C})=\frac{1}{\left(2\pi\right)^{\left|C\right|}}\left|\Theta_{C}\right|^{{n/2}}\exp\left\\{-\mathop{\mathrm{tr}}\left(\Theta_{C}D_{C}\right)/2\right\\},$

where $D$ is the sample covariance matrix, as
$D=n^{{-1}}\sum_{i=1}^{n}y_{i}y_{i}^{\top}$, $\left|\Theta\right|$ is the
determinant of $\Theta$, $\left|C\right|$ is the cardinality of $C$, and
$y_{C}$ is the subvector $y$ indexes by $C$, similarly for all other elements.
Similarly for $f(y_{S}{\,|\,}\Theta_{S})$, where $S\in{\mathcal{S}}$.
Following Dawid and Lauritzen (1993), a conjugate prior for $\Theta$ is a
hyper-Wishart distribution that is also hyper Markov with respect to $G$,
having the form

$\eta\left(\Theta{\,|\,}G,\varphi\right)=\frac{\prod_{C\in\mathcal{C}}\eta\left(\Theta_{C}{\,|\,}\varphi_{C}\right)}{\prod_{S\in{\mathcal{S}}}\eta\left(\Theta_{S}{\,|\,}\varphi_{S}\right)}.$

Each prior marginal $\eta\left(\Theta_{C}{\,|\,}\varphi_{C}\right)$ is of the
form

$\eta\left(\Theta_{C}{\,|\,}\varphi_{C}\right)\propto{\left|\Theta_{C}\right|^{\beta_{C}}\exp\left\\{-\mathop{\mathrm{tr}}\left(\Theta_{C}Q_{C}\right)/2\right\\}},$

with normalization constant
$\int\eta\left(\Theta_{C}{\,|\,}\varphi_{C}\right)d\Theta_{C}=2^{\delta\left|C\right|/2}{\Gamma_{\left|C\right|}(\beta_{C})}/{\left|Q_{C}\right|^{\beta_{C}}}$,
where $\Gamma_{k}$ denotes the multivariate gamma function. Here
$\beta_{C}=(\delta+\left|C\right|-1)$, $\varphi_{C}=(\delta,Q_{C})$,
$Q\in\mathcal{M}_{p}^{+}$ a scale positive definite matrix and
$\delta>\max_{C\in\mathcal{C}}\left\\{\left|C\right|-1\right\\}$ the number of
degrees of freedom. By Dawid and Lauritzen (1993, Thm. 3.9), the collection of
priors ${\left\\{\eta\left(\Theta_{C}{\,|\,}\varphi_{C}\right)\right\\}}$ are
pairwise hyper consistent and there exists a unique hyper Wishart joint
distribution of the form

$f(y_{C}\mid\Theta_{C})\eta\left(\Theta_{C}{\,|\,}\varphi_{C}\right)\propto\frac{1}{\left(2\pi\right)^{\left|C\right|}}\left|\Theta_{C}\right|^{\alpha_{Q}}\exp\left\\{-\mathop{\mathrm{tr}}{\left(\Theta_{C}\left(D_{C}+Q_{c}\right)\right)/2}\right\\},$

where $\alpha_{C}=\left(\delta+n+\left|C\right|-1\right)/2$. By conjugacy, it
is possible to integrate out $\Theta$ of the term above. With a junction prior
of the form (5.4), the junction posterior is

$\pi(T{\,|\,}t)\propto\frac{\prod_{C\in\mathcal{C}}\phi(C)\rho(C)}{\prod_{S\in{\mathcal{S}}}\psi(S)\rho(S)},$
(A.1)

with

$\rho(C)=\frac{\left|Q_{C}\right|^{\alpha_{C}}}{\left|Q_{C}+D_{C}\right|^{\beta_{C}}}\frac{\Gamma_{\left|C\right|}\left(\alpha_{C}\right)}{\Gamma_{\left|C\right|}\left(\beta_{C}\right)},$

and $\rho(S),S\in{\mathcal{S}}$ is defined similarly.

##  Appendix B Diagnostic plots of the single-move and parallel samplers

Figure 11: Diagnostic plots of the parallel and single-move samplers when
initiated ten times arbitrary on a single dataset. All plots are computed on
the last 200,000 samples of the chain. except the log-likelihood plot.

Figure 12: Cumulative acceptance rates of the parallel, serial, and single-
move samplers, when initiated ten times arbitrary on a single dataset.
Correlation plots are computed on the last 200,000 samples of a chain. Figure
13: Comparison between uniform prior and the clique exponential family prior,
of the parallel sampler, where
$\phi(C)=\exp\left(\alpha\left\\{\left|C\right|-1\right\\}\right)$ and
$\psi(S)=\exp\left(\beta\left|S\right|\right)$ for parameters $\alpha=2$ and
$\beta=4$. Figure 14: The average time in seconds it takes for 1000 iterations
for the Gaussian decomposable graphical model, with setup in Section 7.2,
under the the parallel sampler, its serial variate, and the single-move
sampler. The number of simulated data samples $n$ is equal to the number of
vertices.

##  Appendix C Diagnostic plots of the parallel sampler

Figure 15: Traceplots of the number of edges of the underlying decomposable
graph over the 10 datasets in Section 7.2, of Green and Thomas (2013) and the
parallel sampler.