Linear Time Computation of Moments in
Sum-Product Networks
Geoff Gordon
Machine Learning Department
Carnegie Mellon University
Pittsburgh, PA 15213
ggordon@cs.cmu.edu

Han Zhao
Machine Learning Department
Carnegie Mellon University
Pittsburgh, PA 15213
han.zhao@cs.cmu.edu

Abstract
Bayesian online algorithms for Sum-Product Networks (SPNs) need to update
their posterior distribution after seeing one single additional instance. To do so,they must compute moments of the model parameters under this distribution. The
best existing method for computing such moments scales quadratically in thesize of the SPN, although it scales linearly for trees. This unfortunate scaling
makes Bayesian online algorithms prohibitively expensive, except for small or
tree-structured SPNs. We propose an optimal linear-time algorithm that works
even when the SPN is a general directed acyclic graph (DAG), which significantly
broadens the applicability of Bayesian online algorithms for SPNs. There are three
key ingredients in the design and analysis of our algorithm: 1). For each edgein the graph, we construct a linear time reduction from the moment computation
problem to a joint inference problem in SPNs. 2). Using the property that each SPN
computes a multilinear polynomial, we give an efficient procedure for polynomial
evaluation by differentiation without expanding the network that may containexponentially many monomials. 3). We propose a dynamic programming method
to further reduce the computation of the moments of all the edges in the graph from
quadratic to linear. We demonstrate the usefulness of our linear time algorithm byapplying it to develop a linear time assume density filter (ADF) for SPNs.

1

Introduction

Sum-Product Networks (SPNs) have recently attracted some interest because of their flexibility inmodeling complex distributions as well as the tractability of performing exact marginal inference [11,5, 6, 9, 16–18, 10]. They are general-purpose inference machines over which one can perform exact
joint, marginal and conditional queries in linear time in the size of the network. It has been shown that
discrete SPNs are equivalent to arithmetic circuits (ACs) [3, 8] in the sense that one can transform
each SPN into an equivalent AC and vice versa in linear time and space with respect to the network
size [13]. SPNs are also closely connected to probabilistic graphical models: by interpreting each sum
node in the network as a hidden variable and each product node as a rule encoding context-specificconditional independence [1], every SPN can be equivalently converted into a Bayesian network
where compact data structures are used to represent the local probability distributions [16]. This
relationship characterizes the probabilistic semantics encoded by the network structure and allows
practitioners to design principled and efficient parameter learning algorithms for SPNs [17, 18].
Most existing batch learning algorithms for SPNs can be straightforwardly adapted to the online
setting, where the network updates its parameters after it receives one instance at each time step.
This online learning setting makes SPNs more widely applicable in various real-world scenarios.
This includes the case where either the data set is too large to store at once, or the network needs
to adapt to the change of external data distributions. Recently Rashwan et al. [12] proposed an
31st Conference on Neural Information Processing Systems (NIPS 2017), Long Beach, CA, USA.

online Bayesian Moment Matching (BMM) algorithm to learn the probability distribution of themodel parameters of SPNs based on the method of moments. Later Jaini et al. [7] extended thisalgorithm to the continuous case where the leaf nodes in the network are assumed to be Gaussian
distributions. At a high level BMM can be understood as an instance of the general assumed density
filtering framework [14] where the algorithm finds an approximate posterior distribution within a
tractable family of distributions by the method of moments. Specifically, BMM for SPNs works
by matching the first and second order moments of the approximate tractable posterior distribution
to the exact but intractable posterior. An essential sub-routine of the above two algorithms [12, 7]
is to efficiently compute the exact first and second order moments of the one-step update posterior
distribution (cf. 3.2). Rashwan et al. [12] designed a recursive algorithm to achieve this goal in linear
time when the underlying network structure is a tree, and this algorithm is also used by Jaini et al. [7]in the continuous case. However, the algorithm only works when the underlying network structure isa tree, and a naive computation of such moments in a DAG will scale quadratically w.r.t. the network
size. Often this quadratic computation is prohibitively expensive even for SPNs with moderate sizes.
In this paper we propose a linear time (and space) algorithm that is able to compute any moments of all
the network parameters simultaneously even when the underlying network structure is a DAG. There
are three key ingredients in the design and analysis of our algorithm: 1). A linear time reduction fromthe moment computation problem to the joint inference problem in SPNs, 2). A succinct evaluation
procedure of polynomial by differentiation without expanding it, and 3). A dynamic programming
method to further reduce the quadratic computation to linear. The differential approach [3] used forpolynomial evaluation can also be applied for exact inference in Bayesian networks. This technique
has also been implicitly used in the recent development of a concave-convex procedure (CCCP) for
optimizing the weights of SPNs [18]. Essentially, by reducing the moment computation problem to a
joint inference problem in SPNs, we are able to exploit the fact that the network polynomial of an
SPN computes a multilinear function in the model parameters, so we can efficiently evaluate thispolynomial by differentiation even if the polynomial may contain exponentially many monomials,
provided that the polynomial admits a tractable circuit complexity. Dynamic programming can befurther used to trade off a constant factor in space complexity (using two additional copies of the
network) to reduce the quadratic time complexity to linear so that all the edge moments can becomputed simultaneously in two passes of the network. To demonstrate the usefulness of our linear
time sub-routine for computing moments, we apply it to design an efficient assumed density filter [14]to learn the parameters of SPNs in an online fashion. ADF runs in linear time and space due to ourefficient sub-routine. As an additional contribution, we also show that ADF and BMM can both beunderstood under a general framework of moment matching, where the only difference lies in the
moments chosen to be matched and how to match the chosen moments.

2

Preliminaries

We use [n] to abbreviate {1, 2, . . . , n}, and we reserve S to represent an SPN, and use |S| to mean
the size of an SPN, i.e., the number of edges plus the number of nodes in the graph.
2.1

Sum-Product Networks

A sum-product network S is a computational circuit over a set of random variables X ={X1 , . . . , Xn }. It is a rooted directed acyclic graph. The internal nodes of S are sums or products and the leaves are univariate distributions over Xi . In its simplest form, the leaves of S are
indicator variables IX=x , which can also be understood as categorical distributions whose entire
probability mass is on a single value. Edges from sum nodes are parameterized with positive weights.
Sum node computes a weighted sum of its children and product node computes the product of its
children. If we interpret each node in an SPN as a function of leaf nodes, then the scope of a node
in SPN is defined as the set of variables that appear in this function. More formally, for any node
v in an SPN, if v is a terminal node, say, an indicator variable over X, then scope(v) = {X}, else
scope(v) = ∪ṽ∈Ch(v) scope(ṽ). An SPN is complete iff each sum node has children with the same
scope, and is decomposable iff for every product node v, scope(vi ) ∩ scope(vj ) = ∅ for every pair(vi , vj ) of children of v. It has been shown that every valid SPN can be converted into a complete anddecomposable SPN with at most a quadratic increase in size [16] without changing the underlying
distribution. As a result, in this work we assume that all the SPNs we discuss are complete anddecomposable.
2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:= xi )
if k is a leaf node over Xi
p(XQ iif k is a product node
Vk (x; w) =(1)
j∈Ch(k) Vj (x; w)

P
w
V
(x;
w)
if
k
is
a
sum
node
j∈Ch(k) k,j jwhere Ch(k) is the child list of node k in the graph and wk,j is the edge weight associated withsum node k and its child node j. The probability of a joint assignment X = x is computed by
the value at the root of S with input x divided by a normalization constant Vroot (1; w): p(x) =Vroot (x; w)/Vroot (1; w), where Vroot (1; w) is the value of the root node when all the values of leaf
nodes are set to be 1. This essentially corresponds to marginalizing out the random vector X, which
will ensure p(x) defines a proper probability distribution. Remarkably, all queries w.r.t. x, including
joint, marginal, and conditional, can be answered in linear time in the size of the network.
2.2

Bayesian Networks and Mixture Models

We provide two alternative interpretations of SPNs that will be useful later to design our linear
time moment computation algorithm. The first one relates SPNs with Bayesian networks (BNs).
Informally, any complete and decomposable SPN S over X = {X1 , . . . , Xn } can be converted into
a bipartite BN with O(n|S|) size [16]. In this construction, each internal sum node in S corresponds
to one latent variable in the constructed BN, and each leaf distribution node corresponds to one
observable variable in the BN. Furthermore, the constructed BN will be a simple bipartite graph with
one layer of local latent variables pointing to one layer of observable variables X. An observablevariable is a child of a local latent variable if and only if the observable variable appears as adescendant of the latent variable (sum node) in the original SPN. This means that the SPN S can beunderstood as a BN where the number of latent variables per instance is O(|S|).
The second perspective is to view an SPN S as a mixture model with exponentially many mixture
components [4, 18]. More specifically, we can decompose each complete and decomposable SPN S
into a sum of induced trees, where each tree corresponds to a product of univariate distributions. To
proceed, we first formally define what we called induced trees:
Definition 1 (Induced tree SPN). Given a complete and decomposable SPN S over X ={X1 , . . . , Xn }, T = (TV , TE ) is called an induced tree SPN from S if 1). Root(S) ∈ TV ; 2).
If v ∈ 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 ∈ TV is a product node, then all the children of v in S are in TV , and the correspondingedges are in TE .
It has been shown that Def. 1 produces subgraphs of S that are trees as long as the original SPN S iscomplete and decomposable [4, 18]. One useful result based on the concept of induced trees is:
Theorem 1 ([18]). Let τS = P
Vroot (1;1). τS counts Q
the number of unique induced trees in S, andτS QnVroot (x; w) can be written as t=1
w(k,j)∈TtE k,j
i=1 pt (Xi = xi ), where Tt is the tth uniqueinduced tree of S and pt (Xi ) is a univariate distribution over Xi in Tt as a leaf node.
Thm. 1 shows that τS = Vroot (1; 1) can also be computed efficiently by setting all the edge weights
to be 1. In general counting problems are in the #P complexity class [15], and the fact that bothprobabilistic inference and counting problem are tractable in SPNs also implies that SPNs work on
subsets of distributions that have succinct/efficient circuit representation. Without loss of generality
assuming that sum layers alternate with product layers in S, we have τS = Ω(2H(S) ), where H(S)is the height of S. Hence the mixture model represented by S has number of mixture componentsthat is exponential in the height of S. Thm. 1 characterizes both the number of components and theform of each component in the mixture model, as well as their mixture weights. For the convenienceof later discussion, we call Vroot (x; w) the network polynomial of S.
Corollary 1. The network polynomial Vroot (x; w) is a multilinear function of w with positivecoefficients on each monomial.
Corollary 1 holds since each monomial corresponds to an induced tree and each edge appears at mostonce in the tree. This property will be crucial and useful in our derivation of a linear time algorithm
for moment computation in SPNs.
3

3

Linear Time Exact Moment Computation

3.1

Exact Posterior Has Exponentially Many Modes

Let m be theQ
number of sum nodes in S. Suppose we are given a fully factorized prior distributionmp0 (w; α ) = k=1 p0 (wk ; αk ) over w. It is worth pointing out the fully factorized prior distribution
is well justified by the bipartite graph structure of the equivalent BN we introduced in section 2.2. We
are interested in computing the moments of the posterior distribution after we receive one observationfrom the world. Essentially, this is the Bayesian online learning setting where we update the belief
about the distribution of model parameters as we observe data from the world sequentially. Notethat wk corresponds to the weight vector associated with sum node k, so wk is a vector that satisfies
wk > 0 and 1T wk = 1. Let us assume that the prior distribution for each wk is Dirichlet, i.e.,
P
m
m
Y
Y
Γ( j αk,j ) Y αk,j −1Q
p0 (w; α ) =Dir(wk ; αk ) =wk,jj Γ(αk,j ) j
k=1

k=1

After observing one instance x, we have the exact posterior distribution to be: p(w | x) =p0 (w; α )p(x | w)/p(x). Let Zx , p(x) and realize that the network polynomial also computes
the likelihood p(x | w). Plugging the expression for the prior distribution as well as the network
polynomial into the above Bayes formula, we have
τS Y
n
m
Y
Y
1 X
wk,j
pt (xi )p(w | x) =Dir(wk ; αk )Zx t=1i=1
k=1(k,j)∈TtE

Since Dirichlet is a conjugate distribution to the multinomial, each term in the summation is an
updated Dirichlet with a multiplicative constant. So, the above equation suggests that the exact
posterior distribution becomes a mixture of τS Dirichlets after one observation. In a data set of
D instances, the exact posterior will become a mixture of τSD components, which is intractable tomaintain since τS = Ω(2H(S) ).
The hardness of maintaining the exact posterior distribution appeals for an approximate scheme
where we can sequentially update our belief about the distribution while at the same time efficiently
maintain the approximation. Assumed density filtering [14] is such a framework: the algorithmchooses an approximate distribution from a tractable family of distributions after observing each
instance. A typical choice is to match the moments of an approximation to the exact posterior.
3.2

The Hardness of Computing Moments

In order to find an approximate distribution to match the moments of the exact posterior, we need to
be able to compute those moments under the exact posterior. This is not a problem for traditional
mixture models including mixture of Gaussians, latent Dirichlet allocation, etc., since the number of
mixture components in those models are assumed to be small constants. However, this is not the casefor SPNs, where the effective number of mixture components is τS = Ω(2H(S) ), which also depends
on the input network S.
Qn
1To
i=1 pt (xi ) and ut ,R simplify
Q the notation, for each t ∈ [τS ], we define ct ,p (w) (k,j)∈TtE wk,j dw. That is, ct corresponds to the product of leaf distributions in the tth
w 0Qinduced tree Tt , and ut is the moment of (k,j)∈TtE wk,j , i.e., the product of tree edges, under theprior distribution p0 (w). Realizing that the posterior distribution needs to satisfy the normalization
constraint, we have:
Z
τS
τS
X
Y
X
ct
p0 (w)
wk,j dw =ct ut = Zx
(2)
t=1

w

t=1

(k,j)∈TtE

Note that the prior distribution for a sum node is a Dirichlet distribution. In this case we can computea closed form expression for ut as:
Y Z
Y
Y
α
P k,j
ut =
p0 (wk )wk,j dwk =Ep0 (wk ) [wk,j ] =(3)0wk
j 0 αk,j(k,j)∈TtE

1(k,j)∈TtE(k,j)∈TtE

For ease of notation, we omit the explicit dependency of ct on the instance x .

4More generally, let f (·) be a function applied to each edge weight in an SPN. We use the notationMp (f ) to mean the moment of function f evaluated under distribution p. We are interested incomputing Mp (f ) where p = p(w | x), which we call the one-step update posterior distribution.
More specifically, for each edge weight wk,j , we would like to compute the following quantity:
Z
Z
τS
Y
1 X
ct
wk0 ,j 0 dw (4)Mp (f (wk,j )) =f (wk,j )p(w | x) dw =p0 (w)f (wk,j )Zx t=1w
w
0 0(k ,j )∈TtE

We note that (4) is not trivial to compute as it involves τS = Ω(2H(S) ) terms. Furthermore, in order
to conduct moment matching, we need to compute the above moment for each edge (k, j) from asum node. A naive computation will lead to a total time complexity Ω(|S| · 2H(S) ). A linear timealgorithm to compute these moments has been designed by Rashwan et al. [12] when the underlying
structure of S is a tree. This algorithm recursively computes the moments in a top-down fashion
along the tree. However, this algorithm breaks down when the graph is a DAG.
In what follows we will present a O(|S|) time and space algorithm that is able to compute all the
moments simultaneously for general SPNs with DAG structures. We will first show a linear time
reduction from the moment computation in (4) to a joint inference problem in S, and then proceed touse the differential trick to efficiently compute (4) for each edge in the graph. The final component
will be a dynamic program to simultaneously compute (4) for all edges wk,j in the graph by trading
constant factors of space complexity to reduce time complexity.
3.3

Linear Time Reduction from Moment Computation to Joint Inference

Let us first compute (4) for a fixed edge (k, j). Our strategy is to partition all the induced trees basedon whether they contain the tree edge (k, j) or not. Define TF = {Tt | (k, j) 6∈ Tt , t ∈ [τS ]} andTT = {Tt | (k, j) ∈ Tt , t ∈ [τS ]}. In other words, TF corresponds to the set of trees that do notcontain edge (k, j) and TT corresponds to the set of trees that contain edge (k, j). Then,Z
Y
1 X
Mp (f (wk,j )) =ct
p0 (w)f (wk,j )wk0 ,j 0 dw
Zxw
Tt ∈TT(k0 ,j 0 )∈TtE
Z
Y
1 X
+
ct
p0 (w)f (wk,j )wk0 ,j 0 dw
(5)
Zx
w
0 0
Tt ∈TF(k ,j )∈TtE

For the induced trees that contain edge (k, j), we have
Z
Y
1 X
1 X
wk0 ,j 0 dw =ct
p0 (w)f (wk,j )ct ut Mp00,k (f (wk,j ))
Zx
Zx
w
0 0
Tt ∈TT(6)

Tt ∈TT(k ,j )∈TtE

where p00,k is the one-step update posterior Dirichlet distribution for sum node k after absorbing theterm wk,j . Similarly, for the induced trees that do not contain the edge (k, j):
Z
Y
1 X
1 X
ct
p0 (w)f (wk,j )wk0 ,j 0 dw =ct ut Mp0,k (f (wk,j ))
(7)
Zx
Zx
w
0 0
Tt ∈TF

Tt ∈TF(k ,j )∈TtE

where p0,k is the prior Dirichlet distribution for sum node k. The aboveQ
equation holds by changing
the order of integration and realize that since (k, j) is not in tree Tt , (k0 ,j 0 )∈TtE wk0 ,j 0 does notcontain the term wk,j . Note that both Mp0,k (f (wk,j )) and Mp00,k (f (wk,j )) are independent ofspecific induced trees, so we can combine the above two parts to express Mp (f (wk,j )) as:
!
!
1 X
1 X
Mp (f (wk,j )) =ct ut Mp0,k (f (wk,j )) +ct ut Mp00,k (f (wk,j )) (8)Zx
Zx
Tt ∈TF

Tt ∈TT

From (2) we haveτS
1 X
ct ut = 1 andZx t=1

τS
X

ct ut =X
Tt ∈TT

t=1

5

ct ut +X
Tt ∈TF

ct ut

This implies that Mp (f ) is in fact a convex combination of Mp0,k (f ) and Mp00,k (f ). In other words,
since both Mp0,k (f ) and Mp00,k (f ) can be computed in closed form for each edge (k, j), so in order
to compute (4), we only need to be able to compute the two
Recall that for
Qcoefficients efficiently.
P
each induced tree Tt , we have the expression of ut as ut = (k,j)∈TtE αk,j / j 0 αk,j 0 . So the termPτS
t=1 ct ut can thus be expressed as:
τSX

ct ut =

t=1

τS
X

Y

t=1 (k,j)∈TtE

n
Y
α
P k,j
pt (xi )0j 0 αk,j i=1(9)

The key observation that allows us to find the linear time reduction lies in the fact that (9) sharesexactly the same functional form as the network polynomial, with the only difference being thespecification of edge weights in the network. The following lemma formalizes our argument.
PτS
Lemma 1.
t=1 ct ut can be computed in O(|S|) time and space in a bottom-up evaluation of S.
Proof. Compare the form of (9) to the network polynomial:p(x | w) = Vroot (x; w) =

τS
X

Y

wk,j

t=1 (k,j)∈TtEnYpt (xi )

(10)

i=1

Clearly (9) and (10) share the Psame functional form and the only difference lies in that the edge weightused in (9) is given by αk,j / j 0 αk,j 0 while the edge weight used in (10) is given by wk,j , both ofwhich are constrained to be positive and locally normalized. This means
P that in order to compute thevalue of (9), we can replace all the edge weights wk,j with αk,j / j 0 αk,j 0 , and a bottom-up passevaluation of S will give us the desired result at the root of the network. The linear time and spacecomplexity then follows from the linear time and space inference complexity of SPNs.

In other words, we reduce the original moment computation problem for edge (k, j) to a joint
inference problem in S with a set of weights determined by α .
3.4

Efficient Polynomial Evaluation by Differentiation
P
To evaluate (8), we also need to compute Tt ∈TT ct ut efficiently, where the sum is over a subset of
induced trees that contain edge (k, j). Again, due to the exponential lower bound on the number of
unique induced trees, a brute force computation is infeasible in the worst case. The key observation
is
PτS
that we can use the differential trick
P to solve this problem by realizing the fact that Zx = t=1 ct ut
is a multilinear function in αk,j / j 0 αk,j 0 , ∀k, j and it has a tractable circuit representation since itshares the same network structure with S.
PPτSLemma 2.
Tt ∈TT ct ut = wk,j (∂t=1 ct ut /∂wk,j ), and it can be computed in O(|S|) time andspace in a top-down differentiation of S.
Proof. Define wk,j , αk,j /X
Tt ∈TTct ut =PX

j 0 αk,j , then0Y

wk0 ,j 0nYpt (xi )Tt ∈TT (k0 ,j 0 )∈TtE

i=1X

wk0 ,j 0= wk,jY

Tt ∈TT (k0 ,j 0 )∈TtE(k0 ,j 0 )6=(k,j)= wk,jnY
i=1

pt (xi ) + 0 ·XX∂∂
ct ut +
ct ut
∂wk,j
∂wk,j
Tt ∈TT

Tt ∈TFX

ct ut

Tt ∈TF!
= wk,j

τS
∂ X
ct ut
∂wk,j t=1!

where the second equality is by Corollary 1 that the network polynomial is a multilinear functionof wk,j and the third equality holds because TF is the set of trees that do not contain wk,j . The lastequality follows by simple algebraic transformations. In summary, the above lemma holds becauseof the fact that differential operator applied to a multilinear function acts as a selector for all the
6

P
P τS
P
ct ut − Tt ∈TT ct ut canmonomials containing a specific variable. Hence, Tt ∈TF ct ut = t=1
also be computed. To show the linear time and space complexity, recall that the differentiation
w.r.t.wk,j canPτSbe efficiently computed by back-propagation in a top-down pass of S once we havecomputed t=1
ct ut in a bottom-up pass of S.
Remark. The fact that we can compute the differentiation w.r.t. wk,j using the original circuit
without expanding it underlies many recent advances in the algorithmic design of SPNs. Zhao et al.
[18, 17] used the above differential trick to design linear time collapsed variational algorithm and theconcave-convex produce for parameter estimation in SPNs. A different but related approach, where
the differential operator is taken w.r.t. input indicators, not model parameters, is applied in computing
the marginal probability in Bayesian networks and junction trees [3, 8]. We finish this discussion
by concluding that when the polynomial computed by the network is a multilinear function in termsof model parameters or input indicators (such as in SPNs), then the differential operator w.r.t. avariable can be used as an efficient way to compute the sum of the subset of monomials that containthe specific variable.
3.5

Dynamic Programming: from Quadratic to Linear

PτS
Define Dk (x; w) = ∂Vroot (x; w)/∂Vk (x; w). Then the differentiation term ∂ t=1
ct ut /∂wk,j inLemma 2 can be computed via back-propagation in a top-down pass of the network as follows:
PτS
∂ t=1
ct ut
∂Vroot (x; w) ∂Vk (x; w)== Dk (x; w)Vj (x; w)
(11)
∂wk,j
∂Vk (x; w)∂wk,jLet λk,j = (wk,j Vj (x; w)Dk (x; w)) /Vroot (x; w) and fk,j = f (wk,j ), then the final formula forcomputing the moment of edge weight wk,j under the one-step update posterior p is given by
Mp (fk,j ) = (1 − λk,j ) Mp0 (fk,j ) + λk,j Mp00 (fk,j )(12)

Corollary 2. For each edges (k, j), (8) can be computed in O(|S|) time and space.
The corollary simply follows from Lemma 1 andLemma 2 with the assumption that the moments
under the prior has closed formPsolution. By definition, we also have λk,j = Tt ∈TT ct ut /Zx ,Dk (x; w)+
hence 0 ≤ λk,j ≤ 1, ∀(k, j). This formula
wk,jshows that λk,j computes the ratio of all the
induced trees that contain edge (k, j) to the network. Roughly speaking, this measures how
important the contribution of a specific edge isVj (x; w) ×××
to the whole network polynomial. As a result,
we can interpret (12) as follows: the more important the edge is, the more portion of the moment
comes from the new observation. We visualize
our moment computation method for a single
edge (k, j) in Fig. 1.
Figure 1: The moment computation only needs
Remark. CCCP for SPNs was originally de- three quantities: the forward evaluation value atrived using a sequential convex relaxation tech- node j, the backward differentiation value node k,nique, where in each iteration a concave surro- and the weight of edge (k, j).
gate function is constructed and optimized. The key update in each iteration of CCCP ([18], (7)) is0
given as follows: wk,j∝ wk,j Vj (x; w)Dk (x; w)/Vroot (x; w), where the R.H.S. is exactly the sameas λk,j defined above. From this perspective, CCCP can also be understood as implicitly applying thedifferential trick to compute λk,j , i.e., the relative importance of edge (k, j), and then take updates
according to this importance measure.
In order to compute the moments of all the edge weights wk,j , a naive computation would scale
O(|S|2 ) because there are O(|S|) edges in the graph and from Cor. 2 each such computation takes
O(|S|) time. The key observation that allows us to further reduce the complexity to linear comesfrom the structure of λk,j : λk,j only depends on three terms, i.e., the forward evaluation value7Vj (x; w), the backward differentiation value Dk (x; w) and the original weight of the edge wk,j .
This implies that we can use dynamic programming to cache both Vj (x; w) and Dk (x; w) in abottom-up evaluation pass and a top-down differentiation pass, respectively. At a high level, we trade
off a constant factor in space complexity (using two additional copies of the network) to reduce thequadratic time complexity to linear.
Theorem 2. For all edges (k, j), (8) can be computed in O(|S|) time and space.
Proof. During the bottom-up evaluation pass, in order to compute the value Vroot (x; w) at the root ofS, we will also obtain all the values Vj (x; w) at each node j in the graph. So instead of discardingthese intermediate Vj (x; w), we cache them by allocating additional space at each node j. So after
one bottom-up evaluation pass of the network, we will also have all the Vj (x; w) for each node j, atthe cost of one additional copy of the network. Similarly, during the top-down differentiation pass ofthe network, because of the chain rule, we will also obtain all the intermediate Dk (x; w) at each nodek. Again, we cache them. Once we have both Vj (x; w) and Dk (x; w) for each edge (k, j), from(12), we can get all the moments for all the weighted edges in S simultaneously. Because the wholeprocess only requires one bottom-up evaluation pass and one top-down differentiation pass of S, thetime complexity is 2|S|. Since we use two additional copies of S, the space complexity is 3|S|. 
We summarize the linear time algorithm for moment computation in Alg. 1.
Algorithm 1 Linear Time Exact Moment Computation
Input: Prior p0 (w | α ), moment f , SPN S and input x.
Output: Mp (f (wk,jP )), ∀(k, j).
1: wk,j ← αk,j / j 0 αk,j 0 , ∀(k, j).
2: Compute Mp0 (f (wk,j )) and Mp00 (f (wk,j )), ∀(k, j).
3: Bottom-up evaluation pass of S with input x. Record Vk (x; w) at each node k.
4: Top-down differentiation pass of S with input x. Record Dk (x; w) at each node k.
5: Compute the exact moment for each (k, j): Mp (fk,j ) = (1 − λk,j ) Mp0 (fk,j ) + λk,j Mp00 (fk,j ).

4Applications in Online Moment Matching

In this section we use Alg. 1 as a sub-routine to develop a new Bayesian online learning algorithm
for SPNs based on assumed density filtering [14]. To do so, we find an approximate distribution byminimizing the KLQ
divergence between the one-step update posterior and the approximate distribution.
m
Let P = {q | q = k=1 Dir(wk ; βk )}, i.e., P is the space of product of Dirichlet densities that aredecomposable over all the sum nodes in S. Note that since p0 (w; α ) is fully decomposable, we havep0 ∈ P. One natural choice is to try to find an approximate distribution q ∈ P such that q minimizes
the KL-divergence between p(w|x) and q, i.e.,p̂ = arg min KL(p(w | x) || q)q∈P

It is not hard to show that when q is an exponential family distribution, which is the case in our
setting, the minimization problem corresponds to solving the following moment matching equation:
Ep(w|x) [T (wk )] = Eq(w) [T (wk )](13)
where T (wk ) is the vector of sufficient statistics of q(wk ). When q(·) is a Dirichlet, we have T (wk ) =log wk , where the log is understood to be taken elementwise. This principle of finding an approximate
distribution is also known as reverse information projection in the literature of information theory [2].
As a comparison, information projection corresponds to minimizing KL(q || p(w | x)) within thesame family of distributions q ∈ P. By utilizing our efficient linear time algorithm for exact moment
computation, we propose a Bayesian online learning algorithm for SPNs based on the above moment
matching principle, called assumed density filtering (ADF). The pseudocode is shown in Alg. 2.
In the ADF algorithm, for each edge wk,j the above moment matching equation amounts to solving
the following equation:Xψ(βk,j ) − ψ(βk,j 0 ) = Ep(w|x) [log wk,j ]j0

8

where ψ(·) is the digamma function. This is a system of nonlinear equations about β where the
R.H.S. of the above equation can be computed using Alg. 1 in O(|S|) time for all the edges (k, j). To
efficiently solve it, we take exp(·) at both sides of the equation and approximate the L.H.S. using thefact that exp(ψ(βk,j )) ≈ βk,j − 12 for βk,j > 1. Expanding the R.H.S. of the above equation using
the identity from (12), we have:X
exp ψ(βk,j ) − ψ(βw,j 0 ) = exp Ep(w|x) [log wk,j ]j0!(1−λk,j )!λk,jαk,j − 12
αk,j + 12
βk,j − 12
P
× P
(14)
⇔P
1 =
1
1
0
0
0
j 0 βk,j − 2
j 0 αk,j − 2
j 0 αk,j + 2
P
Note that (αk,j − 0.5)/( j 0 αk,j 0 − 0.5) is approximately the mean of the prior Dirichlet under p0
P
and (αk,j + 0.5)/( j 0 αk,j 0 + 0.5) is approximately the mean of p00 , where p00 is the posterior byadding one pseudo-count to wk,j . So (14) is essentially finding a posterior with hyperparameter βsuch that the posterior mean is approximately the weighted geometric mean of the means given by p0
and p00 , weighted by λk,j .
Instead of matching the moments given by the sufficient statistics, also known as the natural moments,
BMM tries to find an approximate distribution q by matching the first order moments, i.e., the meanof the prior and the one-step update posterior. Using the same notation, we want q to match thefollowing equation:
βk,j
αk,j
αk,j + 1
Eq(w) [wk ] = Ep(w|x) [wk ] ⇔ P= (1 − λk,j ) P+ λk,j P
(15)
0
0
0
β
α
j 0 k,jj 0 k,jj 0 αk,j + 1
Again, we can interpret the above equation as to find the posterior hyperparameter β such that theposterior mean is given by the weighted arithmetic mean of the means given by p0 and p00 , weighted
by λk,j . Notice that due to the normalization constraint, we cannot solve for β directly from theabove equations, and in order to solve for β we will need one more equation to be added into thesystem. However, from line 1 of Alg. 1, what we need in the next iteration of the algorithm is notβ, but only its normalized version. So we can get rid of the additional equation and use (15) as the
update formula directly in our algorithm.
Using Alg. 1 as a sub-routine, both ADF and BMM enjoy linear running time, sharing the same orderof time complexity as CCCP. However, since CCCP directly optimizes over the data log-likelihood,in practice we observe that CCCP often outperforms ADF and BMM in log-likelihood scores.
Algorithm 2 Assumed Density Filtering for SPN
Input: Prior p0 (w | α ), SPN S and input {xi }∞i=1 .
1: p(w) ← p0 (w | α )2: for i = 1, . . . , ∞ do3:
Apply Alg. 1 to compute Ep(w|xi ) [log wk,j ] for all edges (k, j).
4:
Find p̂ = arg minq∈P KL(p(w | xi ) || q) by solving the moment matching equation (13).
5:p(w) ← p̂(w).
6: end for5

Conclusion

We propose an optimal linear time algorithm to efficiently compute the moments of model parameters
in SPNs under online settings. The key techniques used in the design of our algorithm include
the liner time reduction from moment computation to joint inference, the differential trick that isable to efficiently evaluate a multilinear function, and the dynamic programming to further reduce
redundant computations. Using the proposed algorithm as a sub-routine, we are able to improve thetime complexity of BMM from quadratic to linear on general SPNs with DAG structures. We also use
the proposed algorithm as a sub-routine to design a new online algorithm, ADF. As a future direction,
we hope to apply the proposed moment computation algorithm in the design of efficient structurelearning algorithms for SPNs. We also expect that the analysis techniques we develop might findother uses for learning SPNs.
9

Acknowledgements
HZ thanks Pascal Poupart for providing insightful comments. HZ and GG are supported in part by
ONR award N000141512365.

References[1] C. Boutilier, N. Friedman, M. Goldszmidt, and D. Koller. Context-specific independence inBayesian networks. In Proceedings of the Twelfth international conference on Uncertainty inartificial intelligence, pages 115–123. Morgan Kaufmann Publishers Inc., 1996.
[2] I. Csiszár and F. Matus. Information projections revisited. IEEE Transactions on Information
Theory, 49(6):1474–1490, 2003.
[3] A. Darwiche. A differential approach to inference in Bayesian networks. Journal of the ACM(JACM), 50(3):280–305, 2003.
[4] A. Dennis and D. Ventura. Greedy structure search for sum-product networks. In International
Joint Conference on Artificial Intelligence, volume 24, 2015.
[5] R. Gens and P. Domingos. Discriminative learning of sum-product networks. In Advances in
Neural Information Processing Systems, pages 3248–3256, 2012.
[6] R. Gens and P. Domingos. Learning the structure of sum-product networks. In Proceedings of
The 30th International Conference on Machine Learning, pages 873–880, 2013.
[7] P. Jaini, A. Rashwan, H. Zhao, Y. Liu, E. Banijamali, Z. Chen, and P. Poupart. Online algorithms
for sum-product networks with continuous variables. In Proceedings of the Eighth International
Conference on Probabilistic Graphical Models, pages 228–239, 2016.
[8] J. D. Park and A. Darwiche. A differential semantics for jointree algorithms. Artificial
Intelligence, 156(2):197–216, 2004.
[9] R. Peharz, S. Tschiatschek, F. Pernkopf, and P. Domingos. On theoretical properties of sumproduct networks. In AISTATS, 2015.
[10] R. Peharz, R. Gens, F. Pernkopf, and P. Domingos. On the latent variable interpretation insum-product networks. IEEE Transactions on Pattern Analysis and Machine Intelligence, 39(10):2030–2044, 2017.
[11] H. Poon and P. Domingos. Sum-product networks: A new deep architecture. In Proc. 12th Conf.
on Uncertainty in Artificial Intelligence, pages 2551–2558, 2011.
[12] A. Rashwan, H. Zhao, and P. Poupart. Online and distributed bayesian moment matchingfor parameter learning in sum-product networks. In Proceedings of the 19th International
Conference on Artificial Intelligence and Statistics, pages 1469–1477, 2016.
[13] A. Rooshenas and D. Lowd. Learning sum-product networks with direct and indirect variableinteractions. In ICML, 2014.
[14] H. W. Sorenson and A. R. Stubberud. Non-linear filtering by approximation of the a posteriori
density. International Journal of Control, 8(1):33–51, 1968.
[15] L. G. Valiant. The complexity of computing the permanent. Theoretical Computer Science, 8(2):189–201, 1979.
[16] H. Zhao, M. Melibari, and P. Poupart. On the relationship between sum-product networks andbayesian networks. In ICML, 2015.
[17] H. Zhao, T. Adel, G. Gordon, and B. Amos. Collapsed variational inference for sum-productnetworks. In ICML, 2016.
[18] H. Zhao, P. Poupart, and G. Gordon. A unified approach for learning the parameters of
sum-product networks. NIPS, 2016.
10