# Likelihood Based Inference in Fully and Partially Observed Exponential
Family Graphical Models with Intractable Normalizing Constants

Yujie Chen  Anindya Bhadra and Antik Chakraborty
Department of Statistics Correspondence: 150 N. University St., West
Lafayette, IN 47907, USA. Email<EMAIL_ADDRESS> Purdue University

###### Abstract

Probabilistic graphical models that encode an underlying Markov random field
are fundamental building blocks of _generative modeling_ to learn _latent
representations_ in modern multivariate data sets with complex dependency
structures. Among these, the exponential family graphical models are
especially popular, given their fairly well-understood statistical properties
and computational scalability to high-dimensional data based on pseudo-
likelihood methods. These models have been successfully applied in many
fields, such as the Ising model in statistical physics and count graphical
models in genomics. Another strand of models allows some nodes to be _latent_
, so as to allow the marginal distribution of the observable nodes to depart
from exponential family to capture more complex dependence. These approaches
form the basis of generative models in artificial intelligence, such as the
Boltzmann machines and their restricted versions. A fundamental barrier to
likelihood-based (i.e., both maximum likelihood and fully Bayesian) inference
in both fully and partially observed cases is the intractability of the
likelihood. The usual workaround is via adopting pseudo-likelihood based
approaches, following the pioneering work of Besag (1974). The goal of this
paper is to demonstrate that full likelihood based analysis of these models is
feasible in a computationally efficient manner. The chief innovation lies in
using a technique of Geyer (1991) to estimate the intractable normalizing
constant, as well as its gradient, for intractable graphical models. Extensive
numerical results, supporting theory and comparisons with pseudo-likelihood
based approaches demonstrate the applicability of the proposed method.

Keywords: Boltzmann Machines, Contrastive Divergence, Generative models, Ising
Model, Poisson Graphical Model.

##

###  1 Introduction

Exponential family graphical models provide a coherent framework for
describing dependence in multivariate data. However, likelihood-based
statistical inference is typically hindered by the presence of an intractable
normalizing constant, and due to the lack of scalability of doubly intractable
procedures (e.g., Murray et al., 2012; Møller et al., 2006; Liang, 2010),
existing methods rely on pseudo-likelihood or approximate Bayesian methods to
facilitate estimation (Ravikumar et al., 2010; Besag, 1974). Although these
methods could be consistent under certain assumptions, it is a classical
result that full likelihood methods are statistically efficient, at least in
an asymptotic sense and for correctly specified models (Fisher, 1922; Rao,
1945). This apparent bottleneck also applies to the _partially observed
exponential family graphical models_ , where latent variables are introduced
in such a way that the joint distribution of latent and observed variables
belongs to an exponential family, although the marginal distribution of the
observed variables may not. This leads to the so called _product of experts_
family of models for the observed variables (Hinton et al., 1998; Hinton,
2002), which allows complicated dependence structure not belonging in the
exponential family to be modeled, and is considered a foundational building
block for _generative artificial intelligence_ , through an architecture known
as the _Hopfield Network_ or _Boltzmann Machine_ (Hopfield, 1982), whose
purpose is to model Hebbian learning in the brain neural circuitry (Hebb,
1949). The intractable normalizing constant in this case has precluded a
practically successful deployment of a full Boltzmann Machine so far, and has
led to _reduced form architectures_ where training is possible. For example,
in Restricted Boltzman Machines (RBM) or Deep Boltzmann Machines (DBM)
(Hinton, 2007; Fischer and Igel, 2012; Salakhutdinov and Hinton, 2009), where
observed multivariate data is modeled as the marginal distribution over
unobserved hidden variables; a key assumption is the independence of the
hidden variables given the visible variables, and vice versa, which leads to
training by _approximate likelihood_ based methods, such as the contrastive
divergence (CD) algorithm of Hinton (2002). Yet, it is known that with high
probability the CD solution differs from the maximum likelihood solution
(Sutskever and Tieleman, 2010), and the justification for the artificial
bipartite grouping of latent and observed variables in an RBM, or its deep
extensions such as the DBM (Salakhutdinov and Hinton, 2009), is primarily
computational, rather than substantive.

To address these issues, this article develops an array of computationally
feasible methods that enable full-likelihood based analysis (both maximum
likelihood and fully Bayesian) for both fully and partially observed
exponential family graphical models. This is achieved without resorting to
computationally expensive _doubly intractable_ Metropolis-Hastings approaches
(Murray et al., 2012; Møller et al., 2006), which allows our method to scale
to hundreds of variables in both fully and partially observed cases. We also
propose a technique, which we believe to be the first viable method, for
training a _full_ Boltzmann machine, without any restrictions on its
architecture, such as a bipartite graph. The key to our scalability is
sampling under the _independence model_ , underlying any exponential family
graphical model, which accomplishes two things. First, it allows us to use
vanilla Monte Carlo instead of running a prohibitive MCMC scheme after each
parameter update as in CD. Second, we show that the variance of the importance
sampling estimate is well controlled with our choice. We also develop
theoretical results for likelihood based inference in both fully and partially
observed cases. These have received little attention so far, due to the
computational intractability of likelihood-based inference, in the first
place.

The remainder of this article is organized as follows. In Section 2 we provide
the necessary background material on fully and partially observed exponential
family graphical models, where for the latter we focus in particular on
Boltzmann machines and their variants. The main methodological innovation is
discussed in Section 3, where we outline the use of Geyer (1991) for computing
the intractable likelihood and its gradient. Applications to fully and
partially observed models are discussed in Section 4 and 5, followed by
theoretical properties of likelihood-based inference in Section 6. Numerical
experiments and data analysis results are in 7 and 8, followed by concluding
remarks in Section 9.

###  2 Preliminaries

####  2.1 Pairwise _fully observed_ exponential family graphical models

Consider a graph $G=(V,E)$ with $|V|=p$ nodes. Associate each node with a
random variable $X_{j}$ with sample space $\mathcal{X}_{j}$, for
$j=1,\ldots,p$. The edge set consists of pairs $(j,k)$ such that the nodes $j$
and $k$ are connected. We assume the connection is undirected, i.e. if
$(j,k)\in E$ then $(k,j)\in E$. A clique $C\subset V$ of a graph is a fully-
connected subset of vertices. Define $\mathcal{C}$ to be the set of all
cliques of the graph. A probability distribution on the graph can be defined
using a product of clique-wise compatibility functions $\psi_{C}$, i.e.,

$p(x_{1},\ldots,x_{p})=\frac{1}{z}\prod_{C\in\mathcal{C}}\psi_{C}(x_{C}),$ (1)

where $z$ is such that $p(x_{1},\ldots,x_{p})$ is a valid density or mass
function. Valid distributions such as (1) under certain conditions can be
arrived at from another perspective. Indeed, consider the node-conditional
distribution $X_{j}\mid X_{-j}$, i.e. the conditional distribution of $X_{j}$
given the other variables. Suppose this distribution belongs to the univariate
exponential family of models. Then, by the Hammersley–Clifford theorem
(Lauritzen, 1996), the joint distribution produced by these node-conditional
distribution is of the form (1). This particular class of models is known as
exponential family graphical models (Yang et al., 2015), for which the
compatibility functions for each clique are the dot product of the respective
sufficient statistics and the parameter of the distribution. When the clique
set is the set of all nodes and the set of all edges, one obtains what is
known as the pairwise exponential family graphical models (PEGM), with pdf
given by:

$p_{\theta}(x_{1},\ldots,x_{p})=\frac{1}{z(\theta)}\exp\left\\{\sum_{j\in
V}\theta_{j}T_{j}(x_{j})+\sum_{(j,k)\in
E}\theta_{jk}T_{jk}(x_{j},x_{k})+\sum_{j\in V}C(x_{j})\right\\},$ (2)

where $T_{j}(x_{j})$ and $T_{jk}(x_{j},x_{k})$ are the sufficient statistics,
and $\theta\in\Theta\subset\mathbb{R}^{p\times p}$ is the parameter matrix
with $\theta_{jk}=\theta_{j}$ if $j=k$ and $\theta_{jk}=\theta_{kj}$ if $j\neq
k$. Further generalizations of these models when the cliques are not
restricted to be the set of vertices and edges are discussed in Yang et al.
(2015). Conditional independence structures in (2) can also be defined in
terms of the pairwise parameters $\theta_{jk}$, since $\theta_{jk}=0$ if and
only if $(j,k)\notin E$. In other words, $X_{j}\perp\\!\\!\\!\perp X_{k}\mid
X_{-(j,k)}$ if and only if $\theta_{jk}=0$. It is straightforward to observe
that if $X$ has density given by (2), then the node-conditional distribution
of $X_{j}\mid X_{-j}$ is:

$\displaystyle p_{\theta}(x_{j}\mid
x_{-j})=\frac{1}{z(\theta;x_{-j})}\exp\left\\{\theta_{j}T_{j}(x_{j})+C(x_{j})+2\underset{k\in
N(j)}{\sum}\theta_{jk}T_{jk}(x_{j},x_{k})\right\\},$

where $z(\theta;x_{-j})$ is the normalizing constant belonging to a univariate
exponential family distribution for $(X_{j}\mid X_{-j})$, and $N(j)$ denotes
all neighbors of $j$, as encoded by $G$. The parameter space $\Theta$ for
which (2) is a valid probability model is convex (Barndorff-Nielsen, 1978;
Wainwright and Jordan, 2008) and the sample space is
$\mathcal{X}=\mathcal{X}_{1}\times\ldots\times\mathcal{X}_{p}$. When $\Theta$
is open, the model is said to be regular. Some specific PEGM examples are as
follows.

######  Example 2.1.

(Ising). The Ising model (Ising, 1924) is a PEGM with $C(x_{j})=0$,
$T_{j}(x_{j})=x_{j}$ and $T_{jk}(x_{j},x_{k})=x_{j}x_{k}$. Moreover,
$X_{j}\mid X_{-j}$ has a Bernoulli distribution with natural parameter
$\theta_{j}+2\sum_{k\in N(j)}\theta_{jk}x_{k}$. The sample space for this
model is $\\{0,1\\}^{p}$ and $\Theta=\mathbb{R}^{p\times p}$.

######  Example 2.2.

(PGM). The Poisson graphical model (PGM), proposed by Besag (1974), and
further developed by Yang et al. (2013) is a PEGM with $C(x_{j})=\log x_{j}!$,
$T_{j}(x_{j})=x_{j}$ and $T_{jk}(x_{j},x_{k})=x_{j}x_{k}$. Here $X_{j}\mid
X_{-j}\sim\text{Poi}(\theta_{j}+2\sum_{k\in N(j)}\theta_{jk}x_{k})$. The
sample space is $\\{\mathbb{N}\cup\\{0\\}\\}^{p}$, i.e. the set of all
$p$-dimensional count vectors, and $\Theta$ is the set of all $p\times p$ real
matrices with non-positive off-diagonal elements.

Other examples include the Potts model (Potts, 1952), or the truncated Poisson
graphical model (Yang et al., 2013). In all of the examples above, likelihood-
based inference is generally infeasible due to the intractability of
$z(\theta)$ and available methods approximate the likelihood $\ell(\theta)$ by
$\prod_{j}\ell(X_{j}\mid X_{-j})$, the _pseudo-likelihood_. There are two key
reasons for this: (1) $X_{j}\mid X_{-j}$ is a member of a univariate
exponential family and (2) due to the product form of the pseudo-likelihood,
inference can be carried out in parallel for the $p$ nodes.

####  2.2 Pairwise _partially observed_ exponential family graphical models:
product of experts and _approximate_ maximum likelihood via contrastive
divergence

Consider a binary vector $\mathbf{v}\in\\{0,1\\}^{p}$ which might encode the
on and off pixels of a black-and-white image. In generative artificial
intelligence (AI), the goal is to learn the underlying features of the
distribution of $\mathbf{v}$ so that high-quality synthetic data from the
trained AI can be produced for new content. One potential modeling choice
could be an exponential family model, e.g., the Ising. For modern complex data
such as images, this restriction to exponential family could severely limit
the model’s ability to capture more complex dependence not belonging to the
exponential family. As an alternative, Boltzmann Machines have emerged as one
of the fundamental tools for generative AI. These can be thought of as
stochastic neural networks represented by a graph that encodes symmetric
connections between observed variables and hidden or latent variables (Ackley
et al., 1985). The joint distribution of the $p$ observed and $m$ hidden
variables are assumed to belong to a PEGM (Ising in this case):
$p_{\theta}(\mathbf{v},\mathbf{h})=z(\theta)^{-1}\exp({-E_{\theta}(\mathbf{v},\mathbf{h})})$
where $E_{\theta}(\mathbf{v},\mathbf{h})$ is known as the energy function in
the literature and
$z(\theta)=\sum_{\mathbf{v},\mathbf{h}}\exp({-E_{\theta}(\mathbf{v},\mathbf{h})})$.
Specifically,

$\displaystyle\log
p_{\theta}(\mathbf{v},\mathbf{h})=\sum_{j}\theta_{jj}v_{j}+\sum_{k}\theta_{kk}h_{k}+\sum_{j\neq
j^{\prime}}\theta_{jj^{\prime}}v_{j}v_{j^{\prime}}+\sum_{k\neq
k^{\prime}}\theta_{kk^{\prime}}h_{k}h_{k^{\prime}}+\sum_{j,k}\theta_{jk}v_{j}h_{k}-\log
z(\theta),$

for $j,j^{\prime}=1,\ldots,p;\;k,k^{\prime}=p+1,\ldots,p+m$ and
$\theta\in\mathbb{R}^{(p+m)\times(p+m)}$. The marginal distribution of the
visible variables,
$p_{\theta}(\mathbf{v})=\sum_{\mathbf{h}}p_{\theta}(\mathbf{v},\mathbf{h})$,
does not necessarily belong to the exponential family. Introduction of the
hidden variables allows one to describe complex dependency structures in the
visible variables using simple conditional distributions. Moreover, the hidden
variables act as non-linear feature detectors (Hinton, 2007). It can be seen
that for this model, the gradient of the log-likelihood for a visible sample
$\mathbf{v}$ is:

$\displaystyle\dfrac{\partial\ell(\theta)}{\partial\theta}=\dfrac{\partial\log\sum_{\mathbf{h}}p_{\theta}(\mathbf{v},\mathbf{h})}{\partial\theta}=\sum_{v,h}\dfrac{\partial
E_{\theta}(\mathbf{v},\mathbf{h})}{\partial\theta}p_{\theta}(\mathbf{v},\mathbf{h})-\sum_{v,h}\dfrac{\partial
E_{\theta}(\mathbf{v},\mathbf{h})}{\partial\theta}p_{\theta}(\mathbf{h}\mid\mathbf{v}),$
(3)

where (3) can be interpreted as the difference of two expectations: expected
value of the gradient of the energy function under the joint distribution of
$(\mathbf{v},\mathbf{h})$ and under the conditional distribution of
$\mathbf{h}\mid\mathbf{v}$. For most general versions of this model, existing
learning algorithms approximate either both or at least one of these two terms
using MCMC sampling, which is too slow for practical deployment, and lacks
scalability to larger dimensions.

Equation (3) also illustrates that the most general version of the Boltzmann
machine is hard to train and computationally demanding, which has paved the
way for simplified versions of Boltzmann machines, specifically the restricted
Boltzmann machine (RBM) where the underlying graph is bipartite. RBM
(Smolensky, 1986) is inspired by the biological connectivity network among
neurons. RBM’s multi-layer extension, or deep Boltzmann machine (DBM)
(Salakhutdinov et al., 2007) has proven to be a popular machine learning tool
for unsupervised learning (Zhang et al., 2019) and is a _universal
approximator_ for any distribution over $\\{0,1\\}^{p}$ (Le Roux and Bengio,
2008; Montufar and Ay, 2011). Adjacent layers of a DBM have the same
architecture as an RBM. Figure 1 shows the graphical model for a BM, RBM and a
DBM with three hidden layers.

$v_{1}$$v_{2}$$v_{3}$$h_{1}$$h_{2}$$h_{3}$ Fig. 1(a): BM

$v_{1}$$v_{2}$$v_{3}$$h_{1}$$h_{2}$$h_{3}$ Fig. 1.(b): RBM

$v_{1}$$v_{2}$$v_{3}$$h_{1}^{(1)}$$h_{2}^{(1)}$$h_{3}^{(1)}$$h_{1}^{(2)}$$h_{2}^{(2)}$$h_{3}^{(2)}$$h_{1}^{(3)}$$h_{2}^{(3)}$$h_{3}^{(3)}$
Fig. 1.(c): DBM

$v_{1}$$v_{2}$$v_{3}$$h_{1}^{(2)}$$h_{2}^{(2)}$$h_{3}^{(2)}$$h_{1}^{(1)}$$h_{2}^{(1)}$$h_{3}^{(1)}$$h_{1}^{(3)}$$h_{2}^{(3)}$$h_{3}^{(3)}$
Fig. 1.(d): DBM rearranged as an RBM

Figure 1: From left to right: BM, RBM, DBM with three visible nodes. DBM has
three layers of hidden variables. The notation is: hidden nodes
($h_{j}\in\\{0,1\\}$, shaded in gray), visible nodes ($v_{k}\in\\{0,1\\}$,
transparent). Deep hidden nodes in layer $l$ are denoted by $h^{(l)}$.

Learning in the RBM is facilitated by the conditional independence of the
hidden variables given the visible variables and vice versa. Specifically, for
RBM:
$-E_{\theta}(\mathbf{v},\mathbf{h})=\sum_{k}\theta_{kk}h_{k}+\sum_{j}\theta_{jj}v_{j}+\sum_{j}\sum_{k}\theta_{jk}v_{j}h_{k}$,
which implies
$\mathbb{P}(h_{k}=1\mid\mathbf{v})=\sigma(\theta_{kk}+\sum_{j}\theta_{jk}v_{j})$
and
$\mathbb{P}(v_{j}=1\mid\mathbf{h})=\sigma(\theta_{jj}+\sum_{k}\theta_{jk}h_{k})$
where $\sigma(x)=\\{1+\exp(-x)\\}^{-1},\;x\in\mathbb{R}$ is the logistic
sigmoid function. The bipartite dependence structure gives a _product of
experts_ marginal model for the visible variables (Fischer and Igel, 2012):

$p_{\theta}(\mathbf{v})=z(\theta)^{-1}\prod_{j=1}^{p}\exp({\theta_{jj}v_{j}})\prod_{k=p+1}^{p+m}\left(1+\exp({{\theta_{kk}+\sum_{j=1}^{p}\theta_{jk}v_{j}}})\right).$
(4)

The contrastive divergence (CD) algorithm (Hinton, 2002) is an approximate
procedure for finding the MLE under this model via gradient ascent. Suppose
$\theta=\theta^{(t)}$ at iteration $t$. Using (3) for an RBM one obtains:

$\displaystyle\nabla_{\theta^{(t)}}\log
p_{\theta_{t}}(\mathbf{h},\mathbf{v})=\langle\mathbf{h}\mathbf{v}^{T}\rangle_{\mathrm{data}}-\langle\mathbf{h}\mathbf{v}^{T}\rangle_{\mathrm{model}},$

where the angle brackets with subscript “data” denote an expectation with
respect to $p_{\theta^{(t)}}(\mathbf{h}\mid\mathbf{v})$ at the observed
$\mathbf{v}$, which is analytic; and those with subscript “model” denote an
expectation with respect to $p_{\theta^{(t)}}(\mathbf{h},\mathbf{v})$; which
is typically not available in closed form and is evaluated using MCMC.
However, the difficulty here is that direct and exact simulation from
$p_{\theta^{(t)}}(\mathbf{h},\mathbf{v})$ is infeasible. Thus, the CD
algorithm uses a Gibbs sampler to iteratively sample
$(\mathbf{h}\mid\mathbf{v},\theta^{(t)})$ and
$(\mathbf{v}\mid\mathbf{h},\theta^{(t)})$. The key issues are:

  1. 1.

The procedure requires running a new Gibbs sampler every time $\theta^{(t)}$
is updated, until convergence to the target distribution
$p_{\theta^{(t)}}(\mathbf{h},\mathbf{v})$ is achieved. This is reminiscent of
the double MH procedure (Liang, 2010), which is computationally prohibitive
and is one of the primary motivations behind the bipartite graph of RBM that
allows _batch sampling_ of $\mathbf{h}\mid\mathbf{v},\theta^{(t)}$ and
$\mathbf{v}\mid\mathbf{h},\theta^{(t)}$ by drawing independent Bernoulli
vectors. The conditional independence structure is lost if one allows within
layer connections, precluding the use of more general Boltzmann machines.

  2. 2.

Even for an RBM with batch sampling, convergence to the stationary target
$p_{\theta^{(t)}}(\mathbf{h},\mathbf{v})$ is only achieved if the chain is run
for $K\to\infty$ MCMC iterations at every iteration $t$, under standard MCMC
theory. In practice, the CD algorithm is only run for a small $K$ (CD-$K$) or
even $K=1$ Gibbs iterations for every $t$. In this case, it can be shown that
CD is not an estimate of the gradient of _any_ function (Section 4, Sutskever
and Tieleman, 2010). Thus, gradient ascent cannot converge to the MLE,
although in practice it might have reasonable practical performance.

###  3 A Monte Carlo estimate of $z(\theta)$ and its Gradient for Exponential
Family Graphical Models

Let $p_{\theta}(x)={q_{\theta}(x)}/{z(\theta)},\,x\in\mathbb{R}^{p}$ be a
probability density function such that $z(\theta)=\int
q_{\theta}(x)dx=\exp(A(\theta))$ and $\theta\in\Theta$. Consider the situation
where the explicit form of $q_{\theta}(x)$ is known for any $\theta$ in the
parameter space and is computationally cheap to evaluate but the analytical
form of $z(\theta)$ is intractable for every $\theta\in\Theta$. Given $n$
independent observations from such a model, form the matrix
$\mathbf{X}\in\mathbb{R}^{n\times p}$. For notational convenience, we later
use $\mathbf{X}=(X_{1},\ldots,X_{p})$, where each $X_{j}$ is understood to be
a vector of the $j$th variable over $n$ samples for $j=\\{1,\ldots,p\\}$.
Otherwise, when the indexing is over the $i$th row of $\mathbf{X}$, it is made
explicit by writing $X_{i\bullet}$ for $i=\\{1,\ldots,n\\}$. Let the log-
likelihood at two parameter values $\theta,\phi\in\Theta$ be $\ell(\theta)$
and $\ell(\phi)$, respectively. Then,

$\displaystyle\ell(\theta)-\ell(\phi)=\log\frac{p_{\theta}(\mathbf{X})}{p_{\phi}(\mathbf{X})}=\sum_{i=1}^{n}\log\frac{q_{\theta}(X_{i\bullet})}{q_{\phi}(X_{i\bullet})}-n\log\frac{z(\theta)}{z(\phi)},$

where the ratio ${z(\theta)}/{z(\phi)}$ cannot be explicitly computed.
However, Geyer (1991) noted:

$\displaystyle\dfrac{z(\theta)}{z(\phi)}$ $\displaystyle=\frac{1}{z(\phi)}\int
q_{\theta}(x)dx=\dfrac{1}{z(\phi)}\int\dfrac{q_{\theta}(x)}{q_{\phi}(x)}q_{\phi}(x)dx=\int\dfrac{q_{\theta}(x)}{q_{\phi}(x)}p_{\phi}(x)dx=\mathbb{E}_{Y\sim
p_{\phi}}\left[\dfrac{q_{\theta}(Y)}{q_{\phi}(Y)}\right],$ (5)

motivating the Monte Carlo estimate
$\frac{1}{N}\sum_{i=1}^{N}\frac{q_{\theta}(Y_{i\bullet})}{q_{\phi}(Y_{i\bullet})}$,
where $Y_{i\bullet}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}p_{\phi}$. For any
generic $\phi\in\Theta$, the node-conditional distributions can be used to
devise a Markov chain Monte Carlo (MCMC) scheme to sample from $p_{\phi}$.
However, if specific choices of $\phi$ are available such that independent
samples are possible to simulate cheaply, then these samples could also be
utilized to approximate the expectation. Moreover, if for such a $\phi$, the
corresponding normalizing constant $z(\phi)$ is available in closed form, one
has an estimate of $z(\theta)$, for all $\theta$. The key to the current work
is our choice: $\boxed{\phi=\text{diag}(\theta)}$. With this choice: (a) a
sample $Y\sim p_{\phi}$ can be obtained by sampling $Y_{j}\sim
p_{\theta_{jj}}$ independently and setting $Y=(Y_{1},\ldots,Y_{p})$, and (b)
$z(\phi)$ is analytically tractable. This is in general true for any PEGM
since the independent model is simply a distribution over $p$-univariate
exponential families, and a few examples are provided next. For the Ising
model,
$z(\phi)=\prod_{j=1}^{p}z(\theta_{jj})=\prod_{j=1}^{p}(1+\exp({\theta_{jj}}))$.
For the Poisson Graphical model,
$z(\phi)=\prod_{j=1}^{p}z(\theta_{jj})=\prod_{j=1}^{p}\exp(\exp({\theta_{jj}}))$.
For observed samples $\mathbf{X}$, the log-likelihood is
$\ell(\theta)=\sum_{i=1}^{n}\log q_{\theta}(X_{i\bullet})-n\log z(\theta)$. To
obtain the maximum likelihood estimator of $\theta$, one can use the gradient
ascent method leading to updates of the form
$\theta^{(t+1)}=\theta^{(t)}+\gamma\nabla_{\theta}\ell(\theta)=\theta^{(t)}+\gamma[\sum_{i=1}^{n}\nabla_{\theta}\log
q_{\theta}(X_{i\bullet})-n\nabla_{\theta}\log z(\theta)]$, where $\gamma$ is a
suitable step-size. Clearly, for such an algorithm to succeed, one needs
access to the gradient $\nabla_{\theta}\log
z(\theta)=\nabla_{\theta}z(\theta)/z(\theta)$. In the next proposition, an
unbiased estimate for $\nabla_{\theta}z(\theta)$ is proposed, similar in
spirit to the estimate of $z(\theta)$.

######  Proposition 3.1.

$\dfrac{\nabla_{\theta}z(\theta)}{z(\phi)}=\mathbb{E}_{Y\sim
p_{\phi}}\left[\dfrac{\nabla_{\theta}q_{\theta}(Y)}{q_{\phi}(Y)}\right].$

###### Proof.

In (5) replace $z(\theta)$ with $\nabla_{\theta}z(\theta)$ and differentiate
under the integral sign.∎

While the Monte Carlo version of estimates of $z(\theta)$ and
$\nabla_{\theta}z(\theta)$ are unbiased by construction, their practical
utility is determined by their variance, which need not be bounded for an
arbitrary $\phi$ (Geyer and Thompson, 1992). However, for our specific choice
of $\phi$, the following key proposition establishes that the variances of
Monte Carlo estimates of $z(\theta)$ and $\nabla_{\theta}z(\theta)$ are finite
under mild conditions.

######  Proposition 3.2.

Let $u=\theta-\mathrm{diag}(\theta)$, and $\phi=\mathrm{diag}(\theta)$. Then,
(a) Monte Carlo estimate for $z(\theta)$ has bounded variance if
$(2u+\phi)\in\Theta$ and, moreover, (b) the Monte Carlo estimate of
$\nabla_{\theta}z(\theta)$ also has bounded variance if
$(2(1+\delta)u+\phi)\in\Theta$ for some $\delta>0$.

Refer to Supplementary Section S.1.1 for a proof of Proposition 3.2. We
proceed to make several remarks on the implications of this proposition.

######  Remark 3.1.

Note that,
$\frac{1}{1+\delta}(2(1+\delta)u+\phi)+\frac{\delta}{1+\delta}\phi=2u+\phi.$
Thus, by the convexity of exponential family parameter space (see, e.g.,
Theorem 9.1 of Barndorff-Nielsen, 1978), the result of Proposition 3.2, Part
(b) implies that of Part (a), for all $\delta>0$.

######  Remark 3.2.

Note that if $\phi=\phi_{0}$ is an arbitrary reference point not depending on
$\theta$, the Monte Carlo variance of the importance ratio
$q_{\theta}(X)/q_{\phi_{0}}(X)$ can easily become unbounded, even if
$\phi_{0}\in\Theta$. This is observed, for example, in a Poisson graphical
model, by choosing a $\phi_{0}$, such that
$\mathrm{offdiag}(\theta-\phi_{0})>0$. This shows the importance of setting
$\phi=\mathrm{diag}(\theta)$, and the technique of Geyer (1991) is not very
useful under an arbitrary $\phi$, possibly explaining why the technique has
not been widely attempted in this model. Setting $\phi=\mathrm{diag}(\theta)$
removes the first order $T_{j}(X_{j})$ terms from the ratio
$q_{\theta}(X)/q_{\phi}(X)$, which is helpful for bounding the variance of
what remains.

######  Remark 3.3.

Another benefit of setting $\phi=\mathrm{diag}(\theta)$, is that one is
(exactly) sampling from the independence model $p_{\phi}(\cdot)$ for the
purpose of Monte Carlo, which allows batch sampling, and is computationally
far more efficient than (approximately) sampling from $p_{\theta}(\cdot)$ at
every iteration, for which an iterative technique such as MCMC must be
deployed. Variations of this latter approach, i.e., sampling auxiliary data
$Y_{1},\ldots,Y_{N}$ under $p_{\theta}(\cdot)$ given the current $\theta$ to
approximate $z(\theta)$ and $\nabla_{\theta}\log z(\theta)$, are used in the
double Metropolis–Hastings of Liang (2010), the contrastive divergence
approach of Hinton (2002) that is popular for fitting a restricted Boltzmann
machine (Salakhutdinov et al., 2007), and the Hamiltonian Monte Carlo approach
under intractable likelihood as in Stoehr et al. (2019). All of these are
based on the identity $\nabla_{\theta}\log z(\theta)=\mathbb{E}_{Y\sim
p_{\theta}}\\{\nabla_{\theta}(\log q_{\theta}(Y))\\}$, for which the RHS can
be approximated by $N^{-1}\sum_{i=1}^{N}\nabla_{\theta}(\log
q_{\theta}(Y_{i\bullet}))$, where
$Y_{1\bullet},\ldots,Y_{N\bullet}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}p_{\theta}$. The need to sample iteratively from
$p_{\theta}(\cdot)$ lies at the heart of the lack of scalability for MCMC
based methods, in both fully Bayesian as well as approximate maximum
likelihood calculations for doubly intractable models. The choice of
$\phi=\mathrm{diag}(\theta)$, and subsequent sampling from $p_{\phi}$,
eliminates this difficulty, resulting in far greater scalability.

######  Remark 3.4.

Although not our main focus, the same construction may potentially be used to
design a simple MCMC-free sampler for $p_{\theta}(\cdot)$. See Supplementary
Section S.2 for details.

The following result provides further support on the well-behaved nature of
our Monte Carlo estimates, providing exponential concentration bounds.

######  Proposition 3.3.

Consider the class of PEGMs in (2) with $T_{j}(x_{j})=x_{j}$ and
$T_{jk}(x_{j},x_{k})=x_{j}x_{k}$ and suppose the $j$th sample space is bounded
in the interval $[a_{j},b_{j}]$ for $j=1,\ldots,p$. Set
$U=q_{\theta}(Y)/q_{\phi}(Y)$ and
$V=\nabla_{\theta_{jk}}q_{\theta}(Y)/q_{\phi}(Y)$ with $Y\sim p_{\phi}$. Then
for every $\theta$, both $U$ and $V$ as a function of $Y$ have bounded
differences in the sense of Boucheron et al. (2003, Corollary 1). Moreover,
there exists $C>0$ such that for every $t>0$,
$\mathbb{P}[|U-\mathbb{E}(U)|>t]\leq 2\exp\\{-t^{2}/4C\\}$, and
$\mathbb{P}[|V-\mathbb{E}(V)|>t]\leq 2\exp\\{-t^{2}/4C\\}$.

The proof is given in Supplementary Section S.1.2. Proposition 3.3 establishes
that for PEGMs with finite support, the specific choice of
$\phi=\text{diag}(\theta)$ not only leads to finite variance Monte Carlo
estimates of $z(\theta)$ and $\nabla_{\theta}z(\theta)$, but for this specific
class, the estimates are sub-Gaussian, which, in turn, implies exponential
concentration of the sample mean of $U,V$. The class of PEGMs with finite
support is quite rich. Indeed, the Ising model, Potts model, and the truncated
Poisson graphical model are all members of this family.

###  4 Applications to Fully Observed Intractable Pairwise Exponential Family
Graphical Models

####  4.1 MLE and Bayesian inference in low dimensions

Suppose we observe i.i.d. data $\mathbf{X}=(X_{1\bullet},\ldots,X_{n\bullet})$
where $X_{i\bullet}$’s are observations from a PEGM with density
$p_{\theta}(x)=q_{\theta}(x)/z(\theta)$. The maximum likelihood estimate
satisfies
$\hat{\theta}=\operatorname*{arg\,max}_{\theta\in\Theta}\ell(\theta)$ or
$\nabla_{\theta}\ell(\hat{\theta})=0$ where $\ell(\theta)=\sum_{i=1}^{n}\log
q_{\theta}(X_{i\bullet})-\log z(\theta)$. Denote $\log
q_{\theta}(\mathbf{X})=\sum_{i=1}^{n}\log q_{\theta}(X_{i\bullet})$. The
estimate can be obtained by implementing a gradient ascent algorithm subject
to the constraint $\theta\in\Theta$. We achieve this by the method of
projected gradient ascent, which is a special case of proximal methods under
box type constraints such as ours (p. 149, Parikh and Boyd, 2014). The
projection operation is defined as
$\mathcal{P}_{\Theta}(\theta^{*})=\operatorname*{arg\,min}_{\theta\in\Theta}\left\lVert\theta-\theta^{*}\right\rVert_{2}^{2}$.
Thus, the updates of the gradient ascent algorithm have the form:

$\displaystyle\theta^{(t+1)}=\mathcal{P}_{\Theta}(\tilde{\theta}^{(t+1)})$
$\displaystyle=\mathcal{P}_{\Theta}\left(\theta^{(t)}+\gamma\nabla_{\theta}\ell(\theta^{(t)})\right)=\mathcal{P}_{\Theta}\left(\theta^{(t)}+\gamma\dfrac{\nabla_{\theta}\,q_{\theta^{(t)}}(\mathbf{X})}{q_{\theta^{(t)}}(\mathbf{X})}-\gamma\dfrac{\nabla_{\theta}z(\theta^{(t)})}{z(\theta^{(t)})}\right),$
(6)

where $\gamma>0$ is the step size. For the Ising model,
$\Theta=\mathbb{R}^{p\times p}$, hence $\mathcal{P}_{\Theta}$ is simply the
identity map. However, for the Poisson graphical model,
$\Theta=\\{\theta:\theta_{jj}\in\mathbb{R},\theta_{jk}\leq 0,\,j\neq
k=1,\ldots,p\\},$ for which $\mathcal{P}_{\Theta}$ simply amounts to
thresholding any positive $\\{\tilde{\theta}^{(t+1)}\\}_{jk},j\neq k$,
computed via the gradient ascent step, to zero.

For Bayesian inference, consider proper priors $\pi(\theta)$ such that
$\int_{\Theta}\pi(d\theta)=1$ and assume $\nabla_{\theta}\pi(\theta)$ exists.
Then,
$U(\theta)=-\log\pi(\theta\mid\mathbf{X})=-\ell(\theta)-\log\pi(\theta)+\log
m(\mathbf{X})$ where
$m(\mathbf{X})=\int_{\Theta}\ell(\theta)\pi(\theta)d\theta$. Posterior
sampling is usually achieved by constructing a random walk Markov kernel that
transitions from the current state $\theta$ to a new candidate state
$\theta^{\prime}$. Existing algorithms for doing this (Liang, 2010; Møller et
al., 2006; Murray et al., 2012) are based on generating auxiliary samples from
the candidate $\theta^{\prime}$. However, it is well known that the
Hamiltonian Monte Carlo sampler could be more efficient in exploring a
complicated likelihood surface compared to random walk Metropolis kernels,
since the former uses gradient information for the MCMC updates (Neal et al.,
2011). This often leads to better exploration of the target density, and the
advantages could be substantial in high dimensions. One potential drawback is
the gradient may not always be available. Since we have an estimate of
$\nabla_{\theta}\ell(\theta)$ following Proposition 3.1, and hence of
$\nabla_{\theta}U(\theta)$, we consider a candidate $\theta^{\prime}$
constructed using reflected Hamiltonian dynamics (Betancourt, 2011), so that
$\theta^{\prime}\in\Theta$. This proposal is then finally accepted with a
Metropolis-Hastings correction (Neal et al., 2011, Equation 3.6). To complete
the algorithmic specifications, in low-dimensional situations, we use a flat
Gaussian prior (mean 0 and variance 100) on the parameters for Ising model.
For PGM, we use Gaussian priors on the diagonal elements of $\theta$ (same
parameters as Ising) and a flat exponential prior on the off-diagonal elements
(rate 0.01).

####  4.2 Penalized MLE and Bayesian inference in high dimensions

The more practically useful and statistically challenging situation arises
when $p$ is moderate to high-dimensional. We consider the problem of structure
learning under the assumption of a sparse underlying graph $G$. This problem
was considered by Meinshausen and Bühlmann (2006) within the context of
Gaussian graphical models, Ravikumar et al. (2010) addressed Ising models, and
Yang et al. (2012) for general PEGMs. All these approaches are based on the
idea of neighborhood selection, wherein neighbors of each node are estimated
separately. Essentially, these procedures rely on the pseudo-likelihood
$\prod_{j=1}^{p}p_{\theta}(X_{j}\mid X_{-j})$ with a suitable penalty. To our
knowledge, joint structure learning for general PEGMs has not been considered
before. Let $\left\lVert M\right\rVert_{1}=\sum_{j}\sum_{k}|M_{jk}|$ denote
the $\ell_{1}$-norm of the matrix $M$. We propose the following
$\ell_{1}$-penalized estimator for $\theta$:

$\hat{\theta}_{\lambda}=\operatorname*{arg\,min}_{\theta\in\Theta}-\ell(\theta)+\lambda\left\lVert\theta\right\rVert_{1}=\operatorname*{arg\,min}_{\theta\in\Theta}Q(\theta).$
(7)

For computing the estimator, a proximal optimization algorithm (Parikh and
Boyd, 2014) can be used with minor modifications from the updates in (6).
Indeed, the objective function is equivalent to
$Q(\theta)+\mathbbm{1}_{\theta\in\Theta}$. Since both the $\ell_{1}$-penalty
and the parameter space $\Theta$ are convex, proximal maps of the sum of the
two functions amount to the composition of the proximal maps of the functions
(Yu, 2013). This leads to batch updates of the form:

$\theta^{(t+1)}=\mathcal{P}_{\Theta}\,[\text{prox}_{\ell_{1}}(\theta^{(t)}+\gamma\nabla_{\theta}\ell(\theta^{(t)}))],$
(8)

where $\text{prox}_{\ell_{1}}(x)$ is the well-known soft-thresholding operator
applied to each element of $x$ (Parikh and Boyd, 2014). The tuning parameter
$\lambda$ is selected via cross-validation wherein the out-of-sample log-
likelihood is maximized. This usually results in better Frobenius norms but
poor structure learning. When structure learning is the main focus, we
estimate the underlying graph by averaging estimates $\hat{G}_{\lambda}$
obtained over a grid of values for $\lambda$. A more detailed discussion is
deferred to Supplementary Section S.3.

Fully Bayesian inference is also possible in the high-dimensional case. Prior
choices for sparse parameters is a well-developed field. One popular choice is
the family of scale mixtures of normal prior which is elicited as
$\theta_{jk}\mid\rho^{2}_{jk},\tau^{2}\overset{iid}{\sim}{\mathrm{N}}(0,\rho^{2}_{jk}\tau^{2}),\,\rho_{jk}\overset{iid}{\sim}f,\,\tau\sim
g$ for some choice of mixing distribution $f$ and $g$. These are
computationally beneficial (Bhattacharya et al., 2016) compared to the
discrete mixture priors (Mitchell and Beauchamp, 1988; George and McCulloch,
1993), without compromising statistical accuracy (Bhattacharya et al., 2015;
Chakraborty et al., 2020; Bhadra et al., 2017). Furthermore, posterior modes
under these priors can be thought of as a penalized estimator of $\theta$
under the penalty $-\log\pi(\theta)$. In particular, when $\tau=1$ and
$\rho_{jk}^{2}\mid\lambda\sim\text{Exp}(1/2\lambda^{2})$ (Park and Casella,
2008), then $-\log\pi(\theta)=\lambda\left\lVert\theta\right\rVert_{1}$, the
posterior mode of which is $\hat{\theta}_{\lambda}$ in (7). In our work, we
choose this prior augmented by the hyperprior
$\lambda\sim\text{Gamma}(a_{\lambda},b_{\lambda})$, i.e. if $\pi_{L}(\theta)$
denotes the joint prior over $\theta$, it can be specified hierarchically as
$\theta_{jk}\mid\rho^{2}_{jk},\lambda\overset{iid}{\sim}{\mathrm{N}}(0,\rho_{jk}^{2}\lambda^{-1}),\,\rho_{jk}^{2}\overset{iid}{\sim}\text{Exp}(1/2)$
and $\lambda\sim\text{Gamma}(a_{\lambda},b_{\lambda})$. Additionally, to
respect the restriction that $\theta\in\Theta$, we truncate the prior
$\pi(\theta)$ to $\\{\theta:\theta\in\Theta\\}$ so that our working prior is
$\pi(\theta)\propto\pi_{L}(\theta)\mathbb{I}_{\theta\in\Theta}$. An HMC
sampler can be designed aided by the scale mixture representation of
$\pi(\theta)$. Indeed, the conditional log-posterior
$U(\theta\mid\\{\rho_{jk}\\})$ is differentiable, and the HMC sampler can be
implemented conditional on the latent $\\{\rho_{jk}\\}$. The latent
$\\{\rho_{jk}\\}$’s can be updated given $\theta_{jk}$ using an inverse-
Gaussian distribution. Details are given in Supplementary Section S.3.

###  5 Applications to Partially Observed Graphical Models: The Restricted and
Full Boltzmann Machines

First consider the training of RBMs. The complete data model is Ising, which
is a PEGM, and thus the proposed methodology can be used without much effort.
To see this, consider the Expectation Maximization (EM) algorithm for
computing the maximum likelihood estimator. Standard calculations yield the
following updates for a gradient ascent algorithm for the weights between the
hidden and visible variables:

$\displaystyle\theta^{(t+1)}$
$\displaystyle=\theta^{(t)}+\gamma\left\\{\mathbb{P}(\mathbf{h}=\mathbf{1}\mid\mathbf{v},\theta=\theta^{(t)})\mathbf{v}^{T}-\nabla_{\theta}\log
z(\theta^{(t)})\right\\}$ (9) $\displaystyle=\theta^{(t)}+$
$\displaystyle\gamma\left\\{\mathbb{P}(\mathbf{h}=\mathbf{1}\mid\mathbf{v},\theta=\theta^{(t)})\mathbf{v}^{T}-\mathbb{E}_{(\mathbf{h},\mathbf{v})\sim
p_{\phi}}\left[\dfrac{\nabla_{\theta}e^{-E_{\theta}(\mathbf{v},\mathbf{h})}}{e^{-E_{\phi}(\mathbf{v},\mathbf{h})}}\right]\middle/\mathbb{E}_{(\mathbf{h},\mathbf{v})\sim
p_{\phi}}\left[\dfrac{e^{-E_{\theta}(\mathbf{v},\mathbf{h})}}{e^{-E_{\phi}(\mathbf{v},\mathbf{h})}}\right]\right\\},$

where the second equality follows from Proposition 3.1. For RBMs,
$\mathbb{P}(\mathbf{h}=\mathbf{1}\mid\mathbf{v},\theta=\theta^{(t)})=\prod_{k=p+1}^{p+m}\mathbb{P}(h_{k}=1\mid\mathbf{v},\theta=\theta^{(t)})$,
which is a product of sigmoid functions. Updates for the biases are similar.
Similar to the CD algorithm described in Section 2.2, the proposed algorithm
here also makes use of Monte Carlo sampling to estimate $\nabla_{\theta}\log
z(\theta)$. But, unlike CD, in this case samples $(\mathbf{h},\mathbf{v})\sim
p_{\phi}$ are independent, and estimates the correct gradient. A Gibbs sampler
is not needed.

Remarkably, the proposed framework allows a feasible training method for a
full BM as well; a problem long considered intractable. Define the parameter
matrix $\theta\in\mathbb{R}^{(p+m)\times(p+m)}$ as
$\theta=\begin{pmatrix}\theta_{\mathbf{pp}}&\theta_{\mathbf{pm}}\\\
\theta_{\mathbf{pm}}^{\mathrm{\scriptscriptstyle{T}}}&\theta_{\mathbf{mm}}\end{pmatrix}$,
where $\theta_{\mathbf{pp}}$ contains bias and interaction terms between the
visible variables, $\theta_{\mathbf{pm}}$ contains the interaction terms
between the hidden variables, and $\theta_{\mathbf{mm}}$ contains bias and
interaction terms between the hidden variables. At $\theta=\theta^{(t)}$, the
expected value of the complete data log-likelihood given the observed
$\mathbf{v}$ is:

$\displaystyle\mathbb{E}_{\theta^{(t)},\mathbf{v}}\\{\log
p_{\theta}(\mathbf{v},\mathbf{h})\\}=$
$\displaystyle\sum_{j}\theta_{jj}v_{j}+\sum_{k}\theta_{kk}\mathbb{E}_{\theta^{(t)},\mathbf{v}}(h_{k})+\sum_{j\neq
j^{\prime}}\theta_{jj^{\prime}}v_{j}v_{j^{\prime}}$ $\displaystyle+\sum_{k\neq
k^{\prime}}\theta_{kk^{\prime}}\mathbb{E}_{\theta^{(t)},\mathbf{v}}(h_{k}h_{k^{\prime}})+\sum_{j,k}\theta_{jk}v_{j}\mathbb{E}_{\theta^{(t)},\mathbf{v}}(h_{k})-\log
z(\theta),$

for $j,j^{\prime}=1,\ldots,p;\;k,k^{\prime}=p+1,\ldots,p+m$, where
$\mathbb{E}_{\theta^{(t)},\mathbf{v}}$ stands for the expectation taken with
respect to the density $p_{\theta_{t}}(\cdot\mid\mathbf{v})$. Thus, to train
this model using EM, one needs $\\{\mathbb{E}_{\theta_{t},\mathbf{v}}(h_{k})$,
$\mathbb{E}_{\theta^{(t)},\mathbf{v}}(h_{k}h_{k^{\prime}})\\}$, or more
generally, $\mathbb{E}_{\theta^{(t)},\mathbf{v}}\\{g(\mathbf{h})\\}$, where
$g(\mathbf{h})=g(h_{1},\ldots,h_{m})$ is a function of the hidden variables.
Fortunately, the conditional distribution of $\mathbf{h}\mid\mathbf{v}$, i.e.
$p_{\theta}(\cdot\mid\mathbf{v})$ here is also Ising. Indeed,

$p_{\theta_{\mathbf{h}\mid\mathbf{v}}}(\mathbf{h}\mid\mathbf{v})\propto\exp\left\\{\sum_{k}(\theta_{kk}+\sum_{j}\theta_{jk}v_{j})h_{k}+\sum\sum_{k\neq
k^{\prime}}\theta_{kk^{\prime}}h_{k}h_{k^{\prime}}\right\\}=\exp({-E_{\theta_{\mathbf{h}\mid\mathbf{v}}}(\mathbf{v},\mathbf{h})}),$

where the notation $\theta_{\mathbf{h}\mid\mathbf{v}}$ stands for the
parameter of the conditional Ising model for $\mathbf{h}\mid\mathbf{v}$. A
self-normalized importance sampling estimate of
$\mathbb{E}_{\theta,\mathbf{v}}\\{g(\mathbf{h})\\}$ is now presented.

######  Proposition 5.1.

Define
$w(\mathbf{v},\mathbf{h})=\exp\\{-E_{\theta_{\mathbf{h}\mid\mathbf{v}}}(\mathbf{v},\mathbf{h})+E_{\phi_{\mathbf{h}\mid\mathbf{v}}}(\mathbf{v},\mathbf{h})\\}$.
Then,

$\displaystyle\mathbb{E}_{\theta,\mathbf{v}}\\{g(h)\\}=\mathbb{E}_{(\mathbf{h}\mid\mathbf{v})\sim\phi_{\mathbf{h}\mid\mathbf{v}}}\left[g(\mathbf{h})\dfrac{p_{\theta_{\mathbf{h}\mid\mathbf{v}}}(\mathbf{h}\mid\mathbf{v})}{p_{\phi_{\mathbf{h}\mid\mathbf{v}}}(\mathbf{h}\mid\mathbf{v})}\right]=\dfrac{\mathbb{E}_{(\mathbf{h}\mid\mathbf{v})\sim\phi_{\mathbf{h}\mid\mathbf{v}}}\left[g(\mathbf{h})w(\mathbf{v},\mathbf{h})\right]}{\mathbb{E}_{(\mathbf{h}\mid\mathbf{v})\sim\phi_{\mathbf{h}\mid\mathbf{v}}}\left[w(\mathbf{v},\mathbf{h})\right]},$

where a natural choice for $\phi_{\mathbf{h}\mid\mathbf{v}}$ is
$\text{diag}(\theta_{\mathbf{h}\mid\mathbf{v}})$.

###### Proof.

Clearly,

$\displaystyle\mathbb{E}_{(\mathbf{h}\mid\mathbf{v})\sim\phi_{\mathbf{h}\mid\mathbf{v}}}\left[g(\mathbf{h})\dfrac{p_{\theta_{\mathbf{h}\mid\mathbf{v}}}(\mathbf{h}\mid\mathbf{v})}{p_{\phi_{\mathbf{h}\mid\mathbf{v}}}(\mathbf{h}\mid\mathbf{v})}\right]=\mathbb{E}_{(\mathbf{h}\mid\mathbf{v})\sim\phi_{\mathbf{h}\mid\mathbf{v}}}\left[g(\mathbf{h})w(\mathbf{v},\mathbf{h})\right]\left[\dfrac{z(\theta_{\mathbf{h}\mid\mathbf{v}})}{z(\phi_{\mathbf{h}\mid\mathbf{v}})}\right]^{-1}.$

The proof follows by applying (5) to the second term on the right hand side of
the above. ∎

######  Remark 5.1.

Variance of this self-normalized importance sampling estimate is finite when
$\phi_{\mathbf{h}\mid\mathbf{v}}=\text{diag}(\theta_{\mathbf{h}\mid\mathbf{v}})$.
Indeed, the estimator $Z$ is of the form $Z=X/Y$. Thus, the Cauchy-Schwarz
inequality gives
$\text{Var}_{\phi_{\mathbf{h}\mid\mathbf{v}}}(Z)\leq\mathbb{E}_{\phi_{\mathbf{h}\mid
v}}(X^{2}/Y^{2})\leq\mathbb{E}_{\phi_{\mathbf{h}\mid
v}}(X^{2})\mathbb{E}_{\phi_{\mathbf{h}\mid v}}(1/Y^{2})$. Since
$g(\mathbf{h})$ is bounded for $g(\mathbf{h})=h_{k}$ and
$g(\mathbf{h})=h_{k}h_{k^{\prime}};\,k,k^{\prime}=p+1,\ldots,p+m$; another
application of Cauchy–Schwarz inequality along with Proposition 3.2 proves
that $\mathbb{E}_{\phi_{\mathbf{h}\mid v}}(X^{2})<\infty$. The random
denominator $1/Y$ is of the form
$\exp\\{E_{\theta_{\mathbf{h}\mid\mathbf{v}}}(\mathbf{v},\mathbf{h})-E_{\phi_{\mathbf{h}\mid\mathbf{v}}}(\mathbf{v},\mathbf{h})\\}$.
That the variance of this quantity is finite follows from Proposition 3.2.

######  Remark 5.2.

The key difference between training a BM and an RBM is that when using (9) for
RBMs, one Monte Carlo estimate of $z(\theta)$ is sufficient. However, for BMs,
Monte Carlo estimates of $\mathbb{E}_{\theta,\mathbf{v}}\\{g(\mathbf{h})\\}$
are needed. Moreover, the parameters of $\mathbf{h}\mid\mathbf{v}$ are
dependent on $\mathbf{v}$. A CD-type algorithm can in principle be
implemented, but would be immensely expensive since one needs to run an
$m$-dimensional Gibbs sampler for each $\mathbf{v}$. On the other hand, with
the proposed method, batch sampling is possible using the specific choice
$\phi_{\mathbf{h}\mid\mathbf{v}}=\text{diag}(\theta_{\mathbf{h}\mid\mathbf{v}})$.

######  Remark 5.3.

The training of BM also allows one to compute an estimate of the marginal
likelihood of the visible variables. Indeed, for all
$\mathbf{h}\in\\{0,1\\}^{m}$, the following relation holds: $\log
p_{\hat{\theta}}(\mathbf{v})=\log p_{\hat{\theta}}(\mathbf{v},\mathbf{h})-\log
p_{\hat{\theta}}(\mathbf{h}\mid\mathbf{v})$. Since both
$(\mathbf{v},\mathbf{h})$ and $(\mathbf{h}\mid\mathbf{v})$ are Ising, this
quantity is computable using techniques developed here. The marginal
likelihood may then be used for model selection or hyperparameter tuning
purposes.

The proposed framework also encompasses the DBM where there are multiple
layers of hidden variables but the underlying graph can be rearranged in a
two-layer graph (Aggarwal, 2018), with all odd layer variables constituting
one half of the graph and even layers constituting the other half, with the
restriction that connections between variables in the same layer are not
allowed; see Figure 1. Thus, the DBM can be thought of as a relaxation of the
BM and a similar training algorithm can be implemented in principle to train
the model. We omit the details. In the context of DBM, Salakhutdinov and
Murray (2008) proposed a related annealed importance sampling scheme (Neal,
2001) where a similar choice of $\phi=\text{diag}(\theta)$ is proposed for
computing marginal probabilities. But they do not discuss full likelihood
inference, and still train their model via CD.

###  6 Theoretical Properties

In this section we investigate the properties of the likelihood-based
estimates, including the Bayesian posterior, specifically in high-dimensional
settings, since the properties of these procedures in low-dimensional problems
are fairly well-understood. We first focus on the $\ell_{1}$-penalized
estimator in (7). Since its introduction in the regression context
(Tibshirani, 1996), $\ell_{1}$-penalized estimators have been used in many
other settings; for instance in sparse covariance estimation problems (Rothman
et al., 2010), reduced-rank problems (Yuan et al., 2007), generalized linear
models (Ravikumar et al., 2010; Yang et al., 2012), among others. These
estimators have been shown to be consistent with the optimum rate of
convergence for sparse problems. We establish these properties of the
estimator in (7) for PEGMs under fairly general assumptions, drawing parallels
with the existing literature. In the sequel, we write $\theta_{0}\in\Theta$ as
the true data-generating parameter, and $\mathbf{S}=\\{(j,k):\theta_{0,jk}\neq
0,j,k=1,\ldots,p\\}$, $|\mathbf{S}|=s$. We start by stating the assumptions.

######  Assumption 6.1 (Dimension).

We assume that $\log p=o(n)$, $s=o(p)$ and $s\log p=o(n)$ where
$a_{n}=o(b_{n})$ denotes $a_{n}/b_{n}\to 0$ as $n\to\infty$.

######  Assumption 6.2 (Regularity).

The space $\Theta$ of feasible parameters is open and $\theta_{0}$ is in the
interior of $\Theta$. We also assume the PEGM is sufficiently regular, i.e.,
for some $\beta,\rho>0,$

$\lambda_{\min}\left(\dfrac{\partial^{2}}{\partial\theta^{2}}A(\theta_{0})\right)\geq\beta,\quad\left[\sum_{j,k,l=1}^{p}\left(\dfrac{\partial^{3}}{\partial\theta^{3}}A(\theta_{0})\right)_{jkl}^{2}\right]^{1/2}\leq\rho,$

where $\lambda_{\min}$ denotes the minimum eigenvalue.

######  Assumption 6.3 (Bounded moments).

For all $j,k=1,\ldots,p,$ assume,
$\mathbb{E}_{\theta_{0}}[|T(X_{j})|]\leq\kappa_{1}$,
$\mathbb{E}_{\theta_{0}}[T^{2}(X_{j})]\leq\kappa_{2}$, and
$\mathbb{E}_{\theta_{0}}[|T(X_{j},X_{k})|]\leq\kappa_{3}$. Also, assume that,

$\max_{u:|u|<1}\dfrac{\partial^{2}}{\partial\theta_{j}^{2}}A(\theta_{0}+ue_{j})\leq\kappa_{4},$

where $e_{j}\in\mathbb{R}^{p^{2}}$ is the unit vector corresponding to the
index $j$. Assume, furthermore,

$\max_{\eta:|\eta|<1}\dfrac{\partial^{2}}{\partial\eta^{2}}\bar{A}_{jk}(\theta_{0};\eta)\leq\kappa_{5},$

where,

$\bar{A}_{jk}(\theta;\eta)\coloneqq\log\int_{\mathcal{X}}\exp\left\\{\eta
T(X_{j},X_{k})+\sum_{j=1}^{p}\theta_{j}T(X_{j})+\sum_{\begin{subarray}{c}j\neq
k=1\end{subarray}}^{p}\theta_{jk}T(X_{j},X_{k})+\sum_{j=1}^{p}C(X_{j})\right\\}\mathrm{d}x.$

Assumption 6.2 imposes the condition that the covariance matrix of the model
is not singular and can be reasonably well approximated within a small
neighborhood of $\theta_{0}$. Assumption 6.3 controls the size of the
deviations of the sample moments from the population moments. Then we have the
following result, with a proof in Supplementary Section S.1.4.

######  Theorem 6.4.

Suppose we observe $n$ i.i.d copies
$\mathbf{X}=(X_{1\bullet},\ldots,X_{n\bullet})\sim PEGM(\theta)$. Let
$r_{n}=\sqrt{\frac{s\log p}{n}}$ and $\hat{\theta}_{\lambda}$ be a solution to
(7). Under Assumptions 6.1–6.3, for $\lambda\asymp\sqrt{\frac{\log p}{n}}$,
there exists sufficiently large $M>0$ such that as $n\to\infty$, one has:

$\mathbb{P}_{\theta_{0}}\left[\left\lVert\hat{\theta_{\lambda}}-\theta_{0}\right\rVert_{F}^{2}>Mr_{n}\right]\to
0,$

where $\left\lVert\cdot\right\rVert_{F}$ denotes the Frobenius norm.

The rate of convergence of the $\ell_{1}$-penalized estimator matches with
existing rates in sparse parameter estimation theory; see for example Rothman
et al. (2010).

Frequentist evaluation of Bayesian posteriors is typically studied through the
lens of posterior consistency and contraction rates (Ghosal and Van der Vaart,
2017). In models with increasing number of parameters, a standard condition
for posterior consistency is positive prior support on shrinking
Kullback–Liebler balls around the true distribution. Suppose $p_{\theta_{0}}$
is the true distribution of the data and $p_{\theta}$ is the fitted
distribution. Let $\pi(\theta)$ be the prior distribution on $\theta$ and
$\text{KL}(p_{\theta_{0}},p_{\theta})$ denote the Kullback–Liebler divergence
between $p_{\theta_{0}}$ and $p_{\theta}$. In Proposition S.1.2, we show when
$p_{\theta_{0}}$ and $p_{\theta}$ are PEGMs, then the set
$\\{\theta:\text{KL}(p_{\theta_{0}},p_{\theta})<\epsilon\\}$ for any
$\epsilon>0$ can be characterized in terms of
$\left\lVert\text{vech}(\theta_{0})-\text{vech}(\theta)\right\rVert_{2}^{2}$,
where $\text{vech}(\theta)$ stands for the half-vectorization of a symmetric
matrix $\theta$. As a result, any prior distribution that has positive mass on
the set
$\\{\theta:\left\lVert\text{vech}(\theta_{0})-\text{vech}(\theta)\right\rVert_{2}^{2}<\epsilon\\}$
for $\epsilon>0$ yields a consistent posterior. Properties of global-local
shrinkage priors around small Euclidean neighborhoods of a sparse parameter
are well studied; see, e.g., Bhattacharya et al. (2016). We next establish
posterior consistency for PEGMs under prior
$\pi(\theta)\propto\pi_{L}(\theta)\mathbb{I}_{\theta\in\Theta}$ where
$\pi_{L}(\theta)$ is defined at the end of Section 4.2.

######  Theorem 6.5.

Suppose we observe $n$ i.i.d copies
$\mathbf{X}=(X_{1\bullet},\ldots,X_{n\bullet})\sim PEGM(\theta)$ and $\theta$
is endowed with prior
$\pi(\theta)\propto\pi_{L}(\theta)\mathbb{I}_{\theta\in\Theta}$ where
$\pi_{L}(\theta)$ is the marginal prior obtained under hierarchy
$\theta_{j}\mid\lambda\overset{iid}{\sim}\text{Laplace}(\lambda^{-1})$ and
$\lambda\sim\text{Gamma}(a_{\lambda},b_{\lambda})$. If, $a_{\lambda}=O(s\log
p)$ and $b_{\lambda}=O\left(s\sqrt{\frac{s\log p}{np}}\right)$, then, under
Assumptions 6.1–6.3, for every sufficiently small $\epsilon>0$, we have in
$\mathbb{P}_{\theta_{0}}-$probability:

$\pi\left[\theta:\left\lVert\theta-\theta_{0}\right\rVert_{F}>\epsilon\mid\mathbf{X}\right]\to
0.$

See Supplementary Section S.1.5 for a proof. We next turn our attention to
partially observed models, specifically the RBM. Suppose we observe $n$ i.i.d
variables
$\mathbf{V}=(\mathbf{v}_{1\bullet},\ldots,\mathbf{v}_{n\bullet})^{\mathrm{\scriptscriptstyle{T}}}$
with support on $\\{0,1\\}^{p}$. When an RBM is fitted to the observed data,
maximum likelihood inference involves optimizing the marginal likelihood (4).
The iterative scheme (9) essentially performs Monte Carlo EM algorithm on the
complete-data log-likelihood. Let
$Q(\theta;\psi)=\mathbb{E}_{\psi}[-E_{\theta}(\mathbf{v},\mathbf{h})]-\log
z(\theta)=g(\mathbf{v};\theta,\psi)-\log z(\theta)$, and
$Q^{\star}(\theta;\psi)=g(\mathbf{v};\theta,\psi)-\widehat{\log z(\theta)}$,
its Monte Carlo estimate. Our next result shows EM sequences based on
$Q(\theta;\psi)$ and $Q^{\star}(\theta;\psi)$ are not far apart. The proof is
given in Supplementary Section S.1.6.

######  Theorem 6.6.

Let the parameter space $\Theta$ be compact. Suppose $\\{\theta^{(t)}\\}$ is a
generalized EM sequence obtained with initial point $\theta^{(0)}$, step sizes
$\gamma_{t}$, and let this sequence converge to a stationary point
$\underline{\theta}$ of (4), i.e. $\nabla\log
p_{\theta}(\mathbf{v})|_{\theta=\underline{\theta}}=0$. Also, assume that
$\gamma_{t}\lambda_{\max}[\nabla^{2}_{\theta}g(\mathbf{v};\theta;\psi)]+\gamma_{t}\lambda_{\max}[\nabla^{2}\log
z(\theta)]\leq 1$ for all $t\geq 1$ and $\theta\in\Theta$, and
$\sum_{t=1}^{\infty}\gamma_{t}=0$. Let $\theta_{\star}^{(t)}$ denote a
sequence with the same initial point and step-sizes obtained using
$Q^{\star}(\theta;\psi)$. Then for sufficiently large $t$,

$\left\lVert\theta^{(t)}-\theta_{\star}^{(t)}\right\rVert_{F}=O\left(\frac{1}{t^{1+\eta}}\right),$

for some $\eta>0$,
$\mathbb{P}_{\theta_{\star}^{(1)}}\otimes\cdots\otimes\mathbb{P}_{\theta_{\star}^{(t)}}$-almost
surely.

###  7 Numerical Experiments

We compare the performances of (maximum) likelihood based methods (MLE),
including a full Bayes method, to that of (maximum) pseudo-likelihood (MPLE)
based methods for parameter inference in PEGMs. We consider the Ising model
and the PGM as examples of PEGMs. We design our experiment in a low-
dimensional (LD), high-dimensional (HD), and ultra high-dimensional (UHD)
setting for both models. In the low-dimensional setting, our goal is to
demonstrate the efficiency in terms of uncertainty quantification of a full
likelihood based analysis since both MLE and PMLE are consistent. The focus in
the (ultra) high-dimensional setting is to understand the performance of
_penalized_ versions of these methods, PMLE and PMPLE, respectively, in terms
of graph structure learning. We set the dimension $p=3,5$ (low-dimensional)
and $p=20,30,40,50$ (high-dimensional) and $p=100,200,300$ (ultra high-
dimensional). For both low and high-dimensional settings, the sample size is
set as $n=100$. In the ultra high-dimensional setting we set $n=2p$. The true
parameter $\theta_{0}$ for all the cases is generated randomly, first by
sampling $z_{jk}\sim\text{Bernoulli}(\omega)$ and then setting
$\theta_{0,jk}\mid[z_{jk}=0]=0$ and $\theta_{0,jk}\mid[z_{jk}=1]=\eta$. In the
low-dimensional setting we set $\omega=0.9$, $\eta=-0.8$, whereas for the
other two cases we set $\omega=0.05$, $\eta=-3$.

In the LD setting, we compare the methods using the Frobenius norm
$||\hat{\theta}-\theta_{0}||_{F}^{2}$, average coverage probabilities of 95%
confidence intervals, and their corresponding widths. For the Bayes method, we
choose $\hat{\theta}$ as the posterior mean using the Laplace scale mixture
prior described previously. In the HD setting, we focus on the Matthews
correlation coefficient (MCC) (Chicco et al., 2021) between the estimated and
the true graph and $||\hat{\theta}-\theta_{0}||_{F}^{2}$. The MCC values
ranges between $-1$ and $+1$, with $-1$ indicating the worst performance and
$+1$ the best. The advantage of MCC as a measure is that its value is close to
$+1$ if and only if the estimated graph estimates zero and non-zero edges
between nodes correctly at the same time. We only consider PMLE and PMPLE
here, since the Bayesian posterior mean with a continuous shrinkage prior is
not sparse. For both of these methods, the estimated graph is produced
following the procedure described in Supplementary Section S.3. Additionally,
we report $||\hat{\theta}-\theta_{0}||_{F}^{2}$, where for PMLE, the estimate
$\hat{\theta}$ correspond to the one obtained using cross-validation. Finally,
in the UHD setting, we only consider MCC of PMLE and PMPLE. We summarize the
Frobenius norm and MCC for the Ising model and PGM in Table 7 for 10
replications. The average coverage probabilities of $95\%$ confidence
intervals and their corresponding widths for both models in the LD settings
are in Supplementary Table S.4.1. Full likelihood-based methods achieve lower
Frobenius norms in all the settings considered, for both Ising and PGM,
compared to their PMPLE counterparts. In terms of graph reconstruction, the
PMLE estimator (7) outperforms the PMPLE-based methods with significantly
higher MCC values. Figure 2 shows the heatmap comparing the data generating
parameter $\theta_{0}$ and the estimates obtained from PMPLE and PMLE for the
Ising model with $p=20$; the PMLE estimate corresponds to $\lambda=7$.
Clearly, $\hat{\theta}_{PMLE}$ effectively captures the underlying structure
and also yields estimates that are closer to the truth. Finally, we note that
penalization was not used in LD, and the fully Bayesian procedure was
computationally prohibitive in UHD. These cases are marked with a ‘–’ in Table
7. Additional simulation results on the partially observed case, comparing the
representational powers of BM and RBM, are presented in Supplementary Section
S.4.2.

ccc——ccc—cc——ccc—cc Ising  PGM
Setting $p$ $n$ $||\theta_{0}-\hat{\theta}||_{F}^{2}$ MCC
$||\theta_{0}-\hat{\theta}||_{F}^{2}$ MCC
PMLE PMPLE Bayes PMLE PMPLE PMLE PMPLE Bayes PMLE PMPLE
LD 3 100 1.24 (0.52) 6.48 (6.05) 2.09 (1.62) – – 0.38 (0.09) 5.83 (5.65) 0.70
(0.02) – –
5 2.07 (0.85) 16.23 (6.42) 1.68 (0.28) – – 0.98 (0.27) 19.57 (5.21) 1.37
(0.30) – –
HD 20 100 9.18 (0.61) 9.56 (0.57) 10.58 (0.27) 0.67 (0.05) 0.48 (0.01) 11.60
(0.43) 13.62 (0.17) 13.91 (0.05) 0.83 (0.07) 0.24 (0.02)
30 12.85 (0.99) 13.94 (0.61) 16.24 (0.38) 0.78 (0.02) 0.49 (0.02) 15.27 (1.25)
20.43 (0.12) 20.94 (0.01) 0.71 (0.05) 0.22 (0.02)
40 18.61 (1.31) 27.38 (0.11) 24.33 (0.29) 0.71 (0.05) 0.46 (0.03) 21.00 (1.03)
27.06 (0.07) 27.26 (0.01) 0.73 (0.04) 0.25 (0.02)
50 26.89 (0.50) 35.38 (0.09) 32.10 (0.21) 0.71 (0.03) 0.39 (0.03)30.57 (1.32)
34.36 (0.15) 33.85 (0.01) 0.61 (0.04) 0.22 (0.02)
UHD 100 200 54.35 (0.93) 69.71 (0.04) – 0.80 (0.02) 0.44 (0.01) 59.26 (0.62)
67.43 (0.07) – 0.52 (0.02) 0.27 (0.01)
200 400 128.97 (0.14) 137.16 (0.02) – 0.70 (0.01) 0.41 (0.01) 131.84 (0.03)
134.70 (0.05) – 0.47 (0.01) 0.34 (0.01)
300 600 195.97 (0.06) 202.21 (0.02) – 0.56 (0.01) 0.24 (0.01)201.11 (0.03)
203.73 (0.05) – 0.50 (0.01) 0.33 (0.01)

Table 1: Mean (sd) of MSE = $||\hat{\theta}-\theta_{0}||_{F}^{2}$ and Matthews
correlation coefficient (MCC) across 10 datasets for penalized maximum
likelihood estimation (PMLE), penalized maximum pseudo-likelihood estimation
(PMPLE) and the Bayes methods for Ising and PGM.

Figure 2: True parameter $\theta_{0}$ and the corresponding PMLE and PMPLE
estimates $\hat{\theta}$ for the Ising model with $(n,p)=(100,20)$ and
$\lambda=7$.

###  8 Data Analysis

We provide demonstrations on data sets arising from disparate applications,
intentionally so chosen, to demonstrate the scope of the methodology.

####  8.1 Movie ratings network via Ising model

We demonstrate the merits of a likelihood-based analysis for the Ising Model
by analyzing a real-world movie rating dataset. The original $Movielens\ 25M$
dataset collected by Grouplens consists of movie ratings, ranging from 0 to 5
with increments of 0.5, for 62,000 movies by 162,000 users
(https://grouplens.org/datasets/movielens/). We select the most frequently
rated 50 movies ($p=50$) from the same group of $n=303$ viewers. We code
observations that received a rating of 5 as ‘1’, signifying complete
satisfaction with the film, and assign ‘0’ to ratings of 4.5 or below (Roach
and Gao, 2023). Supplementary Table S.3 provides the legend of the movie
titles, genres, and years of release. The binary variables $X_{ij}$ indicates
the preference of the $i$th viewer for the $j$th movie, where $i=1,\ldots,303$
and $j=1,\ldots,50$. Specifically, $X_{ij}=X_{ik}=1$, where $j\neq k$, implies
that the $i$th viewer has strong preference for both movies $j$ and $k$. Other
combinations can also be interpreted similarly. We model the distribution of
$X$ using $\text{Ising}(\theta)$ where $\theta\in\mathbb{R}^{p\times p}$.
Diagonal elements of $\theta$ model the propensity of a movie being rated 1,
whereas off-diagonal elements capture the direction of association between
movies $(j,k)$ conditional on other movies.

We estimate $\theta$ using (7) and build the graph structure using the stable
selection algorithm described in the Supplementary Section S.3. For stable
selection, we select a grid of values of $\lambda=5,10,\dots,100$ and set
threshold $\pi_{thr}=0.6$. The resulting movie rating network with
abbreviations of the movie titles is plotted in Supplementary Figure S.2,
where the edge widths are proportional to the stable selection probability,
and the node sizes are proportional to the degree. Our analysis reveals that
Star Wars IV and Star Wars V have the most connections, partly due to their
appeal to a broader audience, see Supplementary Table S.4. Lord of the Rings:
The Fellowship of the Ring has the second most connections. These two movie
franchise also form two separate cliques within the graph, see Supplementary
Figure S.2. A more detailed summary of the results, in terms of the most
positive and negative connections, is provided in Table 2. Movies within the
Lord of the Rings and Star Wars franchises have strong positive connections,
whereas Inception (2010) and Batman (1989) have the strongest negative
connection. We remark here that Batman (1989) is directed by Tim Burton,
unlike the popular Batman movie franchise of the 2000s directed by Christopher
Nolan, who also directed Inception (2010). A heatmap of
$\hat{\theta}_{\lambda}$ for $\lambda=10$ is in Supplementary Figure S.3.

Positive Edge | $\hat{\theta_{jk}}$ | Negative Edge | $\hat{\theta_{jk}}$
---|---|---|---
Lord of the Rings (2003) - Lord of the Rings (2001) | 0.77 | Inception (2010) - Batman (1989) | -0.95
Lord of the Rings (2002) - Lord of the Rings (2001) | 0.72 | Terminator (1984) - Dances with Wolves (1990) | -0.80
Lord of the Rings (2003) - Lord of the Rings (2002) | 0.61 | True Lies (1994) - Star Wars: Episode IV (1977) | -0.53
Star Wars: Episode V (1980) - Star Wars: Episode IV (1977) | 0.56 | Memento (2000) - Fugitive (1993) | -0.49
Terminator (1984) - Terminator 2: Judgment Day (1991) | 0.55 | Inception (2010) - Terminator (1984) | -0.48
Star Wars: Episode VI (1983) - Star Wars: Episode IV (1977) | 0.53 | Dances with Wolves (1990) - Twelve Monkeys (1995) | -0.44
Star Wars: Episode VI (1983) - Star Wars: Episode V (1980) | 0.48 | Independence Day (1996) - Batman (1989) | -0.44
Raiders of the Lost Ark (1981) - Star Wars: Episode V (1980) | 0.40 | Godfather (1972) - Independence Day (1996) | -0.42
Godfather (1972) - Schindler’s List (1993) | 0.25 | Inception (2010)-Independence Day (1996) | -0.41
Raiders of the Lost Ark (1981) - Star Wars: Episode IV (1977) | 0.24 | Braveheart (1995) - Toy Story (1995) | -0.39

Table 2: Top 10 positive and negative interactions in the Movie Ratings
Network.

####  8.2 Breast cancer microRNA network via Poisson graphical model

We next apply a full-likelihood ($\hat{\theta}_{PMLE}$) and pseudo-likelihood
($\hat{\theta}_{PMPLE}$) analysis to the task of inferring dependency among
count-valued variables that represent microRNA (miRNA) expression profiles of
the $353$ genes. This level III RNA-Seq data is collected from 445 breast
invasive carcinoma (BRCA) patients ($n=445$) from the Cancer Genome Atlas
(TCGA) project (TCGA, 2012) and is available in the R package XMRF. The
observed values are counts of sequencing reads mapped back to a reference
genome. Raw counts are highly skewed and zero-inflated. We use the processSeq
function in the XMRF package to adjust for sequencing depth, mitigate
overdispersion, and normalize the counts. Since the PGM can only model
negative conditional associations, we select a subset of $p=101$ genes with
negative Pearson correlation between their counts.

We select a grid of $\lambda_{PMLE}=[200,240,280,\dots,1000]$, and set the
threshold $\pi_{thr}=0.6$ to implement the stable selection algorithm. For
pseudo-likelihood, we perform node-wise $\ell_{1}$ penalized Poisson
generalized linear model (GLM) using glmnet (Friedman et al., 2010) in R with
a negative constraint on coefficients and consider a grid of 20 evenly spaced
values in the interval $[5\sqrt{{\log p}/{n}},\;7\sqrt{{\log p}/{n}}]$ for the
penalty parameter $\lambda_{PMPLE}$. We then construct the final estimator
$\hat{\theta}_{\lambda,PMPLE}$ for a given $\lambda$ by setting
$\hat{\theta}_{\lambda,PMPLE}=(\hat{\theta}_{\lambda}+\hat{\theta}^{\mathrm{\scriptscriptstyle{T}}}_{\lambda})/2$,
where $\hat{\theta}_{\lambda}$ represents the matrix formed by coefficients
from Poisson GLMs. The sequence of estimates $\hat{\theta}_{\lambda,PMPLE}$
are then used for model selection following the stable selection procedure
with $\pi_{thr}=0.6$.

Estimated networks from a full-likelihood and pseudo-likelihood based analysis
are comparable. Supplementary Figures S.4 and S.5 show the estimated networks
from a likelihood and pseudo-likelihood analysis, respectively, where blue
edges are common to both networks. Let
$S_{PMLE}:=\\{(j,k):\hat{\theta}_{PMLE,jk}\neq 0,1\leq j\leq k\leq p\\}$ and
$S_{PMPLE}:=\\{(j,k):\hat{\theta}_{PMPLE,jk}\neq 0,1\leq j\leq k\leq p\\}$.
The cardinality of the intersection $|S_{PMLE}\cap S_{PMPLE}|=372$ is about
$70\%$ of $|S_{PMLE}|$. Moreover, the genes with the top numbers of
connections in the PMLE-based network also have high degrees of rank in the
PMPLE-based network, as shown in Table 3. Both networks successfully identify
the edge between the CDH1 (located on chromosome 16q) and the BCL9 (located on
chromosome 1q). Existing studies reported that the CDH1 is often
underexpressed, while the BCL9 is typically overexpressed, caused by frequent
chromosomal aberrations involved in breast cancer pathogenesis (Privitera et
al., 2021). Additionally, the top connected genes in Table 3, such as WHSC1L1,
TMPRSS2, CDH1, SLC34A2 and EP300, are all known key players in breast cancer
(Xiao et al., 2022; Kim et al., 2021).

Gene | WHSC1L1 | TMPRSS2 | CDH1 | SLC34A2 | FGFR2 | EP300 | RNF43 | NUP98 | ACSL3 | MAP2K1
---|---|---|---|---|---|---|---|---|---|---
$d_{PMLE}$ (Ranking) | 52 (1) | 36 (2) | 24 (3) | 21 (4) | 21 (5) | 20 (6) | 19 (7) | 19 (8) | 19 (9) | 19 (10)
$d_{PMPLE}$ (Ranking) | 33 (3) | 15 (43) | 34 (2) | 26 (7) | 19 (25) | 17 (32) | 20 (21) | 17 (32) | 14 (49) | 7 (84)

Table 3: Top 10 genes ranked by degrees from the full-likelihood based
($d_{PMLE}$) and pseudo-likelihood based ($d_{PMPLE}$) estimates for Breast
Cancer data.

####  8.3 Analysis of MNIST data via restricted and full Boltzmann machines

We analyze the MNIST data (Deng, 2012) using RBM and BM where we train the
model with the proposed algorithm. The data consists of images of handwritten
digits ($0$–$9$). The number of training samples for the $10$ digits vary
within the range $n=3000$ to $n=4000$. We consider images of 15$\times$15
resolution. Each image is then converted to a binary vector of length $p=225$,
with on pixels are coded as 1 and off pixels are coded as 0. A fitted RBM or
BM model on these images approximates the distribution of this binary vector
over $\\{0,1\\}^{225}$. Furthermore, since both RBM and BM are generative
models, the trained models can be used to produce synthetic data based on held
out test images of these digits. This is known as reconstruction. A well-
trained model will generally capture the global distribution of the pixels for
each digit, and be able to reproduce subject-specific variations within this
global distribution. In Figure 3, we show the reconstructed images of test
images for all the digits $0$–$9$ with $m=50$ hidden variables for both RBM
and BM. Additionally, we compare the performance of reconstructed test images
in terms of Brier loss. The Brier loss measures how well-calibrated the
reconstructed probabilities are. This is done on test images of all the
digits. Given a test image $\mathbf{v}$, and reconstructed probabilities
$\hat{\mathbf{v}}$, we compute the Brier loss as

$B=(n_{t}p)^{-1}\sum_{i=1}^{n_{t}}\sum_{j=1}^{p}(\mathbf{v}_{ij}-\hat{\mathbf{v}}_{ij})^{2},$

where $n_{t}$ is the number of test images. The results are reported in Table
4. Performance of the likelihood trained RBM (RBM-MLE) is better than CD
trained RBM for all digits. Also, with the same number of hidden variables,
BM-MLE achieves a (marginally) lower reconstruction error compared to RBM-MLE
in most cases.

Figure 3: Observed and reconstructed images of handwritten digits. The first
row corresponds to observed of test images of digits $0$–$9$. The second and
third rows show the reconstructed images from the trained RBM-MLE and BM-MLE
model. All models are fitted with $m=50$ hidden variables.

Digit | $\mathbf{0}$ | $\mathbf{1}$ | $\mathbf{2}$ | $\mathbf{3}$ | $\mathbf{4}$ | $\mathbf{5}$ | $\mathbf{6}$ | $\mathbf{7}$ | $\mathbf{8}$ | $\mathbf{9}$
---|---|---|---|---|---|---|---|---|---|---
RBM-MLE | 0.39 | 0.17 | 0.31 | 0.45 | 0.32 | 0.40 | 0.43 | 0.35 | 0.41 | 0.32
RBM-CD | 0.44 | 0.36 | 0.42 | 0.48 | 0.43 | 0.43 | 0.47 | 0.39 | 0.41 | 0.40
BM-MLE | 0.42 | 0.15 | 0.29 | 0.46 | 0.30 | 0.41 | 0.43 | 0.32 | 0.39 | 0.28

Table 4: Brier loss for RBM and BM trained with proposed algorithm and RBM
trained with CD. The results demonstrate the superiority of full-likelihood
based training over CD-based training.

###  9 Conclusions and Future Works

The overarching focus of the current paper has been on the development of
computationally feasible likelihood-based inference procedures in intractable
graphical models, and in particular for models that were long considered
impermeable to such inference, e.g., Boltzmann machine-based architectures.
Our theoretical results on likelihood-based inference complement the
methodological contributions. A number of future investigations can be
considered. First, we only considered first order gradient-based methods for
optimization. However, since we have a Monte Carlo estimate of the gradient, a
quasi-Newton approach such as BFGS (see, e.g., Broyden, 1970) appears
straightforward to develop and may have better convergence speed in practice
compared to a first order gradient descent scheme. Second, the intractability
of the normalizing constant also affects the computation of evidence or
marginal likelihood, and makes model selection and comparison difficult.
Although not explicitly considered in this work, it is a natural future step
is to consider model selection problems in exponential family graphical models
and to consider multiple connected graphs, where we envisage that the
techniques developed in this paper will prove fruitful.

SUPPLEMENTARY MATERIAL

(a) _Supplementary Text:_ contains proofs, algorithmic details, and additional
results from simulations and data analysis.
(b) _Supplementary Code:_ an R package to implement the methodology can be
downloaded from https://github.com/chenyujie1104/ExponentialGM along with a
README file.

### References

  * Ackley et al. (1985) D. H. Ackley, G. E. Hinton, and T. J. Sejnowski (1985).  A learning algorithm for Boltzmann machines.  Cognitive Science 9, 147–169.
  * Aggarwal (2018) C. C. Aggarwal (2018).  Restricted Boltzmann Machines, pages 235–270.  Springer International Publishing, Cham.
  * Ashford and Sowden (1970) J. Ashford and R. Sowden (1970).  Multi-variate probit analysis.  Biometrics pages 535–546.
  * Barndorff-Nielsen (1978) O. Barndorff-Nielsen (1978).  Information and exponential families in statistical theory.  John Wiley & Sons.
  * Besag (1974) J. Besag (1974).  Spatial interaction and the statistical analysis of lattice systems.  Journal of the Royal Statistical Society: Series B (Methodological) 36, 192–225.
  * Betancourt (2011) M. Betancourt (2011).  Nested sampling with constrained Hamiltonian Monte Carlo.  In AIP Conference Proceedings, volume 1305, pages 165–172. American Institute of Physics.
  * Bhadra et al. (2017) A. Bhadra, J. Datta, N. G. Polson, and B. Willard (2017).  The Horseshoe+ Estimator of Ultra-Sparse Signals.  Bayesian Analysis 12, 1105 – 1131.
  * Bhattacharya et al. (2016) A. Bhattacharya, A. Chakraborty, and B. K. Mallick (2016).  Fast sampling with Gaussian scale mixture priors in high-dimensional regression.  Biometrika 103, 985–991.
  * Bhattacharya et al. (2016) A. Bhattacharya, D. B. Dunson, D. Pati, and N. S. Pillai (2016).  Sub-optimality of some continuous shrinkage priors.  Stochastic Processes and their Applications 126, 3828 – 3842.  In Memoriam: Evarist Giné.
  * Bhattacharya et al. (2015) A. Bhattacharya, D. Pati, N. S. Pillai, and D. B. Dunson (2015).  Dirichlet–Laplace priors for optimal shrinkage.  Journal of the American Statistical Association 110, 1479–1490.
  * Boucheron et al. (2003) S. Boucheron, G. Lugosi, and O. Bousquet (2003).  Concentration inequalities.  In Summer school on machine learning, pages 208–240. Springer.
  * Broyden (1970) C. G. Broyden (1970).  The convergence of a class of double-rank minimization algorithms 1. general considerations.  IMA Journal of Applied Mathematics 6, 76–90.
  * Bühlmann and Van De Geer (2011) P. Bühlmann and S. Van De Geer (2011).  Statistics for high-dimensional data: methods, theory and applications.  Springer Science & Business Media.
  * Chakraborty et al. (2020) A. Chakraborty, A. Bhattacharya, and B. K. Mallick (2020).  Bayesian sparse multiple regression for simultaneous rank reduction and variable selection.  Biometrika 107, 205–221.
  * Chakraborty et al. (2023) A. Chakraborty, R. Ou, and D. B. Dunson (2023).  Bayesian inference on high-dimensional multivariate binary responses.  Journal of the American Statistical Association 0, 1–12.
  * Chicco et al. (2021) D. Chicco, N. Tötsch, and G. Jurman (2021).  The matthews correlation coefficient (mcc) is more reliable than balanced accuracy, bookmaker informedness, and markedness in two-class confusion matrix evaluation.  BioData mining 14, 1–22.
  * Deng (2012) L. Deng (2012).  The MNIST database of handwritten digit images for machine learning research.  IEEE Signal Processing Magazine 29, 141–142.
  * Fischer and Igel (2012) A. Fischer and C. Igel (2012).  An introduction to restricted Boltzmann machines.  In Progress in Pattern Recognition, Image Analysis, Computer Vision, and Applications: 17th Iberoamerican Congress, CIARP 2012, Buenos Aires, Argentina, September 3-6, 2012. Proceedings 17, pages 14–36. Springer.
  * Fisher (1922) R. A. Fisher (1922).  On the mathematical foundations of theoretical statistics.  Philosophical transactions of the Royal Society of London, Series A 222, 309–368.
  * Friedman et al. (2010) J. H. Friedman, T. Hastie, and R. Tibshirani (2010).  Regularization paths for generalized linear models via coordinate descent.  Journal of Statistical Software 33, 1–22.
  * George and McCulloch (1993) E. I. George and R. E. McCulloch (1993).  Variable selection via gibbs sampling.  Journal of the American Statistical Association 88, 881–889.
  * Geyer (1991) C. J. Geyer (1991).  Markov chain Monte Carlo maximum likelihood.  In Computing Science and Statistics: Proc. 23rd Symposium of the Interface Foundation, pages 156–163.
  * Geyer and Thompson (1992) C. J. Geyer and E. A. Thompson (1992).  Constrained Monte Carlo maximum likelihood for dependent data.  Journal of the Royal Statistical Society: Series B 54, 657–683.
  * Ghosal and Van der Vaart (2017) S. Ghosal and A. Van der Vaart (2017).  Fundamentals of nonparametric Bayesian inference, volume 44.  Cambridge University Press.
  * Hebb (1949) D. O. Hebb (1949).  The organization of behavior: A neuropsychological theory.  Wiley.
  * Hinton (2002) G. E. Hinton (2002).  Training products of experts by minimizing contrastive divergence.  Neural Computation 14, 1771–1800.
  * Hinton (2007) G. E. Hinton (2007).  Boltzmann machine.  Scholarpedia 2, 1668\.
  * Hinton et al. (1998) G. E. Hinton, B. Sallans, and Z. Ghahramani (1998).  A hierarchical community of experts.  In Learning in graphical models, pages 479–494. Springer.
  * Hopfield (1982) J. J. Hopfield (1982).  Neural networks and physical systems with emergent collective computational abilities.  Proceedings of the National Academy of Sciences 79, 2554–2558.
  * Ising (1924) E. Ising (1924).  Beitrag zur theorie des ferro-und paramagnetismus.  PhD thesis, Grefe & Tiedemann Hamburg.
  * Kim et al. (2021) H.-S. Kim et al. (2021).  High whsc1l1 expression reduces survival rates in operated breast cancer patients with decreased cd8+ t cells: machine learning approach.  Journal of Personalized Medicine 11, 636\.
  * Lauritzen (1996) S. L. Lauritzen (1996).  Graphical models, volume 17.  Clarendon Press.
  * Le Roux and Bengio (2008) N. Le Roux and Y. Bengio (2008).  Representational power of restricted Boltzmann machines and deep belief networks.  Neural computation 20, 1631–1649.
  * Liang (2010) F. Liang (2010).  A double Metropolis–Hastings sampler for spatial models with intractable normalizing constants.  Journal of Statistical Computation and Simulation 80, 1007–1022.
  * Meinshausen and Bühlmann (2006) N. Meinshausen and P. Bühlmann (2006).  High-dimensional graphs and variable selection with the Lasso.  The Annals of Statistics 34, 1436 – 1462.
  * Meinshausen and Bühlmann (2010) N. Meinshausen and P. Bühlmann (2010).  Stability selection.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72, 417–473.
  * Mitchell and Beauchamp (1988) T. J. Mitchell and J. J. Beauchamp (1988).  Bayesian variable selection in linear regression.  Journal of the American Statistical Association 83, 1023–1032.
  * Møller et al. (2006) J. Møller, A. N. Pettitt, R. Reeves, and K. K. Berthelsen (2006).  An efficient Markov chain Monte Carlo method for distributions with intractable normalising constants.  Biometrika 93, 451–458.
  * Montufar and Ay (2011) G. Montufar and N. Ay (2011).  Refinements of universal approximation results for deep belief networks and restricted Boltzmann machines.  Neural computation 23, 1306–1319.
  * Murray et al. (2012) I. Murray, Z. Ghahramani, and D. MacKay (2012).  MCMC for doubly-intractable distributions.  arXiv preprint arXiv:1206.6848 .
  * Neal (2001) R. M. Neal (2001).  Annealed importance sampling.  Statistics and Computing 11, 125–139.
  * Neal et al. (2011) R. M. Neal et al. (2011).  MCMC using Hamiltonian dynamics.  Handbook of Markov chain Monte Carlo 2, 2\.
  * Newey (1991) W. K. Newey (1991).  Uniform convergence in probability and stochastic equicontinuity.  Econometrica: Journal of the Econometric Society pages 1161–1167.
  * Parikh and Boyd (2014) N. Parikh and S. Boyd (2014).  Proximal algorithms.  Foundations and trends® in Optimization 1, 127–239.
  * Park and Casella (2008) T. Park and G. Casella (2008).  The Bayesian lasso.  Journal of the American Statistical Association 103, 681–686.
  * Pinelis (2020) I. Pinelis (2020).  Exact lower and upper bounds on the incomplete gamma function.  arXiv preprint arXiv:2005.06384 .
  * Potts (1952) R. B. Potts (1952).  Some generalized order-disorder transformations.  In Mathematical proceedings of the cambridge philosophical society, volume 48, pages 106–109. Cambridge University Press.
  * Privitera et al. (2021) A. P. Privitera, V. Barresi, and D. F. Condorelli (2021).  Aberrations of chromosomes 1 and 16 in breast cancer: a framework for cooperation of transcriptionally dysregulated genes.  Cancers 13, 1585\.
  * Rao (1945) C. R. Rao (1945).  Information and accuracy attainable in the estimation of statistical parameters. kotz s & johnson nl (eds.), breakthroughs in statistics volume i: Foundations and basic theory, 235–248.
  * Ravikumar et al. (2010) P. Ravikumar, M. J. Wainwright, and J. D. Lafferty (2010).  High-dimensional Ising model selection using $\ell_{1}$-regularized logistic regression.  The Annals of Statistics 38, 1287–1319.
  * Roach and Gao (2023) L. Roach and X. Gao (2023).  Graphical local genetic algorithm for high-dimensional log-linear models.  Mathematics 11, 2514\.
  * Rothman et al. (2010) A. J. Rothman, E. Levina, and J. Zhu (2010).  Sparse multivariate regression with covariance estimation.  Journal of Computational and Graphical Statistics 19, 947–962.
  * Salakhutdinov and Hinton (2009) R. Salakhutdinov and G. Hinton (2009).  Deep Boltzmann machines.  In Artificial intelligence and statistics, pages 448–455. PMLR.
  * Salakhutdinov et al. (2007) R. Salakhutdinov, A. Mnih, and G. Hinton (2007).  Restricted Boltzmann machines for collaborative filtering.  In Proceedings of the 24th international conference on Machine learning, pages 791–798.
  * Salakhutdinov and Murray (2008) R. Salakhutdinov and I. Murray (2008).  On the quantitative analysis of deep belief networks.  In Proceedings of the 25th International Conference on Machine Learning, pages 872–879.
  * Smolensky (1986) P. Smolensky (1986).  Information processing in dynamical systems: Foundations of harmony theory.  In Parallel distributed processing: Explorations in the microstructure of cognition, pages 194–281–. MIT Press, Cambridge, MA.
  * Stoehr et al. (2019) J. Stoehr, A. Benson, and N. Friel (2019).  Noisy Hamiltonian Monte Carlo for doubly intractable distributions.  Journal of Computational and Graphical Statistics 28, 220–232.
  * Sutskever and Tieleman (2010) I. Sutskever and T. Tieleman (2010).  On the convergence properties of contrastive divergence.  In Proc. AISTATS, pages 789–795.
  * TCGA (2012) TCGA (2012).  Comprehensive molecular portraits of human breast tumours.  Nature 490, 61–70.
  * Tibshirani (1996) R. Tibshirani (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society. Series B (Methodological) 58, 267–288.
  * Varga (2010) R. S. Varga (2010).  Gervsgorin and his circles, volume 36.  Springer Science & Business Media.
  * Wainwright and Jordan (2008) M. J. Wainwright and M. I. Jordan (2008).  Graphical models, exponential families, and variational inference.  Foundations and Trends® in Machine Learning 1, 1–305.
  * Wu (1983) C. F. J. Wu (1983).  On the Convergence Properties of the EM Algorithm.  The Annals of Statistics 11, 95 – 103.
  * Xiao et al. (2022) X. Xiao et al. (2022).  Tmprss2 serves as a prognostic biomarker and correlated with immune infiltrates in breast invasive cancer and lung adenocarcinoma.  Frontiers in Molecular Biosciences 9, 647826\.
  * Yang et al. (2012) E. Yang, G. Allen, Z. Liu, and P. Ravikumar (2012).  Graphical models via generalized linear models.  In F. Pereira, C. Burges, L. Bottou, and K. Weinberger, editors, Advances in Neural Information Processing Systems, volume 25. Curran Associates, Inc.
  * Yang et al. (2015) E. Yang, P. Ravikumar, G. I. Allen, and Z. Liu (2015).  Graphical models via univariate exponential family distributions.  Journal of Machine Learning Research 16, 3813–3847.
  * Yang et al. (2013) E. Yang, P. K. Ravikumar, G. I. Allen, and Z. Liu (2013).  On Poisson graphical models.  Advances in neural information processing systems 26,.
  * Yu (2013) Y.-L. Yu (2013).  On decomposing the proximal map.  Advances in neural information processing systems 26,.
  * Yuan et al. (2007) M. Yuan, A. Ekici, Z. Lu, and R. Monteiro (2007).  Dimension reduction and coefficient estimation in multivariate linear regression.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 69, 329–346.
  * Zhang et al. (2019) S. Zhang, L. Yao, A. Sun, and Y. Tay (2019).  Deep learning based recommender system: A survey and new perspectives.  ACM computing surveys (CSUR) 52, 1–38.

Supplementary Material to
Likelihood Based Inference in Fully and Partially Observed Exponential Family
Graphical Models with Intractable Normalizing Constants

##

###  S.1 Proofs

####  S.1.1 Proof of Proposition 3.2

(a). (Proof for $z_{\theta}$). For brevity, define $T_{j}=T_{j}(X_{j})$ and
$T_{jk}=T_{jk}(X_{j},X_{k})$. Recall $\phi=\mathrm{diag}(\theta)$. Thus, from
Equation (2), $q_{\theta}(X)/q_{\phi}(X)=\exp\left\\{\underset{(j,k)\in
E}{\sum}\theta_{jk}T_{jk}\right\\}$. However, the model of Equation (2) is of
the form $p_{\eta}(X)=\exp\\{\eta T(X)+H(X)-A(\eta)\\}$ with
$T(X)=(T_{j},T_{jk})$ as the sufficient statistic. That is, $T(X)$ is a
$p\times p$ symmetric matrix with $T_{j}$ as the $j$th diagonal entry and
$T_{jk}$ as the $(j,k)$th off-diagonal entry. The moment generating function
of exponential family distribution satisfies
$M_{T,\eta}(u)=\mathbb{E}_{\eta}(\exp(u^{T}T(X))=\exp(A(u+\eta)-A(\eta)).$
Setting $u=\theta-\mathrm{diag}(\theta)$, and $\phi=\mathrm{diag}(\theta)$,
one has, $u+\phi=\theta\in\Theta$. Similarly, $\phi\in\Theta$, because this
corresponds to the special case where the $p$ variables are independent. Using
the fact that for any valid parameter $\gamma\in\Theta$ we have
$A(\gamma)<\infty$ in an exponential family model, shows that:

$\displaystyle\mathbb{E}_{\phi}(q_{\theta}(X)/q_{\phi}(X))$
$\displaystyle=M_{T,\phi}(u)=\exp(A(u+\phi)-A(\phi))=\exp(A(\theta)-A(\phi))<\infty.$

Similarly, if, $2u+\phi\in\Theta$, then,
$\mathbb{E}_{\phi}\\{(q_{\theta}(X)/q_{\phi}(X))^{2}\\}=M_{T,\phi}(2u)=\exp(A(2u+\phi)-A(\phi))<\infty.$
Thus, $\mathrm{Var}_{\phi}(q_{\theta}(X)/q_{\phi}(X))<\infty.$

(b). (Proof for $\nabla_{\theta}z(\theta)$). Next consider
$\nabla_{\theta}z(\theta)$. For $(j,k)\in E$,
$\frac{\nabla_{\theta_{jk}}q_{\theta}(X)}{q_{\phi}(X)}=\frac{q_{\theta}(X)}{q_{\phi}(X)}T_{jk}$.
Thus,

$\mathbb{E}_{\phi}\left\\{\left(\frac{\nabla_{\theta_{jk}}q_{\theta}(X)}{q_{\phi}(X)}\right)^{2}\right\\}=\mathbb{E}_{\phi}\left\\{\left(\frac{q_{\theta}(X)}{q_{\phi}(X)}\right)^{2}T_{jk}^{2}\right\\}\leq\left\\{\mathbb{E}_{\phi}\left(\frac{q^{2a}_{\theta}(X)}{q^{2a}_{\phi}(X)}\right)\right\\}^{1/a}\left[\mathbb{E}_{\phi}\left(T_{jk}^{2b}\right)\right]^{1/b}<\infty,$

where $a>1,\;b>1,\;(1/a+1/b)=1$ and the claim in the previous display is a
consequence of the Hölder inequality. Similar to Part (a), and taking
$a=1+\delta$, we have:
$\mathbb{E}_{\phi}\left(\frac{q^{2a}_{\theta}(X)}{q^{2a}_{\phi}(X)}\right)=M_{T,\phi}(2au)=M_{T,\phi}(2(1+\delta)u),$
which is bounded if $(2(1+\delta)u+\phi)\in\Theta$, for some $\delta>0$. The
boundedness of ${E}_{\phi}\left(T_{jk}^{2b}\right)$ follows from the fact that
for an exponential family distribution, the sufficient statistic also belongs
to an exponential family, for which moments of all orders exist (Theorem 8.1,
Barndorff-Nielsen, 1978). Repeating the calculation for all $(j,k)\in E$, the
result follows. If $(j,k)\notin E$, the variance is zero and the result is
trivial.

####  S.1.2 Proof of Proposition 3.3

Since $\phi=\text{diag}(\theta)$ and $X$ is distributed according to a PEGM,
we have $U=q_{\theta}(X)/q_{\phi}(X)=\exp\\{\sum_{(j,k)\in
E}\theta_{jk}X_{j}X_{k}\\}$. Next, recall that a function
$g:\mathbb{R}^{p}\to\mathbb{R}$ is said to have the bounded difference
property if for some non-negative constants $c_{1},\ldots,c_{p}$,

$\displaystyle\sup_{x_{1},\ldots,x_{p},x_{j^{\prime}}\in\mathbb{R}^{p}}|g(x_{1},\ldots,x_{j},\ldots,x_{p})-g(x_{1},\ldots,x_{j^{\prime}},\ldots,x_{p})|\leq
c_{j}.$ (S.1)

In the current context, let $U=U(X_{1},\ldots,X_{j},\ldots,X_{p})$ and
$U^{\prime}=U(X_{1},\ldots,X_{j^{\prime}},\ldots,X_{p})$. By assumption, each
$X_{j}$ has support in $[a_{j},b_{j}]$. Without loss of generality, let
$|b_{j}|\geq|a_{j}|$ for all $j=1,\ldots,p$. Clearly,
$U\leq\exp\\{\sum_{(j,k)\in
E}\theta_{jk}|X_{j}X_{k}|\\}\leq\exp\\{\sum_{(j,k)\in
E}\theta_{jk}|b_{j}b_{k}|\\}$. Similarly, $U^{\prime}$ is also bounded. Then
by the triangle inequality, $|U-U^{\prime}|=|U|+|U^{\prime}|\leq
2\exp\\{\sum_{(j,k)\in E}\theta_{jk}|b_{j}b_{k}|\\}$. Thus, $U$ is of bounded
difference according to the definition (S.1). The second assertion then
follows from Boucheron et al. (2003, Theorem 12). Noting that $V=X_{j}X_{k}U$,
the same argument for $V$ then applies.

####  S.1.3 Auxiliary results for Theorems 6.4 and 6.5

In this subsection, we record a few auxiliary results which will be used to
prove Theorems in Section 6 of the main draft. The first result is an
exponential concentration result, which is used to prove Theorem 6.4. The
second result is needed for the proof of Theorem 6.5, and characterizes the
prior Kullback-Leibler support set in terms of the Euclidean distance between
the true parameter and prior realizations.

######  Proposition S.1.1.

Suppose we observe
$X_{1},\ldots,X_{n}\overset{iid}{\sim}\text{PEGM}(\theta_{0})$ where the model
satisfies Assumption 2.2 and Assumption 2.3. Let $\mu_{0}=(\mu_{0,jk})$ be the
matrix of moments, i.e. $\mu_{0,jk}=\mathbb{E}_{\theta_{0}}[T(X_{j},X_{k})]$
when $j\neq k$, and $\mu_{0,jj}=\mathbb{E}_{\theta_{0}}[T(X_{j})]$. Define
$S_{n}=\sum_{i=1}^{n}T(X_{ij},X_{ik})$. Then for some $c\geq 0$

$\mathbb{P}[|S_{n}-n\mu_{0,jk}|>n\nu]\leq 2\exp(-cn\nu^{2}),$

for $0\leq\nu\leq\kappa_{5}$.

###### Proof.

Define $S_{n}=\sum_{i=1}^{n}T(X_{ij},X_{ik})$. We have for any $s$,

$\displaystyle\mathbb{P}(S_{n}-n\mu_{0,jk}>n\nu)$
$\displaystyle\leq\mathbb{E}[e^{s(S_{n}-n\mu_{0,jk})}]e^{-sn\nu}$
$\displaystyle=\prod_{i=1}^{n}\mathbb{E}[e^{sT(X_{ij},X_{ik})}]e^{-ns(\mu_{0,jk}+\nu)}.$

Write $T_{j}=T(X_{j})$ and $T_{jk}=T(X_{j},X_{k})$. Then,

$\displaystyle\log\mathbb{E}[\exp(\alpha T_{jk})]$
$\displaystyle=\log\int_{\mathcal{X}}\exp\left\\{\alpha
T_{jk}+\sum_{j=1}^{p}\theta_{0j}T_{j}+\sum_{\begin{subarray}{c}j,k=1\\\ j\neq
k\end{subarray}}^{p}\theta_{0,jk}T_{jk}+\sum_{j=1}^{p}C(X_{j})-A(\theta_{0})\right\\}\mathrm{d}x$
$\displaystyle=\bar{A}_{jk}(\theta_{0};\alpha)-\bar{A}_{jk}(\theta_{0};0).$

Letting $\alpha\leq 1$ and using a Taylor expansion we get for some
$\gamma\in(0,1)$,

$\displaystyle\bar{A}_{jk}(\theta_{0};\alpha)-\bar{A}_{jk}(\theta_{0};0)$
$\displaystyle=\alpha\dfrac{\partial}{\partial\eta}\bar{A}_{jk}(\theta_{0};\eta)\rvert_{\eta=0}+\frac{\alpha^{2}}{2}\dfrac{\partial^{2}}{\partial\eta^{2}}\bar{A}_{jk}(\theta_{0};\eta)\rvert_{\eta=\gamma\alpha}$
$\displaystyle\leq\alpha\kappa_{3}+\frac{\alpha^{2}}{2}\kappa_{5}.$

Thus,
$\mathbb{E}[\exp(sT_{jk})]\leq\exp(s\kappa_{3}+\frac{s^{2}}{2}\kappa_{5})$ for
$s\leq 1$. Hence, for $0\leq s\leq 1$, and then minimizing over $s$ and since
$\mu_{0,jk}\geq-\kappa_{3}$ we get,

$\displaystyle\mathbb{P}(S_{n}-n\mu_{0,jk}>n\nu)$ $\displaystyle\leq
e^{-ns\nu+ns^{2}\frac{\kappa_{5}}{2}}\leq e^{-n\frac{\nu^{2}}{2\kappa_{5}}},$
$\displaystyle\leq
e^{-n\frac{(\nu+\mu_{0,jk}-\kappa_{1})^{2}}{2\kappa_{5}}},\quad
s=(\nu+\mu_{0,jk}-\kappa_{1})/\kappa_{5}$ $\displaystyle\leq
e^{-cn\nu^{2}},\quad\nu\leq(\kappa_{1}+\kappa_{3})/2.$

when $\nu\leq\kappa_{5}$. Since $\mathbb{P}(|X|>t)\leq 2\mathbb{P}(X>t)$, the
result follows. ∎

A similar analysis leads to $\mathbb{P}(|S_{n}-\theta_{0,jj}|>n\nu)\leq
2\exp(-cn\nu^{2})$ for $0\leq\nu\leq\kappa_{4}$.

For notational convenience, we write
$\text{vech}(\theta)=\theta\in\mathbb{R}^{p(p+1)/2}$. Let
$\mu\in\mathbb{R}^{p\times p}$ denote the corresponding symmetric matrix of
mean parameters, i.e. $\mu_{jk}=\mathbb{E}[T(X_{j},X_{k})]$ for $j\neq k$ and
$\mu_{jj}=\mathbb{E}[T(X_{j})]$ for $j=1,\ldots,p$. Similarly, $\mu_{0}$ is
the mean matrix when the parameter is $\theta_{0}$. Again, we write
$\text{vech}(\mu)=\mu$.

######  Proposition S.1.2.

Suppose $\pi$ is some distribution on $\Theta$. Then for every $\epsilon>0$
the exists $\delta=\delta(\epsilon)>0$ such that,

$\pi\left\\{\theta\in\Theta:\text{KL}(p_{\theta_{0}},p_{\theta})<\epsilon\right\\}\geq\pi\left\\{\theta\in\Theta:\left\lVert\theta-\theta_{0}\right\rVert_{2}<\delta\right\\}.$

###### Proof.

It is straightforward to establish that
$KL(p_{\theta_{0}},p_{\theta})=\langle\theta_{0}-\theta,\mu_{0}\rangle+\log
z(\theta)-\log z(\theta_{0})$. Indeed,

$\displaystyle KL(p_{\theta_{0}},p_{\theta})=$ $\displaystyle\log
z(\theta)-\log z(\theta_{0})+\sum_{j}(\theta_{0,jj}-\theta_{jj})\int
T(x_{j})p_{\theta_{0}}(x)dx$ $\displaystyle+\sum_{j\neq
k}(\theta_{0,jk}-\theta_{jk})\int T(x_{j},x_{k})p_{\theta_{0}}(x)$
$\displaystyle=$ $\displaystyle\langle\theta-\theta_{0},\mu_{0}\rangle.$

Write $A(\theta)=\log z(\theta)$ and $A(\theta_{0})=\log z(\theta_{0})$. It is
well known that $A(\theta)$ is convex and has derivatives of all orders
(Wainwright and Jordan, 2008, Proposition 3.1). Furthermore, $\nabla
A(\theta)=\mu$ for all $\theta$. Also, Taylor’s theorem implies that for some
$\gamma\in(0,1)$,

$\displaystyle A(\theta)$ $\displaystyle=A(\theta_{0})+\langle\nabla
A(\theta_{0}),\theta-\theta_{0}\rangle+(\theta-\theta_{0})^{\mathrm{\scriptscriptstyle{T}}}\left[\int_{0}^{1}(1-\gamma)\\{\nabla^{2}A(\theta_{0}+\gamma(\theta-\theta_{0}))\\}d\gamma\right](\theta-\theta_{0})$
$\displaystyle\leq A(\theta_{0})+\langle\nabla
A(\theta_{0}),\theta-\theta_{0}\rangle+(1-\gamma)\left\lVert\theta-\theta_{0}\right\rVert_{2}^{2}\lambda_{\max}[\nabla^{2}A(\theta_{0}+\gamma(\theta-\theta_{0}))],$

where $\lambda_{\max}(P)$ is the largest eigenvalue of a matrix $P$. Taking
the limit $\gamma\to 0$ we obtain

$\displaystyle KL(p_{\theta_{0}},p_{\theta})$
$\displaystyle=\langle\theta_{0}-\theta,\mu_{0}\rangle+A(\theta)-A(\theta_{0})$
$\displaystyle\leq\langle\theta_{0}-\theta,\mu_{0}\rangle+\langle-\nabla
A(\theta_{0}),\theta_{0}-\theta\rangle+\lambda_{\max}[\nabla^{2}A(\theta_{0})]\left\lVert\theta-\theta_{0}\right\rVert_{2}^{2}$
$\displaystyle=\lambda_{\max}[\nabla^{2}A(\theta_{0})]\left\lVert\theta-\theta_{0}\right\rVert_{2}^{2}.$

Hence
$\left\\{\theta:\left\lVert\theta-\theta_{0}\right\rVert_{2}<\epsilon\right\\}\subset\left\\{\theta:\text{KL}(p_{\theta_{0}},p_{\theta})<\lambda_{\max}[\nabla^{2}A(\theta_{0})]\epsilon\right\\}$
which implies that

$\pi\left\\{\theta\in\Theta:\text{KL}(p_{\theta_{0}},p_{\theta})<\lambda_{\max}[\nabla^{2}A(\theta_{0})]\delta\right\\}\geq\pi\left\\{\theta\in\Theta:\left\lVert\theta-\theta_{0}\right\rVert_{2}<\delta\right\\}.$

∎

####  S.1.4 Proof of Theorem 6.4

Define the matrix $S=(S_{jk})$ where $S_{jk}=\sum_{i=1}^{n}T(X_{ij},X_{ik})$
if $j\neq k$ and $S_{jj}=\sum_{i=1}^{n}T(X_{ij})$. Then
$\ell(\theta)=\sum_{i=1}^{n}q_{\theta}(X_{i})-nA(\theta)$ $=\mathrm{tr}(\theta
S)-nA(\theta)$. The $\ell_{1}$-penalized estimator is defined as the minimizer
of
$-\ell(\theta)+\lambda\left\lVert\theta\right\rVert_{1}=nA(\theta)-\mathrm{tr}(\theta
S)+\lambda\left\lVert\theta\right\rVert_{1}$. Let

$\displaystyle Q(\theta)$ $\displaystyle=nA(\theta)-\mathrm{tr}(\theta
S)+\lambda\left\lVert\theta\right\rVert_{1}-nA(\theta_{0})+\mathrm{tr}(\theta_{0}S)-\lambda\left\lVert\theta_{0}\right\rVert_{1}$
$\displaystyle=n[A(\theta)-A(\theta_{0})]-\mathrm{tr}[(\theta-\theta_{0})(S-\mu_{0})]-\mathrm{tr}[(\theta-\theta_{0})\mu_{0}]+\lambda[\left\lVert\theta\right\rVert_{1}-\left\lVert\theta_{0}\right\rVert_{1}].$

Then, the $\ell_{1}$ penalized estimate minimizes $Q(\theta)$. Define
$\Delta=\theta-\theta_{0}$ and $G(\Delta)=Q(\theta_{0}+\Delta)$. Equivalently,
we want to minimize $G(\Delta)$ with respect to $\Delta$. First, note that for
any $\alpha\in(0,1)$ and $\theta_{1},\theta_{2}\in\Theta$,

$\displaystyle\small\alpha Q(\theta_{1})=$ $\displaystyle
n\alpha[A(\theta_{1})-A(\theta_{0})]-\mathrm{tr}[\alpha(\theta_{1}-\theta_{0})(S-\mu_{0})]-\mathrm{tr}[\alpha(\theta_{1}-\theta_{0})\mu_{0}]$
$\displaystyle+\lambda\alpha[\left\lVert\theta_{1}\right\rVert_{1}-\left\lVert\theta_{0}\right\rVert_{1}],$
$\displaystyle(1-\alpha)Q(\theta_{2})=$ $\displaystyle
n(1-\alpha)[A(\theta_{2})-A(\theta_{0})]-\mathrm{tr}[(1-\alpha)(\theta_{2}-\theta_{0})S]-\mathrm{tr}[(1-\alpha)(\theta_{2}-\theta_{0})\mu_{0}]$
$\displaystyle+\lambda(1-\alpha)[\left\lVert\theta_{1}\right\rVert_{1}-\left\lVert\theta_{0}\right\rVert_{1}].$

This implies that if $\theta^{\star}=\alpha\theta_{1}+(1-\alpha)\theta_{2}$
then,

$\displaystyle\alpha Q(\theta_{1})+(1-\alpha)Q(\theta_{2})$ $\displaystyle\geq
n[A(\theta^{\star})-A(\theta)]-\mathrm{tr}(\theta^{\star}S)-\mathrm{tr}(\theta^{\star}\mu_{0})+\lambda[\left\lVert\theta^{\star}\right\rVert_{1}-\left\lVert\theta_{0}\right\rVert_{1}]$
$\displaystyle=Q(\theta^{\star}),$

by triangle inequality and since $A(\theta)$ is convex. Similarly, if
$\Delta_{1}=\theta_{1}-\theta_{0}$ and $\Delta_{2}=\theta_{2}-\theta_{0}$,
then $\alpha G(\Delta_{1})+(1-\alpha)G(\Delta_{2})\geq G(\Delta^{\star})$
where $\Delta^{\star}=\theta^{\star}-\theta_{0}$. Thus, $G(\Delta)$ is convex
and $G(0)=0$. Hence, if $\hat{\Delta}=\inf_{\Delta}G(\Delta)$, then
$G(\hat{\Delta})\leq 0$.

Now,
$G(\Delta)=n[A(\theta_{0}+\Delta)-A(\theta_{0})]-\mathrm{tr}[\Delta(S-\mu_{0})]-\mathrm{tr}(\Delta\mu_{0})+\lambda[\left\lVert\theta_{0}+\Delta\right\rVert_{1}-\left\lVert\theta_{0}\right\rVert_{1}]$.
A Taylor expansion of $A(\theta_{0}+\Delta)-A(\theta_{0})$ yields

$\displaystyle A(\theta_{0}+\Delta)-A(\theta_{0})=\mathrm{tr}(\nabla
A(\theta_{0})\Delta)+\tilde{\Delta}^{\mathrm{\scriptscriptstyle{T}}}\left[\int_{0}^{1}(1-v)\frac{\partial^{2}}{\partial
v^{2}}A(\theta_{0}+v\Delta)dv\right]\tilde{\Delta},$

where $\tilde{\Delta}=\text{vec}(\Delta)$. Thus
$G(\Delta)=\tilde{\Delta}^{\mathrm{\scriptscriptstyle{T}}}\left[\int_{0}^{1}(1-v)\partial^{2}A(\theta_{0}+v\Delta)dv\right]\tilde{\Delta}-\mathrm{tr}[\Delta(S-\mu_{0})]+\lambda[\left\lVert\theta_{0}+\Delta\right\rVert_{1}-\left\lVert\theta_{0}\right\rVert_{1}]$.
That is, $G(\Delta)=T_{1}-T_{2}+T_{3}$. Let $\mathbf{S}=\\{j\leq
k:\theta_{0,jk}\neq 0,j,k=1,\ldots,p\\}$. We next work with each term
separately. We have

$\displaystyle|T_{2}|=|\mathrm{tr}[\Delta(S-\mu_{0})]|$
$\displaystyle\leq\left|\sum\sum_{j\leq
k}(s_{jk}-\mu_{0,jk})\Delta_{jk}\right|$ $\displaystyle\leq\left|\sum_{j\leq
k:}\sum_{\theta_{0,jk}=0}(s_{jk}-\mu_{0,jk})\Delta_{jk}\right|+\left|\sum_{j\leq
k:}\sum_{\theta_{0,jk}\neq 0}(s_{jk}-\mu_{0,jk})\Delta_{jk}\right|$
$\displaystyle\leq\max_{\mathbf{S}^{c}}|s_{jk}-\mu_{0,jk}|\left\lVert\Delta_{\mathbf{S}^{c}}\right\rVert_{1}+\sqrt{s}\max_{\mathbf{S}}|s_{jk}-\mu_{0,jk}|\left\lVert\Delta_{\mathbf{S}}\right\rVert_{F},$

where the second bound comes from Cauchy-Schwarz inequality. From Proposition
S.1.1, and the union bound we have that with probability at least
$1-2e^{-C\log p/n^{2}}$, $\max_{j\leq k}|s_{jk}-\mu_{0,jk}|\leq
C_{1}\sqrt{\log p/n}$ for some $C,C_{1}>0$. Thus we have with probability
tending to 1,

$\displaystyle|T_{2}|\leq C_{1}\sqrt{\frac{\log
p}{n}}\left\lVert\Delta_{\mathbf{S}^{c}}\right\rVert_{1}+C_{1}\sqrt{\frac{s\log
p}{n}}\left\lVert\Delta_{\mathbf{S}}\right\rVert_{F}.$

Next, we work with
$T_{3}=\lambda[\left\lVert\theta_{0}+\Delta\right\rVert_{1}-\left\lVert\theta_{0}\right\rVert_{1}]$.
Clearly,
$\left\lVert\theta_{0}\right\rVert_{1}=\left\lVert\theta_{0\mathbf{S}}\right\rVert_{1}$
and
$\left\lVert\theta_{0}+\Delta\right\rVert_{1}=\left\lVert\theta_{0\mathbf{S}}+\Delta_{S}\right\rVert_{1}+\left\lVert\Delta_{\mathbf{S}^{c}}\right\rVert_{1}$.
Thus, by triangle inequality

$T_{3}\geq\lambda(\left\lVert\Delta_{\mathbf{S}^{c}}\right\rVert_{1}-\left\lVert\Delta_{\mathbf{S}}\right\rVert_{1}).$

Set $\lambda=\frac{C_{1}}{\epsilon}\sqrt{\frac{\log p}{n}}$ for some
$\epsilon>0$. Then since
$\left\lVert\Delta_{\mathbf{S}}\right\rVert_{1}\leq\sqrt{s}\left\lVert\Delta_{\mathbf{S}}\right\rVert_{F}$,
we have,

$\displaystyle T_{3}$
$\displaystyle\geq\dfrac{C_{1}}{\epsilon}\sqrt{\dfrac{\log
p}{n}}\left\lVert\Delta_{\mathbf{S}^{c}}\right\rVert_{1}-\dfrac{C_{1}}{\epsilon}\sqrt{\dfrac{s\log
p}{n}}\left\lVert\Delta_{\mathbf{S}}\right\rVert_{F}$
$\displaystyle\geq\dfrac{C_{1}}{\epsilon}\sqrt{\dfrac{\log
p}{n}}\left\lVert\Delta_{\mathbf{S}^{c}}\right\rVert_{1}-\dfrac{C_{1}}{\epsilon}\sqrt{\dfrac{s\log
p}{n}}\left\lVert\Delta\right\rVert_{F}.$

Finally, for $T_{1}$ we have for
$\Delta\in\Theta_{n}(M)=\\{\Delta:\Delta=\Delta^{\mathrm{\scriptscriptstyle{T}}},\left\lVert\Delta\right\rVert_{F}=Mr_{n}\\}$,

$\displaystyle T_{1}$
$\displaystyle=\tilde{\Delta}^{\mathrm{\scriptscriptstyle{T}}}\left[\int_{0}^{1}(1-v)\frac{\partial^{2}}{\partial
v^{2}}A(\theta_{0}+v\Delta)dv\right]\tilde{\Delta}$
$\displaystyle\geq\lambda_{\min}\left(\int_{0}^{1}(1-v)\frac{\partial^{2}}{\partial
v^{2}}A(\theta_{0}+v\Delta)dv\right)\left\lVert\Delta\right\rVert_{F}^{2}$
$\displaystyle\geq\left\lVert\Delta\right\rVert_{F}^{2}\int_{0}^{1}(1-v)\lambda_{\min}\left(\frac{\partial^{2}}{\partial
v^{2}}A(\theta_{0}+v\Delta)\right)dv$
$\displaystyle\geq\left\lVert\Delta\right\rVert_{F}^{2}\min_{0\leq v\leq
1}\lambda_{\min}\left(\frac{\partial^{2}}{\partial
v^{2}}A(\theta_{0}+v\Delta)\right)\int_{0}^{1}(1-v)dv$
$\displaystyle\geq\frac{1}{2}\left\lVert\Delta\right\rVert_{F}^{2}\min\left[\lambda_{\min}\left(\frac{\partial^{2}}{\partial\theta^{2}}A(\theta_{0}+\bar{\Delta})\right):\left\lVert\bar{\Delta}\right\rVert_{F}\leq
Mr_{n}\right]$
$\displaystyle\geq\frac{1}{2}\left\lVert\Delta\right\rVert_{F}^{2}\min\left[\lambda_{\min}\left(\frac{\partial^{2}}{\partial\theta^{2}}A(\theta_{0})+\rho\left\lVert\bar{\Delta}\right\rVert_{F}:\left\lVert\bar{\Delta}\right\rVert_{F}\leq
Mr_{n}\right)\right],$

where the last inequality follows from a Taylor expansion on
$\frac{\partial^{2}}{\partial\theta^{2}}A(\theta_{0})$ and Assumption 6.2.
Since $\left\lVert\bar{\Delta}\right\rVert_{F}\leq Mr_{n}$ and $Mr_{n}\to 0$
as $n\to\infty$ we have that,

$\lambda_{\min}\left(\frac{\partial^{2}}{\partial\theta^{2}}A(\theta_{0})+\rho\left\lVert\bar{\Delta}\right\rVert_{F}:\left\lVert\Delta\right\rVert_{F}\leq
Mr_{n}\right)\geq\frac{\beta}{2}.$

Putting all the pieces together we obtain,

$\displaystyle G(\Delta)$
$\displaystyle\geq\frac{n\beta}{4}\left\lVert\Delta\right\rVert_{F}^{2}-C_{1}\sqrt{\dfrac{\log
p}{n}}\left\lVert\Delta_{S^{c}}\right\rVert_{1}\left(1-\frac{1}{\epsilon}\right)-C_{1}\sqrt{\dfrac{s\log
p}{n}}\left\lVert\Delta_{S}\right\rVert_{F}\left(1+\frac{1}{\epsilon}\right)$
$\displaystyle>0,$

for sufficiently large $M$ and $n$. Hence, we have argued that if
$\hat{\Delta}$ is the minimizer of $G(\Delta)$ then $G(\hat{\Delta})\leq 0$
and $\inf\\{G(\Delta):\Delta\in\Theta_{n}(M)\\}>0$, then $\hat{\Delta}$ must
be inside the ball $\Theta_{n}(M)$, i.e.
$\left\lVert\hat{\Delta}\right\rVert_{F}\leq Mr_{n}$.

####  S.1.5 Proof of Theorem 6.5

The vectorized data generating parameter
$\text{vech}(\theta_{0})\in\Theta\subset\mathbb{R}^{p(p+1)/2}$. By Assumption
6.2, $\Theta$ is open, hence there exists $\epsilon>0$ such that
$\\{\theta:\left\lVert\theta-\theta_{0}\right\rVert_{2}<\epsilon\\}\subset\Theta$.
Moreover, from Proposition S.1.2, we have

$\pi\left\\{\theta\in\Theta:\left\lVert\theta-\theta_{0}\right\rVert_{2}<\epsilon\right\\}\geq\pi\left\\{\theta\in\Theta:\text{KL}(p_{\theta_{0}},p_{\theta})<\delta(\epsilon)\right\\}.$

Hence, if we can show that
$\pi\left\\{\theta:\left\lVert\theta-\theta_{0}\right\rVert_{2}<\epsilon\right\\}>0$
for every sufficiently small $\epsilon>0$, then that would imply posterior
consistency in weak topology (Ghosal and Van der Vaart, 2017, Example 6.20).
To that end, we have

$\displaystyle\pi\left\\{\theta:\left\lVert\theta-\theta_{0}\right\rVert_{2}<\epsilon\right\\}$
$\displaystyle\geq\dfrac{\pi_{L}\left\\{\theta:\left\lVert\theta-\theta_{0}\right\rVert_{2}<\epsilon\right\\}}{\pi_{L}(\Theta)}$
$\displaystyle\geq\pi_{L}\left\\{\theta:\left\lVert\theta-\theta_{0}\right\rVert_{2}<\epsilon\right\\}$
(S.2)

since $\pi_{L}(\Theta)\leq 1$, and where
$\pi(\theta)\propto\pi_{L}(\theta)\mathbb{I}_{\theta\in\Theta}$. Next, we
obtain a lower bound of
$\pi_{L}\left\\{\theta:\left\lVert\theta-\theta_{0}\right\rVert_{2}<\epsilon\right\\}$.
To simplify notation, we show for $\theta,\theta_{0}\in\mathbb{R}^{p}$, and
$\theta_{0}$ having $s$ non-zero entries

$\pi_{L}\left\\{\theta:\left\lVert\theta-\theta_{0}\right\rVert_{2}<\epsilon\right\\}>0.$

The result then applies to our context where
$\theta_{0},\theta\in\mathbb{R}^{p(p+1)/2}$.

Define $S=\\{j:|\theta_{0j}|\neq 0\\}$. We then have

$\displaystyle\mathbb{P}[\left\lVert\theta-\theta_{0}\right\rVert_{2}<\delta]$
$\displaystyle\geq\mathbb{P}\left[\left\lVert\theta_{S}-\theta_{0S}\right\rVert_{2}<\frac{\delta}{2}\right]\mathbb{P}\left[\left\lVert\theta_{S^{c}}\right\rVert_{2}<\frac{\delta}{2}\right]$
$\displaystyle\geq\prod_{j\in
S}\mathbb{P}\left[|\theta_{j}-\theta_{0j}|<\frac{\delta}{2\sqrt{s}}\right]\int_{I_{\lambda}}\mathbb{P}\left[\left\lVert\theta_{S^{c}}\right\rVert_{2}<\frac{\delta}{2}\mid\lambda\in
I_{\lambda}\right]\pi(\lambda)d\lambda,$

where $I_{\lambda}=[\lambda_{1},\lambda_{2}]$ with
$\lambda_{1}=\sqrt{\frac{np\log p}{s}}$ and $\lambda_{2}=2\lambda_{1}$. We
next focus on the conditional probability:

$\displaystyle\mathbb{P}\left[\left\lVert\theta_{0S^{c}}\right\rVert_{2}\leq\frac{\delta}{2}\mid\lambda\in
I_{\lambda}\right]\geq\prod_{j\notin
S}\mathbb{P}\left[|\theta_{j}|<\frac{\delta}{2\sqrt{p}}\mid\lambda\in
I_{\lambda}\right].$

First, we handle the term
$\mathbb{P}\left[|\theta_{j}|<\frac{\delta}{2\sqrt{p}}\mid\lambda\in
I_{\lambda}\right]$. Recall, if the random variable $X$ is sub-exponential
distribution with parameters $(\nu^{2},\rho),$ then,

$\mathbb{P}\left[|X-\mu|\geq
t\right]\leq\begin{cases}&2e^{-\frac{t^{2}}{2\nu^{2}}},\quad\text{if
}0<t\leq\frac{\nu^{2}}{\rho},\\\ &2e^{-\frac{t}{2\rho}},\quad\text{if
}t>\frac{\nu^{2}}{\rho}.\end{cases}$ (S.3)

Here $\theta_{j}\mid\lambda\sim\text{Laplace}(1/\lambda)$, and the Laplace
distribution is sub-exponential with $\nu=2/\lambda$ and $\rho=2/\lambda$. Set
$t=\frac{\delta}{2\sqrt{p}}=\frac{1}{2}\sqrt{\frac{s\log p}{np}}$ in (S.3) and
note that $\nu^{2}/\rho=2/\lambda$. Thus, for any
$\lambda>4\sqrt{\frac{np}{s\log p}}$ we have:

$\displaystyle\mathbb{P}\left[|\theta_{j}|<\frac{\delta}{2\sqrt{p}}\mid\lambda\right]$
$\displaystyle\geq 1-2e^{-\frac{\delta\lambda}{4\sqrt{p}}}.$

Specifically, for any $\lambda\in I_{\lambda}$ we get,

$\displaystyle\mathbb{P}\left[|\theta_{j}|<\frac{\delta}{2\sqrt{p}}\mid\lambda\in
I_{\lambda}\right]$ $\displaystyle\geq
1-2e^{-\frac{\delta\lambda}{4\sqrt{p}}}\geq
1-2e^{\frac{\delta\lambda_{1}}{4\sqrt{p}}}=1-2e^{-\frac{1}{4}\log p}.$

Thus for some $C>0$,
$\mathbb{P}\left[|\theta_{j}|<\frac{\delta}{2\sqrt{p}}\mid\lambda\right]\geq
e^{-C\log p}$ for some $C>0$. We also have that $\mathbb{P}[\lambda\in
I_{\lambda}]=\mathbb{P}[\lambda>\lambda_{1}]-\mathbb{P}[\lambda>\lambda_{2}]$.
Next, we record the following upper and lower bounds for the incomplete Gamma
function (Pinelis, 2020),

$\displaystyle\dfrac{x^{a}e^{-x}(x-a-1)}{(x-a)^{2}+a}\leq\Gamma(a,x)=\int_{x}^{\infty}t^{a-1}e^{-t}dt\leq\dfrac{x^{a}e^{-x}}{x-a}.$

Also, we have $(a/e)^{a-1}\leq\Gamma(a)\leq(a/2)^{a-1}$. Set
$\alpha=a_{\lambda}$ and $\beta=b_{\lambda}$. Thus,

$\displaystyle\mathbb{P}(\lambda\in I_{\lambda})$
$\displaystyle=\frac{\beta^{\alpha}}{\Gamma(\alpha)}\Gamma(\alpha,\beta\lambda_{1})-\frac{\beta^{\alpha}}{\Gamma(\alpha)}\Gamma(\alpha,2\beta\lambda_{1})$
$\displaystyle=e^{\log\beta}e^{-(\alpha-1)\log\alpha}e^{(\alpha-1)\log
2}\Gamma(\alpha,\beta\lambda_{1})-e^{\log\beta}e^{-(\alpha-1)\log\alpha}e^{(\alpha-1)}\Gamma(\alpha,2\beta\lambda_{1})$
$\displaystyle\geq
e^{\log\beta}e^{-(\alpha-1)\log\alpha}e^{K(\alpha-1)}(\beta\lambda_{1})^{\alpha}e^{-C\beta\lambda_{1}},\quad\text{for
}\alpha=O(\beta\lambda_{1})$
$\displaystyle=e^{\log\beta}e^{-(\alpha-1)\log\alpha}e^{(\alpha-1)\log
2}e^{\alpha\log\beta\lambda_{1}}e^{-C\beta\lambda_{1}}$ $\displaystyle\geq
e^{-C\beta\lambda_{1}}\geq e^{-Cs\log p},$

for $\beta=O\left(s\sqrt{\frac{s\log p}{np}}\right)$. Hence, we have that,

$\displaystyle\int_{I_{\lambda}}\mathbb{P}\left[\left\lVert\theta_{S^{c}}\right\rVert_{2}<\frac{\delta}{2}\mid\lambda\in
I_{\lambda}\right]\pi(\lambda)d\lambda\geq e^{-Cs\log p},$

for some $C>0$. We now handle the term
$\mathbb{P}[|\theta_{j}-\theta_{0j}|<\frac{\delta}{2\sqrt{s}}]$. Since
$\lambda\sim\text{Gamma}(\alpha,\beta)$, marginal density of $\theta_{j}$ is:

$\displaystyle
f(\theta_{j})=\dfrac{\beta^{\alpha}}{2\Gamma(\alpha)}\dfrac{\Gamma(\alpha+1)}{(|\theta_{j}|+\beta)^{\alpha}}=\dfrac{(\alpha+1)\beta^{\alpha}}{2(|\theta_{j}|+\beta)^{\alpha}}.$

Let $A=\\{\theta_{j}:|\theta_{j}-{\theta_{0j}}|\leq\delta^{\star},j\in S\\}$
where $\delta^{\star}=\frac{\delta}{2\sqrt{s}}$. We have for $\alpha,\beta$ as
chosen above

$\displaystyle\mathbb{P}\left[\theta_{j}\in A\right]$
$\displaystyle=\int_{A}\dfrac{(\alpha+1)\beta^{\alpha}}{2(|\theta_{j}|+\beta)^{\alpha}}d\theta_{j}$
$\displaystyle\geq e^{-\log
2}e^{\alpha\log\beta}e^{\log(\alpha+1)}e^{\log\delta^{\star}}e^{-\alpha(|\theta_{0j}|+\delta^{\star}+\beta)}$
$\displaystyle\geq e^{-C\alpha\log(\frac{1}{\beta})},$

which establishes that
$\pi_{L}\left\\{\theta:\left\lVert\theta-\theta_{0}\right\rVert_{2}<\epsilon\right\\}>0$.
Hence, the theorem follows.

####  S.1.6 Proof of Theorem 6.6

For simplicity, consider one observed sample, i.e. $n=1$. The same arguments
apply for a general sample size. The complete data log-likelihood is:

$\log
p_{\theta}(\mathbf{v},\mathbf{h})=\sum_{j=1}^{p}\theta_{jj}v_{j}+\sum_{k=p+1}^{p+m}\theta_{kk}h_{k}+\sum_{j=1}^{p}\sum_{k=p+1}^{p+m}\theta_{jk}v_{j}h_{k}-\log
z(\theta).$

Given the current iteration $\theta^{(t)}$, let,

$\displaystyle Q(\theta;\theta^{(t)})$
$\displaystyle=\mathbb{E}_{\theta^{(t)}}[\ell_{c}(\theta)]$
$\displaystyle=\sum_{j=1}^{p}\theta_{jj}v_{j}+\sum_{k=p+1}^{p+m}\theta_{kk}\mathbb{E}_{\theta^{(t)}}[h_{k}\mid\mathbf{v}]+\sum_{j=1}^{p}\sum_{k=p+1}^{p+m}\theta_{jk}v_{j}\mathbb{E}_{\theta^{(t)}}[h_{k}\mid\mathbf{v}]-\log
z(\theta)$ $\displaystyle=g(\mathbf{v};\theta,\theta^{(t)})-\log z(\theta),$

where
$\mathbb{E}_{\theta^{(t)}}[h_{k}\mid\mathbf{v}]=\mathbb{P}_{\theta^{(t)}}[h_{k}\mid\mathbf{v}]=\sigma(\theta_{kk}^{(t)}+\sum_{j=1}^{p}\theta_{jk}^{(t)}v_{j})$,
and $g(\mathbf{v};\theta,\theta^{(t)})=Q(\theta;\theta^{(t)})-\log z(\theta)$.
Then the update (9), is obtained by setting $\dfrac{\partial
Q(\theta;\theta^{(t)})}{\partial\theta}=0$. Also, there exists a suitable step
size $\gamma$ such that $\theta^{(t+1)}$ will satisfy
$Q(\theta^{(t+1)};\theta^{(t)})\geq Q(\theta^{(t)};\theta^{(t)})$.
Furthermore, $Q(\theta_{1};\theta_{2})$ is continuous for all $\theta_{1}$ and
$\theta_{2}$. Then from (Wu, 1983, Theorem 2) any generalized EM sequence
$\theta^{(t)}$ converges to a stationary point of $\log
p_{\theta}(\mathbf{v})$.

Consider any such sequence with starting value $\theta^{(0)}$. In what
follows, we write $\left\lVert M\right\rVert_{2}$ to be the Euclidean norm of
a symmetric matrix $M$, i.e. $\left\lVert
M\right\rVert_{F}=\left\lVert\text{vech}(M)\right\rVert_{2}$. Let
$\theta^{(t)}$ and $\theta_{\star}^{(t)}$ denote the sequences obtained by an
exact EM algorithm and a Monte Carlo EM algorithm (which is our case). Then
the EM and the Monte Carlo EM are of the following form:

$\displaystyle\theta^{(t+1)}$ $\displaystyle=\theta^{(t)}+\gamma_{t}\nabla
Q(\theta^{(t)};\theta^{(t)}),$ $\displaystyle\theta_{\star}^{(t+1)}$
$\displaystyle=\theta^{(t)}+\gamma_{t}\nabla
Q^{\star}(\theta_{\star}^{(t)};\theta_{\star}^{(t)}),$

where $\nabla Q^{\star}(\theta;\psi)$ is the Monte Carlo estimate of $\nabla
Q(\theta;\psi)$ for any $\theta,\psi\in\Theta$. Let $T_{N}^{(1)}(\theta)$ and
$T_{N}^{(2)}(\theta)$ denote the Monte Carlo estimates of $\log z(\theta)$ and
$\nabla\log z(\theta)$, respectively. Clearly, $T_{N}^{(1)}(\theta)$ converges
to $\log z(\theta)$, pointwise and almost surely. That this is true for
$T_{N}^{(2)}(\theta)$ as an estimator of $\nabla\log z(\theta)$ follows from
the continuous mapping theorem. Hence, we have,

$\displaystyle\left\lVert\theta^{(1)}-\theta_{\star}^{(1)}\right\rVert_{2}$
$\displaystyle=\left\lVert\gamma_{1}\nabla
Q(\theta^{(0)};\theta^{(0)})-\gamma_{1}Q^{\star}(\theta^{(0)};\theta^{(0)})\right\rVert_{2}$
$\displaystyle\leq\gamma_{1}\left\lVert\nabla\log
z(\theta^{(0)})-T_{N}^{(2)}(\theta^{(0)})\right\rVert_{2}$
$\displaystyle\leq\gamma_{1}\sup_{\theta\in\Theta}\left\lVert\nabla\log
z(\theta)-T_{N}^{(2)}(\theta)\right\rVert_{2}\leq\gamma_{1}\epsilon,$

where taking supremum is possible as the convergence is uniform as we shall
prove later. More generally, if
$\left\lVert\theta^{(t)}-\theta^{(t)}_{\star}\right\rVert_{2}\leq\delta$, for
some $\delta>0$, then,

$\displaystyle\left\lVert\theta^{(t+1)}-\theta_{\star}^{(t+1)}\right\rVert_{2}$
$\displaystyle\leq\left\lVert\theta^{(t)}-\theta_{\star}^{(t)}\right\rVert_{2}+\left\lVert\gamma_{t}\nabla
Q(\theta^{(t)};\theta^{(t)})-\gamma_{t}\nabla
Q^{\star}(\theta^{(t)}_{\star};\theta^{(t)}_{\star})\right\rVert_{2}$
$\displaystyle\leq\delta+\gamma_{t}\left\lVert\nabla
g(\mathbf{v};\theta,\theta^{(t)})-\nabla
g(\mathbf{v};\theta,\theta_{\star}^{(t)})\right\rVert_{2}+\gamma_{t}\left\lVert\nabla\log
z(\theta^{(t)})-T_{N}^{(2)}(\theta_{\star}^{(t)})\right\rVert_{2}$
$\displaystyle\leq\delta+\gamma_{t}\lambda_{\max}[\nabla^{2}g(\mathbf{v};\theta,\tilde{\theta}_{1})]\left\lVert\theta^{(t)}-\theta_{\star}^{(t)}\right\rVert_{2}+\gamma_{t}\left\lVert\nabla\log
z(\theta^{(t)})-T_{N}^{(2)}(\theta_{\star}^{(t)})\right\rVert_{2}$
$\displaystyle\leq\delta+\gamma_{t}C_{1}\delta+$
$\displaystyle\gamma_{t}\left\lVert\nabla\log z(\theta^{(t)})-\nabla\log
z(\theta_{\star}^{(t)})\right\rVert_{2}+\gamma_{t}\left\lVert\nabla\log
z(\theta_{\star}^{(t)})-T_{N}^{(2)}(\theta_{\star}^{(t)})\right\rVert_{2}$
$\displaystyle\leq\delta+\gamma_{t}C_{1}\delta+\gamma_{t}\lambda_{\max}[\nabla^{2}\log
z(\tilde{\theta}_{2})]+\gamma_{t}\epsilon$
$\displaystyle\leq\delta+\gamma_{t}C_{1}\delta+\gamma_{t}C_{2}\delta+\gamma_{t}\epsilon\leq
C\delta\gamma_{t},$

where
$\tilde{\theta}_{i}=\theta^{(t)}+\rho_{i}(\theta^{(1)}-\theta_{\star}^{(t)})$
for $\rho_{i}\in(0,1)$, and $i=1,2$. That the Hessian of
$g(\mathbf{v};\theta,\phi)$ is bounded by a global constant follows from the
fact that $\nabla^{2}g(\mathbf{v};\theta,\phi)$ is continuous on the compact
support of $\phi$, and $\mathbf{v}$ is bounded from the Gershgorin circle
theorem (Varga, 2010). Similarly, the Hessian of $\log z(\theta)$ is also
bounded on the compact support of $\theta$. Hence, if $C\gamma_{t}\leq 1$, for
all $t\geq 1$, then,

$\displaystyle\sum_{t=1}^{\infty}\left\lVert\theta^{(t)}-\theta_{\star}^{(t)}\right\rVert_{2}\leq\epsilon\sum_{t=1}^{\infty}{\gamma_{t}}=0.$

In other words, for sufficiently large $t$,

$\left\lVert\theta^{(t)}-\theta^{(t)}_{\star}\right\rVert_{2}=O\left(\frac{1}{t^{1+\eta}}\right),$

for some $\eta>0$. It then remains to show that the convergence of
$T_{N}^{(2)}(\theta)$ to $\nabla\log z(\theta)$ is uniform over $\Theta$.
Observe that:

$T_{N}^{(2)}(\theta)=\dfrac{\sum_{i=1}^{N}\dfrac{\nabla
q_{\theta}(Y_{i})}{q_{\phi}(Y_{i})}}{\sum_{i=1}^{N}\dfrac{q_{\theta}(Y_{i})}{q_{\phi}(Y_{i})}}=\dfrac{\sum_{i=1}^{N}D_{1}(Y_{i};\theta)}{\sum_{i=1}^{N}D_{2}(Y_{i};\theta)},\,\,(\text{say})$

where $Y_{i}\in\\{0,1\\}^{p+m}$ and $\phi=\text{diag}(\theta)$. Then,

$\displaystyle\nabla T_{N}^{(2)}(\theta)=\dfrac{\sum_{i=1}^{N}\nabla
D_{1}(Y_{i};\theta)}{\sum_{i=1}^{N}D_{2}(Y_{i};\theta)}-\dfrac{\sum_{i=1}^{N}D_{1}(Y_{i};\theta)}{\left(\sum_{i=1}^{N}D_{2}(Y_{i};\theta)\right)^{2}}\sum_{i=1}^{N}\nabla
D_{2}(Y_{i};\theta).$

Therefore,

$\displaystyle\left\lVert\nabla T_{N}^{(2)}(\theta)\right\rVert_{2}$
$\displaystyle\leq\dfrac{1}{[\sum_{i=1}^{N}D_{2}(Y_{i};\theta)]^{1/2}}\left\lVert\sum_{i=1}^{N}\nabla
D_{1}(Y_{i};\theta)\right\rVert_{2}+\left[\dfrac{\sum_{i=1}^{N}D_{1}(Y_{i};\theta)}{\left(\sum_{i=1}^{N}D_{2}(Y_{i};\theta)\right)^{2}}\right]^{1/2}\left\lVert\sum_{i=1}^{N}\nabla
D_{2}(Y_{i};\theta)\right\rVert_{2}$ $\displaystyle\leq C,$

for some $C>0$, due to the fact that $D_{1},D_{2},\nabla D_{1},\nabla D_{2}$
are all continuous with respect to $\theta$ and $Y_{i}$ is a bounded random
variable, and since $\Theta$ is compact. Since by definition $\Theta$ is also
convex, and $T_{N}^{(2)}(\theta)$ is infinitely differentiable, uniform
convergence follows from (Newey, 1991, Corollary 2.2).

###  S.2 A Simple Accept–Reject Sampler for $p_{\theta}(\cdot)$

Generating samples from a Boltzmann distribution, of which PEGMs are special
cases, remains a challenge; especially without resorting to computationally
expensive MCMC techniques. However, our approach also accomplishes this goal.
Specifically, for a given $\theta$, suppose we simulate $Y\sim
p_{\phi}(\cdot),$ where $\phi=\mathrm{diag}(\theta)$. This is cheap, and does
not require MCMC. Moreover, $z(\phi)$ is analytically available. Next, recall,
$p_{\theta}(x)=q_{\theta}(x)/z(\theta)$. Thus, an accept–reject sampler for
$p_{\theta}(\cdot)$ consists of the following two steps:

  1. 1.

Simulate $Y\sim p_{\phi}(\cdot)$

  2. 2.

Accept $Y$ with probability $r=q_{\theta}(Y)/\\{Mp_{\phi}(Y)\\}$ where $M$ is
such that $q_{\theta}(Y)<Mp_{\phi}(Y)$ for all $Y$.

Thus, $M$ must be chosen to satisfy the criterion:

$\displaystyle\exp\left(\sum_{j}\theta_{j}Y_{j}+\sum_{j\neq
k}\theta_{jk}Y_{j}Y_{k}\right)$
$\displaystyle<M\exp\left(\sum_{j}\theta_{j}Y_{j}\right)/z(\phi),$
$\displaystyle\text{i.e., }\exp\left(\sum_{j\neq
k}\theta_{jk}Y_{j}Y_{k}\right)$ $\displaystyle<Mz(\phi)^{-1}.$

However, a simple choice of $M$ is available for several PEGMs. For example,
in a Poisson graphical model, all off-diagonal elements of a valid $\theta$
are negative, and one may take $Mz(\phi)^{-1}=1$, i.e., $M=z(\phi)$.
Similarly, for Ising, where $Y_{j}=\\{0,1\\}$ for all $j$, a simple choice is
$\\{M:Mz(\phi)^{-1}>\exp(p^{2}\max_{jk}\theta_{jk})\\}$. Other, less obvious
choices of $M$ may result in better average acceptance probabilities. We do
not explore the details in this work.

###  S.3 Algorithmic Details for Section 4

####  S.3.1 Tuning $\lambda$ in (8) and model selection

Here we describe how the value of the tuning parameter $\lambda$ is selected
in (8). Generally, pseudo-likelihood methods have a regression structure in
that the likelihood of $X_{j}\mid X_{-j}$ is a generalized linear model for
PEGMs. As a result, out-of-sample prediction is typically used for tuning
parameter selection. This is, however, not the case when a full-likelihood
analysis is implemented. In our work, we instead consider the out-of-sample
log-likelihood to select the value of $\lambda$ given estimated
$\hat{\theta}_{\lambda}$ is available for a grid of values of $\lambda$:
$[\lambda_{l},\lambda_{u}]$. Specifically, we divide the observed data
randomly into $K$-folds, of which $K-1$ folds are used to obtain
$\hat{\theta}_{\lambda}$ for each value of $\lambda$. The log-likelihood
$\ell(\hat{\theta}_{\lambda})$ is computed on the remaining fold. The process
is repeated for each fold, and we choose $\hat{\theta}_{\lambda}$ to be our
estimator for which we obtain the maximum out-of-sample log-likelihood. This
generally results in improved Frobenius norm estimates of $\theta$ as reported
in Table 7. However, it is well known that cross-validated estimates in
penalized estimation procedures may provide sub-optimal performance in terms
of model selection (Bühlmann and Van De Geer, 2011). To address this, we
construct an estimator of $\theta$ that selects edges between nodes which are
most often selected. Suppose, for $R$ choices of $\lambda$, we obtain
estimates
$\hat{\mathbf{S}}_{\lambda_{r}}=\\{(j,k):\hat{\theta}_{\lambda_{r},jk}\neq
0\\}$. Define
$I_{jk}^{r}=\mathbb{I}\\{(j,k)\in\hat{\mathbf{S}}_{\lambda_{r}}\\}$. Then, we
estimate the structure of the graph $G$ by:

$\hat{\mathbf{S}}=\left\\{(j,k):\frac{1}{R}\sum_{r=1}^{R}I^{r}_{jk}>\pi_{thr},\right\\}$
(S.4)

where $0<\pi_{thr}<1$. The procedure closely resembles stability selection
(Meinshausen and Bühlmann, 2010), except here we do not perform subsampling.
Our experiments with and without subsampling resulted in similar performance.
Hence, we omit that step. For the threshold $\pi_{thr}$, we use $0.6$. Results
from this choice are reported under the MCC column in Table 7. For fair
comparison, we also implemented this structure learning method with pseudo-
likelihood based approaches.

####  S.3.2 HMC for PEGMs

Suppose we observe iid data $\mathbf{X}=(X_{1\bullet},\ldots,X_{n\bullet})$
from some $\text{PEGM}(\theta)$. Consider the following prior structure over
$\theta$: $\pi(\theta)=\int\left[\prod_{j\leq
k}\pi(\theta_{jk}\mid\tau)\right]g(\tau)d\tau$, i.e. the elements
$\theta_{jk}$ are conditionally independent given a global shrinkage parameter
$\tau$. The posterior is then:

$\pi(\theta\mid\mathbf{X})\propto\left[\prod_{i=1}^{n}\dfrac{q_{\theta}(X_{i\bullet})}{z(\theta)}\right]\pi(\theta).$
(S.5)

We make two simplifying assumptions. First, we fix $\tau=1$, noting that our
procedure can be easily incorporated in the general case when $\tau\sim g$
with an additional layer of sampling. Second, we assume that the convexity of
$\Theta$ can be represented as $C(\theta)\geq 0$ for a suitable choice of a
differentiable $C(\theta)$.

Suppose $\pi(\theta_{jk})$ admits the marginal representation
$\pi(\theta_{jk})=\int\pi(\theta_{jk}\mid\rho_{jk})\pi(\rho_{jk})d\rho_{jk},$
with respect to some latent variables $\rho_{jk}$, and
$\pi(\theta_{jk}\mid\rho_{jk})$ is differentiable. This is indeed the case for
any prior that is a scale mixture of Gaussian. Then samples from the augmented
posterior $\pi(\theta,\\{\rho_{jk}\\}_{j\leq k}\mid\mathbf{X})$ can be
obtained by employing a Gibbs sampler that iterates between sampling from
$\pi(\theta\mid\\{\rho_{jk}\\}_{j\leq k},\mathbf{X})$ and
$\pi(\\{\rho_{jk}\\}_{j\leq k}\mid\theta,\mathbf{X})$. We use a HMC sampler to
sample from $\pi(\theta\mid\\{\rho_{jk}\\}_{j\leq k},\mathbf{X})$ and sample
$\rho_{jk}$’s in a block since these are conditionally independent. To ensure
that the HMC sampler lives within $C(\theta)$, we adopt the constrained HMC
sampler from Betancourt (2011). A full description of the transition operation
of the sampler is provided in Algorithm S.1 where we write
$\theta=\text{vech}(\theta)$ with a slight abuse of notation. Let $U(\theta)$
denote the negative of the logarithm of $\pi(\theta\mid\\{\rho_{jk}\\}_{j\leq
k},\mathbf{X})$ up to proportionality constants, i.e.
$U(\theta)=-\sum_{i=1}^{n}\log q_{\theta}(X_{i\bullet})+n\log
z(\theta)-\sum_{j\leq k}\log\pi(\theta_{jk}\mid\rho_{jk})$, and
$\nabla_{\theta}U(\theta)$ denote its gradient which is available following
Proposition 3.1 and due to the assumption that $\pi(\theta_{jk}\mid\rho_{jk})$
is differentiable. The new state generated according to Algorithm S.1 is
accepted with the standard Metropolis correction (Neal et al., 2011, Section
5.3.2.1).

Algorithm S.1 Constrained HMC sampler for updating
$\theta\mid\\{\rho_{jk}\\}_{j\leq k},\mathbf{X}$

Input: $\theta$ [current state], $\epsilon$ [step size], $L$ [number of
leapfrog steps]

Output: $\theta$ [next update]

Sample $p\sim{\mathrm{N}}(0,I)$ and set
$p=p-\frac{1}{2}\epsilon\nabla_{\theta}U(\theta)$,

while $l\leq L$ do

$\theta=\theta+\epsilon p$

if $C(\theta)\geq 0$ then

$p\leftarrow p-\epsilon\nabla_{\theta}U(\theta)$

else

$r\leftarrow\dfrac{\nabla C(\theta)}{\left\lVert\nabla C(\theta)\right\rVert}$

$p\leftarrow p-2(r^{\prime}p)r$

end if

end while

$\theta\leftarrow\theta+\epsilon p$

$p\leftarrow p-\frac{1}{2}\epsilon\nabla_{\theta}U(\theta)$

###  S.4 Additional Numerical Experiments

####  S.4.1 Likelihood vs. pseudo-likelihood based CIs in low dimensions

Results in Table S.4.1 substantiate that likelihood-based inference provides
the best results in low-dimensions. For this, we generate data according to
the same mechanism as in Section 7 of the main manuscript and consider
$p=3,5$. We report the average coverage and width of 95% confidence intervals
using bootstrap, over 100 replications. Results in Table S.4.1 indicate that
while both maximum likelihood estimate (MLE) and pseudo maximum likelihood
estimate (MPLE) have comparable coverage probabilities, widths of confidence
intervals higher for PMLE, often by a factor as large as 10. No penalization
was used in low dimensions, either for full likelihood or for
pseudolikelihood.

ccc——cc—cc——cc—cc Ising  PGM
Setting $p$ $n$ Coverage  Width  Coverage  Width
MLE MPLE MLE MPLE MLE MPLE MLE MPLE
LD 3 100 0.88 (0.08) 0.87 (0.07) 1.33 (0.38) 4.02 (2.86) 0.86 (0.35) 0.87
(0.34) 0.36 (0.16) 3.66 (3.88)
5 0.80 (0.12) 0.82 (0.13) 1.34 (0.22) 4.98 (2.52) 0.66 (0.47) 0.67 (0.47) 0.37
(0.22) 4.17 (4.06)

Table S.1: Coverage and width of confidence intervals of maximum likelihood
(MLE) and maximum pseudo-likelihood (MPLE) procedures for two PEGMs: Ising and
PGM.

####  S.4.2 A comparison of RBM and BM representational power

Existing literature on training performance of the CD algorithm and its other
versions focus mostly on benchmark datasets. Nonetheless, in this section we
report the results of CD training versus the proposed approach in controlled
simulation experiments. We compare the quality of the corresponding density
estimate in terms of the total variation distance:
$\text{TV}(p_{\theta_{0}},p_{\theta})=\sum_{\mathbf{v}\in\\{0,1\\}^{p}}|p_{\theta_{0}}(\mathbf{v})-p_{\hat{\theta}}(\mathbf{v})|$
where the dimension of $\theta_{0}$ is the data generating value of the
parameter and $\hat{\theta}$ is the estimate obtained from the two methods,
noting that the dimensions of $\theta_{0}$ and $\hat{\theta}$ might be
different given that the number of hidden variables is unknown in reality. We
consider two cases $(p,m)=(2,2)$ and $(p,m)=(3,4)$ and set the sample size
$n=1000$. Elements of the true data generating $\theta_{0}$ are set to 0 with
probability 0.5 and are drawn from $\text{Uniform}[-1,1]$ with probability
0.5, respecting the bipartite structure of the graph. When implementing CD we
set $k=1$, although we did not notice any significant difference in results
for $k$ up to 10. We report in Table S.4.2 results averaged over 50
replications where we abbreviate the proposed method here by FL (for _full
likelihood_). All hyperprameters for training for the two algorithms are kept
same, including the step size and convergence criteria. Two key observations
can be made from the results. First, for increasing values of the number of
hidden variables $m$, the fit is better, which is expected. Second, the
performance of FL is better compared to CD in almost all cases, albeit by a
small margin.

We next consider experiments where we study the representational power of RBMs
to that of BMs. Since both method learn distributions of visible variables, we
want to investigate whether allowing connections within layers as in BM (see
also Figure 1), leads to similar or better learning power for a smaller number
of hidden variables. For this we consider two situations $p=3,5$ which is the
dimension of the visible variable. We generate the visible variables from a
multivariate probit model where $y_{j}=\mathbb{I}(z_{j}>0)$ and