arXiv:1508.03606v2 [math.PR] 7 Mar 2016

Hierarchical Models as Marginals of
Hierarchical Models
Guido Montúfar

Johannes Rauh

Max Planck Institute for

Department of Mathematics and Statistics

Mathematics in the Sciences

York University

montufar@mis.mpg.de

jarauh@yorku.ca
March 8, 2016

Abstract
We investigate the representation of hierarchical models in terms of marginals of
other hierarchical models with smaller interactions. We focus on binary variables and
marginals of pairwise interaction models whose hidden variables are conditionally independent given the visible variables. In this case the problem is equivalent to the
representation of linear subspaces of polynomials by feedforward neural networks with
soft-plus computational units. We show that every hidden variable can freely model
multiple interactions among the visible variables, which allows us to generalize and
improve previous results. In particular, we show that a restricted Boltzmann machine
with less than [2(log(v) + 1)/(v + 1)]2v − 1 hidden binary variables can approximate
every distribution of v visible binary variables arbitrarily well, compared to 2v−1 − 1
from the best previously known result.
Keywords: hierarchical model, restricted Boltzmann machine, interaction model, connectionism, graphical model

1

Introduction

Consider a finite set V of random variables. A hierarchical log-linear model is a set of
joint probability
distributions that can be written as products of interaction potentials, as
Q
p(x) = Λ ψΛ (x), where ψΛ (x) = ψΛ (xΛ ) only depends on the subset Λ of variables and
where the product runs over a fixed family of such subsets. By introducing hidden variables,
it is possible to express the same probability
in terms of potentials which involve
Pdistributions
Q
only small sets of variables, as p(x) =
ψ
(x,
y),
with small sets λ. Using small
λ
y
λ
interactions is a central idea in the context of connectionistic models, where the sets λ are
often restricted to have cardinality two. Due to the simplicity of their local characteristics,
these models are particularly well suited for Gibbs sampling [4]. The representation, or
explanation, of complex interactions among observed variables in terms of hidden variables
is also related to the study of common ancestors [16].
We are interested in sufficient and necessary conditions on the number of hidden variables, their values, and the interaction structures, under which the visible marginals are
flexible enough to represent any distribution from a given hierarchical model. Many problems can be formulated as special cases of this general problem. For example, the problem

of calculating the smallest number of layers of variables that a deep Boltzmann machine
needs in order to represent any probability distribution [9].
In this article, we focus on the case that all variables are binary. For the hierarchical
models with hidden variables, we restrict our attention to models involving only pairwise
interactions and whose hidden variables are conditionally independent given the visible variables (no direct interactions between the hidden variables). A prominent example of this
type of models is the restricted Boltzmann machine, which has full bipartite interactions between the visible and hidden variables. The representational power of restricted Boltzmann
machines has been studied assiduously; see, e.g., [3, 17, 7, 10]. The free energy
function of
P
such a model is a sum of soft-plus computational units x 7→ log(1 + exp( i∈V wi xi + c)).
On the other hand, the energy function of a fully observable hierarchical model with binary
variables is a polynomial, with monomials corresponding to pure interactions. Since any
function of binary variables can be expressed as a polynomial, the task is then to characterize the polynomials computable by soft-plus units.
Younes [17] showed that a hierarchical model with N binary variables and a total of
M pure higher order interactions (among three or more variables) can be represented as
the visible marginal of a pairwise interaction model with M hidden binary variables. In
Younes’ construction, each pure interaction is modeled by one hidden binary variable that
interacts pairwise with each of the involved visible variables. In fact, he shows that this
replacement can be accomplished without increasing the number of model parameters, by
imposing linear constraints on the coupling strengths of the hidden variable. In this work
we investigate ways of squeezing more degrees of freedom out of each hidden variable. An
indication thatthis should be possible is the fact that the full interaction model, for which
M = 2N − N2 − N − 1, can be modeled by a pairwise interaction model with 2N −1 − 1
hidden variables [10]. Indeed, by controlling groups of polynomial coefficients at the time,
we show that in general less than M hidden variables are sufficient.
A special case of hierarchical models with hidden variables are mixtures of hierarchical
models. The smallest mixtures of hierarchical models that contain other hierarchical models
have been studied in [8]. The approach followed there is different and complementary to
our analysis of soft-plus polynomials. For the necessary conditions, the idea there is to
compare the possible support sets of the limit distributions of both models. For the sufficient
conditions, the idea is to find a small S-set covering of the set of elementary events. An S-set
of a probability model is a set of elementary events such that every distribution supported
in that set is a limit distribution from the model.
Another type of hierarchical models with hidden variables are tree models. The geometry
of binary tree models has been studied in [18] in terms of moments and cumulants. That
analysis bears some relation to ours in that it also elaborates on Möbius inversions.
This paper is organized as follows. Section 2 introduces hierarchical models and formalizes our problem in the light of previous results. Section 3 pursues a characterization of the
polynomials that can be represented by soft-plus units. Section 4 applies this characterization to study the representation of hierarchical models in terms of pairwise interaction
models with hidden variables. This section addresses principally restricted Boltzmann machines. Section 5 offers our conclusions and outlook.

2

2

Preliminaries

This section introduces hierarchical models, with and without hidden variables, formalizes
the problem that we address in this paper, and presents motivating prior results.

2.1

Hierarchical Models

Consider a finite set V of variables with finitely many joint states x = (xi )i∈V ∈ X =
×i∈V Xi . We write v = |V | for the cardinality of V . For a given set S ⊆ 2V of subsets of
V let
(
)
X
VX,S := g(x) =
gΛ (x) : gΛ (x) = gΛ (xΛ ) .
Λ∈S

This is the linear subspace of R spanned by functions gΛ that only depend on sets of
variables Λ ∈ S. The hierarchical model of probability distributions on X with interactions
S is the set


1
exp(g(x)) : g ∈ VX,S ,
(1)
EX,S := p(x) =
Z(g)
P
where Z(g) = x0 ∈X exp(g(x0 )) is a normalizing factor. We call
X

E(x) = g(x) =

X

gΛ (x)

(2)

Λ∈S

the energy function of the corresponding probability distribution.
For convenience, in all what follows we assume that S is a simplicial complex, meaning
that A ∈ S implies B ∈ S for all B ⊆ A. Furthermore, we assume that the union of elements
of S equals V . In the case of binary variables, Xi = {0, 1} for all i ∈ V , the energy can be
written as a polynomial, as
X
Y
E(x) =
JΛ
xi .
Λ∈S

i∈Λ

Here, JΛ ∈ R, Λ ∈ S, are the interaction weights that parametrize the model.

2.2

Hierarchical Models with Hidden Variables

Consider an additional set H of variables with finitely many joint states y = (yj )j∈H ∈ Y =
V ∪H
×j∈H Yj . We write h = |H| for the cardinality of H. For a simplicial complex
,
P T ⊆2
X×Y
let VX×Y,T ⊆ R
be the linear subspace of functions of the form g(x, y) = λ∈T gλ (x, y),
gλ (x, y) = gλ ((x, y)λ ). The marginal on X of the hierarchical model EX×Y,T is the set




1 X
exp(g(x, y)) : g ∈ VX×Y,T ,
(3)
MX×Y,T := p(x) =


Z(g)
y∈Y

P
where Z(g) = x0 ∈X,y0 ∈Y exp(g(x0 , y 0 )) is again a normalizing factor. The free energy of a
probability distribution from MX×Y,T is given by

X
X
F (x) = log
exp
gλ (x, y) .
(4)
y∈Y

λ∈T

3

Here and throughout “log” denotes the natural logarithm.
In the case of binary visible variables, Xi = {0, 1} for all i ∈ V , the free energy (4) can
be written as a polynomial, as
X
Y
F (x) =
KB
xi ,
i∈B

B⊆V

where the coefficients can be computed from Möbius’ inversion formula as
X

X
X
KB =
(−1)|B\C| log
exp
gλ ((1C , 0V \C ), y) , B ∈ 2V .
C⊆B

y∈Y

(5)

λ∈T

Here (1C , 0V \C ) ∈ {0, 1}V is the vector with value 1 in the entries i ∈ C and value 0 in the
entries i 6∈ C.
If there are no direct interactions between hidden variables, i.e., if |λ ∩ H| ≤ 1 for all
λ ∈ T , then the sum over y factorizes and the free energy (4) can be written as

 X
X
X
X
gλ (x, yj ) .
(6)
F (x) =
gλ (x) +
log
exp
λ:λ∩H=∅

yj ∈Yj

j∈H

λ∈T :j∈λ

Particularly interesting are the models with full bipartite interactions between the set of
visible variables and the set of hidden variables, i.e., models with T = {λ ⊆ V ∪H : |λ∩V | ≤
1, |λ ∩ H| ≤ 1}, called restricted Boltzmann machines (with discrete variables).

2.3

Problem and Previous Results

In general the marginal of a hierarchical model is not a hierarchical model. However, one
may ask which hierarchical models are contained in the marginal of another hierarchical
model.
To represent a hierarchical model in terms of the marginal of another hierarchical model,
we need to represent (1) in terms of (3). Equivalently, we need to represent all possible
energy functions in terms of free energies. Given a set of visible variables V and a simplicial
complex S ⊆ 2V , what conditions on the set of hidden variables H and the simplicial
complex T ⊆ 2V ∪H are sufficient and necessary in order for any function E of the form (2)
to be representable in terms of some function F of the form (4)? We would like to arrive at
a result that generalizes the following.
• A restricted Boltzmann machine with h hidden binary variables can approximate any
probability distribution from a binary hierarchical model ES with |{Λ ∈ S : |Λ| > 1}| ≤
h arbitrarily well [17].
• The restricted Boltzmann machine with h = 2v−1 − 1 hidden binary variables can
approximate any probability distribution of v binary variables arbitrarily well [10].
Our Theorem 11 in Section 4 improves and generalizes these statements. The basis of this
result are soft-plus polynomials, which we discuss in the following section.

4

x1
x2

f (s) = log(1 + exp(s))

x3

φ

..
.
s

xv

Figure 1: Illustration of a soft-plus computational unit. The possible inputs X = {0, 1}V ,
corresponding to the vertices of the unit V -cube, are mapped to the real line by an affine
map x 7→ w> x + c, and then the soft-plus non-linearity f : s 7→ log(1 + exp(s)) is applied.

3

Soft-plus Polynomials

Consider a function of the form
φ : {0, 1}V → R; x 7→ log(1 + exp(w> x + c)),

(7)

parametrized by w = (wi )i∈V ∈ RV and c ∈ R. We regard φ as a soft-plus computational
unit, which integrates each input vector x ∈ {0, 1}V into a scalar via x 7→ w> x + c and
applies the soft-plus non-linearity f : R → R+ ; s 7→ log(1 + exp(s)). See Figure 1 for an
illustration of this function. In view of Equation (6), the function φ corresponds to the
free energy added by one hidden binary variable interacting pairwise with V visible binary
variables. The parameters wi , i ∈ V , correspond to the pair interaction weights and c to
the bias of the hidden variable.
What kinds of polynomials on {0, 1}V can be represented by soft-plus units? Following
Equation (5), the polynomial coefficients of φ are given by
!
X

X
|B\C|
KB (w, c) =
(−1)
log 1 + exp
wi + c
, B ∈ 2V .
(8)
i∈C

C⊆B

For each B ∈ 2V this is an alternating sum of the values φ(x) of the soft-plus unit on
the input vectors x ∈ {0, 1}V with supp(x) ⊆ B. In particular, KB is independent of the
parameters wi , i 6∈ B. We will use the shorthand notation wB for (wi )i∈B .
Note that, if wi = 0 for some i ∈ V , then KC = 0 for all C ∈ 2V with i ∈ C. In the
following we focus on the description of the possible values of the highest degree coefficients.
For example, Younes [17] showed that a soft-plus unit can represent a polynomial with an
arbitrary leading coefficient:
Lemma 1 (Lemma 1 in [17]). Let B ⊆ V and wi = 0 for i 6∈ B. Then, for any JB ∈ R,
there is a choice of wB ∈ RB and c ∈ R such that KB = JB .
The idea of Younes’ proof of Lemma 1 is to choose all non-zero wi of equal magnitude.
This simplifies the calculations and reduces the number of free parameters to one. Our goal
is to show that a soft-plus unit can actually freely model several polynomial coefficients at
5

K{1,2,3}

K{1,2,3,4}

{1, 2, 3, 4}

K{1,2}

{1, 2, 3}

{1, 2}

K{1,2}

{1, 3}

{1}

{1, 2, 4}

{1, 3, 4}

{2, 3}

{1, 4}

{2}

{3}

K{1,2,3}

{2, 3, 4}

{2, 4}

{3, 4}

K{1}

{4}

K{1}

K∅
∅

0
Figure 2: Illustration of Lemma 2. Depicted is for each
P edge pair
Q (B, B ) the set of
2
coefficient pairs (KB , KB 0 ) ∈ R of the polynomials
C⊆V KC
i∈C xi expressible as
log(1 + exp(w> x + c)). Shown is also the set of monomials of partial degree one and degree
at most 4, partially ordered by variable inclusion.

the same time. Our approach to simplify the Möbius inversion formula (8) is to choose the
parameters w and c in such a way that the function φ has many zeros. Clearly this can only
be done in an approximate way, since the soft-plus function is strictly positive. Nevertheless,
these approximations can be made arbitrarily accurate, since log(1 + exp(s)) ≤ exp(s) is
arbitrarily close to zero for sufficiently large negative values of s.
We call a pair of sets (B, B 0 ) an edge pair or a covering pair when B ) B 0 and there
is no set C with B ) C ) B 0 . The next lemma shows that a soft-plus unit can jointly
model the coefficients of an edge pair, at least in part. When the maximum degree |B| is at
most 3, the two coefficients are restricted by an inequality, but when |B| ≥ 4, there are no
such restrictions. The result is illustrated in Figure 2.
Lemma 2. Consider an edge pair (B, B 0 ). Depending on |B|, for any  > 0 there is a
choice of wB ∈ RB and c ∈ R such that k(KB , KB 0 ) − (JB , JB 0 )k ≤  if and only if
JB 0
JB 0

JB 0 ≥ 0, −JB ,
≥ 0, −JB or JB 0 ≤ 0, −JB ,
≥ 0, −JB or JB 0 ≤ 0, −JB ,
(JB , JB 0 ) ∈ R2 ,

for
for
for
for

|B| = 1
|B| = 2
|B| = 3
|B| ≥ 4.

Proof. This proof is deferred to Appendix A.
Remark 3. If (B, B 0 ) is an edge pair with |B| = 3, then, despite having |B| + 1 = 4
parameters to vary (wi , i ∈ B, and c), we can only determine the polynomial coefficients
KB and KB 0 up to a certain inequality. We expect that the same is true in general: If we
6

want to freely control k polynomial coefficients, we need strictly more than k parameters.
Otherwise, the coefficients are restricted by some inequalities. This situation is common
in models with hidden variables. In particular, mixture models often require many more
parameters to eliminate such inequalities than expected from naı̈ve parameter counting [8].
It is natural to ask whether it is possible to control other pairs of coefficients or even larger
groups of coefficients. We discuss a simple example before proceeding with the analysis of
this problem.
Example 4. Consider a soft-plus unit with two binary inputs. Write f : s 7→ log(1+exp(s))
for the soft-plus non-linearity and f0 = f (c), f1 = f (w1 + c), f2 = f (w2 + c), f12 =
f (w1 + w2 + c) for the values of the soft-plus unit on {0, 1}2 . From Equation (8) it is easy
to see that
K ∅ = f0
≥0
K{1} = f1 − f0 ≥ −K∅
K{2} = f2 − f0 ≥ −K∅ .
Now let us investigate the quadratic coefficient K{1,2} = f12 − f1 − f2 + f0 . Using the
convexity of f we find
0 ≤ K{1,2} ,
0 ≤ K{1,2} ≤ −K{1} , −K{2} ,
−K{1} ≤ K{1,2} ≤ 0,
−K{2} ≤ K{1,2} ≤ 0,

if
K{1} , K{2}
if −K{1} , −K{2}
if K{1} , −K{2}
if −K{1} , K{2}

≥0
≥0
≥0
≥ 0.

Hence the computable polynomials have coefficient triples (K{1} , K{2} , K{1,2} ) enclosed in
a polyhedral region of R3 as depicted in Figure 3. However, any pair (K{1} , K{2} ) ∈ R2 is
possible (for K∅ large enough).
The next lemma shows that a soft-plus unit can jointly model certain tuples of polynomial
coefficients corresponding to v − k + 1 monomials of degree k. We call star tuple a set of
the form {B ∪ {j} : j ∈ B 0 }, where B, B 0 ⊆ V satisfy B ∩ B 0 = ∅. Each element of the star
tuple covers the set B. In the Hasse diagram of the power set 2V , the sets B ∪ {j}, j ∈ B 0 ,
are the leaves of a star with root B.
Lemma 5. Consider any B, B 0 ⊆ V with B ∩ B 0 = ∅. Let wi = 0 for i 6∈ B ∪ B 0 . Then,
0
for any JB∪{j} ∈ R, j ∈ B 0 , and  > 0, there is a choice of wB∪B 0 ∈ RB∪B and c ∈ R such
that |KB∪{j} − JB∪{j} | ≤  for all j ∈ B 0 , and |KC | ≤  for all C 6= B, B ∪ {j}, j ∈ B 0 .

Proof of Lemma 5. Since wi = 0 for i 6∈ B ∪ B 0 , we have that KC = 0 for all C 6⊆ B ∪ B 0 .
)ω, wi = ω for all i ∈ B, and wj = JB∪{j} for j ∈ B 0 . Choosing
We choose
c = −(|B| − 21P
P
ω  j∈B 0 |wj | yields f ( i∈C wi + c) ≈ 0 for all C 6⊇ B. In this case,
KC ≈ 0,

for all B 6⊆ C ⊆ B ∪ B 0 .

Furthermore, for all j ∈ B 0 we have
X

X

KB∪{j} ≈f
wi + wj + c − f
wi + c
i∈B

=f (JB∪{j} +

i∈B
1
− f ( 2 ω)
− ( 12 ω) = JB∪{j} .

1
2 ω)

≈(JB∪{j} + 21 ω)

Similarly, KB∪C ≈ 0 for all C ⊆ B 0 with |C| ≥ 2. Note that KB ≈ 12 ω.
7

K{1,2}

K{2}
K{1}

Figure 3: Illustration of Example 4. Depicted is a region of R3 , clipped to [−1, 1]3 , which
contains the coefficient triples (K{1} , K{2} , K{1,2} ) ∈ R3 of the polynomials computable by
a soft-plus unit with two binary inputs. This region consists of 4 solid convex cones.
P
The intuition behind Lemma 5 is simple. When i∈B wi + c  1, the values w> x + c,
for x with xi = 1, i ∈ B, fall in a region where the soft-plus function is nearly linear. In
turn, the soft-plus unit is nearly a linear function of xj , j ∈ B 0 , with coefficients wj , j ∈ B 0 .
Remark 6. Closely related to soft-plus units are rectified linear units, which compute
functions of the form
ϕ : {0, 1}V → R; x 7→ max{0, w> x + c}.
In this case the non-linearity is s 7→ {0, s}. This reflects precisely the zero/linear behavior of
the soft-plus activation for large negative or positive values of s. Our polynomial descriptions
are based on this behavior and hence they apply both to soft-plus and rectified linear units.
We close this section with a brief discussion of dependencies among coefficients. The
next proposition gives a perspective on the possible values of the coefficient KB , depending
on wm , once wB\{m} and c have been fixed.
Proposition 7. Let (B, B 0 ) be an edge pair with B 0 = B \ {m} and let JB ∈ R. For fixed
0
wB 0 ∈ RB and c ∈ R, there is some wm ∈ R such that KB = JB if and only if a certain
0
degree-2|B |−1 polynomial in one real variable has a positive root.
Proof of Proposition 7. Observe that
KB (w, c) = KB 0 (wB 0 , c + wm ) − KB 0 (wB 0 , c).
Hence KB = JB if and only if KB 0 (wB 0 , c + wm ) = KB 0 (wB 0 , c) + JB =: r. We use the

8

abbreviation t̃ = et , which implies positivity. We have
KB 0 (wB 0 , c + wm ) =

X

(−1)|B

0

\C|

log 1 + exp

X

C⊆B 0

wi + c + wm



!

i∈C


= log 

Y 

1 + w̃m c̃

C⊆B 0

Y

w̃i

(−1)|B0 \C|


.

i∈C

Now, KB 0 (wB 0 , c + wm ) = r if and only if
Y 

1 + w˜m c̃

C⊆B 0

Y

w̃i

(−1)|B0 \C|

= r̃,

i∈C

or, equivalently,
Y



1 + w̃m c̃

C⊆B 0 :
B 0 \C even

Y

w̃i



− r̃

Y 

1 + w̃m c̃

C⊆B 0 :
B 0 \C odd

i∈C

This is a polynomial of degree at most 2|B

0

|−1

Y


w̃i = 0.

i∈C

in w̃m = ewm .

This description implies various kinds of constraints. For example, by Descartes’ rule of
signs, a polynomial can only have positive roots if the sequence of polynomial coefficients,
ordered by degree, has sign changes.

4

Conditionally Independent Hidden Variables

In the case of a bipartite graph between V and H with all variables binary, the hierarchical
model (or rather its visible marginal) is a restricted Boltzmann machine, denoted RBMV,H .
This model is illustrated in Figure 4. The free energy takes the form
!
X

X
X
F (x) =
bi x i +
log 1 + exp
wji xi + cj
.
i∈V

j∈H

i∈V

This is the sum of an arbitrary degree-one polynomial, with coefficients bi ∈ R, i ∈ V ,
and h = |H| independent soft-plus units, with parameters wji ∈ R, j ∈ H, i ∈ V and
cj ∈ R, j ∈ H. The free energy contributed by the hidden variables can be thought of as a
feedforward network with soft-plus computational units.
We can use each soft-plus unit to model a group of coefficients of any given polynomial,
starting at the highest degrees. Using the results from Section 3 we arrive at the following
representation result:
Theorem 8. Every distribution from a hierarchical model ES on {0, 1}V can be approximated arbitrarily well by distributions from RBMV,H whenever there exist h sets B1 , . . . , Bh ⊆
2V which cover {Λ ∈ S : |Λ| ≥ 2} in reverse inclusion order, where each Bj is a star tuple
or an edge pair of sets of cardinality at least 3.

9

c1

c2

c3

c4

y1

y2

y3

y4

ch
yh

···

w11

whv
x1

x2

x3

b1

b2

b3

···

xv
bv

Figure 4: A restricted Boltzmann machine. The free energy contributed by the hidden units
is a sum of independent soft-plus units.
Proof of Theorem 8. We need to express the possible energy functions of the hierarchical
model as sums of independent soft-plus units plus linear terms. This problem can be reduced
to covering the appearing monomials of degree two or more by groups of coefficients that
can be jointly controlled by soft-plus units. In view of Lemmas 2 and 5, edge pairs with
sets of cardinality 3 or more and star tuples can be jointly controlled. We start with the
highest degrees and cover monomials downwards, because setting the coefficients of a given
group may produce uncontrolled values for the coefficients of smaller monomials. Since S is
a simplicial complex, we only need to cover the elements of S.
Finding a minimal covering is in general a hard combinatorial problem. In the following
we derive upper bounds for the k-interaction model, which is the hierarchical model ES with
S = Sk := {Λ ⊆V : |Λ| ≤ k}. We will focus on star tuples and consider individual coverings
of the layers Vj = {Λ ⊆ V : |Λ| = j}. Let v = |V |. Denote D(v, j) the smallest number

of star tuples that cover Vj . We use the following notion from the theory of combinatorial
designs (see [1] for an overview on that subject).
For integers v ≥ k > r denoteC(v, k, r) the

smallest possible number of elements of Vk such that every element from Vr is contained
in at least one of them.

Lemma 9. For 0 < j ≤ v, the minimal number of star tuples that cover Vj is D(v, j) =
C(v, v − j + 1, v − j). Inserting known results for C(v, t + 1, t) we obtain the exact values
D(v, 1)
D(v, 2)

= 1
= vl − 1 l
mm
4
v
v−1
D(v, 3)
= v−2
v−3 · · · 2 · · ·
 
 v−2 
D(v, v − 3) = v4 v−1
(v 6≡ 7 mod 12)
3
 v  v−1  2
D(v, v − 2) = 3 2
 
D(v, v − 1) = v2
D(v, v)
= 1
and the general bound
 
1 + log(v − j + 1) v
D(v, j) ≤
, 0 < j ≤ v.
v−j+1
j

Furthermore, we have the simple bound D(v, j) ≤ v−1
j−1 , 0 < j ≤ v.
10


Proof of Lemma 9. A star tuple covering of Vj is given by a collection B1 , . . . , Bn of ele

V
ments of j−1
such that every element of Vj contains at least one of the Bi . The minimal
possible number of elements in such a collection is precisely D(v, j) = C(v, v − j + 1, v − j).
The equalities follow from corresponding equalities for C(v, v − j + 1, v − j) by several authors, which are listed in [13].
follows
from a result by Erdős and Spencer [2]
h The
. inequality
ih
i
k
v
k
showing that C(v, k, r) ≤ r
1 + log r . The simple bound results from the fact
r


V
\{1}
that each set B from j contains a set B 0 from Vj−1
.
Remark 10. Lemma 9 presents widely applicable bounds on the cardinality of star tuple
coverings, which are naturally not always tight. For v ≤ 28, better individual bounds on
C(v, t + 1, t) can be found in [13, Table III]. See also [14] for a list of known exact values.
In another direction [15] offers optimal asymptotic bounds on C(v, k, r) for fixed k and r.
Lemma 9 allows us to formulate the following more explicit version of Theorem 8:
Theorem 11. Let 1 ≤ k ≤ v. Every distribution from the k-interaction model ESk on
{0, 1}V can be approximated arbitrarily well by distributions from RBMV,H whenever h
Pk
surpasses or equals U (v, k) = j=2 D(v, j), which is bounded above as indicated in Lemma 9.


Pk
log(v−1)+1 Pk
v+1
This is the case, in particular, whenever h ≥ j=2 v−1
j=2
v+1
j .
j−1 or h ≥
Proof of Theorem 11. This follows directly from Theorem 8 and Lemma 9. For
 the last
statement we use the simple bound from the lemma, by which D(v, j) ≤ v−1
j−1 , and the
log(v−j+1)+1 v+1
general bound, by which D(v, j) ≤
j .
v+1
In order to provide a numerical sense of Theorem 11 we give upper bounds on U (v, k),
2 ≤ k ≤ v ≤ 14, in Table 1. For convenience we also provide an Octave [5] script for
computing such bounds in http://personal-homepages.mis.mpg.de/montufar/starcover.m.
In the special case k = v, the k-interaction model ESk is the full interaction model
and contains all (strictly positive) probability distributions on {0, 1}V . Hence Theorem 11
entails the following universal approximation result:
Corollary 12. Every distribution on {0, 1}V can be approximated arbitrarily well by distributions from RBMV,H whenever h surpasses or equals U (v, v), which is bounded above
as indicated in Lemma 9. This is the case, in particular, whenever h ≥ 2v−1 − 1 or
h ≥ 2(log(v−1)+1)
(2v − (v + 1) − 1) + 1.
v+1
Corollary 12 provides a significant and unexpected improvement the best previously
known upper bound 2v−1 − 1 from [10]. Whether the upper bound 2v−1 − 1 was optimal
or not had remained an open problem in [10] and several succeeding papers. In Table 2 we
give upper bounds on U (v, v), 2 ≤ v ≤ 40, and compare these with the previous result.
Remark 13. In general an RBM can represent many more distributions than just the
interaction models described above. For several small examples discussed further below,
our bounds for the representation of interaction models are tight. However, Theorem 11 is
based on upper bounds on a specific type of coverings and we suspect that it can be further
improved, at least in some special cases, even if not reaching the hard lower bounds coming
from parameter counting.

11

k\v
2

2
1

3
2

4
3

5
4

6
5

7
6

8
7

9
8

3

-

3

5

8

11

15

19

24

29

55

39

82

54

117

98

69

147

84

234

124

356

182

251

453

80

172

343

570

908

1385

4
5

1

-

3

-

6
7

-

11
12
15

17
20

27
34

53

6

-

-

-

-

21

38

64

7

-

-

-

-

-

39

68

8

-

-

-

-

-

-

9

-

-

-

-

-

10

-

-

-

-

11

-

-

-

12

-

-

13

-

14

-

31

63

84
85

109
184

121
189

10
9

74

175
373

205
390

12
11

13
12

14
13

35

41

48

55

162

216

125

268

341

520

725

282

427

742

1276

789

1534

808

1591

814

1615

815

1621

348

69

126

-

-

127

227

-

-

-

-

228

-

-

-

-

-

-

-

-

-

-

-

-

-

-

-

-

-

-

-

-

-

-

-

-

-

-

-

-

-

-

127

190
255

222

11
10

395
396
511

395
414
420
421

1023

559
672
729
753

759

1622

750

2107

1014
2705

1259
3078

1416
3156

1494
3182

1520
3189

1002
2068

1473
3389

1944
4652

2452
5583

2823
6105

3053
6196

3144
6229

760

1527

-

-

1528

3184

-

-

-

3185

2047

3190
4095

3177
6236
6237
8191

Table 1: Upper bounds on the minimal number of hidden units for which RBMV,H can
approximate every distribution from the k-interaction model ESk on {0, 1}V arbitrarily well,
following from Theorem 11, for 2 ≤ k ≤ v ≤ 14. Shown are upper bounds on U (v, k) =
Pk
j=2 D(v, j) evaluated using Lemma 9 and some individual bounds on D(v, j) = C(v, v −
j + 1, v − j) from [13, Table III]. Upper scripts indicate values obtained using only Lemma 9.
Lower scripts indicate the previous RBM universal approximation bound 2v−1 − 1 from [10].
Entries with v ≤ 9 or k ≤ 3 are exact values of U (v, k).

12

v
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

U (v, v) ≤

1
3
6
12
21
39
69
127
228
421
760
1528
3185
6642
14,269
30,352
63,431
132,195
272,160
553,195
1115,207
2227,484
4427,830
8760,826
17,265,199
33,951,316
66,656,315
132,084,407
257,962,181
504,141,876
985,875,453
1929,093,753
3776,867,237
7398,516,744
14,500,416,431
28,433,369,622
55,779,952,400
109,476,401,847
214,954,581,277

2v−1 − 1 =

1
3
7
15
31
63
127
255
511
1023
2047
4095
8191
16,383
32,767
65,535
131,071
262,143
524,287
1048,575
2097,151
4194,303
8388,607
16,777,215
33,554,431
67,108,863
134,217,727
268,435,455
536,870,911
1073,741,823
2147,483,647
4294,967,295
8589,934,591
17,179,869,183
34,359,738,367
68,719,476,735
137,438,953,471
274,877,906,943
549,755,813,887

l

2v
v+1

m
−1 =

1
1
3
5
9
15
28
51
93
170
315
585
1092
2047
3855
7281
13,797
26,214
49,932
95,325
182,361
349,525
671,088
1290,555
2485,513
4793,490
9256,395
17,895,697
34,636,833
67,108,863
130,150,524
252,645,135
490,853,405
954,437,176
1857,283,155
3616,814,565
7048,151,460
13,743,895,347
26,817,356,775

Table 2: Bounds on the minimal number of hidden units for which RBMV,H can approximate every distribution on {0, 1}V arbitrarily well, for 2 ≤ v ≤ 40. The first column
Pv gives upper bounds following from Corollary 12. Shown are upper bounds on U (v, v) = j=2 D(v, j)
evaluated using Lemma 9 and some individual bounds on D(v, j) = C(v, v − j + 1, v − j)
from [13, Table III]. The second column giveslthe previous
upper bound 2v−1 − 1 from [10].
m

The last column gives the hard lower bound

2v
v+1

− 1 that results from parameter count-

ing, i.e., from demanding that the model RBMV,H has at least (h + 1)(v + 1) − 1 ≥ 2v − 1
parameters.
13

Besides from RBMs we can also consider models that include interactions among the
visible variables other than biases. In this case we only need to cover the interaction sets
from the simplicial complex S that are not already included in the simplicial complex T . In
Theorem 8 one just replaces {Λ ∈ S : |Λ| ≥ 2} by S \ T . We note the following special case:
Corollary 14. Every distribution from the k-interaction model ESk on {0, 1}V can be approximated
Pk arbitrarily well by the visible marginals of a pairwise interaction model with
h =
j=3 D(v, j) hidden binary variables. The latter is bounded above as indicated in
Lemma 9. In particular, every distribution on {0, 1}V can be approximated arbitrarily well
by the visible marginals of a pairwise interaction model with h = 2v−1 − (v − 1) − 1 or
(2v − (v + 1) − 1 − (v+1)v
) + 1 hidden binary variables.
h = 2(log(v−2)+1)
v+1
4
Corollary 14 improves a previous
 result by Younes [17], which showed that a pairwise
interaction model with h = 2v − v2 − v − 1 hidden binary variables can approximate every
distribution on {0, 1}V arbitrarily well.
We close this section with a few small examples illustrating our results.
Example 15. The model RBM3,1 is the same as the two-mixture of product distributions
of 3 binary variables and is also known as the tripod tree model. It has 7 parameters and
the same dimension. What is the largest hierarchical model contained in the closure of this
model? The closure of a model M is the set of all probability distributions that can be
approximated arbitrarily well by probability distributions from M.
The closure of RBM3,1 contains all 3 hierarchical models on {0, 1}3 with two pairwise
interactions. For example, it contains the model ES with S = {{1, 2}, {1, 3}, {1}, {2}, {3}}.
Indeed, two quadratic coefficients can be jointly modeled by one soft-plus unit (Lemma 5)
and the linear coefficients with the biases of the visible variables. In particular, the closure
of RBM3,1 also contains the 3 hierarchical models with a single pairwise interaction.
It does not contain the hierarchical model with 3 pairwise interactions, ES with S = S2 =
{{1, 2}, {1, 3}, {2, 3}, {1}, {2}, {3}}, which is known as the no-three-way interaction model.
One way of proving this is by comparing the possible support sets of the two models, as
proposed in [8]. The support set of a product distribution is a cylinder set. The support set
of a mixture of two product distributions is a union of two cylinder sets. On the other hand,
the possible support
Q sets of a hierarchical model correspond to the faces of its marginal
polytope, conv{( i∈Λ xi )Λ∈S : x ∈ X} ⊂ RS . The marginal polytope of the no-three-way
interaction model is the cyclic polytope C(N, d) with N = 8 vertices and dimension d = 6
(see, e.g., [8, Lemma 18]). This is a neighborly polytope, meaning that every d/2 = 3
or less vertices form a face. In turn, every subset of {0, 1}3 of cardinality d/2 = 3 is the
support set of a distribution in the closure of the no-three-way interaction model.1 Since the
set {(100), (010), (001)} is not a union of two cylinder sets, the closure of RBM3,1 does not
contain the no-three-way interaction model.
Example 16. The closure of RBM3,2 contains the no-three-way interaction model. Two of
the quadratic coefficients can be jointly modeled with one hidden unit and the third with
the second hidden unit (Lemma 5).
It does not contain the full interaction model. Following the ideas explained in the previous example, this can be shown by analyzing the possible support sets of the distributions
in the closure of RBM3,2 . For details on this we refer the reader to [12].
1 More

2k

generally, in [6] it is shown that if S ⊇ {Λ ⊆ V : |Λ| ≤ k}, then the marginal polytope of ES is
− 1 neighborly, meaning that any 2k − 1 or fewer of its vertices define a face.

14

Example 17. The model RBM3,3 is a universal approximator. This follows immediately
from the universal approximation bound 2v−1 − 1 from [10]. This observation can be recovered from our results as follows. The cubic coefficient can be modeled with one hidden unit
(Lemma 1). Two quadratic coefficients can be jointly modeled with one hidden unit and
the third with another hidden unit (Lemma 5).
Example 18. The model RBM4,6 is a universal approximator. The quartic coefficient can
be modeled with one hidden unit. The 4 cubic coefficients can be modeled with two hidden
units (Lemma 5). The 6 quadratic coefficients can be grouped into 3 pairs with a shared
variable in each pair. These can be modeled with 3 hidden units (Lemma 5).

5

Conclusions and Outlook

We studied the kinds of interactions that appear when marginalizing over a hidden variable
that is connected by pair-interactions with all visible variables. We derived upper bounds on
the minimal number of variables of a hierarchical model whose visible marginal distributions
can approximate any distribution from a given fully observable hierarchical model arbitrarily
well. These results generalize and improve previous results on the representational power of
RBMs from [10] and [17].
Many interesting questions remain open at this point: A full characterization of soft-plus
polynomials and the necessary number of hidden variables is missing.
It would be interesting to look at non-binary hidden variables. This corresponds to
analyzing the hierarchical models that can be represented by mixture models. In the case of
conditionally independent binary hidden variables, the partial factorization leads to soft-plus
units, whereas in the case of higher-valued hidden variables, it leads to shifted logarithms
of denormalized mixtures. Similarly, it would be interesting to take a look at non-binary
visible variables. In this case state vectors cannot be identified in a one-to-one manner
with subsets of variables. This means that the correspondence between function values and
polynomial coefficients is not as direct.
Our analysis could also be extended to cover the representation of conditional probability
distributions from hierarchical models in terms of conditional restricted Boltzmann machines
and to refine the results on this problem reported in [11].
Another interesting direction are models where the hidden variables are not conditionally
independent given the visible variables, such as deep Boltzmann machines, which involve
several layers of hidden variables. This case is more challenging, since the free energy does
not decompose into independent terms.

Acknowledgments
We thank Nihat Ay for helpful remarks with the manuscript.

References
[1] C. J. Colbourn and J. H. Dinitz. Handbook of Combinatorial Designs, Second Edition
(Discrete Mathematics and Its Applications). Chapman & Hall/CRC, 2006.
[2] P. Erdős and J. Spencer. Probabilistic Methods in Combinatorics. Academic Press Inc,
1974.
15

[3] Y. Freund and D. Haussler. Unsupervised learning of distributions on binary vectors
using two layer networks. In J. E. Moody, S. J. Hanson, and R. P. Lippmann, editors, Advances in Neural Information Processing Systems 4, pages 912–919. MorganKaufmann, 1992. http://papers.nips.cc/paper/535-unsupervised-learningof-distributions-on-binary-vectors-using-two-layer-networks.pdf.
[4] S. Geman and D. Geman. Stochastic relaxation, Gibbs distributions, and the Bayesian
restoration of images. Pattern Analysis and Machine Intelligence, IEEE Transactions
on, PAMI-6(6):721–741, 1984. http://dx.doi.org/10.1109/TPAMI.1984.4767596.
[5] S. H. John W. Eaton, David Bateman and R. Wehbring. GNU Octave version 4.0.0
manual: a high-level interactive language for numerical computations. 2015. http:
//www.gnu.org/software/octave/doc/interpreter.
[6] T. Kahle. Neighborliness of marginal polytopes. Beiträge zur Algebra und Geometrie,
51(1):45–56, 2010. http://www.emis.de/journals/BAG/vol.51/no.1/4.html.
[7] N. Le Roux and Y. Bengio. Representational power of restricted Boltzmann machines
and deep belief networks. Neural Computation, 20(6):1631–1649, June 2008. http:
//dx.doi.org/10.1162/neco.2008.04-07-510.
[8] G. Montúfar. Mixture decompositions of exponential families using a decomposition of
their sample spaces. Kybernetika, 49(1):23–39, 2013. http://www.kybernetika.cz/
content/2013/1/23.
[9] G. Montúfar. Deep narrow Boltzmann machines are universal approximators. International Conference on Learning Representations 2015 (ICLR 2015), San Diego, CA,
USA, 2015. http://www.arxiv.org/abs/1411.3784.
[10] G. Montúfar and N. Ay. Refinements of universal approximation results for deep belief
networks and restricted Boltzmann machines. Neural Computation, 23(5):1306–1319,
2011. http://dx.doi.org/10.1162/NECO_a_00113.
[11] G. Montúfar, N. Ay, and K. Ghazi-Zahedi. Geometry and expressive power of conditional restricted Boltzmann machines. Journal of Machine Learning Research, 16:2405–
2436, 2015. http://jmlr.org/papers/v16/montufar15b.html.
[12] G. Montúfar and J. Morton. When does a mixture of products contain a product of
mixtures? SIAM Journal on Discrete Mathematics, 29:321–347, 2015. http://dx.
doi.org/10.1137/140957081.
[13] K. J. Nurmela and P. R. J. Östergård. New coverings of t-sets with (t+1)-sets. Journal
of Combinatorial Designs, 7(3):217–226, 1999. http://dx.doi.org/10.1002/(SICI)
1520-6610(1999)7:3<217::AID-JCD5>3.0.CO;2-W.
[14] OEIS. The on-line encyclopedia of integer sequences, A066010 triangle of covering
numbers T(n,k) = C(n,k,k-1), n >= 2, 2 <= k <= n. Published electronically at
http://oeis.org, 2010.
[15] V. Rödl. On a packing and covering problem. European Journal of Combinatorics,
6(1):69–78, 1985. http://dx.doi.org/10.1016/S0195-6698(85)80023-8.

16

[16] B. Steudel and N. Ay. Information-theoretic inference of common ancestors. Entropy,
17(4):2304, 2015. http://dx.doi.org/10.3390/e17042304.
[17] L. Younes. Synchronous Boltzmann machines can be universal approximators. Applied Mathematics Letters, 9(3):109–113, 1996. http://dx.doi.org/10.1016/08939659(96)00041-9.
[18] P. Zwiernik and J. Q. Smith. Tree cumulants and the geometry of binary tree models.
Bernoulli, 18:290–321, 2012. http://dx.doi.org/10.3150/10-BEJ338.

A

Proofs

Proof of Lemma 2. Let B 0 = B \ {m}. The edge coefficients satisfy
KB 0 (wB 0 , c) =

X

(−1)

|B 0 \C|

C⊆B 0

log 1 + exp

X


wi + c

!

i∈C

and
KB (wB , c) = KB 0 (wB 0 , c + wm ) − KB 0 (wB 0 , c).
Using this structure, we now proceed with the proof of the individual cases.
The case |B 0 | = 0. We omit this simple exercise.
The case |B 0 | = 1. The if statement is as follows. The elements of the set {0, 1}B are the
vertices of the |B|-dimensional unit cube. We call two vectors x, x0 ∈ {0, 1}B adjacent if
they differ in exactly one entry, in which case they are the vertices of an edge of the cube.
The weights wB and c can be chosen such that the affine map {0, 1}B → R; xB 7→
>
wB xB + c maps any chosen pair of adjacent vectors to any arbitrary values and all other
vectors to large negative values. The soft-plus function is monotonically increasing, taking
value zero at minus infinity and plus infinity at plus infinity. Hence, for any s, s0 ∈ R+ , one
finds weights w and c such that

(xB 0 , xm ) = (1, . . . , 1, 1)
 s,
s0 ,
(xB 0 , xm ) = (1, . . . , 1, 0) ,
φ(x) =

≈ 0, otherwise
or, alternatively, such that

(xB 0 , xm ) = (1, . . . , 1, 0, 1)
 s,
s0 ,
(xB 0 , xm ) = (1, . . . , 1, 0, 0) .
φ(x) =

≈ 0, otherwise
This leads to KB ≈ (s − s0 ) and KB 0 ≈ s0 or, alternatively, KB ≈ −(s − s0 ) and KB 0 ≈ −s0 .
The approximation can be made arbitrarily precise.
The only if statement is as follows. Denote the soft-plus function by f : R → R+ ; s 7→
log(1 + exp(s)). Since |B 0 | = 1, C ⊆ B 0 implies C = B 0 or C = ∅. We have that
KB 0 (wB 0 , c) = f (wB 0 + c) − f (c) and KB 0 (wB 0 , c + wm ) = f (wB 0 + c + wm ) − f (c + wm ) are
either both positive or both negative, depending on the sign of wB 0 . If both are positive,
then KB (wB , c) = KB 0 (wB 0 , c + wm ) − KB 0 (wB 0 , c) ≥ −KB 0 (wB 0 , c), and similarly in the
case that both are negative.
17

The case |B 0 | = 2. The if statement follows from the previous case |B 0 | = 1. Indeed,
consider an edge pair (C, C 0 ) with an element more than the edge pair (B, B 0 ), such that
B = C \ {n} and B 0 = C 0 \ {n}. Then, for any wB and c, choosing wn large enough
one obtains an arbitrarily accurate approximation KC ((wB , wn ), c − wn ) ≈ KB (wB , c) and
KC 0 ((wB 0 , wn ), c − wn ) ≈ KB 0 (wB 0 , c).
For the only if statement we use a similar argument as previously. We have KB 0 (wB 0 , c) =
f (w1 + w2 + c) + f (c) − f (c + w1 ) − f (c + w2 ). By convexity of f , this is non-negative if and
only if either w1 , w2 ≥ 0 or w1 , w2 ≤ 0. In other words, this is non-negative if and only if
w1 · w2 ≥ 0. Under either of these conditions, KB 0 (wB 0 , c + wm ) is also non-negative. Similarly, KB 0 (wB 0 , c) is non-positive if and only if w1 · w2 ≤ 0. In this case, KB 0 (wB 0 , c + wm )
is also non-positive. Now the statement follows as in the case |B 0 | = 1.

The case |B 0 | = 3. We need to show that any edge pair coefficients can be represented.
Consider first JB 0 ≥ 0. We choose weights of the form wB 0 = ω1B 0 , where ω ∈ R and 1B 0 is
the vector of |B 0 | ones. Then KB 0 (wB 0 , c) = f (3ω + c) − 3f (2ω + c) + 3f (ω + c) − f (c). We
can choose ω and c such that 3ω + c = f −1 (JB 0 ) while 2ω + c, ωc, c take very large negative
values. This yields KB 0 ≈ JB 0 .
Note that the derivative of the soft-plus function is the logistic function, i.e., f 0 (s) =
1/(1+exp(−s)). Choosing ω large enough from the beginning, the function wm 7→ KB 0 (wB 0 , c+
wm ) is monotonically increasing in the interval wm ∈ [0, ω/2] and surpasses the value
1
5 ω. On the other hand, when wm is large enough, depending on ω and c, we have that
5
5
2ω + c + wm ≥ 12
(3ω + c + wm ) and f (2ω + c + wm ) ≥ 12
f (3ω + c + wm ). In this case
1
1
f (3ω + c + wm ) − 3f (2ω + c + wm ) ≤ − 4 (3ω + c + wm ) ≤ − 4 ω. At the same time, ω + c + wm
1
and c + wm are smaller than − 12
ω and so f (ω + c + wm ) and f (c + wm ) are very small in
absolute value.
By the mean value theorem, depending on wm , KB 0 (wB 0 , c + wm ) takes any value in the
0
interval [− 51 ω, 15 ω], where ω is arbitrarily large. In turn, we can obtain KB = KB 0 (wB
,c +
wm ) − KB 0 (wB 0 , c) ≈ JB for any JB ∈ R.
For JB 0 ≤ 0 the proof is analogous after label switching for one variable.
The case |B 0 | > 3. This follows from the previous case |B 0 | = 3 in the same way that the
if part of the case |B 0 | = 2 follows from the case |B 0 | = 1.

18

