Collapsed Variational Inference for Sum-Product Networks
Han Zhao†
HAN . ZHAO @ CS . CMU . EDU
Tameem Adel§
T.M.A.A.H ESHAM @ UVA . NL
Geoff Gordon†
GGORDON @ CS . CMU . EDU
Brandon Amos†
BAMOS @ CS . CMU . EDU†
School of Computer Science, Carnegie Mellon University, Pittsburgh, PA, USA§
Machine Learning Lab, University of Amsterdam, Amsterdam, the Netherlands; Radboud University

Abstract
Sum-Product Networks (SPNs) are probabilistic
inference machines that admit exact inference inlinear time in the size of the network. Existing parameter learning approaches for SPNs arelargely based on the maximum likelihood principle and are subject to overfitting compared tomore Bayesian approaches. Exact Bayesian posterior inference for SPNs is computationally intractable. Even approximation techniques suchas standard variational inference and posterior
sampling for SPNs are computationally infeasible even for networks of moderate size due tothe large number of local latent variables per
instance. In this work, we propose a novel
deterministic collapsed variational inference algorithm for SPNs that is computationally efficient, easy to implement and at the same timeallows us to incorporate prior information into
the optimization formulation. Extensive experiments show a significant improvement in accuracy compared with a maximum likelihood based
approach.

1. Introduction
Parameter learning is an important research problem in
Sum-Product Networks (SPNs) (Poon & Domingos, 2011;
Gens & Domingos, 2012). Due to recent progress, we
can locally maximize the likelihood of SPNs, but those approaches overfit for complex models because they do not
retain posterior uncertainty or allow flexible incorporationof prior information. SPNs are essentially mixture models
with potentially exponentially many components (Zhao &Poupart, 2015). Exact posterior inference for SPNs is comProceedings of the 33 rd International Conference on Machine
Learning, New York, NY, USA, 2016. JMLR: W&CP volume
48. Copyright 2016 by the author(s).

putationally intractable due to the exponential growth of
the number of mixture components in the posterior distribution. Approximate Bayesian inference for SPNs is challenging, mostly because of the vast number of local latentvariables and their intricate connections. Unlike traditional
mixture models, e.g., LDA (Blei et al., 2003), Gaussian
mixture models, etc., where the number of local latent variables is a small constant per observation, the number of
local latent variables in SPNs is linearly proportional to
the size of the network for each observable instance, often around thousands to millions. This property makes
standard approximate Bayesian inference approaches, including mean-field variational inference (VI) (Wainwright& Jordan, 2008) and posterior sampling (Casella & Berger,
2002), computationally infeasible even for SPNs with moderate size. In SPNs, the memory overhead incurred by the
number of local variational parameters or the number of local sampling points increases linearly with both the numberof data points, D, and the size of the network, |S|, and so
can easily become prohibitively large.
In this work we propose a collapsed variational inference
algorithm for SPNs that is robust to overfitting and canbe naturally extended into a stochastic variant to scale tolarge data sets. We call our algorithm CVB-SPN. CVBSPN is a deterministic approximate Bayesian inference algorithm that is computationally efficient and easy to implement while at the same time allowing us to incorporate
prior information into the design. Like other variational
techniques, CVB-SPN is based on optimization. Unlike
other variational techniques, the number of parameters to
be optimized is only linear in the network size (O(|S|)) andis independent of the size of the training set, as opposed toO(D|S|) in ordinary VI. It is worth noting that, different from traditional collapsed variational approaches in theliterature that marginalize out global latent variables (Tehet al., 2006; 2007), here we consider a complementary approach: instead of marginalizing out the global latent variables in order to spread out the interactions among many
local latent variables, CVB-SPN takes advantage of the fastCollapsed Variational Inference for Sum-Product Networks

exact inference in SPNs and marginalizes out the local latent variables in order to maintain a marginal variational
posterior distribution directly on the global latent variables,i.e., the model parameters. The posterior mean of the variational distribution over model parameters can then be used
as a Bayesian estimator. To the best of our knowledge, thisis the first general Bayesian approach to learn the parameters of SPNs efficiently in both batch and stochastic settings.

graph where the interior nodes are sums or products andthe leaves are univariate distributions over Xi . Edges emanating from sum nodes are parameterized with positive
weights. Each node in an SPN encodes an unnormalized
marginal distribution over X.

At first glance, marginalizing out all the local latent variables in graphical models seems to be a bad idea for two
reasons. First, by marginalizing out local latent variables
we naively appear to incur computation exponential in thetree-width of the graph (Jordan et al., 1999; Wainwright& Jordan, 2008). Except for graphical models with special structures, for example, LDA, Gaussian mixture models, thin junction trees, etc., such exact computation is intractable by itself. Second, marginalization will in general invalidate the conjugacy between the prior distribution
and the joint distribution over global and local latent variables, which further makes the expectation over the variational posterior intractable. Fortunately, as we will showin Sec. 3, the ability of SPNs to model context-specificindependence helps to solve the first problem, and by using a change of variables CVB-SPN handles the second
problem efficiently. Besides the reduced space complexity, we also show that the objective of CVB-SPN forms astrictly better lower bound to be optimized than the evidence lower bound (ELBO) in standard VI. To show thevalidity of CVB-SPN, we conduct extensive experiments
and compare it with maximum likelihood based methods
in both batch and stochastic settings.

8>= xi )<p(X
Q i
Vk (x | w) =j2Ch(k) Vj (x | w)>
:P
j2Ch(k) wkj Vj (x | w)

To tackle the second problem described above, there is a
strand of recent work on extending standard VI to general
nonconjugate settings using stochastic sampling from thevariational posterior (Blei et al., 2012; Kingma & Welling,2013; Ranganath et al., 2014; Mnih & Gregor, 2014; Titsias& Lázaro-Gredilla, 2014; Titsias, 2015). However, control
variates need to be designed in order to reduce the highvariance incurred by insufficient samples. Furthermore, foreach such sample from the variational posterior, those algorithms need to go through the whole training data set tocompute the stochastic gradient, leading to a total computational cost O(N D|S|), where N is the sample size. This
is often prohibitively expensive for SPNs.

2. Background2.1. Sum-Product Networks
A sum-product network (SPN) is a graphical representationof a joint probability distribution over a set of random variables X = {X1 , . . . , Xn }. It is a rooted directed acyclic

In more detail, let x be an instantiation of the random vector X. We associate an unnormalized probability Vk (x |w) with each node k when the input to the network is x
with network weights set to be w:
k is a leaf node over Xi
k is a product node
k is a sum node(1)

where Ch(k) is the child list of node k in the graph andwkj is the edge weight associated with sum node k and its
child node j. The joint distribution encoded by an SPN
is then defined by the graphical structure and the weights.
The probability/density of a joint assignment X = x is
proportional to the value at the root of the SPN with input
x divided by a normalization constant Vroot (1 | w):
p(x) =

Vroot (x | w)
Vroot (1 | w)(2)

where the normalization constant Vroot (1 | w) is obtainedby setting V at all the leaf nodes to be 1 and then propagating toward the root according to Eq. 1. Intuitively, setting V at all the leaf nodes to be 1 corresponds to integrating/marginalizing out the random vector X, which will ensure Eq. 2 defines a proper probability distribution. Eq. 2
can also be used to compute the marginal probability ofa partial assignment Y = y: simply set V at leaf nodes
whose corresponding random variable is not in Y to be 1
and other leaf nodes based on the assignment Y = y. Intuitively, this corresponds to integrating out variables outside of the partial assignment. We can compute conditional
probabilities by evaluating two partial assignments:p(Y = y | Z = z) =p(Y = y, Z = z)
Vroot (y, z | w)=p(Z = z)Vroot (z | w)

Since joint, marginal and conditional queries can all becomputed by two network passes, exact inference takes linear time with respect to the network size. Without loss ofgenerality, in this paper we focus our attention on the casewhere all the random variables Xi are Boolean. However,direct extension to the continuous case is not hard; one cansimply replace each univariate distribution at the leaves of
the graph with the desired continuous distribution.
2.2. SPNs as Mixture Models
It has been recently shown that any complete and decomposable SPN S over X = {X1 , . . . , Xn } is equivalent toCollapsed Variational Inference for Sum-Product Networks
+

w1

W1

W2

H1

X1

W3

H2

X2

H3

X3

···

Wm

···

Hm

···

X1

+

w3

w2

⇥

⇥

X1

⇥

X2

= w1

X2

+

+w2

⇥

X1

X2

+

+w3

⇥

X1

X2

⇥

X1

X2

Figure 2. Induced tree decomposition of a complete and decomposable SPN into a mixture model. Each component is characterized by a product of univariate distributions over X1 and X2 .
Xn
D

Figure 1. Graphical model representation of SPN S. The box
is a plate that represents replication over D training instances.
m = O(|S|) corresponds to the number of sum nodes and n is
the number of observable variables in S. Typically, mn.
W, H, X correspond to global latent variables, local latent variables and observable variables, respectively.

a Bayes net with O(n|S|) size (Zhao et al., 2015). The insight behind the construction (Fig. 1) is that each internal
sum node in S corresponds to a latent variable in the constructed Bayes net, and the Bayes net will be a bipartite
graph (excluding the nodes that correspond to global latentvariables) with one layer of local latent variables pointing
to one layer of observable variables X. An observable variable is a child of a local latent variable if and only if theobservable variable appears as a descendant of the latentvariable (sum node) in the original SPN. Equivalently, thisshows that the SPN S can be interpreted as a Bayes net
where the number of latent variables per instance is O(|S|).
Although the Bayes net perspective provides an interesting
connection between the two models, we can hardly work on
the Bayes net in practice due to the large number of latentvariables and high treewidth of the graph.
Another representation of an SPN is to decompose it into asum of induced trees (Zhao & Poupart, 2015). This decomposition works as follows:
Definition 1 (Induced tree SPN). Given a complete anddecomposable SPN S over X = {X1 , . . . , Xn }, T =(TV , TE ) is called an induced tree SPN from S if

1. Root(S) 2 TV .
2. If v 2 TV is a sum node, then exactly one child of v
in S is in TV , and the corresponding edge is in TE .
3. If v 2 TV is a product node, then all the children of v
in S are in TV , and the corresponding edges are in TE .

Here TV is the node set of T and TE is the edge set of T .

The name “induced tree SPN” is justified because it has
been shown that Def. 1 produces mixture components thatare trees whenever the original SPN is complete and decomposable. A key result based on the notion of induced trees that exhibits the distribution modeled by SPNs

is (Zhao & Poupart, 2015):
Theorem 2. Let
P⌧⌧SS =
QVroot (1 | 1). Then
Qn Vroot (x | w) canbe written asw
t=1(k,j)2TtE kj
i=1 pt (Xi = xi ),where Tt is the tth unique induced tree of S and pt (Xi )is a univariate distribution over Xi in Tt as a leaf node.
In Thm. 2, ⌧S equals the number of induced trees in S.
Fig. 2 gives an illustrative example of Thm. 2. Thm. 2 characterizes both the number of components and the form ofeach component in the mixture model, as well as their mixture weights.
2.3. Variational Bayes InferenceStandard variational Bayes inference is an optimizationbased approach to approximate the full posterior distribution of a probabilistic model (Jordan et al., 1999). It
works by constructing and maximizing an evidence lower
bound (ELBO) of the log marginal likelihood function
log p(x). Equivalently, one can minimize the KullbackLeibler (KL) divergence between the true posterior distribution and the variational posterior distribution. To review,let ⇥ = {H, W} represent the set of latent variables in aprobabilistic model (including both the local, H, and theglobal, W, latent variables) and let X represent the data.
For example, in mixture models such as LDA, H corresponds to the topic assignments of the words and W corresponds to the topic-word distribution matrix. The joint
probability distribution of ⇥ and X is p(X, ⇥ | ↵ ) where ↵is the set of hyperparameters of the model. Standard VI approximates the true posterior p(⇥ | X, ↵ ) with a variational
distribution q(⇥ | ) with a set of variational parameters. This reduces an inference problem into an optimization problem in which the objective function is given by
the ELBO:
log p(x | ↵ )b )↵)] + H[q(⇥| )] =: L(
Eq [log p(x, ⇥|↵The variational distribution q(⇥ | ) is typically assumed to
be fully factorized, with each variational parameter i governing one latent variable ✓i . The variational parameters
are then optimized to maximize the ELBO, and the optimal variational posterior will be used as a surrogate to thetrue posterior distribution p(⇥ | X, ↵ ). It is worth notingthat in standard VI the number of variational parameters
is linearly proportional to the number of latent variables,
including both the global and local latent variables.

Collapsed Variational Inference for Sum-Product Networks

3. Collapsed Variational Bayesian Inference
3.1. Motivation
Standard VI methods are computationally infeasible for
SPNs due to the large number of local latent variables foreach training instance, as shown in Fig. 1. Let W denotethe set of global latent variables (model parameters) and
Hd denote the set of local latent variables, where d indexes
over the training instances. In standard VI we need to maintain a set of variational parameters for each of the latentvariables, i.e., W and {Hd }D
d=1 . In the case of SPNs, thenumber of local latent variables is exactly the number ofinternal sum nodes in the network, which can be linearly
proportional to the size of the network, |S|. Together withthe global latent variables, this leads to a total number ofO(D|S| + |S|) variational parameters to be maintained instandard VI. As we will see in the experiments, this is prohibitively expensive for SPNs that range from tens of thousands to millions of nodes.
To deal with the expensive computation and storage in standard VI, we develop CVB-SPN, a collapsed variational algorithm, which, instead of assuming independence among
the local latent variables, models the dependence among
them in an exact fashion. More specifically, we marginalize all the local latent variables out of the joint posterior
distribution and maintain a marginal posterior distribution
over the global latent variables directly. This approach
would not be an option for many graphical models, but
as we will see later, we can take advantage of fast exact
inference in SPNs. On the other hand, we still variationally approximate the posterior distribution of the global
variables (model parameters) using a mean-field approach.
Note this basic assumption made in CVB-SPN is different from typical collapsed variational approaches developed for LDA and HDP (Teh et al., 2006; 2007), where
global latent variables are integrated out and local latentvariables are assumed to be mutually independent. As a result, CVB-SPN models the marginal posterior distribution
over global latent variables only and hence the number ofvariational parameters to be optimized is O(|S|) comparedwith O(D|S|) in the standard case.
Intuitively, this is not an unreasonable assumption to makesince compared with local latent variables, the global latent variables in Fig. 1 are further away from the influenceof observations of X and they are mutually independenta priori. Besides the computational consideration, anothermotivation to collapse out the local latent variables is thatno interpretation associated with the local latent variables
in an SPN has so far proven to be useful in practice. Unlike other mixture models such as LDA, where local latentvariables correspond to topic assignment of words, the sum
nodes in SPNs do not share typical statistical interpretations that may be useful in real-world applications.

CVB-SPN makes fewer assumptions about the independence among random variables and has fewer variationalvariables to be optimized, but to achieve these benefits, wemust overcome the following difficulties: the cost of exact marginalization over local latent variables and the nonconjugacy between the prior distribution and the likelihood
function introduced by the marginalization. The first problem is elegantly handled by the property of SPN that exact
inference over X is always tractable. In what follows weproceed to derive the CVB-SPN algorithm that efficiently
solves the second problem.
3.2. Collapsed Variational Inference
Throughout the derivation we will assume that the weights
wkj associatedwith a sum node k are locally normalized,P
i.e.,j2Ch(k) wkj = 1, 8k. This can be achieved bya bottom-up pass of the network in time O(|S|) withoutchanging the joint probability distribution over X; see Peharz et al. (2015) and Zhao et al. (2015) for more details.
For SPNs with locally normalized weights w, it can be verified that Vroot (1 | w) = 1, hence the marginal likelihood
function p(x | w) that marginalizes out local latent variables h is given by Vroot (x | w).
Since the weights associated with each sum node k are
locally normalized, we can interpret each sum node as amultinomial random variable with one possible value foreach child of the sum node. It follows that we can specify
a Dirichlet prior Dir(wk |↵k ) for each sum node k. Since
all the global latent variables are d-separated before we get
the observation x, a prior distribution over all the weights
can be factorized as:
↵) =
p(w|↵

m
Y

k=1

p(wk |↵k ) =

m
Y

k=1

Dir(wk |↵k )(3)

Based on Thm. 2, the true posterior distribution after a sequence of observations {xd }D
d=1 is:↵)p(w|{xd }D
d=1 , ↵ ) / p(w|↵
=

m
Y

k=1

Dir(wk |↵k )

⌧S
D X
Y

D
Y

d=1

Y

Vroot (xd |w)

d=1 t=1 (k,j)2TtE

wkj

n
Y

pt (xdi ) (4)i=1Intuitively, Eq. 4 indicates that the true posterior distribution of the model parameters is a mixture model where thenumber of components scales as ⌧SD and each component
is a product of m Dirichlets, one for each sum node. Here⌧S corresponds to the number of induced trees in S (cf.
Thm. 2). For discrete SPNs the leaf univariate distributions
are simply point mass distributions, i.e., indicator variables.
For continuous distributions such as Gaussians, we alsoneed to specify the priors for the parameters of those leaf

Collapsed Variational Inference for Sum-Product Networks

distributions, but the analysis goes the same as the discrete
case. Eq. 4 is intractable to compute exactly, hence weresort to collapsed variational inference. To simplify notation, we will assume that there is only one instance in thedata set, i.e., D = 1. Extension of the following derivation
to multiple training instances is straightforward. Consider
the log marginal probability over observable variables thatupper bounds the new ELBO L( ):log p(x | ↵ )↵)] + H[q(w| )] (5)Eq(w| ) [log p(x, w|↵
X↵)] + H[q(w| )]= Eq(w| ) [logp(x, h, w|↵h=: L( )
Qm
where q(w| ) = k=1 Dir(wk | k ) is the factorized variational distribution over w to approximate
the true marginalP
posterior distribution p(w|x, ↵ ) = h p(w, h|x, ↵ ). Notethat here we are using q(w) as opposed to q(w)q(h) instandard VI. We argue that L( ) gives us a better lower
bound to optimize than the one given by standard VI. Tob ) be the ELBO given by standard VI, i.e.,
see this, let L(b ) = Eq(w)q(h) [log p(x, w, h|↵↵)] + H[q(w)q(h)]L(

we have:by standard VI. This conclusion is also consistent with theone obtained in collapsed variational inference for LDA
and HDP (Teh et al., 2006; 2007) where the authors collapsed out the global latent variables instead of the local latent variables. It is straightforward to verify that↵) and L( ) is given by
the difference between log p(x|↵KL(q(w| ) k p(w|x, ↵ )), i.e., the KL-divergence between
the variational marginal posterior and the exact marginal
posterior distribution. Substituting q(w| ) into the KLdivergence objective and simplifying it, we reach the following optimization objective:↵))
minimize KL(q(w| ) k p(w|↵Eq(w| ) [log p(x|w)](7)
where the first part can be interpreted as a regularization
term that penalizes a variational posterior that is too far
away from the prior, and the second part is a data-fitting
term that requires the variational posterior to have a good
fit to the training data set. The first part in (7) can be efficiently computed due to the factorization assumption ofboth the prior and the variational posterior. However, thesecond part does not admit an analytic form because after
we marginalize out all the local latent variables, q(w| )is no longer conjugate to the likelihood function p(x|w) =Vroot (x|w). We address this problem in the next subsection.

↵)] + H[q(w)]L( ) = Eq(w) [log p(x, w|↵

3.3. Upper Bound by Logarithmic Transformationp(h, w, x | ↵ )] + H[q(w)]= Eq(w) Ep(h|w,x,↵↵) [logp(h|w, x, ↵ )The hardness of computing Eq(w| ) [log p(x|w)] makes thedirect optimization of (7) infeasible in practice. Whilep(h, w, x | ↵ )= max Eq(w) Eq(h|w) [log] + H[q(w)]nonconjugate inference with little analytic work is possiq(h|w)
q(h|w)
ble with recent innovations (Blei et al., 2012; Kingma &p(h, w, x | ↵ )max Eq(w) Eq(h) [log] + H[q(w)] Welling, 2013; Ranganath et al., 2014; Mnih & Gregor,
q(h)
q(h)
2014; Titsias & Lázaro-Gredilla, 2014; Titsias, 2015), they↵)] + H[q(w)q(h)]Eq(w)q(h) [log p(x, w, h|↵fail to make use of a key analytic property of SPNs, i.e.,easy marginalization, that allows for the optimal variationalb )= L(
distributions of local variables to be used. In this section we
The first equality holds by the definition of L( ) and theshow how to use a logarithmic transformation trick to de↵)velop an upper bound of the objective in (7) which leads to↵) = log p(h,w,x|↵
second equality holds since log p(x, w|↵↵)p(h|w,x,↵the efficient CVB-SPN algorithm.
is constant w.r.t. Ep(h|w,x,↵↵) [·]. The third equality holds
because the inner expectation is the maximum that can
The key observation is that the likelihood of w is
be achieved by the negative KL divergence between
a posynomial function (Boyd et al., 2007). To see
the approximate posterior q(h|w) and the true posterior
this,
as shown in Q
Thm. 2, we have p(x|w) =P⌧S Qnp(h|w, x, ↵ ). The inequality in the fourth line is due tow
kj
t=1(k,j)2TtE
i=1 pt (Xi = xi ). The productQ
the fact that the maximization over q(h|w) is free of connterm with respect to x, i=1 pt (Xi = xi ), is guaranteed
straint hence the optimal posterior is given by the true
to be nonnegative and can be treated as a constant w.r.t.
conditional posterior q ⇤ (h|w) = p(h|w, x, ↵ ), while thew. Hence each component in p(x|w), in terms of w,
maximization over q(h) in the fourth line needs to satis a monomial function with positive multiplicative conisfy the independence assumption made in standard VI, i.e.,
stant, and it follows that p(x|w) is a posynomial funcq(h|w) = q(h). Combining the inequality above with (5),tion of w. A natural implication of this observation is that
we havewe can do a change of variables transformation to w suchb )log p(x | ↵ ) L( ) L(
(6)
that log p(x|w) becomes a convex function in terms of thewhich shows that the ELBO given by the collapsed variatransformed variables. More specifically, let w0 = log w,
tional inference is a better lower bound than the one given
where log(·) is taken elementwise. The log likelihood, nowCollapsed Variational Inference for Sum-Product Networks

expressed in terms of w0 , is0⌧SXYlog p(x | w) = log @wkj

t=1 (k,j)2TtE

0
0
⌧S
X
= log @
exp @ct +
t=1

i=1

pt (Xi = xi )AX(k,j)2TtE

0=: log p(x | w )

n
Y

1

11

0 AA
wkj

(8)

Qn

where ct is defined as ctP= log i=1 pt (Xi = xi ). For
0each unique tree Tt , ct + (k,j)2TtE wkj
is an affine func0
tion of w . Since the log-sum-exp function is convex in its
argument, log p(x | w0 ) is convex in w0 . Such a change ofvariables trick is frequently applied in the geometric programming literature to transform a non-convex posynomial
optimization problem into a convex programming problem (Boyd et al., 2007; Chiang, 2005).
The convexity of log p(x | w0 ) helps us to develop a lower
bound of Eq(w| ) [log p(x|w)] that is efficient to compute.
Let q 0 (w0 ) be the corresponding variational distribution
over w0 induced from q(w) by the bijective transformation
between w and w0 . We haveZ
Eq(w) [log p(x | w)] = q(w) log p(x | w) dw
Z= q 0 (w0 ) log p(x | w0 ) dw0
log p(x | Eq0 (w0 ) [w0 ])where log p(x | Eq0 (w0 ) [w0 ]) means the log-likelihood ofx with theQedge weights set to be exp(Eq0 (w0 ) [w0 ])). Sincemq(w) = k=1 Dir(wk | k ) is
⇣ a product of⌘Dirichlets, we
0
can compute the weight exp Eq0 (w0 ) [wkj
] for each edge(k, j) as

Z
Z
0
0
0
Eq0 (w0 ) [wkj] = q 0 (w0 )wkj
dw0 = q 0 (wk0 )wkj
dwk0
ZX
= q(wk ) log wkj dwk = ( kj )((9)
kj 0 )
j0

where (·) is the digamma function. The equation above
then implies the new edge weight can be computed by
0
1
X
0
A
exp Eq0 (w0 ) [wkj] = exp @ ( kj )(kj 0 )
j0

⇡P

1
2

kj

j0

kj 0

1
2

(10)

where the approximation is by exp( (x)) ⇡ x 12 whenx > 1. Note that the mean of the variational
posteriorP
is given by w̄kj = Eq(w) [wkj ] = kj / j 0 kj 0 , which

Algorithm 1 CVB-SPN
Input: Initial , prior hyperparameter ↵ , training instances {xd }D
d=1 .
Output: Locally optimal ⇤ .
1: while not converged do2:Update w = exp(Eq0 (w0 | ) [w0 ]) with Eq. 10.
3:
Set r = 0.
4:
for d = 1 to D do
5:
Bottom-up evaluation of log p(xd |w).
@6:
Top-down differentiation of @w
log p(xd |w).
7:
Update r based on xd .
8:
end for↵)).
9:
Update r based on KL(q(w| ) k p(w|↵10:
Update with projected GD.
11: end while⇣⌘0
is close to exp Eq0 (w0 ) [wkj
] when kj > 1. Roughly
speaking, this shows that the lower bound we obtained for
Eq(w) [log p(x | w)] by utilizing the logarithmic transformation is trying to optimize the variational parameterssuch that the variational posterior mean has a good fit to
the training data. This is exactly what we hope for sinceat the end we need to use the variational posterior mean asa Bayesian estimator of our model parameter. Combining
all the analysis above, we formulate the following objective function to be minimized that is an upper bound of theobjective function given in (7):↵))
KL(q(w| ) k p(w|↵log p(x | Eq0 (w0 | ) [w0 ]) (11)Note that we will use the approximation given by Eq. 10
in the above optimization formulation and in Alg. 1. Thisis a non-convex optimization problem due to the digammafunction involved. Nevertheless we can still achieve a local optimum with projected gradient descent in the experiments. We summarize our algorithm, CVB-SPN, in Alg. 1.
For each instance, the gradient of the objective function in(11) with respect to can be computed by a bottom-up andtop-down pass of the network in time O(|S|). Please refer
to (Zhao & Poupart, 2015; Peharz et al., 2016) for moredetails.
3.4. Parallel and Stochastic Variants
Note that Alg. 1 is easily parallelizable by splitting training
instances in the loop (Line 4–Line 8) across threads. It
can also be extended to the stochastic setting where at eachround we sample one instance or a mini-batch of instances{xik }sk=1 (s is the size of the mini-batch) from the training
set. The stochastic formulation will be helpful when thesize of the full training data set is too large to be stored in
the main memory as a whole, or when the training instances
are streaming so that at each time we only have access to
one instance (Hoffman et al., 2013).

Collapsed Variational Inference for Sum-Product Networks
Table 1. Statistics of data sets and models. n is the number of observable random variables modeled by the network, |S| is the size of the
network and p is the number of parameters to be estimated. n ⇥ D/p means the ratio of training instances times the number of variables
to the number of parameters.
Data set
n
|S|
p
Train
Valid
Test n ⇥ D/p
NLTCS
1613,733
1,716
16,181
2,157
3,236
150.871MSNBC
1754,83924,452 291,326 38,843 58,265
202.541
KDD 2k64
48,27914,292 180,092 19,907 34,955
806.457
Plants69
132,959
58,853
17,412
2,321
3,482
20.414
Audio
100
739,525
196,103
15,000
2,000
3,000
7.649
Jester100
314,013
180,750
9,000
1,000
4,116
4.979
Netflix
100
161,655
51,601
15,0002,000
3,000
29.069Accidents
111
204,501
74,804
12,7581,700
2,551
18.931Retail135
56,931
22,11322,041
2,938
4,408
134.560
Pumsb-star
163
140,339
63,173
12,262
1,635
2,452
31.638DNA
180
108,02152,121
1,600
400
1,186
5.526
Kosarak
190
203,321
53,20433,3754,450
6,675
119.187
MSWeb
294
68,853
20,34629,441
3,270
5,000
425.423Book
500
190,625
41,122
8,700
1,159
1,739
105.783
EachMovie
500
522,753
188,387
4,524
1,002
591
12.007
WebKB
839
1,439,751
879,893
2,803
558
838
2.673
Reuters-52
889
2,210,325 1,453,390
6,532
1,028
1,540
3.995
20 Newsgrp
910 14,561,965 8,295,407
11,293
3,764
3,764
1.239
BBC
1058
1,879,921 1,222,536
1,670
225
330
1.445
Ad
1556
4,133,421 1,380,676
2,461
327
491
2.774

4. Experiments
4.1. Experimental Setting
We conduct experiments on a set of 20 benchmark data sets
to compare the performance of the proposed collapsed variational inference method with maximum likelihood estimation (Gens & Domingos, 2012). The 20 real-world data sets
used in the experiments have been widely used (Rooshenas& Lowd, 2014) to assess the modeling performance of
SPNs. All the features are binary. The 20 data sets also
cover both low dimensional and high dimensional statistical estimation, hence, they enable a thorough experimental
comparison. Detailed information about each data set aswell as the SPN models is shown in Table 1. All the SPNs
are built using LearnSPN (Gens & Domingos, 2013).
We evaluate and compare the performance of CVB-SPN
with other parameter learning algorithms in both the batchand online learning settings. For a baseline, we compare to the state-of-the-art parameter learning algorithm for
SPNs where the algorithm optimizes the training set loglikelihood directly in order to find a maximum likelihood
estimator, which we will denote as MLE-SPN. We compare
CVB-SPN and MLE-SPN in both batch and online cases.
To have a fair comparison, we apply the same optimization method, i.e., projected gradient descent, to optimize
the objective functions in CVB-SPN and MLE-SPN. CVBSPN optimizes over the variational parameters of the posterior distribution based on (11) while MLE-SPN optimizes
directly over the model parameters to maximize the loglikelihood of the training set. Since the optimization vari-ables are constrained to be positive in both CVB-SPN and
MLE-SPN, we need to project the parameters back onto thepositive orthant after every iteration. We fix the projection
margin ✏ to 0.01, i.e., w = max{w, 0.01} to avoid numerical issues. We implement both methods with backtracking
line search to automatically adjust the learning rate at each
iteration. In all experiments, the maximum number of iterations is fixed to 50 for both methods. We discard the model
parameters returned by LearnSPN and use random weights
as initial model parameters. CVB-SPN is more flexible toincorporate those model weights returned by LearnSPN as
the hyperparameters for prior distributions. In the experiments we multiply the weights returned by LearnSPN by apositive scalar and treat them as the hyperparameters of theprior Dirichlet distributions. This is slightly better than using randomly initialized priors, but the differences are negligible on most data sets. MLE-SPN can also incorporate
the model weights returned by LearnSPN by treating themas the hyperparameters of fixed prior Dirichlets. This corresponds to an MAP formulation. However in practice wefind the MAP formulation performs no better than MLESPN, and on small data sets, MAP-SPN gives consistently
worse results than MLE-SPN; so, we report only resultsfor MLE-SPN. For both MLE-SPN and CVB-SPN we use
a held-out validation set to pick the best solution during the
optimization process. We report average log-likelihood oneach data set for each method.
Both CVB-SPN and MLE-SPN are easily extended to theonline setting where training instances are coming in a
streaming fashion. In this case we also compare CVB-SPN
and MLE-SPN to an online Bayesian moment matchingCollapsed Variational Inference for Sum-Product Networks
Table 2. Average log-likelihoods on test data. Highest average log-likelihoods are highlighted in bold. " / # are used to represent
statistically better/worse results than (O)CVB-SPN respectively.
Data set
MLE-SPN CVB-SPN
OBMM OMLE-SPN OCVB-SPN
NLTCS# -6.44-6.08"-6.07#-7.78-6.12
MSNBC#-7.02-6.29-6.35#-6.94-6.34
KDD 2k
#-4.24-2.14-2.14#-27.99-2.16
Plants#-28.78-12.86"-15.14#-30.23-16.03Audio#-46.42-40.36#-40.70#-48.90-40.58Jester#-59.55-54.26-53.86#-63.67-53.84
Netflix#-64.88-60.69-57.99#-65.72-57.96
Accidents#-50.14-29.55#-42.66#-58.63-38.07Retail#-15.53-10.91#-11.42#-82.42-11.31
Pumsb-star#-80.61-25.93#-45.27#-80.19-37.05
DNA#-102.62-86.73#-99.61#-96.84-91.52
Kosarak#-47.16-10.70-11.22#-111.95-11.12MSWeb#-19.69-9.89#-11.33#-140.86-10.73
Book#-88.16-34.44#-35.55#-299.02-34.77
EachMovie#-97.15-52.63"-59.50#-284.92-64.75
WebKB#-199.15-161.46"-165.57#-413.94-169.31
Reuters-52#-218.97-85.45-108.01#-513.97-108.04
20 Newsgrp#-260.69-155.61#-158.01#-728.11-156.63BBC#-372.45-251.23#-275.43#-517.36-272.56Ad#-311.87-19.00#-63.81#-572.01-57.56

method (OBMM) (Rashwan et al., 2016) that is designed
for learning SPNs. OBMM is a purely online algorithm
in the sense that it can only process one instance in each
update in order to avoid the exponential blow-up in thenumber of mixture components. OBMM constructs an approximate posterior distribution after seeing each instanceby matching the first and second moments of the approximate posterior to the exact posterior. In this experiment,
we compare both CVB-SPN (OCVB-SPN) and MLE-SPN(OMLE-SPN) to OBMM where only one pass over the
training set is allowed for learning and at each round only
one instance is available to each of the algorithms.
4.2. Results
All experiments are run on a server with Intel Xeon CPU
E5 2.00GHz. The running time ranges from 2 min toaround 5 days depending on the size of the data set andthe size of the network. All algorithms have roughly thesame running time on the same data set as they scale linearly in the size of the training set and the size of the network. Table 2 shows the average joint log-likelihood scoresof different parameter learning algorithms on 20 data sets.
For each configuration, we use bold numbers to highlightthe best score among the offline methods and among theonline methods. We also use " / # to indicate whether
the competitor methods achieve statistically significant better/worse results than CVB-SPN on the corresponding test
data set under the Wilcoxon signed-rank test (Wilcoxon,1950) with p-value  0.05. In the batch learning experiment, CVB-SPN consistently dominates MLE-SPN on every data set with a large margin. The same results canbe observed in the online learning scenarios: both OCVBSPN and OBMM significantly outperform OMLE-SPN on
all the 20 experiments. These experiments demonstrate theeffectiveness and robustness of Bayesian inference methods in parameter estimation on large graphical models like
SPNs. In the online learning scenario, OCVB-SPN beats
OBMM on 14 out of the 20 data sets, and achieves statistically better results on 10 out of the 14. On the other hand,
OBMM obtains statistically better results than OCVB-SPN
on 4 of the data sets. This is partly due to the fact that
OCVB-SPN explicitly optimizes over an objective function that includes the evaluation of the variational posteriormean on the likelihood function, while OBMM only tries
to match the first and second moments of the distribution
after each update.

5. Conclusion
We develop a collapsed variational inference method,
CVB-SPN, to learn the parameters of SPNs. CVB-SPN
directly maintains a variational posterior distribution overthe global latent variables by marginalizing out all the local latent variables. As a result, CVB-SPN is more memory efficient than standard VI. We also show that the collapsed ELBO in CVB-SPN is a better lower bound than
the standard ELBO. We construct a logarithmic transformation trick to avoid the intractable computation of a highdimensional expectation. We conduct experiments on 20
data sets to compare the proposed algorithm with state-ofthe-art learning algorithms in both batch and online settings. The results demonstrate the effectiveness of CVBSPN in both cases.

Collapsed Variational Inference for Sum-Product Networks

Acknowledgements
HZ and GG gratefully acknowledge support from ONR
contract N000141512365. TA is partially funded by
the Netherlands NWO, project 612.001.119. BA acknowledges support from the National Science Foundation(NSF) grant number CNS-1518865, the Intel Corporation,
Google, Vodafone, NVIDIA, and the Conklin Kistler family fund.

References
Blei, David M, Ng, Andrew Y, and Jordan, Michael I. Latent Dirichlet Allocation. the Journal of machine Learning Research, 3:993–1022, 2003.
Blei, David M, Jordan, Michael I, and Paisley, John W.
Variational Bayesian Inference with Stochastic Search.
In ICML, 2012.
Boyd, Stephen, Kim, Seung-Jean, Vandenberghe, Lieven,and Hassibi, Arash. A Tutorial on Geometric Programming. Optimization and Engineering, 8(1):67–127,
2007.
Casella, George and Berger, Roger L. Statistical Inference,
volume 2. Duxbury Pacific Grove, CA, 2002.
Chiang, Mung. Geometric Programming for Communication Systems. Now Publishers Inc, 2005.
Gens, Robert and Domingos, Pedro. Discriminative Learning of Sum-Product Networks. pp. 3248–3256, 2012.
Gens, Robert and Domingos, Pedro. Learning the Structureof Sum-Product Networks. In ICML, 2013.
Hoffman, Matthew D, Blei, David M, Wang, Chong, andPaisley, John. Stochastic Variational Inference. The
Journal of Machine Learning Research, 14(1):1303–
1347, 2013.
Jordan, Michael I, Ghahramani, Zoubin, Jaakkola,
Tommi S, and Saul, Lawrence K. An Introduction to
Variational Methods for Graphical Models. Machine
Learning, 37(2):183–233, 1999.
Kingma, Diederik P and Welling, Max. Auto-encoding
Variational Bayes. arXiv preprint arXiv:1312.6114,
2013.

Peharz, Robert, Gens, Robert, Pernkopf, Franz, andDomingos, Pedro.
On the Latent Variable Interpretation in Sum-Product Networks. arXiv preprint
arXiv:1601.06180, 2016.
Poon, Hoifung and Domingos, Pedro. Sum-Product Networks: A New Deep Architecture. In Proc. 12th Conf.
on Uncertainty in Artificial Intelligence, 2011.
Ranganath, Rajesh, Gerrish, Sean, and Blei, David. BlackBox Variational Inference. In AISTATS, pp. 814–822,
2014.
Rashwan, Abdullah, Zhao, Han, and Poupart, Pascal. Online and Distributed Bayesian Moment Matching for Parameter Learning in Sum-Product Networks. In International Conference on Artificial Intelligence and Statistics, 2016.
Rooshenas, Amirmohammad and Lowd, Daniel. Learning
Sum-Product Networks with Direct and Indirect Variable
Interactions. In Proceedings of the 31st International
Conference on Machine Learning, pp. 710–718, 2014.
Teh, Yee W, Newman, David, and Welling, Max. A Collapsed Variational Bayesian Inference Algorithm for Latent Dirichlet Allocation. In Advances in neural information processing systems, 2006.
Teh, Yee W, Kurihara, Kenichi, and Welling, Max. Collapsed Variational Inference for HDP. In Advances inneural information processing systems, pp. 1481–1488,
2007.
Titsias, Michalis and Lázaro-Gredilla, Miguel. Doubly
Stochastic Variational Bayes for Non-conjugate Inference. In Proceedings of the 31st International Conference on Machine Learning, 2014.
Titsias, Michalis K. Local Expectation Gradients for Doubly Stochastic Variational Inference. arXiv preprint
arXiv:1503.01494, 2015.
Wainwright, Martin J and Jordan, Michael I. Graphical Models, Exponential Families, and Variational Inference. Foundations and Trends R in Machine Learning,
1(1-2):1–305, 2008.
Wilcoxon, Frank. Some Rapid Approximate Statistical
Procedures. Annals of the New York Academy of Sciences, pp. 808–814, 1950.

Mnih, Andriy and Gregor, Karol. Neural Variational Inference and Learning in Belief Networks. In ICML, 2014.

Zhao, Han and Poupart, Pascal. A Unified Approachfor Learning the Parameters of Sum-Product Networks.
arXiv preprint arXiv:1601.00318, 2015.

Peharz, Robert, Tschiatschek, Sebastian, Pernkopf, Franz,
and Domingos, Pedro. On Theoretical Properties of
Sum-Product Networks. In AISTATS, 2015.

Zhao, Han, Melibari, Mazen, and Poupart, Pascal. On
the Relationship between Sum-Product Networks andBayesian Networks. In ICML, 2015.