APPROXIMATING FACES OF MARGINAL POLYTOPES
IN DISCRETE HIERARCHICAL MODELS
By Nanwei Wang

arXiv:1603.04843v2 [] 18 Jan 2018

York University, Toronto, Canada
and
By Johannes Rauh
York University, Toronto, Canada
Max Planck Institute for Mathematics in the Sciences, Germany
and
By Hélène Massam
York University, Toronto, Canada
The existence of the maximum likelihood estimate in hierarchical loglinear models is crucial to the reliability of inference for this
model. Determining whether the estimate exists is equivalent to finding whether the sufficient statistics vector t belongs to the boundary
of the marginal polytope of the model. The dimension of the smallest face Ft containing t determines the dimension of the reduced
model which should be considered for correct inference. For higherdimensional problems, it is not possible to compute Ft exactly. Massam and Wang (2015) found an outer approximation to Ft using
a collection of sub-models of the original model. This paper refines
the methodology to find an outer approximation and devises a new
methodology to find an inner approximation. The inner approximation is given not in terms of a face of the marginal polytope, but in
terms of a subset of the vertices of Ft .
Knowing Ft exactly indicates which cell probabilities have maximum likelihood estimates equal to 0. When Ft cannot be obtained
exactly, we can use, first, the outer approximation F2 to reduce the
dimension of the problem and, then, the inner approximation F1 to
obtain correct estimates of cell probabilities corresponding to elements of F1 and improve the estimates of the remaining probabilities
corresponding to elements in F2 \ F1 . Using both real-world and simulated data, we illustrate our results, and show that our methodology
scales to high dimensions.
Keywords: existence of the maximum likelihood estimate, marginal
polytope, faces, facial sets, extended maximum likelihood estimate.

∗

H. Massam gratefully acknowledges support from an NSERC Discovery Grant

1

2

WANG & RAUH & MASSAM

1. Introduction. Discrete hierarchical models are an essential tool for
the analysis of categorical data given under the form of a contingency table.
The study of these models goes back more than a century, and a detailed
history of their development is given in Fienberg and Rinaldo (2007). Nowadays, discrete hierarchical models are used for the analysis of large sparse
contingency tables where many, if not most, of the entries are small or zero
counts. It is well-known that in such cases, the maximum likelihood estimate
(henceforth abbreviated MLE) of the parameters may not exist. The non existence of the MLE has problematic consequences for inference, clearly for
estimation, but also for testing and model selection. Fienberg and Rinaldo
(2012) list the statistical implications of the non existence of the MLE, such
as the unreliability of the estimates of some of the parameters or the usage of the wrong degrees of freedom for testing one model against another.
Geyer (2009) describes the problems attached to the nonexistence of the
MLE and presents an R program that yields meaningful confidence intervals and tests. Letac and Massam (2012) study the statistical implications
of the nonexistence of the MLE on model selection in Bayesian inference.
Fienberg and Rinaldo (2012) also give necessary and sufficient conditions for the existence of the MLE, which are not restricted to hierarchical models, but apply to all discrete exponential families (loglinear models).
These conditions are extensions of results given earlier by Haberman (1974),
Barndorff-Nielsen (1978) and Eriksson et al. (2006). They are essentially as
follows. Denote by I the set of outcomes of a statistical experiment, that
is the set of cells of the contingency table where the data is classified. Let
I+ = {i1 , . . . , iN } be the outcome of N independent repetitions of the experiment within either a multinomial or Poisson setting. The data I+ is
summarized
by the vector t of sufficient statistics, which is of the form
P
f
t= N
j=1 ij for some vectors fi , i ∈ I, determined by the given hierarchical loglinear model. Under these assumptions, the distribution of the data
belongs to a natural exponential family with density
(1)

f (i1 , . . . , iN ; θ) = exp{hθ, ti − N k(θ)}

with respect to the counting measure, where θ is a loglinear parameter.
To each exponential family is associated a polytope P, called the marginal
polytope, which is the convex hull of the vectors fi , i ∈ I. Furthermore,
P contains all possible realizations of Nt for arbitrary repetitions of the
statistical experiment. For given data and a given hierarchical model, the
MLE then exists if and only if Nt belongs to the relative interior of P. If
the MLE does not exist, then Nt belongs to the relative interior of a face
denoted Ft . It is the smallest face of P containing Nt , and it is proper

APPROXIMATING FACES OF MARGINAL POLYTOPES

3

(i.e. Ft 6= P). Thus, determining whether, for a given data set, the MLE of
the parameter of a discrete hierarchical loglinear model exists is equivalent
to determining whether Nt belongs to a proper face of P. The parameter
may be the loglinear parameter θ, or the cell probabilities p = (p(i), i ∈ I)
obeying the constraints of the model and in 1-1 correspondence with θ. The
MLE can thus be thought of in terms of θ, or in terms of p.
If the MLE does not exist, it is still possible to compute the extended
MLE (EMLE) (Barndorff-Nielsen, 1978; Lauritzen, 1996; Csiszár and Matúš,
2008), which is a probability distribution that maximizes the likelihood over
the closure of the hierarchical model (that is, the EMLE can be approximated arbitrarily well by distributions from the hierarchical model). The
support of the EMLE is given by the set Ft = {i ∈ I : fi ∈ Ft }, called
the facial set of Ft . When this support is known, computing the EMLE
is equivalent to an ordinary MLE computation on a smaller exponential
family EFt , of dimension dim(Ft ), generated by a measure with support Ft
(Geyer, 2009). Therefore, precise knowledge of Ft and Ft yields which is the
proper dimension of the model to be used in testing and which outcomes
i ∈ I are attributed a probability of 0 by the EMLE, and it allows us to compute the EMLE. One should also note that the usual regularity conditions
used for the asymptotic properties of the MLE, which are not satisfied for
the given model when the MLE does not exist, are satisfied for the reduced
model EFt .
The problem is then to find Ft . This is easy when the face lattice of P is
known or can be computed using a standard discrete geometry toolkit such
as, for example, polymake (Gawrilow and Joswig, 2000). For some classes of
marginal polytopes, the face lattice is known, for example for decomposable
models and no-three-way-interaction models with small variables (Vlach,
1986). For binary variables, the marginal polytope is a cut polytope (Deza
and Laurent, 2009). Other authors have studied convex support polytopes,
which replace marginal polytopes for more general exponential families. Notably, many such polytopes have been described for exponential random
graph models, see, for example, Karwa and Slavkovic (2016) and papers
cited therein. When the face lattice of P cannot be computed, algorithms
to compute Ft that are based on linear programming have been proposed
by Eriksson et al. (2006), by Geyer (2009), and by Fienberg and Rinaldo
(2012). These methods, however, become computationally infeasible in large
dimensions, which happens, in our experience, for hierarchical models when
the set of random variables V contains more than 16 binary variables (or
correspondingly fewer larger variables).
For larger models, Massam and Wang (2015) propose to approximate Ft

4

WANG & RAUH & MASSAM

by relating it to faces of smaller hierarchical models as follows. A hierarchical model for the discrete random variable X = (Xv , v ∈ V ) is determined
by a set of interactions among its components Xv , v ∈ V, that is represented by a simplicial complex ∆. Massam and Wang (2015) consider subsets
Vi , i = 1, . . . , k, of V containing less than 16 variables and the hierarchical
models Markov with respect to the induced simplicial subcomplexes. Linear
programming can be used to compute the smallest faces Fti containing the
corresponding sufficient statistic ti , i = 1, . . . , k. These faces, which a priori
are faces of the marginal polytopes of the submodels, naturally correspond
to faces of the original marginal polytope. Massam and Wang (2015) prove
that the intersection of these is a face F2 of P containing Ft . Thus, if F2 is
a proper face of P, then Ft is necessarily a proper face, and therefore the
MLE does not exist. While Massam and Wang (2015) work with graphical
models, we show that their results are also true for hierarchical loglinear
models.
We call F2 an outer approximation to Ft . This is similar to the notion
of an outer approximation in optimization, which describes a polytope that
contains the original polytope of interest. While the outer approximation
polytope in optimization usually has the same dimension as the original
polytope, the outer approximation face F2 does not necessarily have the
same dimension as Ft .
The purpose of this paper is to add to this outer approximation F2 an
inner approximation F1 that is a subset of Ft . While F2 is derived from
looking at simplicial subcomplexes of ∆, the inner approximation is constructed by enlarging the simplicial complex through added interactions. In
particular, we propose a process of “completing a separator,” which leads to
a decomposable simplicial complex which, in turn, can be studied by looking
at the sub-simplices corresponding to its components with a small number
of vertices in V . Thus, both F1 and F2 can be obtained by computing facial
sets on smaller hierarchical models involving fewer nodes.
The inner and outer approximations to Ft satisfy
F1 ⊂ Ft ⊂ F2 .
Clearly, we want F1 as large as possible and F2 as small as possible to have as
much information about Ft as possible. In our simulations, we observe that
Ft = F2 most of the time and that Ft = F1 quite often. The approximations
F1 and F2 allow to bound the dimension of Ft , and thus knowledge from
F1 and F2 can be taken into account whenever the dimension of Ft plays a
role, for example in hypothesis testing.
When the MLE does not exist, even though the maximum likelihood

APPROXIMATING FACES OF MARGINAL POLYTOPES

5

procedure cannot be used to obtain a point estimate for the parameter
vector θ, some of its components θj may still be finite and well-defined
in this situation. In Section 4 we introduce a loglinear parametrization µ,
different from θ, that allows to say precisely which parameter combinations
have a finite well-defined limit and thus remain meaningful for statistical
inference. Moreover, we demonstrate that even when Ft is unknown, the
parametrization µ can be adjusted to incorporate knowledge that is available
in the form of inner and outer approximations F1 and F2 .
We extend the work of Fienberg and Rinaldo (2012) and that of Geyer
(2009) in several directions: first, we construct approximations to Ft in high
dimensions when a direct computation of Ft is not feasible. Second, we
explicitly identify all parameter combinations that remain finite and meaningful when the MLE does not exist and Ft is known, and we also discuss
what can be said when only approximations to Ft are available.
The remainder of this paper is organized as follows. In Section 2, we give
preliminaries on hierarchical models, and faces and facial sets. Section 3
contains the original methodology to obtain the approximations F1 and F2 .
In Section 4, we show how to use F1 and F2 to identify the parameters of
the hierarchical models that can be estimated and those that cannot. In
Section 5, we present two examples. A simulated data set is used to assess
how often our approximations succeed to identify the true facial set Ft . The
NLTCS data set, studied by Dobra et al. (2003) and Dobra and Lenkoski
(2011), illustrates how the outer approximation F2 improves estimates of
cell probabilities and log-linear parameters. Both of these examples have
16 nodes. In Section 6, we discuss how to apply the methodology to larger
models and how to use for inference the information that it yields. Two
examples illustrate this: simulated data from the graphical model of the
5 × 10 grid, and the real-world data set of voting records in the US Senate.
Appendix A describes the concrete parametrization that we use in the
examples. Appendix B discusses the case of two binary variables to illustrate what happens to the usual parameters when the MLE does not exist.
Appendix C discusses how to further improve the parametrization µL introduced in Section 2. Appendices D and E give further results for the examples
from Section 5. Appendix F briefly summarizes the linear programming algorithm to compute Ft by Fienberg and Rinaldo (2012). Appendix G contains
proofs of the results of Barndorff-Nielsen (1978) that describe the closure of
an exponential family and the EMLE.
Our results apply not only to hierarchical models, but to arbitrary discrete
exponential families. In this paper, the focus is on hierarchical and graphical

6

WANG & RAUH & MASSAM

models, which are the main application, and for which the construction of the
inner and outer approximations can be described in terms of the underlying
simplicial complex or graph.
2. Preliminaries. In the following four subsections, we recall basic
facts about hierarchical models, discrete exponential families, polytopes and
the closure of exponential families, and we define the extended MLE.
2.1. Hierarchical models and discrete exponential families. For details
and proofs on the material in this subsection, we refer to Letac and Massam
(2012) and Rauh et al. (2011). Let X = (Xv , v ∈ V ) be a discrete random
vector with components indexed by a finite set V = {1, . . . , p}. Each variable
Xv takes
Q values in a finite set Iv , v ∈ V . The vector X takes its values in
I = v∈V Iv , the set of cells i = (iv , v ∈ V ) of a p-dimensional contingency
table. Q
For any D ⊆ V , the subvector XD = (Xv )v∈D takes its values in
ID = v∈D Iv . The D-marginal cell of i ∈ I will be denoted by iD = (iv )v∈D .
The corresponding restriction is the coordinate projection map i 7→ iD and
is denoted by πD .
Let ∆ be a simplicial complex on V , that is, ∆ is a set of subsets D ⊂ V
such that D ∈ ∆ and D0 ⊂ D implies D0 ∈ ∆. The joint distribution of X
is hierarchical with underlying simplicial complex ∆ (or generating set ∆)
if the probability p(i) = P (X = i) of a single cell i = (iv , v ∈ V ) is of the
form
X
(2)
log p(i) =
θD (iD )
D∈∆

where θD (iD ) is a function of the marginal cell iD = (iv , v ∈ D) only. To
make precise the dependence on θ, we sometimes write pθ (i) instead of p(i).
The set of all such distributions E∆ := {pθ } is called the hierarchical model
of ∆.
Equation (2) is essentially a linear condition on log p(i). It is possible
to parametrize the hierarchical model using a finite vector of parameters
(θj )j∈J such that
X
(3)
log pθ (i) =
θj aj,i − k(θ),
j∈J

where A∆ = (aj,i )P
is a fixed real matrix (depending only
j∈J,i∈I P
Pon ∆) and
where k(θ) = log i exp( i θj aj,i ) ensures the normalization i∈I pθ (i) =
1. This parametrization is not unique. In the examples, we use an explicit
parametrization that is used, for example, by Letac and Massam (2012). For
convenience, we recall this parametrization in Appendix A.

APPROXIMATING FACES OF MARGINAL POLYTOPES

7

An important subclass of hierarchical models is the class of graphical
models. Let G = (V, E) be an undirected graph with vertex set V and edge
set E. A subset D ⊆ V is a clique of G if for any i, j ∈ V , i 6= j, the edge (i, j)
is in E. The set of cliques of G, denoted by ∆(G), is a simplicial complex. The
graphical model of G is defined as the hierarchical model of ∆(G). Graphical
models are important because of their interpretation in terms of conditional
independence, see Lauritzen (1996).
Hierarchical models are examples of discrete exponential families, see
Barndorff-Nielsen (1978); Fienberg (2007); Rauh et al. (2011). Generalizing (3), let I and J be finite sets and let A ∈ RJ×I be a real matrix. Denote
the columns of A by fi , i ∈ I. The discrete exponential family corresponding
to A, denoted by EA , consists of all probability distributions on I that are
of the form


(4) pθ (i) = exp hθ, fi i − k(θ) = exp (At θ)i − k(θ) ,
θ ∈ RJ ,
P
P
where, as above, k(θ) = log i exp( j θj aj,i ). It is convenient to write Ã

for the (1 + |J|) × I matrix with columns equal to f1i , i ∈ I, and to set
θ0 := −k(θ) and θ̃ = (θ0 , θ) (as a column vector). Then (4) rewrites to

(5)
pθ (i) = exp Ãt θ̃ ,
θ ∈ RJ .
Both A and Ã are called design matrices of the model. The convex hull of
the columns fi , i ∈ I, is called the convex support polytope, denoted by PA .
In the case of a graphical or hierarchical model, P∆ := PA∆ is called a
marginal polytope.
The parametrization θ → pθ is identifiable if and only if Ã has full rank.
If Ã does not have full rank, then one can drop rows of A to obtain a
submatrix A0 such that Ã0 has full rank. This is equivalent to setting certain
parameters to zero until the remaining parameters are identifiable.
Later, the following reparametrization will be useful: select an element
of I, which we will denote by 0. Let A0 be the matrix with columns fi − f0 ,
i ∈ I \{0}. It is not difficult to see that A and A0 define the same exponential
family (since Ã and Ã0 have the same row span). Let h0 = rank(A0 ) =
rank(Ã0 ) − 1, and select a set L of h0 linearly independent vectors among
the columns of A0 . For i ∈ L, let µi = µi (θ) := hθ, fi − f0 i, and let µL =
(µi , i ∈ L). Then the µL are identifiable parameters on EA : in fact, their
number is equal to h0 , and they are independent by construction.
It is possible to extend the definition of µi (θ) to all i ∈ I. Note that only
the parameters µi with i ∈ L are free parameters, while the parameters

8

WANG & RAUH & MASSAM

µi with i ∈ I \ L are linear functions of µL . The µi can be interpreted as
log-likelihood ratios:
µi (θ) = log

pθ (i)
,
pθ (0)

µ0 (θ) = 0.

Let n = (n(i), i ∈ I) be an I-dimensional column vector of cell counts
summarizing the outcome of a statistical experiment. Then


N
(6)
Ãn =
and An = t,
t
P
where N = i∈I n(i) is the total cell counts and t is the column vector of
Q
sufficient statistic. The likelihood i∈I pθ (i)n(i) can be written under the
form of a natural exponential family. Indeed,
Y

X
pθ (i)n(i) = exp Ãn, θ̃ = exp
θj tj − N k(θ) .
i∈I

j∈J

The log-likelihood function for the loglinear parameters θ of EA is therefore
X
(7)
l(θ) =
θj tj − N k(θ).
j∈J

It is well-known that l(θ) is concave. If the parameters are identifiable, then
it is strictly concave. We can also express the log-likelihood as a function of
µ = (µi , i ∈ I):
X
X
X
(8)
l(µ) =
n(i) log p(i) =
n(i)µi − N log(
exp µi ).
i∈I

i∈I

i∈I

As stated before, only a subset µL of the parameters µ are independent, and
the remaining µi , i ∈
/ L, can be expressed as linear functions of µL .
2.2. The convex support and its facial sets. We next recall some facts
about facial sets. We refer to Ziegler (1998) for a general introduction to
polytopes and their face lattices.
The convex support polytope PA is defined as the convex hull of a finite number of points fi , i ∈ I. It is of interest to know which subsets of
{fi }i∈I lie on a given face F. Thus, we describe a face F by identifying the
corresponding facial set F = {i ∈ I : fi ∈ F}. For any subset S ⊆ I, denote by FA (S) the smallest facial set that contains S. The intersection of
facial sets is again facial, and so FA (S) is well-defined. When PA = P∆ is a
marginal polytope, we abbreviate FA∆ (S) by F∆ (S).

APPROXIMATING FACES OF MARGINAL POLYTOPES

9

As mentioned in the introduction, to derive the inner approximation F1 to
Ft and its outer approximation F2 , we need to consider sub-models of a given
model. When one exponential family EA0 is a subset of another family EA ,
then the convex support polytope PA0 is a linear projection of PA , and the
columns fi0 of A0 are indexed by the same set I as the columns fi of A. Since
inverses of linear projections preserve faces, it follows from basic results
about polytopes that FA (S) ⊆ FA0 (S); see Chapter 1 in Ziegler (1998). For
hierarchical models, these facts are summarized in the following result:
Lemma 2.1. Let ∆ and ∆0 be simplicial complexes on the same vertex
set with ∆0 ⊆ ∆. Then P∆0 is a coordinate projection of P∆ . The inverse
image of any face of P0 is a face of P. Moreover, for any S ⊆ I, we have
F∆ (S) ⊆ F∆0 (S).
Remark 2.2. It is convenient to embed PA in a vector space with one
additional dimension using a map Rh → Rh+1 , t 7→ t̃ := (1, t). This has the
advantage that all defining inequalities are brought into a homogeneous form
with vanishing constant: note that hg, fi i − c = hg̃c , f˜i i, where g̃c := (−c, g).
When a defining inequality of a face F is given, its facial set F can be
obtained by checking whether fi ∈ F for each i ∈ I. In the other direction,
when a facial set F is given, it is much more difficult to compute a defining
inequality of the corresponding face F. However, it is straightforward to
compute the linear equations defining F: the set of such equations 0 =
hg, xi − c = hg̃, x̃i corresponds to the set of vectors g̃ ∈ ker ÃtF , where ÃF is
the matrix obtained from Ã by dropping the columns not in F .
2.3. The closure of an exponential family and existence of the MLE. For
a family EA and cell counts n = (n(i) : i ∈ I) given as above, a parameter
value θ∗ is an MLE if it is a global maximum of l(θ). An MLE need not
exist, since the domain of the parameters θ is unbounded. The likelihood
can also be written as a function of cell probabilities. For any probability
distribution p on I let
Y
˜l(p) = log{
p(i)n(i) }.
i∈I

Then l(θ) = ˜l(pθ ), and θ∗ is an MLE if and only if pθ∗ maximizes ˜l subject
to the constraint that p belongs to EA , i.e. is of the form (4). When ˜l has
no maximum on EA , we can pass to the topological closure EA . It can be
characterized in terms of the convex support polytope PA and its facial sets
as follows:

10

WANG & RAUH & MASSAM

Theorem
S 2.3 (Barndorff-Nielsen (1978)). The topological closure of EA
is EA = F EF,A , where F runs over all facial sets of PA and where EF,A
consists of all probability distributions of the form pF,θ , with
(
exp(hθ, fi i − kF (θ)), if i ∈ F,
(9)
pF,θ (i) =
0,
otherwise,
where kF (θ) = log

P

i∈F

exp(hθ, fi i).

Proof. See Barndorff-Nielsen (1978). For self-containedness we provide
a proof in our notation in Appendix G.1.
Thus, EA is a finite union of sets EF,A that are exponential families themselves with a very similar parametrization, using the same number of parameters. The design matrix of EF,A is the submatrix AF of A consisting of
the columns indexed by F . However, for any proper facial set F 6= I, the
parametrization θ 7→ pF,θ is never identifiable since all columns of AF lie on
a supporting hyperplane defining F and thus ÃF never has full rank.
Although the parameters θ on EA and the parameters θ on EF,A play
similar roles, they are very different in the following sense: if θ(s) is a sequence
(s)
of parameters with pθ(s) → pF,θ for some θ, then, in general, lims→∞ θj 6= θj
for all j ∈ J.
Theorem 2.4 (Barndorff-Nielsen (1978)). For any vector of observed
counts n, there is a unique maximum p∗ of ˜l in EA . This maximum p∗
satisfies: (1) Ap∗ = Nt , where t = An, (2) supp(p∗ ) = Ft , (3) p∗ ∈ EFt ,A .
Proof. See Barndorff-Nielsen (1978). For self-containedness we provide
a proof in our notation in Appendix G.2.
The maximum p∗ in Theorem 2.4 is called the extended maximum likelihood estimate (EMLE). By Theorem 2.4, when Ft is known, the EMLE p∗
can be computed by computing the MLE on EFt ,A . If the MLE θ∗ exists,
then p∗ = pθ∗ .
2.4. Decomposable models. Computing Ft or finding an approximation
is easier when the simplicial complex ∆ of the model is decomposable. We
need the following definitions.
Let V 0 ⊂ V . The restriction or induced subcomplex to V 0 is ∆|V 0 = {S ∈
∆ | S ⊆ V 0 }. The subcomplex ∆|V 0 is complete, if ∆|V 0 contains V 0 (and
thus all subsets of V 0 ). In this case, we also say that V 0 is complete in ∆.

APPROXIMATING FACES OF MARGINAL POLYTOPES

11

A subset S ⊂ V is a separator of ∆ if there exist V1 , V2 ⊂ V with V1 ∩
V2 = S, ∆ = ∆|V1 ∪ ∆|V2 and V1 6= S 6= V2 . A simplicial complex that
has a complete separator is called reducible. By extension, we also call the
hierarchical model reducible.
A hierarchical model is decomposable if its generating set is a union ∆ =
∆1 ∪ ∆2 ∪ . . . ∆r of induced sub-complexes ∆i = ∆|Vi in such a way that
1. each ∆i is a complete simplex: ∆i = {S ⊆ Vi }; and
2. (∆1 ∪ · · · ∪ ∆i ) ∩ ∆i+1 is a complete simplex.
In other words, ∆ arises by iteratively gluing simplices along complete subsimplices.
Lemma 2.5 below states that, if ∆ is reducible, then any facial set for ∆
is the intersection of the preimage of facial sets for its components. It is a
simple reformulation of Lemma 8 in (Eriksson et al., 2006).
Lemma 2.5.

Let ∆ be reducible into two components ∆|V1 and ∆|V2 .

1. If F ⊆ I is facial with respect to ∆, then πV1 (F ) and πV2 (F ) are facial
with respect to ∆|V1 and ∆|V2 .
2. Conversely, if F1 ⊆ IV1 and F2 ⊆ IV2 are facial with respect to ∆|V1
and ∆|V2 , then πV−1
(F1 ) ∩ πV−1
(F2 ) is facial with respect to ∆.
1
2
Thus, for any T ⊆ I, let T1 = πV1 (T ) and T2 = πV2 (T ). Then
F∆ (T ) = πV−1
(F∆|V1 (T1 )) ∩ πV−1
(F∆|V2 (T2 )).
1
2
Lemma 2.5 generalizes to more than one separator and thus to more than
two components. It becomes particularly simple when these components are
complete: in that case, F∆|V1 (T1 ) = T1 . Taking the preimage we obtain

πV−1
(π
(T
))
=
i ∈ I : ∃i0 ∈ T such that πV1 (i) = πV1 (i0 ) ⊇ T.
V
1
1
Thus, for a decomposable complex ∆ = ∆1 ∪ ∆2 ∪ · · · ∪ ∆r , we have
(10)

F∆ (T ) = π1−1 (π1 (T )) ∩ π2−1 (π2 (T )) ∩ · · · ∩ πr−1 (πr (T ))

for any T ⊆ I, where πi = πV (∆i ) .
3. Approximations of facial sets. We consider a hierarchical model
with simplicial complex ∆ and marginal polytope P∆ . In this section, we
develop the details of our methodology to obtain inner and outer approximations to the facial set Ft of the data vector t.

12

WANG & RAUH & MASSAM

3.1. Inner approximations. To obtain an inner approximation, our strategy is to find a separator S of ∆ and to complete it. To be precise, we
augment ∆ by adding all subsets of S. Thus, we obtain a simplicial complex ∆S = ∆ ∪ {M : M ⊆ S} in which S is a complete separator. We
can apply Lemma 2.5 to find the facial set F∆S (I+ ), and this will be our
inner approximation of Ft , because F∆S (I+ ) ⊆ F∆ (I+ ) = Ft according to
Lemma 2.1.
An even simpler approximation is obtained by not only completing the
separator itself, but also the two parts V1 , V2 separated by S: the simplicial
complex ∆V1 ,V2 := {M : M ⊆ V1 } ∪ {M : M ⊆ V2 } is decomposable and
contains ∆. Its facial sets can be computed from (10).
In general, the approximation obtained from a single separator (or, in general, a single super-complex) is not good; that is, Ft = F∆ (I+ ) tends to be
much larger than F∆S (I+ ) or F∆V1 ,V2 (I+ ). Thus we need to combine information from several separators. For example, given two separators S, S 0 ⊆ V ,
we find a chain of approximations
G00 := I+ ,
G1 := F∆S (G00 ), G01 := F∆S0 (G1 ),

G2 := F∆S (G01 ), G02 := F∆S0 (G2 ),
..
.

that satisfy
I+ ⊆ G1 ⊆ G01 ⊆ G2 ⊆ · · · ⊆ Ft ,
where all inclusions except the last one are due to the definition of F∆S (T )
or F∆S0 (T ) as the smallest facial sets containing T in ∆S or ∆S 0 . The last
inclusion is a consequence of Lemma 2.1 since both ∆S and ∆S 0 contain ∆.
This chain of approximations has to stabilize; that is, after a certain number of iterations,
will not improve any more. The limit
S the approximations
S
FS,S 0 (I + ) := i Gi = i G0i can be characterized as the smallest subset
of I that contains I + and is facial both with respect to ∆S and ∆S 0 . The
same iteration can be done replacing ∆S and ∆S 0 by ∆V1 ,V2 and ∆V10 ,V20 . Applying in turn F∆V1 ,V2 and F∆V 0 ,V 0 gives another approximation F̃S,S 0 (I + ),
1

2

namely the smallest subset of I that contains I + and is facial both with
respect to ∆V1 ,V2 and ∆V10 ,V20 . This latter approximation will be used in Section 5.1. Clearly, F̃S,S 0 (I + ) is a worse approximation than FS,S 0 (I + ), since
F̃S,S 0 (I + ) ⊆ FS,S 0 (I + ) ⊆ Ft , but it is easier to compute.
We use the following strategies:
1. if possible, use all separators of a graph.

APPROXIMATING FACES OF MARGINAL POLYTOPES

13

We illustrate this strategy in Section 5.2, using a graphical model associated
with the NLTCS data set.
There are two problems with this strategy: First, if S is such that either
V1 or V2 is large, then it becomes difficult to compute F∆|V1 and F∆|V2 .
Such “bad” separators always exist: namely, each node i ∈ V is separated
by its neighbours from all other nodes. In this case, V1 consists of i and its
neighbours, and V2 consists of V \{i}. For such a “bad” separator we can only
compute F∆V1 ,V2 , but not F∆S . Second, the number of separators may be
large. Thus, when computing the inner approximation, it may take a long
time until the iteration over all separators converges. A faster alternative
strategy is the following:
2. use all separators such that both V1 \ S and V2 \ S are not too small
(for example, min{|V1 \ S|, |V2 \ S|} ≥ 3).
In the case of the grids studied in Sections 5.1 and 6.2, which have a lot
of regularity, we use an adapted strategy:
3. in a grid, use the horizontal, vertical and diagonal separators.
In the case of grids, the vertical separators form a family of pairwise disjoint
separators. In Section 6 we show how to make use of such a family to study
faces of hierarchical models, even if the facial sets are so large that they
become computationally intractable.
3.2. Outer approximations. By Lemma 2.1, the facial set F∆0 (S) for a
simplicial sub-complex ∆0 ⊆ ∆ provides an outer approximation of F∆ (S).
Removing sets from ∆ decreases the dimension of the marginal polytope,
so it is often easier to compute F∆0 (S) than to compute F∆ (S). Our main
strategy is to look at induced sub-complexes.
When comparing ∆ with an induced sub-complex ∆|V 0 for some V 0 ⊂ V ,
we have to be precise about whether we consider ∆|V 0 as a complex on V
or on V 0 . When we consider it on V , then its design matrix A has columns
fi indexed by i ∈ I. When we consider it on V 0 , its design matrix A0 has
columns fi0 indexed by IV 0 . Because we have the same set of interactions
whether we are on V or V 0 , we have for i ∈ I and i0 ∈ IV 0 ,
(11)

fi = fi00 ⇔ i ∈ πV−10 (i0 ).

Therefore the marginal polytopes of the two models are the same since they
are the convex hull of the same set of vectors {fi , i ∈ I} = {fi00 , i0 ∈ IV 0 }.
The relationship between the facial sets on V and V 0 is as follows:

14

WANG & RAUH & MASSAM

Lemma 3.1.

Let V 0 ⊆ V . For K ⊂ I, we have
0
F∆|V 0 (K) = πV−10 (F∆|
(πV 0 (K))).
0
V

0
Here, F∆|
denotes the facial set when ∆|V 0 is considered as a simplicial
V0
complex on V 0 , and F∆|V 0 denotes the facial set when ∆|V 0 is considered as
a simplicial complex on V .

Proof. For K ⊂ I, the two sets A = {ai , i ∈ K} and B = {bi0 , i0 ∈
πV 0 (K)} are identical and therefore the smallest faces of the marginal polytopes for ∆V 0 on V or V 0 containing A and B respectively are the same.
0
By definition of F∆
(πV 0 (K)), the smallest face containing B is defined
V0
0
0
by {bi0 , i ∈ F∆ 0 (πV 0 (K))}. By definition of F∆V 0 (K), the smallest face
V
containing A is {ai , i ∈ F∆V 0 (K)}. Also by (11), we have that {ai , i ∈
0
0
πV−10 (F∆
(πV 0 (K)))} = {bi0 , i0 ∈ F∆
(πV 0 (K))}. Therefore F∆V 0 (K) =
V0
V0
−1
0
(πV 0 (K))).
πV 0 (F∆
0
V

In general, F∆|V 0 (I+ ) is not a good approximation of F∆ (I+ ). This approximation can be improved by considering several subsets of V . To be
precise, if V1 , . . . , Vr ⊆ V , then F∆ (I+ ) ⊆ F∆|V (I+ ) for i = 1, . . . , r, and
i
T
thus F∆ (I+ ) ⊆ ri=1 F∆|V (I+ ) =: FV1 ,...,Vr ;∆ (I+ ). In contrast to the case of
i
the inner approximation, no repeated iteration is needed. Thus, the outer
approximation is faster to compute.
The question is now how to choose the subsets Vi . Clearly, the subsets
Vi should cover V , and, more precisely, they should cover ∆, in the sense
that for any D ∈ ∆ there should be one Vi with D ∈ ∆|Vi . The larger the
sets Vi , the better the approximation becomes, but the more difficult it is
to compute FV1 ,...,Vr ;∆ (I+ ). One generic strategy is the following:
1. use all subsets of V of fixed cardinality k plus all facets D ∈ ∆
with |D| ≥ k.
This choice of subsets indeed covers ∆. The parameter k should be chosen
as large as possible such that computing FV1 ,...,Vr ;∆ (I+ ) is still feasible. Note
that computing F∆|D (I+ ) for D ∈ ∆ is trivial, since P∆|D is a simplex.
Another natural strategy, due to Massam and Wang (2015), is the following:
2. for fixed k, use balls Bk (v) = {w : d(v, w) ≤ k} around the nodes
v ∈ V , where d(·, ·) denotes the edge distance in the graph.
Our general philosophy is that the subsets Vi should be large enough to
preserve some of the structure of ∆. For example, for the grid graphs, we

APPROXIMATING FACES OF MARGINAL POLYTOPES

15

suggest to use 3 × 3 sub-grids. These graphs have two nice properties: First,
they already have the appearance of a small grid. Second, for any vertex
v ∈ V , there is a 3 × 3 sub-grid that contains v and all neighbours of v. We
will compare two different strategies:
3. for a grid, use all 3 × 3 sub-grids;
4. cover a grid by 3 × 3 sub-grids.
In Section 6.2 we compare these two methods, and we observe that, in the
example of the 5 × 10 grid, it suffices to only look at a covering.
In general, it is not enough to look at induced sub-complexes, unless ∆ has
a complete separator (see Section 2.4). However, the approximation tends
to be good and gives the correct facial set in many cases.
3.3. Comparing the two approximations. Suppose that we have computed two approximations F1 , F2 of Ft such that F1 ⊆ Ft ⊆ F2 . If we are in
the lucky case that F1 = F2 , then we know that Ft = F1 = F2 . In general,
the cardinality of F2 \ F1 indicates the quality of our approximations.
The sets F1 , F2 and Ft can also be compared by the ranks of the matrices
ÃF1 , ÃF2 and ÃFt obtained from Ã by keeping only the columns indexed
by F1 , F2 and Ft , respectively. Clearly, rank ÃF1 ≤ rank ÃFt ≤ rank ÃF2 .
Note that rank ÃF2 − 1 equals the dimension of the corresponding face F2
of P, and rank ÃFt − 1 equals the dimension of Ft . Although F1 does not
necessarily correspond to a face of P, we can bound the codimension of Ft
in F2 by
dim F2 − dim Ft ≤ rank ÃF2 − rank ÃF1 .
In particular, if rank ÃF2 = rank ÃF1 , then we know that Ft = F2 . In this
case, our approximations give us a precise answer, even if F1 6= F2 and the
lower approximation F1 is not tight.
4. Parameter Estimation when the MLE does not exist.
4.1. Computing the extended MLE. If the MLE θ∗ exists, then it can
be computed by finding the unique maximum of the log-likelihood function l(θ) given in (7). As mentioned before, l(θ) is concave (or even strictly
concave, if the parameters θ are identifiable), and thus the maximum is, at
least in principle, easy to find (in practice, for larger models, it may be difficult to evaluate the function k(θ), which involves a sum over I; but we will
not discuss this problem here). In general, the maximum cannot be found
symbolically, but there are efficient numerical algorithms to maximize concave functions. Any reasonable hill-climbing algorithm should be capable

16

WANG & RAUH & MASSAM

of finding the MLE. An example of an algorithm commonly used is iterative proportional fitting (IPF), which can be thought of as an algorithm of
Gauss-Seidel type (Csiszár and Shields, 2004).
When the MLE does not exist but the facial set F = Ft of the data is
known, then it is straight forward to compute the extended MLE p∗ . In
this case, we know that p∗ lies in EF,A . To find p∗ , we need to optimize the
log-likelihood ˜l over EF,A = {pF,θ : θ ∈ Rh }. Plugging the parametrization
pF,θ (see Theorem 2.3) into ˜l tells us that we need to optimize the restricted
log-likelihood function
Y
X
(12)
lF (θ) = log(
pF,θ (i)n(i) ) =
θj tj − N kF (θ).
i∈I+

j∈J

This problem is of a similar type as the problem to maximize l in the case
that the MLE exists, and the same algorithms as discussed above can be
used. The problem here is slightly easier, since F is smaller than I. However,
as stated above, the parametrization θ 7→ pF,θ is never identifiable. Of course,
this problem is easy to solve by selecting a set of independent parameters
among the θj . However, depending on the choice of the independent subset,
the values of the parameters change, and in particular, it is meaningless to
compare the values of the parameters θj with parameter values of any other
distribution in EA or in the closure EA .
Before explaining how to find better parameters on EF,A , let us discuss
what happens if the facial set Ft of the data is not known. As mentioned before, whether or not the MLE exists, the log-likelihood function l(θ) is always
strictly concave (assuming that the parametrization is identifiable). When
the MLE does not exist, then the maximum is not at a finite value θ∗ , but lies
“at infinity.” Still, as observed by Geyer (2009, Section 3.15), any reasonable
numerical “hill-climbing” algorithm that tries to maximize the likelihood will
tend towards the right direction. Such a numeric algorithms generates a sequence of parameter values θ(1) , θ(2) , θ(3) , . . . with increasing log-likelihood
values l(θ(1) ) ≤ l(θ(2) ) ≤ . . . . Since l(θ) is concave, our optimization problem is numerically easy (at least in theory), and for any reasonable such
algorithm, the limit lims→∞ l(θ(s) ) will equal supθ l(θ) = maxp∈EA ˜l(p). The
algorithm will stop when the difference l(θ(s+1) ) − l(θ(s) ) becomes negligeably small. The output, θ(s) , then gives a good approximation of the EMLE,
in the sense that p∗ and pθ(s) are close to each other. For many applications,
such as in machine learning, where it is more important to have good values
of the parameters instead of trying to model the “true underlying distribution,” or when doing a likelihood test, where the value of the likelihood is
more important than the parameter values, this may be good enough.

APPROXIMATING FACES OF MARGINAL POLYTOPES

17

However, in this numerical optimization, some of the parameters θj will
tend to ±∞, which may lead to numerical problems. For example, it may
happen that one parameter goes to +∞ and a second parameter to −∞
in such a way that their sum remains finite (see Appendix B for a simple
such example with two variables). This implies that a difference between
two large numbers has to computed, which is numerically unstable. Also,
it is not clear, which parameters tend to infinity numerically. In fact, this
may depend on the chosen algorithm; i.e. different algorithms may yield
approximations of the EMLE that are qualitatively different in the sense
that different parameters diverge.
To avoid such problems, we propose a change of coordinates that allows
us to control which parameters diverge, at least in the case where we know
the facial set Ft . If Ft is unknown, but if we know approximations F1 ⊆ Ft ⊆
F2 , we can use this knowledge to identify some parameters that definitely
remain finite, while some parameters definitely diverge. We cannot control
the behaviour of the remaining parameters, but, as will be illustrated in
Section 5.2, the MLE obtained with the model on F2 lies closer to the EMLE
than the MLE on the original model. The more information we have about
the facial set Ft , the better we can control the above mentioned pathologies.
4.2. An identifiable parametrization. We have seen that when we use
the parametrization θ 7→ pFt ,θ of EA,Ft in the case where Ft 6= I, we have to
expect the following (interrelated) issues:
1. The parametrization is not identifiable, i.e. there are parameters θ, θ0
with pFt ,θ = pFt ,θ0 .
2. While the parametrization θ 7→ pFt ,θ of EFt ,A looks similar to the
parametrization θ 7→ pθ of EA , the values of the parameters in both
parametrizations are not related to each other.
3. When pθ(s) → pFt ,θ as s → ∞ for some parameter values θ(s) , θ, then
some of the parameter values θ(s) diverge to ±∞. When computing
probabilities, there may be linear combinations of these diverging parameters that remain finite.
Next we show that if Ft is known, then, with a convenient choice of L, the
parameters µL (introduced in Section 2.1) solve 1 and 2 and improve 3.
Afterwards, we discuss what can be done if Ft is not known. We briefly
discuss the general solution towards 3 in Appendix C. In any case, the choice
of the parameters will depend on the facial set Ft : it is not possible to define
a single parametrization that works for all facial sets simultaneously.
Suppose that Ft is known. We choose a zero element in I+ and consider

18

WANG & RAUH & MASSAM

the parameters µi as in Section 2. Recall that
µi (θ) = hθ, fi − f0 i = log p(i)/p(0), i ∈ I.
As mentioned in Section 2, the parameters µi are not independent, and we
need to choose an independent subset L. We will do this in two steps.
1. Choose a maximal subset Lt of Ft such that the parameters µi , i ∈ Lt
are independent.
2. Then extend Lt to a maximal subset L ⊆ I such that the parameters µi , i ∈ L are independent by adding elements i ∈ I \ Ft .
It follows from Theorem 2.4 that the following holds.
1. The subset µi , i ∈ Lt , of the parameters µL gives an identifiable
parametrization of EFt ,A .
2. Let µ∗i , i ∈ Lt , be the parameter values that maximize lFt (and thus
give the EMLE). When the likelihood l(µ) in (8) is maximized numerically on I, then in successive iterations of the maximization, the
(s)
estimates µi are such that
(
µ∗i ,
i ∈ Lt ,
(s)
µi →
−∞, otherwise.
In particular, no parameter tends to +∞.
The last property ensures a consistency of the parameters µi on EA and
on EFt ,A . This is important in those cases where the parameters have an interpretation and where it is of interest to know the value of some parameters,
if it is well-defined. For example, in hierarchical models, the parameters correspond to “interactions” of the random variables, and it may be of interest
to know, which of these interactions are important. Thus, it is of interest to
know the size of the corresponding parameter. Usually, it is not the parameter µi , but the original parameters θi that have an interpretation. But when
we understand the parameters µi , we can also tell which of the paramters θi
or which combinations of the parameters θi have finite well-defined values
and can be computed, and which parameters diverge:
Lemma 4.1. Suppose that θ(s) , s ∈ N, are parameter values such that
pθ(s) → p∗ as s → ∞. For any i ∈ Lt , the linear combination
(s)

µi

= hθ(s) , fi i
(s)

has a well-defined finite limit as s → ∞. Any linear combination of the θi
that has a well-defined finite limit (that is, a limit that is independent of the
(s)
choice of the sequence θ(s) ) is a linear-combination of the µi with i ∈ Lt .

19

APPROXIMATING FACES OF MARGINAL POLYTOPES

Proof. The first statement follows from
(s)

µi

= log pθ(s) (i)/pθ(s) (0) → log p∗ (i)/p∗ (0).

For the second statement, note that any linear combination of the θ is also
a linear combination of the µ, since the linear
P map θ 7→ µ(θ) is invertible.
We now show that if a linear combination i ai µi involves some µj with
j∈
/ Lt , then there exist sequences µ(s) , µ0(s) of parameters with
X
X
(s)
0(s)
lim pµ(s) = lim pµ0(s) and lim
ai µi 6= lim
a i µi .
s→∞

s→∞

s→∞

s→∞

i

i

So suppose that µ(s) is a sequence of parameters such that lims→∞ pµ(s)
P
(s)
exists and such that lims→∞ i ai µi is finite. Define
( (s)
µj + 1, if i=j,
0(s)
µi =
(s)
µi ,
otherwise.
An easy computation shows that
lim pµ0(s) = lim pµ(s)

s→∞

s→∞

and

lim

s→∞

X

0(s)

ai µi

= lim

i

s→∞

X

(s)

ai µi + aj .

i

Suppose now that we do not know Ft , but that instead we have approximations F1 , F2 that satisfy
I+ ⊆ F1 ⊆ Ft ⊆ F2 ⊆ I.
In this case, we proceed as follows to obtain an independent subset L among
the parameters µi :
1. Choose a maximal subset L1 of F1 such that the parameters µi , i ∈ L1
are independent.
2. Then extend L1 to a maximal subset L2 ⊆ F2 such that the parameters µi , i ∈ L2 are independent by adding elements i ∈ F2 \ F1 .
3. Finally, extend L2 to a maximal subset L ⊆ I such that the parameters µi , i ∈ L are independent by adding elements i ∈ I \ F2 .
These parameters have the following properties that follow directly from
Lemma 4.1:
Lemma 4.2. Suppose that θ(s) , s ∈ N, are parameter values such that
(s)
pθ(s) → p∗ as s → ∞, and let µi = hθ(s) , fi i.

20

WANG & RAUH & MASSAM

1. For any i ∈ L1 , the linear combination
(s)

µi

= hθ, fi i

has a well-defined finite limit as s → ∞. Thus, any linear combination
(s)
of the µi with i ∈ L1 has a well-defined limit as s → ∞.
P
(s)
2. Any linear combination i ai µi that has a well-defined limit as s →
(s)
∞ is in fact a linear combination of the µi with i ∈ L2 . Thus, a
(s)
linear combination that involves at least one µj with j ∈ L \ L2 does
not have a well-defined limit.
5. Simulation study and applications to real data. In this section,
we illustrate our methodology. In 5.1, we simulate data for the graphical
model of the 4×4 grid and show how to exploit the various types of separators
in order to obtain good inner and outer approximations. We find that our
method gives very accurate result in this model of modest size. In 5.2, we
work with the NLTCS data set, a real-world data set. We compare different
inner approximations F1 and find that most of the time, F1 and F2 are
equal, and thus they are both equal to Ft . We also compute the EMLE
and compare these exact estimates to those obtained when maximizing the
likelihood functions l and lF2 . We find the results given by lF2 better than
those given by l, and extremely close to the finite components of the EMLE.
5.1. 4 × 4 grid graph. We generated random samples of varying sizes for
the graphical model of the 4 × 4 grid graph with binary variables (Fig. 1a).
For each sample, we compute inner and outer approximations F1 and F2 ,
and we compare them to the true facial set Ft , which we can obtain using
linear programming. To obtain an inner approximation, we use two strategies. Either, we iterate over all possible separators, of which there are 106
(Strategy (1) in Section 3.1), or we iterate over the 3 horizontal, 3 vertical
and 8 diagonal separators only (Strategy (3) in Section 3.1). We obtain the
same result with either strategy. Clearly, Strategy (3) is much faster. To
compute the outer approximation, we cover the 4 × 4 grid by four 3 × 3 grids
(Strategy (3) in Section 3.2).
We generate samples from the hierarchical model Pθ , where the vector of
parameters θ is drawn from a multivariate standard normal distribution (for
each sample, new parameters were drawn). The results are given in Table 1.
For each sample size, 1000 samples were obtained. Observe that the squared
length of the parameter vector θ is χ2 -distributed with 40 degrees of freedom
(since the number of parameters is 40). Thus, the expected length of θ is 40,
which is large enough to move the distribution pθ close to the boundary of the

APPROXIMATING FACES OF MARGINAL POLYTOPES

21

Table 1
Facial set approximation of 4 × 4 grid graph(hierarchical log-linear model with parameters
from standard normal distribution)
sample size
10
50
100
150

MLE does not exist

F1 = F t

F2 = Ft

100.0%
89.5%
71.0%
52.0%

97.7%
100.0%
100.0%
100.0%

100.0%
100.0%
100.0%
100.0%

model. Indeed, we observed that when the MLE does not exist, the length
of the numerical estimate of the MLE vector is of the order of magnitude
of 40 (see also the next example in Section 5.2). In all samples that we
generated, Ft = F2 , and F1 = F2 in the vast majority of cases. Thus, for
this graph of relatively modest size, our approximations are very good. We
present additional simulation results in Appendix D.
5.2. NLTCS data set. To illustrate how approximate knowledge of the
facial set allows us to say which parameters can be estimated (as explained in
Section 4), we study the NLTCS data set, which consists of 21 574 observations on 16 binary variables, called ADL1, . . . , ADL6, IADL1, . . . , IADL10.
The reader is referred to Dobra and Lenkoski (2011) for a detailed description of the data set. To associate a hierarchical model to this data, we rely
on the results of Dobra and Lenkoski (2011) who use a Bayesian approach
to estimate the posterior inclusion probabilities of edges. We construct a
graph by saying that (x, y) is an edge if and only if the posterior inclusion
probability of (x, y) is at least 0.40: we obtain Figure 1b. Then we take the
corresponding clique complex of this graph so that our hierarchical model
is a graphical model. There are 314 parameters in this model, including up
to 6-way interactions. In total, the graph has 40 separators.
In order to compare the estimates obtained with or without worrying
about the existence of the MLE and with or without an approximation
to Ft , we maximize the loglikelihood given in terms of µ, rather than θ, as
in (8). First we ignore the fact that the MLE might not exist and numerically optimize the likelihood directly: we call this estimate µ̂MLE . Second, we
find Ft and compute the EMLE with parameters denoted µ̂EMLE . Third, we
obtain an inner and outer approximation to Ft and consider the resulting in0
0
formation on likelihood maximization. We call the resulting estimate µ̂F 2 \F1 .
All estimates are computed using the Matlab function minFunc (Schmidt,
2005).
To compute µ̂EMLE , we first compute the inner approximation F1 that

22

WANG & RAUH & MASSAM

I5
I7

I9
I10

I8
I2
I6

A5
I4
A6

A2

I1

I3

A1

(a)
A3

A4

(b)
Fig 1. (a) 4 × 4 grid graph. (b) Graphical model for NLTCS data set. The label “ An”
abbreviates ADLn, “ In” abbreviates IADLn.

makes use of all the separators in the graph (Strategy 1 in Section 3.1). We
also compute an outer approximation F2 from all 16
5 = 4368 size five local
models and the cliques of size six (Strategy 1 in Section 3.2). We obtain
F1 = F2 and thus deduce that Ft = F1 = F2 . We find |Ft | = 49 536, and so
|Ftc | = 216 − 49 536 = 16 000. Therefore, 16 000 cell probabilities are zero in
the EMLE, a precise estimate of those cells that we could not obtain from
the MLE. We obtain the EMLE by maximizing the loglikelihood function
lFt as in (12). Since rank(ÃFt ) = 303, the dimension of Ft is 302, and there
are only 302 parameters in lF . This information is most important when
testing the present model against another model M2 of smaller dimension.
As pointed out by Geyer (2009) and Fienberg and Rinaldo (2012), the test
statistic, chi-square or loglikelihood, has to be compared to the chi-square
distribution with 302−d2 degrees of freedom, not 314−d2 . Of course, for M2
also, d2 is the dimension of the smallest face of the corresponding polytope
containing the data.
To show how to use the inner and outer approximations when Ft is not
known, we construct coarser inner and outer approximations to Ft , respectively denoted F10 and F20 , and use them to compute another approximation
0
0
µ̂F2 \F1 to the EMLE. To compute F10 , we just use 10 random separators. We

APPROXIMATING FACES OF MARGINAL POLYTOPES

23

find |F10 | = 36 954 and dim F01 = rank ÃF10 − 1 = 300. To compute the outer
approximation F20 , we consider the 4368 local size-five induced models and
select among them the 1000 with the facial sets of smallest cardinality, which
we glue together. We find |F20 | = 50 688 and dim F02 = rank ÃF20 − 1 = 310.
Thus, we know that at least |I \ F20 | = 216 − 50 688 = 14 848 cell probabilities
vanish in the EMLE. Since we pretend not to know Ft , we replace lFt by
X
X
(13)
lF20 (µ) =
µi n(i) − N
exp(µi ).
i∈I+

i∈F20

We know that µi is estimable for i ∈ F10 , that µi goes to negative infinity for
0
i ∈ F2c , and we cannot say anything for µi with i ∈ F20 \ F10 .
As explained in Section 4.2, the components of µ are not functionally
independent. We choose L1 ⊆ F10 , L2 ⊆ F20 and L ⊆ I as in Section 4.2
(we note that the zero cell belongs to I+ ). Then any µi , i ∈ F20 , can be
written as a linear combination of µL2 = (µi , i ∈ L2 ), and we can write
µi = hbi , µL2 i for an appropriate vector bi . Thus, lF20 (µ) only depends on
µL2 = (µi , i ∈ L2 ), and (13) can be rewritten as
X
X
(14)
lF20 (µL2 ) =
hbi , µL2 in(i) − N
exphbi , µL2 i.
i∈I+

i∈F20

Of course, the maximum of lF20 does not exist but, as for the maximization of l, the computer still gives us a numerical approximation, µ̂L2 , and
thus also a numerical estimate µ̂i = hbi , µ̂L2 i, i ∈ F20 . In total, there are
|L2 | = rank(ÃF20 ) − 1 = 310 independent parameters in the loglikelihood
function (14). Among them, there are |L1 | = rank(ÃF10 ) − 1 = 300 estimable
parameters µi , i ∈ L1 . We cannot say anything about the 10 parameters
indexed by L2 \ L1 . If we know Ft , we can identify two more estimable
parameters.
0
0
Table 2 gives the three estimates µ̂MLE
, µ̂EMLE
and µ̂Fi 2 \F1 . The naive
i
i
estimator log nn0i is also listed. Estimates are listed for 19 arbitrarily chosen parameters among the 310 possible ones. The first column of the table
indicates whether the index i belongs to F10 , Ft or F20 . The second column
lists the particular parameters considered. By Theorem 2.4, the only parameters µi with a finite estimate are those for i ∈ Ft . This is illustrated in the
µ̂EMLE column of Table 2, with finite values for µ̂EMLE
, i ∈ Ft (green and
i
EMLE
pink rows), and infinite values for µ̂i
, i ∈ I \ Ft (yellow and blue rows).
When we choose the coarser approximations, F10 and F20 , to Ft , we compute
0
0
0
0
the estimate µ̂F 2 \F1 using (14). The components µ̂Fi 2 \F1 indexed by i ∈ Ft
are excellent. They are finite and close to the corresponding components

24

WANG & RAUH & MASSAM

Table 2
Parameter estimates from 3 methods compared with the relative frequency in the NLTCS
data. Here, each i = (i1 , . .P
. , i16 ) ∈ I = {0, 1}16 is represented by the natural
j−1
∈ {0, . . . , 216 − 1}.
number 16
j=1 ij 2
naive estimate
log ni /n0
−1.2472
−1.7644
−2.3958
−2.5429
−2.8813
−∞
−∞
−∞
−∞
−∞
−∞
−∞
−∞
−∞
−∞
−∞
−∞
−∞
−∞

Parameter
i ∈ F10
µ512
µ65536
µ16
µ528
µ2048
i ∈ Ft \ F10 µ32960
µ34881
i ∈ F20 \ Ft µ36864
µ36880
µ388
µ32769
µ385
µ449
µ32785
µ389
i ∈ I \ F20
µ256
µ320
µ257
µ321

maximum likelihood estimates
0
0
2 \F1
µ̂MLE
µ̂EMLE
µ̂F
i
i
i
−1.2482
−1.2482
−1.2482
−1.7976
−1.7975
−1.7975
−2.3844
−2.3846
−2.3846
−2.6504
−2.6504
−2.6504
−2.7246
−2.7243
−2.7243
−13.8205 −13.8207 −13.8205
−14.3693 −14.3693 −14.3692
−30.8729
−∞
−34.9805
−39.6536
−∞
−45.2229
−28.9090
−∞
−29.4525
−32.3799
−∞
−36.9537
−37.1365
−∞
−35.9399
−38.9673
−∞
−44.9405
−40.1221
−∞
−45.8318
−43.7297
−∞
−40.0158
−35.5482
−∞
−∞
−42.5454
−∞
−∞
−52.9224
−∞
−∞
−60.2208
−∞
−∞

of µ̂EMLE . This can be seen by verifying numerically that the square length
of the projection on Ft of the difference between µ̂MLE and µ̂EMLE is greater
0
0
than that between µ̂F 2 \F1 and µ̂EMLE . Indeed, we have
0

0

LE 2
LE
LE 2
kµ̂FFt2 \F1 − µ̂EM
k ≈ 6.49 < kµ̂M
− µ̂EM
k ≈ 8.52 .
Ft
Ft
Ft
0

0

The components µ̂Fi 2 \F1 indexed by i ∈ F20 \ Ft are finite while the corresponding components of µ̂EM LE are infinite but they are better than those
of µ̂M LE : numerically, we have
X
X
2
0
0 2
µ̂Fi 2 \F1 ≈ 5184 >
µ̂MLE
≈ 4752 .
i
i∈F02 \Ft
0

i∈F02 \Ft
0

The estimates µ̂Fi 2 \F1 , i ∈ F20 \ Ft are better than the corresponding µ̂MLE
i
since they are larger and thus “closer to the truth”. For i ∈ I \ F20 , cor0
0
responding to the blue rows of Table 2, the components µ̂Fi 2 \F1 are better
0
0
than the µMLE
since, by construction, the µ̂Fi 2 \F1 are infinite.
i
For reference, we list the estimates of the top five cell counts obtained
using our method and compare them with those obtained by other methods
in Dobra and Lenkoski (2011) in Appendix E.

APPROXIMATING FACES OF MARGINAL POLYTOPES

25

6. Computing faces for large complexes. If our statistical model
contains many variables and is not reducible, the problem of determining
Ft quickly becomes infeasible. Not only does the marginal polytope become
very complicated, but also the size of the objects that one has to store or
compute grows exponentially. Consider for example a 10 × 10 grid of binary
random variables. This hierarchical model has 280 parameters, and the total
sample space has cardinality |I| = 2100 ≈ 1.27 × 1030 . If Ft is close to I, we
cannot even list the elements of Ft , which consists of approximately 1030
elements. Therefore, we take a local approach and look for separators.
If the simplicial complex ∆ contains a complete separator separating V
into V1 and V2 , we can identify a facial set F implicitly without listing it
explicitly. We only need the two projections FV1 = πV1 (F ) and FV2 = πV2 (F ).
Since F = πV−1
(FV1 ) ∩ πV−1
(FV2 ) (by Lemma 2.5), these two projections
1
2
identify F , and they allow us to do most of the operations that we would
want to do with F . For example, for any i ∈ I, we can check whether i ∈ F by
checking whether πV1 (i) ∈ FV1 and πV2 (i) ∈ FV2 , and we can check whether
F ⊇ I by checking whether FV1 ⊇ IV1 and FV2 ⊇ IV2 . In particular, we can
check whether the MLE exists by looking only at the two subsets V1 and V2 .
Similar ideas apply if ∆ contains a separator that is not complete. Suppose
that S separates V1 from V2 in ∆. We want to use F2 := F∆|V1 (I+ )∩F∆|V2 (I+ )
as an outer approximation and F1 := F∆S (I+ ) as an inner approximation
to Ft . Due to the problems mentioned above, we do not directly compute F1
and F2 , but we compute their projections on V1 and V2 . Instead of F2 , we
compute the facial set F2,V1 := F∆|V1 (πV1 (I+ )) of the V1 -marginal πV1 (I+ )
with respect to ∆|V1 . Similarly we compute F2,V2 := F∆|V2 (πV2 (I+ )). Instead
of F1 , we compute F1,V1 := F∆S |V1 (πV1 (I+ )) and F1,V2 := F∆S |V2 (πV2 (I+ )).
Then we could recover F1 and F2 from the equations
F2 = πV−1
(F2,V1 ) ∩ πV−1
(F2,V2 )
1
2

and

F1 = πV−1
(F1,V1 ) ∩ πV−1
(F1,V2 ).
1
2

For any x ∈ I, we can check whether x ∈ F1 by checking whether πV1 (x) ∈
F1,V1 and πV2 (x) ∈ F1,V2 . More importantly, we can check whether F1 = F2
by checking whether F1,V1 = F2,V1 and F1,V2 = F2,V2 . This idea can be
applied iteratively when either ∆|V1 or ∆|V2 has a separator.
The next two subsections illustrate these ideas. In Section 6.1, we consider
a graph with no particular regularity pattern on 50 nodes and identify convenient separators. In Section 6.2, we consider a grid graph and work with
two families of “parallel” separators that can be used to iteratively improve
the inner approximation.
6.1. US Senate Voting Records Data. We consider the voting record of
all 100 US senators on 309 bills from January 1 to November 19 2015. Similar

26

WANG & RAUH & MASSAM
Reid
Gillibrand

Coons
McCaskill

Markey
Murray

Bennet

Schumer

Brown

Warren
Cantwell

Merkley

Heitkamp

Casey

Menéndez

Wyden

Feinstein

Donnelly

Murphy

Cardin

Blumenthal

Manchin
Heinrich
Peters

Booker

King

Shaheen
Udall
Reed
Carper

Kaine

Whitehouse

Boxer

Mikulski

Warner

Sanders

Durbin

Alexander

Capito

Tester

Schatz

Baldwin

Hatch
Cochran

Leahy

Rounds
Enzi

Thune

Barrasso

Fischer

Klobuchar

McConnell

Cornyn

Tillis

Cassidy

Stabenow

Hirono

Nelson

Sessions

McCain
Roberts

Franken
Ayotte

Wicker
Collins

Shelby

Cruz
Burr

Sasse

Inhofe

Cotton
Isakson

Risch
Graham

Crapo

Scott

Heller
Moran

Perdue

Daines

Gardner

Lankford

Kirk
Lee

Corker

Boozman

Flake

Toomey

Hoeven
Johnson

Ernst
Coats Grassley

Murkowski
Sullivan
Paul
Blunt

Portman
Rubio
Vitter

Fig 2. The graph for the US senate voting records data. Golden nodes are independent
senators, blue nodes are Democratic, and red nodes are Republican.

data for the years 2004–2006 was analyzed by Banerjee et al. (2008). The
votes are recorded as “yea,” “nay” or “not voting.” We transformed the “not
voting” into “nay” and consequently have a 100-dimensional binary data
set. To fit a hierarchical model to this data set, we use the `1 -regularized
logistic regression method proposed by Ravikumar et al. (2010) to identify
the neighbours of each variable andpconstruct an Ising model. We set the
regularization parameter to λ = 32 log p/n ≈ 0.35. The underlying graph
of the Ising model is given in Figure 2. This figure should not be interpreted
as the graph of a graphical model. Rather, the edges in the graph indicate
where the two-way interactions lie. There are 277 parameters in this model
(the number of vertices plus the number of edges). The graph consists of

27

APPROXIMATING FACES OF MARGINAL POLYTOPES
ID

Senator

ID

Senator

22
23
26

Nelson
Reed
Schumer

37
41
47

Cardin
Markey
Udall

ID

Senator

52 Murphy
53 Hirono
56 Gillibrand
Table 3
Numbering of some senators

ID

Senator

61
87

Whitehouse
Warren

two separate large connected components and 14 independent nodes.
There are 309 sample points, and |I+ | = 278. We want to characterize
the face Ft of the data on the marginal polytope. The graph in Figure 2
has many complete separators, and it decomposes as a union of several
small irreducible simplicial subgraphs and two large irreducible subgraphs,
one in each of the large connected components, as shown in Figure 3. By
Lemma 2.5, we can restrict attention to these irreducible subgraphs. For the
small irreducible subgraphs, one easily verifies that the data does not lie on
a proper face of their corresponding marginal polytopes. We are left with
the two large irreducible prime components in Figure 3.
Cornyn

Tillis

Hatch
Rounds

Gillibrand

Coons

Cochran

Markey
Murray

Cassidy

Barrasso

Fischer
Inhofe

Roberts
Sasse

Wicker

Moran

Warren
Cantwell
Feinstein

Merkley

Murphy

Menéndez

Wyden

Risch
Crapo

Schumer

Brown

Cardin

Blumenthal

Heinrich
Peters

Booker

Shaheen

Daines

Udall
Reed

Boozman

Hirono

Nelson
Whitehouse

Hoeven

Ernst

Blunt

(a)

Mikulski
Durbin

(b)

Fig 3. The simplicial complexes after cutting off the small prime components: (a) the
Republican party prime component ∆r . (b) the Democratic party prime component ∆d .
The yellow and pink nodes are the two separator sets we used to compute the facial set.

The Democratic party simplicial complex ∆d consists of 26 variables, and
the model induced from ∆d contains 77 parameters, which is too large to
use linear programming to compute the face of P∆d containing the vector td .
Therefore, we look for separators in order to obtain good inner and outer
approximations. Figure 3b indicates (in yellow and pink) two separators that
separate ∆d into three simplicial complexes denoted, from left to right, by
∆α , ∆β and ∆γ that are small enough for the linear programming method.
∆α has 9 nodes. The corresponding vector tα lies in the relative interior
of P∆α .
∆β has 13 nodes, and tβ lies on a facet Ftβ of P∆β . To simplify the

28

WANG & RAUH & MASSAM

notation, we denote the 100 senators not by their name but by an integer
between 1 and 100. We only need to identify a few and their numbers are
given in Table 3. The inequality of Ftβ is
(15)

t87 − t56,87 ≥ 0,

where t87 denotes the marginal count of senator Warren voting “yea” and
t56,87 denotes the marginal counts of both senators Gillibrand and Warren
voting “yea.”
∆γ has 11 nodes.The data vector tγ lies on the facet of P∆γ with inequality
(16)

t23 − t23,53 ≥ 0.

The intersection of the two facets (15) and (16) gives the outer aproximation
F2,d to Ftd .
To get an inner approximation, we complete each separator, i.e. the yellow
vertices are completed and the pink vertices are completed in Figure 3b. Denote the three simplicial complexes with complete separators as ∆α̃ , ∆β̃ , ∆γ̃
respectively. Then ∆d˜ = ∆α̃ ∪ ∆β̃ ∪ ∆γ̃ is a simplicial complex with two
complete separators. The smallest face Ftd˜ of the marginal polytope P∆d˜
containing the data vector td˜ is our inner approximation. The models of
∆α̃ , ∆β̃ , ∆γ̃ and ∆d˜ include main effects, two-, three- and four-way interactions. The dimension of the model induced by ∆d˜ is 91 (completing the two
separators adds 14 parameters to the original model). We apply the linear
programming method to P∆α̃ , P∆β̃ and P∆γ̃ .
The dimension of the model of ∆α̃ is 27, and Ftα̃ is a facet with equation
(17)

hg1 , tα̃ i = t41 − t22,41 − t41,70 + t22,41,70 = 0.

It follows that {g1 } is a basis of the kernel of AtFα̃ .
The dimension of the model for ∆β̃ is 48. The face Ftβ̃ has codimension 5,
with defining equations

hg2 , tβ̃ i = t87 − t56,87 = 0






hg3 , tβ̃ i = t47,52,61 + t37,52 − t37,52,61 − t37,47,52 = 0
.
(18)
hg4 , tβ̃ i = t37,47,52,61 − t47,52,61 = 0



hg5 , tβ̃ i = t37,52 + t26 − t26,52 − t26,37 = 0



hg6 , tβ̃ i = t41 − t22,41 − t41,70 + t22,41,70 = 0
Again, {g2 , g3 , g4 , g5 , g6 } is a basis of the kernel of AFβ̃ .

APPROXIMATING FACES OF MARGINAL POLYTOPES

29

The dimension of the model for ∆γ̃ is 38. The face Ftγ̃ has codimension 3.
It is defined by the equations


hg7 , tγ̃ i = t47,52,61 + t37,52 − t37,52,61 − t37,47,52 = 0
(19)
.
hg8 , tγ̃ i = t37,47,52,61 − t47,52,61 = 0


hg9 , tγ̃ i = t23 − t23,53 = 0 .
Again, {g7 , g8 , g9 } is a basis of the kernel of AFγ̃ .
From Lemma 2.5, Ftd˜ = Fα̃ ∩ Fβ̃ ∩ Fγ̃ , and the equations for Ftd˜ are


hg10 , td˜i = t41 − t22,41 − t41,70 + t22,41,70 = 0





hg20 , td˜i = t87 − t56,87 = 0



hg 0 , t i = t
47,52,61 + t37,52 − t37,52,61 − t37,47,52 = 0
3
d˜
(20)
0

hg4 , td˜i = t37,47,52,61 − t47,52,61 = 0





hg50 , td˜i = t37,52 + t26 − t26,52 − t26,37 = 0



hg 0 , t i = t − t
23
23,53 = 0,
9
d˜
where the vectors g10 , . . . , g90 are the vectors g1 , . . . , g9 extended to R91 by
adding zeros on the corresponding complementary coordinates (since g10 =
g60 , g30 = g70 , g40 = g80 , only six of the nine equations are needed). Thus,
F1,d := Ftd˜, defined by (20), is a strict subset of the face F2,d defined by (15)
and (16). Next, we refine our argument and show that indeed Ftd = F2,d .
From what we know, it follows that the orthogonal complement of the
subspace generated by Ftd˜ is
n
o
G = g 0 ∈ R91 g 0 = k1 g10 + k2 g20 + k3 g30 + k4 g40 + k5 g50 + k9 g90 .
To describe Ftd , we note that each defining equation of Ftd is of the form
hg, td i = 0, where g is orthogonal to Ftd . For any such g, let g 0 be its
extension to a vector in R91 by adding zero components. Then g 0 ⊥ Ftd˜,
which implies that g 0 ∈ G. Therefore, we can find g by finding all vectors
g 0 ∈ G that vanish on all added components. This yields a system of linear
equations in k1 , . . . , k5 , k9 . We claim that all solutions must satisfy k1 = k3 =
k4 = k5 = 0. Indeed, the coefficient of any triple or quadruple interaction
must vanish (since these do not belong to the original Ising model), which
implies k1 = k3 = k4 = 0, and also the coefficient of t37,52 must vanish,
which implies k5 = 0. On the other hand, the vectors g20 and g90 only contain
interactions that are already present in ∆, and so the coefficients k2 and k9
are free. Thus the equations for Ftd are
(
hg2 , tβ̃ i = t87 − t56,87 = 0,
(21)
hg9 , tγ̃ i = t23 − t23,53 = 0.

30

WANG & RAUH & MASSAM

This is the same as the outer approximation F2,d .
The Republican simplicial complex ∆r consists of 20 variables, and the
model induced from ∆r contains 46 parameters, which is also too large to
directly compute Ftr . The yellow nodes in Figure 3a separate ∆r into two
simplicial complexes denoted (from left to right) by ∆a and ∆b . To compute
the inner approximation, we complete the yellow separator and obtain two
new simplicial complexes ∆ã and ∆b̃ . With linear programming, we find
that the corresponding vectors tã and tb̃ lie in the relative interior of the
polytopes P∆ã and P∆b̃ , respectively. Therefore, F1 = P∆r , from which we
conclude that the corresponding vector tr lies in the relative interior of P∆r .
Thus, the face Ft is now determined: it is characterized by the equalities (21). What insight is there in the knowledge (i) of the non-existence of
the MLE and (ii) of the exact face Ft ? While we have given general remarks
in the introduction, let us illustrate here how knowledge about Ft points to
some issues with the statistical analysis that would possibly be overlooked
if Ft was not known.
First, knowing Ft , and its defining inequalities, for one model also gives
information about other models. It follows from (21) that the MLE does not
exist for any hierarchical model that includes one of the edges (23, 53) or
(56, 87) (to see this, note that the inequality (16) defines a proper face for
any model containing the edge (23, 53), since the corresponding sufficient
statistics vector satisfies the equality in (16)). Thus, if one wants to find
a smaller model, within the realm of hierarchical models, for which the
MLE exists, both edges have to be dropped. However, from the data, here,
evidence for both edges is quite strong, and thus the edges should not be
dropped.
Second, let us consider the computation of the EMLE. As we know Ft ,
instead of running an MLE computation for a model with 277 parameters
and 2100 outcomes, we are left with an MLE computation for a model with
9
277 − 2 = 275 parameters and |Ft | = 16
· 2100 outcomes, those without the
configurations (X23 , X53 ) = (1, 0) or (X56 , X87 ) = (1, 0) that all have counts
zero. These numbers are still too large for a direct computation, even when
taking into account that the EMLE can be computed by restricting to each
of the irreducible components. So, we turn to an approximate method and
compute the maximum composite likelihood estimate. The maximum composite likelihood estimate can be obtained by combining estimates from the
local conditional likelihoods derived from the distribution of each variable
given its neighbours: see for example Liu and Ihler (2012) or Massam and
Wang (2017). Thus the reliability of the maximum composite likelihood estimates depends upon the existence of the maximum in each of the local

APPROXIMATING FACES OF MARGINAL POLYTOPES

31

conditional likelihood. These local conditional likelihoods are derived from
the global model built on the entire cell set I and, certainly in practice,
without worrying about the existence of the global MLE. Let us consider,
for example, the likelihood obtained from the conditional distribution of X23
given its neighbours X19,53,61,78 . For convenience, let 19, 23, 53, 61, 78 be denoted as a, b, c, d, e. This likelihood is the product over all configurations of
iacde in the data set of conditional binomial distributions for the variable Xb
and can be written as
Y
n(i ,i ,i ,i ,i )
(22)
p Xb = ib Xacde = (ia , ic , id , ie ) a b c d e ,
ia ,ib ,ic ,id ,ie

where n(ia , ib , ic , id , ie ) denotes the corresponding marginal cell count. It is
easy to show that the MLE of each p(Xb = 1 | ia , ic , id , ie ) is the empirical estimate n(ia , ib = 1, ic , id , ie )/n(ia , ic , id , ie ). In the data set, n(ia , ib =
1, ic = 0, id , ie ) = 0 for all ia , id , ie ∈ {0, 1}. Thus,
p̂(Xb = 1 | ia = 1, ic = 0, id = 1, ie = 1) =

exp(θ̂b + θ̂ab + θ̂bd + θ̂be )
= 0,
1 + exp(θˆb + θ̂ab + θ̂bd + θ̂be )

so that θ̂b + θ̂ab + θ̂bd + θ̂be = −∞, and the MLE of at least some of these parameters, which are the corresponding parameters of the global model, does
not exist. Now the composite maximum likelihood estimate is obtained by
averaging the estimates obtained from various local conditional likelihoods.
From the b-local conditional model, we also obtain p̂(Xb = 1 | ia = 1, ic =
1, id = 1, ie = 0) = 1/2 and p̂(Xb = 1 | ia = 0, ic = 1, id = 0, ie = 1) = 4/5,
which yield the linear combinations
(23)

θ̂b + θ̂bc + θ̂ab + θ̂bd = 0, θ̂b + θ̂bc + θ̂be = 1.4.

The remarks above are verified numerically. Let θ = (θb , θab , θbc , θbd , θbe ),
and denote by llocal (θ) the local conditional likelihood. Starting at θ0 =
(0, 0, 0, 0, 0) and optimizing (22) in terms of θ in Matlab, we obtain llocal (θ̂) =
3.88830675 for the maximum θ̂ = θ̂1 ≈ (−62.3, −16.8, 35.2, 43.9, 28.5).
If we change the starting point to θ0 = (100, 100, 100, 100, 100), we obtain
θ̂2 = (−162.2, −26.1, 91.2, 97.1, 72.4) and llocal (θ̂2 ) = 3.88830648. Clearly,
the values for θ̂ are unreliable since the MLE in the local conditional model
does not exist. However, both θ̂1 and θ̂2 satisfy equations (23). One can, of
course, obtain estimates of θbc , θab , θbd , θbe from the local conditional models
centered at c, a, d and e respectively but these estimates will not have been
obtained through the method of composite likelihood and it remains to study
their properties.

32

WANG & RAUH & MASSAM
1

6

11 16 21 26 31 36 41 46

5

10 15 20 25 30 35 40 45 50

Fig 4. 5 × 10 grid graph, the red and blue nodes are the set of separators we use to
compute F1 , they are used iteratively to get a better lower approximation

This example shows that our method scales well and makes it possible to
obtain the face Ft for very large examples. It also illustrates how knowing
the face Ft gives us precious information on the reliability of the maximum
composite likelihood estimate.
6.2. The 5 × 10-grid. Let ∆ be the simplicial complex of the 5 × 10 grid
graph. We exploit the regularity of this graph and make use of the vertical
separators in the grid to obtain inner and outer approximations of the facial
sets. The graph has 50 nodes, which is too many to directly compute a facial
set or even to store it. However, the 5 × 10 grid has 8 vertical separators
marked in red and blue in Figure 4, and we can use these to approximate Ft .
Since facial sets for 5 × 3 grids can be computed reasonably fast (3 to 4
seconds on a laptop with 2.50 GHz processor and 12 GB memory), we only
use three of these vertical separators at a time, say the blue separators
S2 = {11, . . . , 15}, S4 = {21, . . . , 25}, S6 = {31, . . . , 35}, S8 = {41, . . . , 45}.
These separate the vertex sets
V1 = {1, . . . , 15}, V3 = {11, . . . , 25}, V5 = {21, . . . , 35},
V7 = {31, . . . , 45}, V9 = {41, . . . , 50}.
Adding the blue separators to ∆ gives a simplicial complex
[
∆S2 ;S4 ;S6 ;S8 := ∆
{F : F ⊆ Sj }
j=2,4,6,8

with five irreducible components supported on the vertex sets V1 , V3 , V5 , V7
and V9 (Figure 6). To compute a facial set with respect to ∆S2 ;S4 ;S6 ;S8 ,

33

APPROXIMATING FACES OF MARGINAL POLYTOPES
1

6

11 11 16 21 21 26 31 31 36 41 41 46

5

10 15 15 20 25 25 30 35 35 40 45 45 50
Fig 5. Five induced sub-grids

according to Lemma 2.5 applied four times, we need to compute
G1,V1 := F∆S2 |V1 (πV1 (I+ )),
G1,V5 := F∆S4 ;S6 |V5 (πV5 (I+ )),

G1,V3 := F∆S2 ;S4 |V3 (πV3 (I+ )),
G1,V7 := F∆S6 ;S8 |V7 (πV7 (I+ )),

G1,V9 := F∆S8 |V9 (πV9 (I+ )).
T
(G1,Vi ) is equal to F∆S2 ;S4 ;S6 ;S8 (I+ ), and thus an inner
Then G1 := i πV−1
i
approximation of Ft . As stated before, we do not need to compute G1 explicitly, but we represent it by means of the G1,Vi .
We can improve the approximations by also considering the red separators
S1 = {6, . . . , 10}, S3 = {16, . . . , 20}, S5 = {26, . . . , 30}, S7 = {36, . . . , 40},
that separate
V0 = {1, . . . , 10}, V2 = {6, . . . , 20}, V4 = {16, . . . , 30},
V6 = {26, . . . , 40}, V8 = {36, . . . , 50}.
(2)

As explained in Section 3.1, we want to compute G1 := F∆S1 ;S3 ;S5 ;S7 (G1 ).
(2)

Again, instead of computing G1 directly, we need only compute the much
(2)
(2)
(2)
(2)
(2)
(2)
smaller sets G1,V0 := πV0 (G1 ), G1,V2 := πV2 (G1 ), . . . , G1,V8 := πV8 (G1 ).
(2)

(2)

(2)

So the question is: is it possible to compute G1,V0 , G1,V2 , . . . , G1,V8 from
G1,V1 , G1,V3 , . . . , G1,V9 , without computing G1 in between?
It turns out that this is indeed possible: By Lemma 2.5, all we need
(2)
to compute G1,Vi is G1,Vj := πVj (G1 ), j = i − 1, i + 1. For i = 0, since
V0 ⊂ V1 , we can compute G1,V0 from πV1 (G1 ) = G1,V1 . For i = 2, 4, 6, 8, since
Vi ⊂ Vi−1 ∪ Vi+1 , we can compute G1,Vi from πVi−1 ∪Vi+1 (G1 ), which itself
can be obtained by “gluing” πVi−1 (G1 ) = G1,Vi−1 and πVi+1 (G1 ) = G1,Vi+1 :




Vi−1 ∪Vi+1 −1
Vi−1 ∪Vi+1 −1
πVi−1 ∪Vi+1 (G1 ) = πVi−1
(G1,Vi−1 ) ∩ πVi+1
(G1,Vi+1 ),

34

WANG & RAUH & MASSAM
0

where πVV 00 for V 00 ⊆ V 0 denotes the marginalization map from IV 0 to IV 00
 0 −1
and where πVV 00
denotes the lifting from IV 00 to IV 0 .
(2)

As explained in Section 3.1, we have to iterate this procedure: From G1
(3)
we want to compute G1 := F∆S2 ;S4 ;S6 ;S8 (G01 ) or, more precisely, we want
(3)

(3)

to compute G1,Vi = πVi (G1 ) for i = 1, 3, . . . , 9. Again, we do this with(2)

out looking at G1 directly by just using the information available through
(3)
(k)
(k)
the G1,Vi . Iterating this procedure, we obtain a sequence of sets G1,Vi , G1,Vj
(with odd i and even j), which stabilizes after a finite number of steps. Let
[ (k)
F1,Vi :=
G1,Vi ,
T
(F1,Vi ). Again, we do
Our best inner approximation is then F1 = 9i=0 πV−1
i
not compute F1 explicitly, but we represent it in terms of the F1,Vi . The
process is visualized in Figure 7.
1

6 11 16 21 26 31 36 41 46

5 10 15 20 25 30 35 40 45 50
(a)

1

6 11 11 16 21 21 26 31 31 36 41 41 46

5 10 15 15 20 25 25 30 35 35 40 45 45 50
(b)
Fig 6. (a) The 5 × 10-grid with the blue separators completed. (b) The five irreducible
subcomplexes after completing the separators.

Let us now consider the outer approximation F2 . We adapt Strategy 3 of
Section 3.2 and cover the graph with 5 × 3 grid subgraphs, since the facial
sets for such graphs can easily be computed. These subgrids are supported

APPROXIMATING FACES OF MARGINAL POLYTOPES

35

Table 4
facial set approximation of 5 × 10 grid graph
sample size
50
100
150
200
300
400
500
1000

F2 6= I

F1 = F 2

100.0%
100.0%
99.9%
99.6%
96.4%
92.9%
84.8%
44.7%

94.3%
82.5%
76.5%
81.2%
87.7%
91.5%
93.9%
99.9%

on the same vertex subsets Vi , i = 1, . . . , 8 as used when computing F1 .
This makes it possible to compare F1 and F2 . For i = 1, 3, . . . , T
8 we compute
(F2,Vi ).
F2,Vi = F∆|V (πVi (I+ )). The outer approximation is then F2 = i πV−1
i
i
Again, we don’t compute F2 explicitly, but we only store F2,Vi in a computer
as a representation of F2 . To compare the two approximations F1 and F2 , we
need only compare their projections F1,Vi and F2,Vi pairwise, i = 1, . . . , 8.
We generated random data of varying sample size. For each fixed sample
size, we generated 100 data samples. The simulation results are shown in
Table 4. For each simulated sample, we compute the sets F1,Vi and F2,Vi as
described above. When computing F1,Vi , we found that 2 iterations actually
suffice. Then we checked whether F2 is a proper subset of I (second column),
and we checked whether F1 = F2 (third column). Both for small and large
sample sizes, we found that the F1 = F2 quite often.
We also investigated what happens when the outer approximation is not
computed using all 3×5-subgrids, but only a cover of four 3×5-subgrids and
one 2 × 5-subgrid (as in Figure 5). In all simulations, this easier approximation gave the same result. The same is not true for the inner approximation:
when using just one of the two families of parallel separators we obtain an
inner approximation that is much too small.
7. Conclusion. As mentioned before, previous work had made it possible to identify Ft for hierarchical models with up to 16 variables. In this
paper, we offer a methodology to approximate, and sometimes, completely
identify the facial set Ft for high-dimensional models. To find an inner and
an outer approximation to Ft , first, we divided the original problem into
subproblems of dimension at most 16 for which we could use linear programming and, second, we combined the facial sets of the subproblems and
related then to Ft . Identifying the subproblems and relating the facial sets
to Ft is numerically easy and the corresponding software can be obtained

36

WANG & RAUH & MASSAM

data I+ on 5 × 10 grid

marginalize

1 = π (I )
I+
V1 +

3 = π (I )
I+
V3 +

LP

LP

G1,V1

marginalize

G1,V0

LP

9 = π (I )
I+
V9 +

LP

G1,V5

LP

G1,V7

G1,V9

marginalize
and glue

marginalize
and glue

marginalize
and glue

marginalize
and glue

G1,V2

G1,V4

G1,V6

G1,V8

LP

G01,V2
..
.

7 = π (I )
I+
V7 +

LP

G1,V3

LP

G01,V0

5 := π (I )
I+
V5 +

LP

G01,V4
..
.

LP

G01,V6
..
.

Fig 7. Flow chart

G01,V8
..
.

..
.

APPROXIMATING FACES OF MARGINAL POLYTOPES

37

upon request, from the authors.
It has long been established that determining the existence of the MLE
is essential to correct statistical inference. In our paper, we have emphasized the problem of parameter estimation and shown how working with
the likelihood lF2 yields much better estimates of the parameter that when
working with l. When testing one model versus another, the correct degrees
of freedom for the asymptotic distribution of the test statistic is the difference between the dimensions of the facial sets for the two models being
compared and not the difference between the dimensions of the two models.
If we only know approximations F1 and F2 , we can use their dimensions to
approximate the correct degrees of freedom.
In high dimensions, when the likelihood functions has so many terms
that the (E)MLE cannot be computed, a popular approach is to compute
the maximum composite likelihood estimate. We have shown through an
example that, when the global MLE does not exist, the local MLE for some
of the same parameters might not exist either. So, combining the values of
the MLE of local likelihoods without being aware that the data lies on a
face of the marginal polytope, one might also obtain misleading estimates
of the parameters through composite likelihood.
We have not addressed the question of how to obtain reliable confidence
intervals for the parameters by exploiting the properties of the inner and
outer approximations to Ft . This subject clearly deserves attention and
should be the subject of further work.
APPENDIX A: PARAMETRIZING HIERARCHICAL MODELS
In this section, we recall the usual parametrization of hierarchical models, see, for example, Letac and Massam (2012). The starting point is the
parametrization (2) of the hierarchical model, which we repeat here for convenience:
X
(A24)
log p(i) =
θD (iD )
D∈∆

This parametrization is not identifiable; that is, for any joint distribution p
from the hierarchical model there are different choices for the functions θD
that satisfy (A24). One way to make the parameters unique is to choose
a special element within each set Iv , which we denote by 0. The choice of
0 is arbitrary, and a different choice of 0 leads to a simple affine change of
parameters. With this choice, the functions θD become unique if one requires
θD (iD ) = 0 whenever iv = 0 for some v ∈ D.

38

WANG & RAUH & MASSAM

A parametrization in terms of real numbers is obtained using the following
definitions: for i ∈ I, we write
S(i) = {v ∈ V ; iv 6= 0},

J = {j ∈ I \ {0}, S(j) ∈ ∆}.

For any j ∈ J, let
θj = θD (iD ) for the unique i ∈ I with S(i) = D, iD = jD .
To simplify the notation, we write j / i whenever S(j) ⊆ S(i) and jS(j) =
iS(j) . It is convenient to introduce the vectors
fi =

X

ej , i ∈ I

j∈J:j/i

where ej , j ∈ J are the unit vectors in RJ . Moreover, let A be the J × I
matrix with columnsfi , i ∈ I, and let Ã be the (1 + |J|) × I matrix with
columns equal to f1i , i ∈ I. Then (2) can be rewritten in the following
equivalent forms
(A25)

log pθ (i) =

X

θj − k(θ) = hθ, fi i − k(θ) = At θ − k(θ) = Ãt θ̃,

j∈J:j/i

where θ̃ = (θ0 , θ) as a column vector and
X
 X

(A26)
− θ0 = k(θ) = log
exp
θj
i∈I

j∈J:j/i

acts as a normalization constant. If n = (n(i), i ∈ I) denotes the I-dimensional column vector of cell counts, then


N
(A27)
Ãn =
and An = t,
t
P
where N =
i∈I n(i) is the total cell counts and t is the column vector
of sufficient statistic with components equal to P
the jS(j) -marginal counts
n(jS(j) ), i.e. t = (tj , j ∈ J) where tj = n(jS(j) ) = i|iS(j) =jS(j) n(i), j ∈ J.
P
It follows from (A27) that Nt = i∈I n(i)
N fi . Therefore, t belongs to the
convex polytope with extreme points fi , i ∈ I. This polytope is the marginal
polytope of the hierarchical model, denoted by P∆ .

39

APPROXIMATING FACES OF MARGINAL POLYTOPES

Example A.1. Let V = {a, b, c}, Ia = {0, 1} = Ib = Ic and ∆ =
{a, b, c, ab, bc}. Then
I = (000, 100, 010, 110, 001, 101, 011, 111),
J = {(100), (010), (001), (110), (011)},
f000

f001

f010

f011

f100

f101

f110

f111

z}|{ z}|{ z}|{ z}|{ z}|{ z}|{ z}|{ z}|{
1
1
1
1
1
1
1
1
 0
1
0
1
0
1
0
1 


 0

0
1
1
0
0
1
1

Ã = 
 0
0
0
0
1
1
1
1 


 0
0
0
1
0
0
0
1 
0
0
0
0
0
0
1
1

θ000
θ100
θ010
θ001
θ110
θ011

APPENDIX B: EXAMPLE: TWO BINARY RANDOM VARIABLES
Consider two binary random variables, and let ∆ = {∅, {1}, {2}, {1, 2}}.
The hierarchical model E∆ is the saturated model ; that is, it contains all
possible probability distributions with full support. Then
f00

f01

f10

f11

z}|{ z}|{ z}|{ z}|{
1
1
1
1
 0
1
0
1 

Ã = 
 0
0
1
1 
0
0
0
1

θ00
θ01
θ10
θ11

The marginal polytope is a 3-simplex (a tetrahedron) with facets
F00 : 1 − t01 − t10 + t11 ≥ 0,
F10 : t10 − t11 ≥ 0,

F01 : t01 − t11 ≥ 0,
F11 : t11 ≥ 0.

Each of the corresponding facets contains three columns of A. In fact, the
facet Fi in the above list does not contain the column fi of A.
The EMLE of the saturated model is just the empirical distribution; that
is, p∗ = N1 n. Suppose that t lies on the facet F00 (i.e. n = (0, n01 , n10 , n11 )
with n(01), n(10), n(11) > 0). If pθ(s) → p∗ , then pθ(s) (00) → 0, while all

40

WANG & RAUH & MASSAM

other probabilities converge to a non-zero value. It follows that
(s)

θ00 = log pθ(s) (00) → −∞,
p (s) (01)
(s)
→ +∞,
θ01 = log θ
pθ(s) (00)
p (s) (10)
(s)
θ10 = log θ
→ +∞,
pθ(s) (00)
p (s) (11)pθ(s) (00)
(s)
θ11 = log θ
→ −∞.
pθ(s) (01)pθ(s) (10)
(s)

(s)

On the other hand, θ01 + θ00 = log pθ(s) (01) converges to a finite value, as
(s)
(s)
(s)
(s)
do θ10 + θ00 = log pθ(s) (10) and θ11 + θ01 = log pθ(s) (11)/pθ(s) (10).
Proceeding similarly for the other facets, one can show for the limits
(s)
θij := lims→∞ θij :

F00
F01
F10
F11

θ00

θ01

θ10

θ11

finite parameter combinations:

−∞
finite
finite
finite

+∞
−∞
finite
finite

+∞
finite
−∞
finite

−∞
+∞
+∞
−∞

θ01 + θ00 , θ10 + θ00 , θ11 + θ01
(s)
(s)
(s)
(s)
θ00 , θ10 , θ01 + θ11
(s)
(s)
(s)
(s)
θ00 , θ01 , θ10 + θ11
(s)
(s)
(s)
θ00 , θ10 , θ01

(s)

(s)

(s)

(s)

(s)

(s)

Each line of the last column contains three combinations of the parameters
(s)
θi that converge to a finite value. Any other parameter combination that
converges is a linear combination of these three. This can be seen by using
the coordinates µi introduced in Section 4.2. For example, on the facet F01 ,
consider the parameters
µ10 = log p(10)/p(00) = θ10 ,

µ11 = log p(11)/p(00) = θ10 + θ01 + θ11 ,

µ01 = log p(01)/p(00) = θ01 .
Then µ10 and µ11 are identifiable parameters on EF01 , and µ01 diverges close
to F01 . By Lemma 4.1, the linear combinations that are well-defined are
µ10 = θ10 and µ11 = θ10 + (θ01 + θ11 ). The above table also lists θ00 , which
is not a linear combination of those but that is fine because it is not free.
We obtain similar results for the facets F01 and F11 . The results are
summarized in the following table:
facet
F01
F10
F11

µ01
−∞
finite
finite

µ10
finite
−∞
finite

µ11
finite
finite
−∞

41

APPROXIMATING FACES OF MARGINAL POLYTOPES

Of course, by definition of the µi s, we cannot consider the facet F00 where
n(00) = 0. To study F00 , we have to choose another zero cell and redefine
the parameters µi .
The situation is more complicated for faces smaller than facets, because
sending a single parameter to plus or minus infinity can be enough to send
the distribution to a face F of higher codimension, as we will see below. The
remaining parameters then determine the position within E∆,F . Thus, in this
case there are more remaining parameters than the dimension of E∆,F .
For example, the data vector n = (n00 , 0, n10 , 0) (with n00 , n10 > 0) lies
on the face F = F01 ∩ F11 of codimension two. If pθ(s) → p∗ , then
(s)

n00
,
N

(s)

pθ(s) (01)
→ −∞,
pθ(s) (00)
n10
p (s) (10)
→ log
= log θ
.
pθ(s) (00)
n00

θ00 = log pθ(s) (00) → log
θ01 = log
(s)

θ10

p

(s)

However, the limit of θ11 = log pθ

(s) (11)pθ (s) (00)

(01)pθ(s) (10) is not determined. The only
θ (s)
(s)
(s)
constraint is that θ11 cannot go to +∞ faster than θ01 goes to −∞, since
(s)
(s)
(s)
(s)
pθ(s) = exp(θ00 + θ01 + θ10 + θ11 ) has to converge to zero.
11

With the same data vector n = (n00 , 0, n10 , 0), suppose we use a numerical algorithm to optimize the likelihood function by optimizing the parameters θj in turn. To be precise, we order the parameters θj in some way. For
simplicity, say that the parameters are θ1 , θ2 , . . . , θh . Then we let
(k+1)

θj

(k+1)

= arg max l(θ1
y∈R

(k+1)

(k)

(k)

, . . . , θj−1 , y, θj+1 , . . . , θh )

(this is called the non-linear Gauss-Seidel method ). Let us choose the ordering θ01 , θ10 , θ11 (note that θ00 = −k(θ) is not a free parameter). We start
(0)
(0)
(0)
at θ01 = θ10 = θ11 = 0. In the first step, we only look at θ01 . That is, we
want to solve
(1)

(B28) 0 =

(1)

(0)

(0)

∂
exp(θ01 ) + exp(θ01 + θ10 + θ11 )
l(θ) = −
(1)
(0)
(1)
(0)
(0)
∂θ01
1 + exp(θ ) + exp(θ ) + exp(θ + θ + θ )
01

10

01

=−
(1)

1

10

11
(1)
2 exp(θ01 )
.
(1)
+ 2 exp(θ01 )

This derivative is negative for any finite value of θ01 , and thus the critical
equation has no finite solution. If we try to solve this equation numerically,

42

WANG & RAUH & MASSAM
(1)

we will find that θ01 will be a large negative number. Next, we look at θ10 .
We fix the other variables and try to solve
(1)

0=

(1)

(1)

(0)

∂
n10
exp(θ10 ) + exp(θ01 + θ10 + θ11 )
l(θ) =
−
(1)
(1)
(1)
(1)
(0)
∂θ10
N
1 + exp(θ01 ) + exp(θ10 ) + exp(θ01 + θ10 + θ11 )
(1)

≈

n10
exp(θ10 )
,
−
(1)
N
1 + exp(θ10 )

(1)

where we have used that θ01 is a large negative number. This equation
always has a unique solution
(1)

θ10 ≈ log

n10
.
N − n10

Finally, we look at θ11 . We have to solve
(1)

0=

(1)

(1)

∂
exp(θ01 + θ10 + θ11 )
≈ 0.
l(θ) = −
(1)
(1)
(1)
(1)
(1)
∂θ11
1 + exp(θ ) + exp(θ ) + exp(θ + θ + θ )
01

10

01

10

11

Actually, this equation again has no solution, and the numerical solution
(1)
(1)
for θ11 should be close to numerical minus infinity. However, since θ01 is
already close to −∞, the equation is already approximately satisfied. Thus,
(1)
there is no need to change θ11 . In simulations, we observed that usually θ11
(1)
will be negative, but not as negative as θ01 . In theory, we would have to
iterate and now optimize θ01 again. But the values will not change much,
since the critical equations are already satisfied to a high numerical precision
after one iteration.
It is not difficult to see that the result is different if we change the order
1 will in any case be a
of the variables. If θ11 is optimized before θ01 , then θ11
large negative number.
For general data, the derivative of l(θ) with respect to θ01 (equation (B28))
takes the form
(1)

(1)

(0)

(0)

∂
t01
exp(θ01 ) + exp(θ01 + θ10 + θ11 )
l(θ) =
−
.
(1)
(0)
(1)
(0)
(0)
∂θ01
N
1 + exp(θ01 ) + exp(θ10 ) + exp(θ01 + θ10 + θ11 )
(1)

Setting this derivative to zero and solving for θ01 leads to a linear equation
(1)
in θ01 with symbolic solution
(1)

θ01 = log

(0)

1 + exp(θ10 )
(0)

(0)

1 + exp(θ10 + θ11 ) 1

t01
N
− tN01

.

APPROXIMATING FACES OF MARGINAL POLYTOPES

43

In fact, for any hierarchical model, the likelihood equation is linear in any
single parameter θj , as long as all other parameters are kept fixed (more
generally this is true when the design matrix A is a 0-1-matrix). Instead of
optimizing the likelihood numerically with respect to one parameter, it is
possible to use these symbolic solutions. This leads to the Iterative Proportional Fitting Procedure (IPFP). In our example, the IPFP would lead to
a division by zero right in the first step, indicating that the MLE does not
exist (unfortunately, IPFP does not always fail that quickly when the MLE
does not exist).
APPENDIX C: PARAMETRIZATIONS ADAPTED TO FACIAL SETS
Let us briefly discuss how to remedy problems 1. (identifiability), 2. (relation between parameters on E and EFt ) and 3. (cancellation of infinities
in linear combinations of diverging parameters) from the beginning of Section 4.2. The idea to remedy 1. and 2. is to define parameters µi , i ∈ L,
of EA such that a subset Lt ⊆ L of the parameters parametrizes EFt ,A in a
consistent way. Denote by Aµ = (aµj,i , j ∈ L, i ∈ I) the design matrix of EA
corresponding to the new parameters µ. Then the necessary conditions are:
(∗) Let AµLt ,Ft := (aµj,i , j ∈ Lt , i ∈ Ft ) be the submatrix of Aµ with rows
indexed by Lt and columns indexed by Lt , and denote by ÃµLt ,Ft the
same matrix with an additional row of ones. The rank of ÃµLt ,Ft is equal
to |Lt | + 1, the number of its rows (and thus, AµLt ,Ft has rank |Lt |).
(∗∗) aµj,i = 0 for all i ∈ Ft and j ∈ L \ Lt .
In fact, (∗∗) implies that AµLt ,Ft is the design matrix of EA,Ft , since the parameters µi with i ∈
/ Lt do not play a role in the parametrization µ 7→ pFt ,µ .
Moreover, (∗) implies that the parametrization µ 7→ pFt ,µ is identifiable. In
this sense, we have remedied problem 1.
Since ÃµLt ,Ft has full row rank, it has a right inverse matrix C̃, such that
ÃµLt ,Ft C̃ = I|Lt |+1 equals the identity matrix of size |Lt | + 1. Recall that
log pFt ,µ (i) = hµt , fiµ i − kF (µ),
log pµ (i) = hµ̃t , fiµ i − k(µ),
for any parameter vector µ and all i ∈ Ft . Since fiµ are the columns of Aµ
and since the components of fiµ corresponding to L \ Lt vanish by (∗∗), we
may apply the matrix C obtained from C̃ by dropping the row corresponding
to kF or k and obtain
(C29)

(log pµ )C = µLt

and

(log pFt ,µ )C = µL .

44

WANG & RAUH & MASSAM

When pµ(s) is a sequence in EA with limit pµ in EFt ,A , then (C29) shows that
(s)

→ µi for i ∈ Lt . In this sense, we have remedied problem 2.
Finally, we solve problem 3. Suppose that we have chosen parameters µL
as in Section 4.2, and let AµL be the design matrix with respect to these
parameters. Then (AµL )j,i = 0 if i ∈ Ft and j ∈
/ Lt . Moreover, for j ∈ Lt ,
the jth column of AµL has a single non-vanishing entry (equal to one) at
position j. Suppose that Ft corresponds to a face Ft of codimension c. Then
there are c facets of P whose intersection is Ft . Thus, following the notation
introduced in Remark 2.2, there exist c inequalities
µi

(C30)

hg̃1 , x̃i ≥ 0,

...,

hg̃c , x̃i ≥ 0

that together define Ft . In this case, the vectors g̃1 , . . . , g̃c are linearly independent and satisfy hg̃j , f˜i i = 0 (thus, they are a basis of the kernel
of (ÃµFLt )t ). It follows that the kth component of gj , denoted by gj,k , vanishes if k ∈ Lt ; that is, the inequalities (C30) do not involve the variables
corresponding to Lt . Let G be the square matrix, indexed by L \ Lt with
entries gj,k , j, k ∈ L \ Lt . Then the square matrix


1 0
G̃ =
0 G
is invertible. We claim that the parameters λ = G̃−1 µL are what we are
looking for.
The design matrix with respect to the parameters λ is Aλ = G̃AµL . For
any j ∈
/ Lt ,
Aλj,i = 0,

if i ∈ Ft ,

and

Aλj,i = hg̃j , f˜i i ≥ 0,

if i ∈
/ Ft .

This implies the following properties:
1. If all parameters λj with j ∈
/ Lt are sent to −∞, then pλ tends towards
a limit distribution with support Ft .
2. The coefficient of λj in any log-probability is non-negative, so there is
no cancellation of ±∞.
So far, we only used the fact that the vectors g̃j define valid inequalities
for the face Ft . Suppose that we choose g̃j in such a way that each inequality
hg̃j , x̃i ≥ 0 defines a facet. The intersection of less than c facets is a face that
strictly contains Ft . This implies that for each j, there exists ij ∈ I \ Ft such
that fij satisfies
hg̃j , f˜ij i > 0,

and

hg̃j 0 , f˜ij i = 0 for all j 0 6= j,

APPROXIMATING FACES OF MARGINAL POLYTOPES

45

and so
Aλj,ij > 0,

and

Aλj0 , ij = 0 for all j 0 6= j.

Hence:
(s)

3. If λj

are sequences of parameters such that pλ(s) tends towards a
(s)

limit distribution with support Ft , then λj → −∞ for all j ∈
/ Lt .
It is not difficult to see that, conversely, any parametrization that satisfies
these three properties comes from facets defining the face Ft .
APPENDIX D: UNIFORM SAMPLING FOR THE 4 × 4-GRID
This section enhances the example in Section 5.1. In a second experiment,
we generated random samples from the uniform distribution, that is from the
probability distribution Pθ in the hierarchical model where all parameters
θj , j ∈ J, are set to zero. For each sample size, 1000 samples were obtained.
The results are given in the following table:
sample size
10
15
20
50

MLE does not exist
98.5%
68.9%
29.0%
0.0%

F1 = Ft

F2 = Ft

96.3%
99.9%
100.0%
100.0%

100.0%
100.0%
100.0%
100.0%

As the table shows, for larger samples the probability that a random sample
lies on a proper face becomes very small. If Ft = I, then clearly Ft = F2 .
But we also found Ft = F2 for all samples with t lying on a proper face,
which shows that F2 is an excellent approximation of Ft in this model. For
the inner approximation, we observed some samples with F1 6= Ft , but they
seem to be very rare.
APPENDIX E: ESTIMATED CELL FREQUENCIES FOR THE NLTCS
DATA
The following table lists the estimates of the top five cell counts obtained
using our method and compares them with those obtained by other methods
in Dobra and Lenkoski (2011).

46

WANG & RAUH & MASSAM

Support of Cell

Observed

GoM

LC

CGGMs

MLE on facial set

∅
{10}
{1 : 16}
{5}
{5, 10}
{12}

3853
1107
660
351
303
216

3269
1010
612
331
273
202

3836.01
1111.51
646.39
360.52
285.27
220.47

3767.76
1145.86
574.76
452.75
350.24
202.12

3647.4
1046.9
604.4
336
257.59
239.24

APPENDIX F: A LINEAR PROGRAMMING ALGORITHM TO
COMPUTE FACIAL SETS
Let A be the design matrix, let A+ be the sub-matrix with columns indexed by the positive cells I+ , and let A0 be the sub-matrix indexed by the
empty cells (thus A = (A+ , A0 ) after reordering the columns).
Lemma F.1.
(F31)

Any solution g ∗ of the non-linear problem
(
A+ g = 0,
max kAgk0
subject to
A0 g > 0,

is perpendicular to the smallest face containing t. The corresponding facial
set is Ft = I \ supp(Ag ∗ ).
The optimization problem (F31) is highly non-linear and non-convex. It
can be solved by repeatedly solving the associated `1 -norm optimization
problem:


A+ g = 0,
(F32)
max kA0 gk1
subject to A0 g ≥ 0,


A0 g ≤ 1.
Due to the constraints, kA0 gk1 = A0 g, and so Problem (F32) is a linear
programming problem. It can be solved repeatedly until the smallest facial
set Ft is obtained. The process is as follows:
The algorithm is introduced in the supplementary material to (Fienberg
and Rinaldo, 2012), where it is also proved that it outputs the correct result.
APPENDIX G: SOME PROOFS
G.1. Proof of Theorem 2.3. Theorem 2.3 goes back to BarndorffNielsen (1978), who studies the closure of much more general exponential
families. The case of a discrete exponential family is much easier.
The theorem follows from the following lemmas:

APPROXIMATING FACES OF MARGINAL POLYTOPES

47

Algorithm 1 Computing Ft Linear Programming
Require: Design matrix A and positive cell index I+
Let A+ := A(I+ , :), A0 := A \ A+
repeat
Solve problem (F32), get the solution g ∗ and the corresponding maximum z ∗
A0 := the submatrix of A0 with those columns fi of A0 that satisfy hfi , g ∗ i = 0.
until A0 = ∅ or z ∗ 6= 0
if A0 = ∅ then
Ft := I+
end if
if z ∗ = 0 then

Ft := I+ ∪ i column i of A belongs to A0
end if

Lemma G.1.

Let p ∈ EA . Then p ∈ EA,supp(p) .

Lemma G.2.

Let p ∈ EA . Then EA,supp(p) ⊆ EA .

Lemma G.3.

Let p ∈ EA . Then supp(p) is facial.

Lemma G.4.

If F is facial, then there exists p ∈ EA with supp(p) = F .
S
Indeed, Lemma G.1 shows that EA ⊆ F EA,F , where the union is over
all support
sets F . Lemma G.2 shows the converse containment, so that
S
EA = F EA,F . It remains to see that a subset F ⊆ I is a support set if and
only if F is facial. This follows from Lemmas G.3 and G.4.
In the proofs of Lemmas G.1 to G.4, we need the following easy lemma
of which we omit the proof:
Lemma G.5.

p ∈ EA if and only if log(p) ⊥ ker A.

Proof of Lemma G.1. Let p = limk→∞ pk , where pk ∈ EA , and let
F = supp(p). Then EA,F is the exponential family EAF , where AF consists
of the columns of A indexed by F . Any v ∈ ker AF can be extended by zeros
to v 0 ∈ ker A. By Lemma G.5,
X
0 = hlog(pk ), v 0 i =
log(pk (i))v(i) → hlog(p), vi.
i∈F

Thus, log(p) ⊥ ker AF , which implies p ∈ EA,F .
Proof of Lemma G.2. Let p = limk→∞ pk , where pk ∈ EA , let F =
supp(p), and let q ∈ EA,F . Then there exists parameters θ with log(q(i)) −
log(p(i)) = hθ, fi i for all i ∈ F . For any k, there exists ck > 0 such that
qk := ck pk exp(hθ, Ai) ∈ EA . Then qk → q as k → ∞, and so q ∈ EA .

48

WANG & RAUH & MASSAM

Proof of Lemma G.3. Let p =Plimk→∞ pk , where pk ∈ EA , and let
1
F = FA (supp(p)). Then x = | supp(p)|
i∈supp(p) fi is an interior point of the
face corresponding
to
F
,
and
thus
there
exist positive coefficients λi > 0,
P
i ∈ F , with x = i∈F λi fi . The vector v = (vi , i ∈ I) defined by
 1

 | supp(p)| − λi , i ∈ supp(p),
vi = −λi ,
i ∈ F \ supp(p),


0,
i∈
/ F,
satisfies Av = x − x = 0. By Lemma G.5, log(pk ) ⊥ v for all k. In particular,
X
X
X
λi log(pk (i)) =
log(pk (i))vi →
log(p(i))vi .
i∈F \supp(p)

i∈supp(p)

i∈supp(p)

On the other hand, note that each coefficient λi for i ∈ F \ supp(p) on the
left hand side is positive, while log(pk (i)) → −∞ for i ∈
/ supp(p). This shows
that F \ supp(p) = ∅.
Proof of Lemma G.4. If F is facial, there exist g ∈ Rh and c ∈ R with
hg, fi i ≥ c for all i ∈ I and hg, fi i = c if and only if i ∈ F . Let θ(s) = −s · g.
Then
X
kF (θ(s) ) + sc = log
exp(−shg, fi i + sc) → log |F |,
i∈I

and so
log pθ(s) (i) = −shg, fi i − kF (θ(s) ) = (sc − shg, fi i) − (kF (θ(s) ) + sc)
(
− log |F |, if i ∈ F,
→
−∞,
if i ∈
/ F,
as s → ∞. Thus, pθ(s) converges to the uniform distribution on F .
G.2. Proof of Theorem 2.4. By definition, any EMLE p∗ belongs
to the closure of the model. According to Theorem 2.3, the support of p∗
is facial. If supp(p) does not contain supp(n), then the log-likelihood goes
to minus infinity, ˜l(p) = −∞, and so p does not maximize the likelihood,
Therefore, supp(p∗ ) is a facial set containing supp(n). Thus, Ft ⊆ supp(p∗ ).
By Lemma G.1, p∗ belongs to E∆,supp(p∗ ) , which is parametrized by a
vector θ, see Theorem 2.3. On E∆,supp(p∗ ) , the log-likelihood function in terms
of this parameter θ is
X
lF (θ) =
θj tj − N kF (θ).
j∈J

APPROXIMATING FACES OF MARGINAL POLYTOPES

49

The critical equations satisfied by p∗ are
Ap∗ =

t
,
N

proving the first property. Note that these equations are independent of the
parameters and the support of p∗ . We now show that any solution to these
equations is supported on the same face of P as Nt .
Let p be a probability distribution on I such that supp(p) does not contain Ft . This means that there is a linear inequality hg, ti ≥ c that is valid
on P and such that (i) hg, fi i = c for all i ∈ Ft ; and (ii) hg, fi i > c for
some i ∈ supp(p). Then
hg, Api =

X
1 X
t
hg, fi ip(i) > c =
n(i)hg, fi i = hg, i,
N
N
i

i

which implies Ap 6= Nt . This shows supp(p∗ ) ⊆ Ft and finishes the proof of
supp(p∗ ) = Ft .
We have shown the two properties, and it remains to argue that the
EMLE is unique. But this follows from the fact that supp(p∗ ) is equal to Ft ,
and lF is strictly convex (up to identifiability), such that there is a unique
distribution on E∆,Ft that maximizes the likelihood.
REFERENCES
O. Banerjee, L. El Ghaoui, and A. d’Aspremont. Model selection through sparse maximum
likelihood estimation for multivariate gaussian or binary data. The Journal of Machine
Learning Research, 9:485–516, 2008.
O. Barndorff-Nielsen. Information and Exponential Families. Wiley, Chichester, first
edition, 1978.
I. Csiszár and F. Matúš. Generalized maximum likelihood estimates for exponential families. Probability Theory and Related Fields, 141:213–246, 2008.
I. Csiszár and P. Shields. Information Theory and Statistics: A Tutorial. Foundations
and Trends in Communications and Information Theory. now Publishers, first edition,
2004.
M. Deza and M. Laurent. Geometry of Cuts and Metrics. Springer Publishing Company,
Incorporated, 1st edition, 2009. ISBN 3642042945, 9783642042942.
A. Dobra and A. Lenkoski. Copula Gaussian graphical models and their application to
modelling functional disability data. Ann. Appl. Stat., 5(2A):969–993, 06 2011.
A. Dobra, E. Erosheva, and S. Fienberg. Disclosure limitation methods based on bounds
for large contingency tables with application to disability data. In Proceedings of conference on the new frontiers of statistical data mining, pages 93–116, 2003.
N. Eriksson, S. Fienberg, A. Rinaldo, and S. Sullivant. Polyhedral conditions for the
nonexistence of the MLE for hierarchical log-linear models. J. Symbolic Comput., 41:
222–233, 2006.
S. Fienberg. The analysis of cross-classified categorical data. Springer-Verlag, New York,
second edition, 2007.

50

WANG & RAUH & MASSAM

S. Fienberg and A. Rinaldo. Three centuries of categorical data analysis: log-linear models
and maximum likelihood estimation. J. of Statistical PLanning and Inference, 137:
3430–3445, 2007.
S. Fienberg and A. Rinaldo. Maximum likelihood estimation in log-linear models. Annals
of Statistics, 40:996–1023, 2012.
E. Gawrilow and M. Joswig. polymake: a framework for analyzing convex polytopes. In Gil
Kalai and Günter M. Ziegler, editors, Polytopes — Combinatorics and Computation,
pages 43–74. Birkhäuser, 2000.
C. Geyer. Likelihood inference in exponential families and directions of recession. Electron.
J. Statist., 3:259–289, 2009.
S. Haberman. The Analysis of Frequency Data. Univ. Chicago Press, Chicago, IL, 1974.
V. Karwa and A. Slavkovic. Inference using noisy degrees. Annals of Statistics, 44:87–112,
2016.
S. Lauritzen. Graphical Models. Oxford Science Publications, first edition, 1996.
G. Letac and H. Massam. Bayes regularization and the geometry of discrete hierarchical
loglinear models. Annals of Statistics, 40:861–890, 2012.
Q. Liu and A. Ihler. Distributed parameter estimation via pseudo-likelihood. International
Conference on Machine Learning ICML, 2012.
H. Massam and N. Wang. A local approach to estimation in discrete loglinear models.
Preprint, 2015. arXiv:1504.05434.
H. Massam and N. Wang. Local conditional and marginal approach to parameter estimation in discrete graphical models. Journal of Multivariate Analysis, 164:1–21, 2017.
J. Rauh, T. Kahle, and N. Ay. Support sets of exponential families and oriented matroids.
International Journal of Approximate Reasoning, 52(5):613–626, 2011.
P. Ravikumar, M. Wainwright, and J. Lafferty. High-dimensional Ising model selection
using `1 -regularized logistic regression. Annals of Statistics, 38(3):1287–1319, 2010.
M. Schmidt. minFunc: unconstrained differentiable multivariate optimization in Matlab.
http://www.cs.ubc.ca/~schmidtm/Software/minFunc.html, 2005.
M. Vlach. Conditions for the existence of solutions of the three-dimensional planar transportation problem. Discrete Applied Mathematics, 13(1):61–78, 1986.
G. Ziegler. Lectures on Polytopes. Springer, second edition, 1998.
E-mail: wangnanw@yorku.ca

E-mail: jrauh@mis.mpg.de
E-mail: massamh@yorku.ca

