Learning General Latent-Variable Graphical Models with
Predictive Belief Propagation and Hilbert Space Embeddings

arXiv:1712.02046v1 [cs.LG] 6 Dec 2017

Borui Wang
Stanford University

Abstract
In this paper, we propose a new algorithm for
learning general latent-variable probabilistic graphical models using the techniques of
predictive state representation, instrumental
variable regression, and reproducing-kernel
Hilbert space embeddings of distributions.
Under this new learning framework, we first
convert latent-variable graphical models into
corresponding latent-variable junction trees,
and then reduce the hard parameter learning
problem into a pipeline of supervised learning problems, whose results will then be used
to perform predictive belief propagation over
the latent junction tree during the actual inference procedure. We then give proofs of our
algorithm’s correctness, and demonstrate its
good performance in experiments on one synthetic dataset and two real-world tasks from
computational biology and computer vision
— classifying DNA splice junctions and recognizing human actions in videos.

1

Introduction

Many important problems in modern machine learning
and artificial intelligence can be gracefully captured by
probabilistic graphical models with latent variables.
However, learning such general latent-variable probabilistic graphical models has remained very difficult,
largely due to the intractability introduced by the
latent variables and the possibly complicated loopy
structures of the graphical models.
Previously, most learning algorithms for latentvariable graphical models would resort to local search
heuristics such as the Expectation Maximization algorithm (Dempster et al., 1977) that turn the parameter
learning problem into a maximization problem over a
non-convex objective function. However, such methods based on iterative local search can lead to bad local
optima and can be very slow to converge.

Geoffrey J. Gordon
Carnegie Mellon University

In order to overcome these shortcomings, recently researchers have been proposing and studying a new
family of non-parametric learning algorithms called
spectral learning that are based on the idea of method
of moments (Hsu et al., 2009; Anandkumar et al.,
2014). These spectral learning algorithms seek to learn
alternative parameterizations of the latent-variable
graphical models that are purely based on the observable variables and can be efficiently recovered from
their low-order observable moments through tensor decompositions. Therefore these algorithms enjoy the
advantages of being local-minimum-free and provably
consistent, and also can often be orders of magnitude
faster than the search-based methods such as the EM
algorithm.
Although the spectral learning algorithms have recently been successfully applied to tackle the learning
problems of more and more types of latent-variable
graphical models, there are still several severe limitations and problems to these spectral algorithms that
need to be solved. First of all, most of the spectral
learning algorithms that have been proposed so far
are only restricted to learning one or several specific
types of latent-variable graphical models (the majority of them are tree-structured), thus making themselves hard to generalize to other types of latentvariable graphical model structures (possibly beyond
trees). Second, most of the current spectral learning
algorithms can only deal with discrete-valued random
variables with moderate cardinalities due to their parameterization assumptions, thus cannot be easily applied to scenarios where the variables in the graphical models are continuous-valued or have large cardinalities, which are quite often encountered in many
important real-world applications. Third, the current spectral algorithms are generally idiosyncratic to
the specific graphical model structures that they are
targeted to learn, and thus cannot provide a flexible
learning framework or template for users to easily incorporate different prior knowledge and probabilistic
assumptions when facing different learning problems.
In order to solve these problems and to overcome these

Learning General Latent-Variable Graphical Models with Predictive Belief Propagation and HSE

2

Predictive State Representations
and Predictive Belief Propagation

In its original formulation, a predictive state representation (PSR) (Singh et al., 2004) is a compact
way to model the state of a controlled dynamical system using a vector consisting of predictions over a set
of observable future experiments (or tests) that can
be performed on the system. Under this representation, the state of the dynamical system can be compactly summarized as a K-dimensional vector of the
predictions over a finite set of K tests, which is called
the “core set of tests”, and the success probabilities of
all other future tests are simply determined by linear
functions of this state vector.
A key insight is that the idea of predictive state representations can actually be carried over to characterize the information flows of message passing and
belief propagation in probabilistic graphical models.
In essence, the state of a system is a time bottleneck that compactly summarizes everything we need
to know about the past in order to predict the future,
as illustrated in Figure 1(a) below.
When we perform belief propagation during inference, at each separator set, the classical approach is
to propagate a message containing the partial sumproduct results from one part of the junction tree to

Intuition: bottleneck
Past

Future
Belief Propagation

predict
separator set

data about past
(many samples)

data about future
(many samples)
state

limitations, in this paper we will develop a new algorithmic framework for learning general latent-variable
probabilistic graphical models based on the ideas of
predictive belief propagation and reproducing-kernel
Hilbert space embeddings of distributions. Our new
learning framework is built on the latent junction tree
graphical representation of latent-variable probabilistic graphical models, and reduces the hard parameter learning problem down to a pipeline of supervised
learning problems named two-stage regression, where
we can flexibly insert any supervised learning algorithms into the procedure as we see fit. When dealing with continuous variables, we will embed their distributions as elements in reproducing-kernel Hilbert
spaces, and then use the kernel trick to propagate and
integrate these embeddings over the latent junction
tree via tensor operations in a non-parametric fashion,
thus seamlessly extending our algorithm from the discrete domain to the continuous domain under the same
general framework. As a result, our learning framework is general enough to apply to all different structures of latent-variable graphical models (both trees
and loopy graphs), can handle both discrete variables
and continuous variables, and can allow users to easily incorporate prior information and probabilistic assumptions during learning.

compress

“Outside Tree”

expand

“Inside Tree”

bottleneck

Can find best rank-k bottleneck via matrix
factorization spectral method

(a)

Geoff Gordon—ICML tutorial—July, 2012

(b)
3

for regularization, don’t just want a state: want a low-d (hopefully minimal) one
so, select k (linear combinations of) predictions which explain as much variance in future as possible

Figure 1: (a) The classical interpretation of state as
a compact time bottleneck in dynamical systems. (b)
The analogy of predictive state in belief propagation
over latent-variable graphical models (and their corresponding latent junction trees). Blue nodes indicate
observable variables and white nodes indicate latent
variables.
the other across the separator set in order to compute
the posterior distribution of downstream variables using the evidence information from the upstream variables. Therefore, when the belief propagation agent
visits each separator set, the message that it is sending over can be essentially viewed as a representation
of the current “state” of the junction tree inference
system at that paticular separator set. This important analogy is illustrated in Figure 1(b) below.
Once this analogy of state is built for junctions trees,
we can further generalize the idea of PSR onto junction
trees. Now instead of using the partial result of sumproduct calculations over the “past” part (i.e., outside tree) of the junction tree as the state as well as
the message to send across the current separator set
(as were done by the classical Shafer-Shenoy algorithm
(Shenoy & Shafer, 1990) and Hugin algorithm (Lauritzen & Spiegelhalter, 1988; Anderson et al., 1989)),
we can use the predictions of the success probabilities
of a core set of future tests consisting of observations
in the “future” part (i.e., inside tree1 ) of the junction
tree to represent the state of the junction tree at the
current separator set, and send this prediction vector as the message into the inside tree during belief
propagation (see Figure 1(b)). More concretely, at a
separator set in a junction tree, among all the observable variables in the inside tree, there always exists a
minimal subset of observable variables whose posterior
joint probabilistic distributions given the evidence information passed in from the outside tree would completely determine the posterior joint distribution of all
the observable variables in the inside tree. This minimal subset of observable variables in the inside tree
would be designated as the core set of variables
1
The inside tree T (S) of a separator set S is defined
to be the subtree rooted at S (including S) in the juntion
tree; the outside tree is defined to be rest of the junction
tree excluding the inside tree T (S).

Borui Wang, Geoffrey J. Gordon

for this particular separator set, and thus all the possible joint value realizations of these variables can be
defined as the core set of tests for a PSR over the latent junction tree. Therefore, the predictive state at
a separator set can be defined as the posterior joint
probability distribution of its core set of variables (or
equivalently a sufficient statistic for this joint distribution).
Now during the inference, instead of passing partial
sum-product results around as the messages, we can
pass the predictive state vectors around as the new
form of predictive messages to calculate the desired
posterior. And during the learning process, we can use
two-stage instrumental regression to learn the operators at each separator set that process incoming predictive messages and distribute outgoing predictive messages, purely based on the observed quantities in the
training dataset. We name this novel message passing
framework predictive belief propagation.

…

…

Cp

β(S)

…

S

outside tree T\T(S)

C

α(SK)

α(S2)

…

k=1

cannot be directly observed in the training data. However, this estimation of W S that we would obtain from

K

can overlap with the variables in ∪ α(Sk ) and thus
k=1

K

the noise term on ⊗ θ

α(S)
T(S2)

inside tree T(S)

Figure 2: Illustration of a typical non-leaf separator
set S in a latent junction tree

3

K

⊗ E[θSk (α(Sk )) | E] and E[θS (α(S)) | E] themselves

θS (α(S)) may be biased, because the variables in α(S)
T(SK)

…

T(S1)

K

pirical samples of ⊗ θSk (α(Sk )) and θS (α(S)), since

k=1

…

…

Therefore, the goal of our new learning algorithm under this PSR representation is mainly to learn these
linear operators W S for each of the non-leaf separator sets S in the latent junction tree from our training data. A first naive attempt for this would be
to perform a linear regression directly using the em-

K

S2

α(S1)

for any outside tree evidence E.

linear regression over samples of ⊗ θSk (α(Sk )) and

SK

S1

K

⊗ E[θSk (α(Sk )) | E] = W S ×S E[θS (α(S)) | E]

k=1

k=1

Root

…

dence information that we observed in the outside tree
T \ T (S) of S, then we have that the outer product of E[θS1 (α(S1 )) | E], E[θS2 (α(S2 )) | E], ... ,
E[θSK (α(SK )) | E] is completely determined by a linear function of E[θS (α(S)) | E], for any E, because
α(S1 ), α(S2 ), ... , α(SK ) are all sets of variables that
are contained in the inside tree below S. That is to
say, there exists a linear operator W S such that2 :

Two-Stage Regression

After determining the form of messages that we are
passing around during the inference process, we now
need to figure out how different messages relate to each
other and how we could learn these inter-message relationships from the training data. Without loss of
generality, let’s consider a non-leaf separator set S in
a latent junction tree T , where S is connected with K
“childen” separator sets {S1 , S2 , ..., SK } below it (see
Figure 2). Let α(S) indicate the core set of observable variables that is associated with a separator set
S. Then according to the definition of PSR on latent junction trees, knowing the posterior joint probability distribution over α(S) would completely determine the posterior joint probability distribution over
all the observable variables in the inside tree T (S)
below S. Let θS [α(S)] denote a sufficient statistics
feature vector for the posterior joint distribution over
α(S) of a separator set S, and let E denote the evi-

k=1

Sk

(α(Sk )) and θS (α(S)) can

be correlated with each other. One powerful method
to overcome this problem of correlated noise in modern statistics and econometrics is instrumental variable regression (Stock & Watson, 2011; Hefny et al.,
2015). In our current setting, a valid instrumental variable must satisfy the criteria that it is correlated with
the input variable θS (α(S)), and that it is not correK

lated with the noise term on ⊗ θSk (α(Sk )). Thus an
k=1

ideal choice of instrumental variable that satisfies the
above criteria would be some features of the observable variables that are contained in the outside tree
T \T (S), because these observable variables in the outside tree are strongly correlated with θS (α(S)) but also
K

do not overlap with the variables in ∪ α(Sk ). Now let
k=1

β(S) indicate the set of observable variables in the outside tree T \ T (S) whose features η S (β(S)) we decide
to use as our instrumental variable. Then we can efficiently obtain an unbiased estimation of the linear operator W S in three straightforward steps: (1) regress
K

θS (α(S)) upon η S (β(S)); (2) regress ⊗ θSk (α(Sk ))
k=1

2
The mathematical notation ×S indicates the modespecific tensor multiplication along mode S.

Learning General Latent-Variable Graphical Models with Predictive Belief Propagation and HSE

upon η S (β(S)); (3) run a final linear regression from
the predictions obtained in step (1) to the predictions
obtained in step (2) to recover an unbiased estimation
of the value of W S . These three steps of supervised
learning constitute the two-stage regression procedure in our new learning algorithm.

Root
A

A

G

A

H

B

J

D

B

Main Algorithm

In this section, we present our main algorithm3 for
learning general latent-variable graphical models in
two consecutive parts: the learning algorithm followed
by the inference algorithm. Here we will describe the
basic version of the algorithm for the case where all the
variables are discrete-valued, and the extension to handle models with continuous-valued variables is given in
section 5. We give the full proof of the consistency of
our algorithm in the appendix.
At a high level, the desired input-output behavior of
our algorithm can be formulated as:
Input: the topology of the latent-variable probabilistic graphical model G, a training set of N i.i.d.
samples {xd1 , ..., xd|O| }N
d=1 , a set of observed evidence
{Xi = xi }i∈E (where E denotes the set of index for
the set of observable variables that are observed as
evidence), and the index Q of the query node XQ .
Output: Estimated posterior probability distribution
of the query node conditioned on the observations of
b Q | {Xi = xi }i∈E ].
the set of evidence nodes: P[X
4.1

The Learning Algorithm:

Model Construction: Convert the latent-variable
graphical model G into an appropriate corresponding
latent-variable junction tree T such that each observable variable in the junction tree T can be associated
with one leaf clique in T . (See Figure 3(a) and 3(b)
for an example.)
Model Specification: For each separator set S in
the junction tree T , among all the observable variables
that are associated with its inside tree T (S), determine a minimal subset α(S) of observable variables4
3
Some mathematical notations that we will use here are:
O denotes the set of observed variables {X1 , ..., X|O| }; H
denotes the set of hidden variables {X|O|+1 , ..., X|O|+|H | };
N( ) denotes the number of elements in a set; N ( ) denotes the number of possible realization values of a discrete
random variable; l( ) denotes the length of a vector; I(a)
denotes the order index of a certain value a among all the
possible realization values of a discrete random variable A;
γ(C) denotes the set of all separator sets that are connected
with a clique node C in a junction tree.
4
For each separator set Sl that is right above a leaf
clique node Cl , we require that α(Sl ) must be the set of
observable variables that are associated with Cl .

B

D

D

B

E

A

D

F

D

G

A

H

C
D

4

D

C

A
B

F

B

B

B

E

D

F

F

J

(a)

(b)

Figure 3: (a) An example latent-variable graphical
model. Blue nodes indicate observable variables and
white nodes indicate latent variables. (b) A corresponding latent junction tree that the graphical model
in (a) is converted to. Pink squares indicate clique
nodes and yellow squares indicate separator sets. Variable nodes with red circles around them are associated
with their current leaf clique nodes.

whose posterior joint probabilistic distributions given
the evidence information passed in from the outside
tree T \ T (S) would completely determine the posterior joint distribution of all the observable variables in
the T (S). This is the core set of variables for S, and
will be denoted as α(S) = {A1 , A2 , ..., AN[α(S)] }. Then
among all the observable variables that are associated
with the outside tree T \ T (S), select a subset of variables β(S) = {B1 , B2 , ..., BN[β(S)] }. Now pick features
for the inside tree φS = θS [α(S)] and the outside tree
ψ S = η S [β(S)], where we require that φS must be a
sufficient statistic for P[α(S) | β(S)]. (See Figure 4 for
an example.)
Stage 1A Regression (S1A): At each non-leaf separator set S in the latent junction tree T , learn a (possibly non-linear) regression model to estimate φ¯S =
E[φS | ψ S ]. The training data for this regression model
n
oN
are
θS (α(S)d ), η S (β(S)d )
across all N i.i.d.
d=1
training samples.
Stage 1B Regression (S1B): At each non-leaf separator set S in the latent junction tree T , where
S is connected with K “childen” separator sets
{S1 , S2 , ..., SK } below it (see Figure 3), learn a (possiK ¯
bly non-linear) regression model to estimate ⊗ φSk =
k=1

K

E[ ⊗ φ

Sk

S

| ψ ]. The training data for this regression
n K
oN
model are
⊗ θSk (α(Sk )d ), η S (β(S)d )
across
k=1

k=1

d=1

all N i.i.d. training samples.
Stage 2 Regression (S2): At each non-leaf separator set S in the latent junction tree T , use the feature
expectations estimated in S1A and S1B to train a lin-

Borui Wang, Geoffrey J. Gordon
K ¯
ear regression model to predict ⊗ φSk = W S ×S φ¯S ,
k=1

where W S is a linear operator associated with S.
The training data for this model are estimates of

K ¯
⊗ φSk , φ¯S for all the training samples that we obk=1

tained from S1A and S1B regressions.

(
~1, if X is not observed
ζ(X) =
ex , if X is observed to have value x
Then the upward message that we send from C to its
parent clique node Cp can be calculated as:


†
mC→Cp = ΦS ×{δ(C)}
⊗ ζ(X)
X∈δ(C)

Root Covariance Estimation: At the root clique
node Cr of the latent junction tree T , estimate the
expectation of the outer product of the inside tree feature vectors of all adjecent separator sets that are connected with Cr by taking average across all the N i.i.d.
training samples:
N
P
⊗ θS (α(S)d )
S∈γ(C
r)
d=1
b
.
T Cr = E[
⊗ θS (α(S))] =
N
S∈γ(Cr )

Message Passing:
 From leaf to root (the upward phase):
For each non-root parent clique node C in T where
it is separated by a separator set S from its parent
node Cp , once it has received the messages from all of
its children nodes C1 , C2 , ..., CK , which are separated
from it by the separator sets S1 , S2 , ..., SK repectively,
compute
the upward message:
m
C→Cp

Root
A

S

A

Outside Tree

C

= WS

AB
A

Inside Tree

B

B

B

B

D

B

E

D

D

A

(S AB )
B

D

D

F

D

G

A

H

F

F

J

↵(S AB )

Figure 4: Illustration of the inside tree, outside tree,
α set and β set for an example separator set S AB .
4.2

×Sj
mCj →C
j∈{1,2,...,K}

and then send this message to its parent node Cp .
B

D

= W S ×S1 mC1 →C ×S2 mC2 →C ×S3 ... ×SK mCK →C

The Inference Algorithm:

Design Tensor Construction: For each leaf separator set Sl of the latent junction tree T , construct an
inside feature design tensor with modes φSl , A1 , A2 ,
..., AN[α(Sl )] , and dimensions l(φSl )×N (A1 )×N (A2 )×
... × N (AN[α(Sl )] ). Each fiber of ΦSl along mode φSl
would be the realization vector of the inside feature
function θSl at the corresponding values for variables
A1 , A2 , ..., AN[α(Sl )] , i.e.,

ΦSl :, I(a1 ), I(a2 ), ..., I(aN[α(Sl )] )
= θSl (a1 , a2 , ..., aN[α(Sl )] )
Initial Leaf Message Generation: At each leaf
clique node C of the latent junction tree T , let δ(C)
indicate the set of observable variables that are contained in C, and let S indicate the separator set right
above C. We define a function ζ(X) of an observable
variable X that evaluates to an all-one vector if X is
not observed in the evidence and evaluates to a onehot-encoding value indicator vector ex if X is observed
to have value x in the evidence. That is to say:

 At the Root:
At the root clique node Cr in T where it is surrounded
by all of its children node C1 , ..., CK (each separated
by the separator sets S1 , ..., SK respectively), for each
k ∈ {1, 2, ..., K}, once Cr has received all the K − 1
upward messages from all of its other K − 1 children
cliques except for Ck , compute the downward message:
mCr →Ck
= T Cr ×S1 mC1 →Cr ×S2 ... ×Sk−1 mCk−1 →Cr
×Sk+1 mCk+1 →Cr ×Sk+2 ... ×SK mCK →Cr
= T Cr

×Sj
mCj →Cr
j∈{1,2,...,K}\k

and then send this message to its childen node Ck .
 From root to leaf (the downward phase):
For each non-root parent clique node C in T where it
is separated from its parent node Cp by a separator set
S and separated from its children nodes C1 , C2 , ..., CK
by the separator sets S1 , S2 , ..., SK respectively, once
it has received the downward message mCp →C from
its parent node Cp , for each k ∈ {1, 2, ..., K}, compute
the downward message:
mC→Ck
= W S ×S mCp →C ×S1 mC1 →C ×S2 ... ×Sk−1 mCk−1 →C
×Sk+1 mCk+1 →C ×Sk+2 ... ×SK mCK →C
= W S ×S mCp →C

×Sj
mCj →C
j∈{1,2,...,K}\k

and then send this message to its children node Ck .

Learning General Latent-Variable Graphical Models with Predictive Belief Propagation and HSE

Computing Query Result: For the query node XQ ,
locate the leaf clique node that XQ is associated with
and call it CQ . Call CQ ’s parent node Cp and the separator set between them SQ . We first use the design ten†

sor ΦSQ and its Moore-Penrose psuedoinverse ΦSQ to
transform the downward incoming message mCp →CQ
and the upward outgoing message mCQ →Cp , respectively, and then compute the Hadamard product of
these transformed versions of the two messages to obtain an estimate of the unnormalized conditional probability of δ(CQ ) given all the evidence {Xi = xi }i∈E :
b
P[δ(C
Q ) | {Xi = xi }i∈E ] ∝

† 
mCp →CQ ×SQ ΦSQ ◦ ΦSQ ×SQ mCQ →Cp
We can then marginalize out the variables in
δ(CQ )\XQ and renormalize to obtain the final query
result - the estimate of the conditional probability distribution of the query variable XQ given all the evib Q | {Xi = xi }i∈E ].
dence: P[X
[Additional Note]: Another important type of queries
that is very commonly encountered in practice and
that we can also easily compute here is the joint probb
ability of all the observed evidence: P[{X
i = xi }i∈E ].
It turns out that in the previous step, before marginalization and renormalization, the Hadamard product is
actually equal to:

† 
mCp →CQ ×SQ ΦSQ ◦ ΦSQ ×SQ mCQ →Cp
b
= P[δ(C
Q ), {Xi = xi }i∈E ]
(See the appendix for the proof.) Now if we marginalize out all the variables in δ(CQ ) from the above probb
ability table, we can easily obtain P[{X
i = xi }i∈E ].
4.3

Proof of Consistency

The learning and inference algorithm that we presented above can be proved to be a consistent estimator of the inquired ground truth conditional probability distribution P[XQ | {Xi = xi }i∈E ]. Here we list
the four key results in the derivation of the proof, and
we present the actual full proof in the appendix.
Theorem 1. At each non-leaf separator set S, the
linear operator W S that is obtained from the Stage 2
Regression converges in probability to:
K
plim W S = P[ ∪ α(Si ) | α(S)] ×{α(S1 )} ΦS1 ×{α(S2 )}

And then the next two theorems directly follows
through mathematical induction:
Theorem 3. During the upward phase of the inference
procudure, all the upward messages that each non-leaf
non-root clique node C sends to its parent node Cp
(which is separated from C by the separator set S) converges in probability to5 :
plim mC→Cp

n→∞

= P[{Xi = xi }Xi ∈ρ(C) | α(S)] ×{α(S)} ΦS

†

Theorem 4. During the downward phase of the inference procudure, all the downward messages that each
non-leaf clique node C sends to any of its children node
Ci (which is separated from C by the separator set Si )
converges in probability to6 :
plim mC→Ci

n→∞

= P[{Xj = xj }Xj ∈µ(Ci ) , α(Si )] ×{α(Si )} ΦSi

5

Extending to Continuous Domain
through Hilbert Space Embeddings
of Distributions

One of the biggest advantages of our new learning
framework is that it can be seamlessly entended from
the discrete domain to the continuous domain in a
very general non-parametric fashion. Previously all
the existing algorithms for learning graphical models
with continuous variables would have to make certain
parametric assumptions (such as multivariate Gaussian) about the function forms that they are using
to model or approximate the continuous distributions
of variables. However, in many real-world data, the
underlying continuous distributions are highly complex and irregular, and such parametric approach
would severely limit the modeling expressiveness of
the graphical models and can often deviate from the
true underlying distributions. In contrast, in our algorithm, we can simply use the reproducing-kernel
Hilbert space (RKHS) embeddings of distributions as
sufficient statistic features and express all the (alternative) parameter learning and belief propagation operations as tensor algebra in the infinite-dimensional
Hilbert space, and then employ the kernel trick to

i=1

n→∞

ΦS2 ×{α(S3 )} ... ×{α(SK )} ΦSK ×{α(S)} ΦS

†

Theorem 2. At the root clique of the latent junction tree, the tensor T Cr that is obtained from the root
covariance estimation step converges in probability to:
plim T Cr = P[ ∪ α(S)] ×{α(S)} ΦS
n→∞

S∈γ(Cr )

S∈γ(Cr )

5
Here ρ(C) denotes the set of all observable variables
that are associated with the leaf clique nodes of the inside
tree T (C) that are actually being observed in the evidence.
6

Here µ(C) denotes the set of all observable variables
that are associated with the leaf clique nodes of the outside tree T \ T (C) that are actually being observed in the
evidence.

Borui Wang, Geoffrey J. Gordon

transform these operations back into tractable finitedimensional linear algebra calculations over Gram matrices. We now explain how to formulate our learning
algorithm over continuous domain with RKHS embeddings7 .
5.1

RKHS Embeddings of Continuous
Distributions

First, for each observable variables X in the latent
graphical model, specify a characteristic kernel:
def
KX (x, x0 ) = hφX (x), φX (x0 )iF

and pick the features for the inside tree to be φS =
φα(S) (α(S)), and the features for the outside tree to
be ψ S = φβ(S) (β(S)). And for each separator set Sl
that is right above a leaf clique node, we require that
φSl = φα(Sl ) (α(Sl )) =
⊗
φXj (xj ).
Xj ∈α(Sl )

n
oN
K
Now given the N i.i.d. α(S)d , β(S)d , ∪ α(Si )d

d=1

tuples in our training dataset, we define the following
linear operators:


Υα(S) = φα(S) (α(S)1 ), ... , φα(S) (α(S)N )


Υβ(S) = φβ(S) (β(S)1 ), ... , φβ(S) (β(S)N )
 K

K
Υext(S) = ⊗ φα(Sk ) (α(Sk )1 )... ⊗ φα(Sk ) (α(Sk )N )
k=1

k=1

Let Υ∗ denote the adjoint (conjugate transpose) of Υ,
we also compute the corresponding Gram matrices:
α(S)∗

−1
= Cbα(S),β(S) × Cbβ(S),β(S)
× φβ(S) (b)
1

1
−1
∗
∗
=
Υα(S) × Υβ(S) ×
Υβ(S) × Υβ(S) + λI
N
N
× φβ(S) (b)

−1
∗
= Υα(S) × Gβ(S),β(S) + λN I
× Υβ(S) × φβ(S) (b)
{z
}
|

ωb

α(S)

=Υ

ωb , for Stage 1A;

and very similarly E[ ⊗ φα(Sk ) (α(Sk )) | β(S) = b] =
k=1

Υext(S) ωb for Stage 1B, where λ is the regularization
parameter.
We now collect all the estimated features for the inside
tree and the extended inside trees that we calculated
from the above S1A and S1B regressions into two operators:


Υµα(S) = Υα(S) ωβ(S)1 , ... , Υα(S) ωβ(S)N


Υµext(S) = Υext(S) ωβ(S)1 , ... , Υext(S) ωβ(S)N
And then for the stage 2 regression, we obtain the
linear operator W S as:
1

∗
WS =
Υµext(S) × Υµα(S)
N
1
−1
∗
×
Υµα(S) × Υµα(S) + λI
N

−1
∗
µext(S)
=Υ
× Gµα(S) ,µα(S) + λN I
× Υµα(S) ,
where the last equality follows from the matrix inversion lemma.
5.3

α(S)

Gα(S),α(S) = Υ
×Υ
∗
Gβ(S),β(S) = Υβ(S) × Υβ(S)

5.2

E[φα(S) (α(S)) | β(S) = b] = µ
bα(S)|β(S)=b

K

Then for each non-leaf separator set S in the latent
junction tree T , specify characteristic kernels for α(S)
and β(S):
def
Kα(S) (α(S), α(S)0 ) = hφα(S) (α(S)), φα(S) (α(S)0 )iG
def
Kβ(S) (β(S), β(S)0 ) = hφβ(S) (β(S)), φβ(S) (β(S)0 )iH

i=1

the Stage 1A and Stage 1B regressions, we would have:

Two-Stage Regression in the Hilbert
Space

After representing the continuous distributions of variables as points embedded in the corresponding RKHS,
we can easily perform the two-stage regression through
linear algebra operations over the Hilbert Space. For
example, if we are using kernel ridge regression to run
7
Here we will follow the notations of (Boots et al., 2013)
in our derivation.

Message Passing Inference in Hilbert
Space

Now we first compute the estimated root covariance
tensor:
N
P
⊗ φα(S) (α(S)d )
d=1 S∈γ(Cr )
Cr
,
T
=
N
and then the initial leaf messages that we send from
each leaf clique node C to its parent clique node Cp
would be:
mC→Cp = ⊗ ξ(X)
X∈δ(C)

where the function ξ(X) is defined to be8 :

φRX (x),
if X is observed to be x
v
ξ(X) =
φ
(x)dx
 u X
, if X is not observed.
v−u
8

Let [u,v] denote the range of value that X can take.

Learning General Latent-Variable Graphical Models with Predictive Belief Propagation and HSE

In the actual calculations, we can use the empirical
average of φX (x)Rover finite number of samples in [u, v]
v
φX (x)dx
.
to approximate u
v−u
The remaining steps of the message-passing inference
procedure follow exactly the same form as in the discrete case, and eventually all the computation involved
will be nicely carried out over the finite-dimensional
Gram matrices. To our knowledge, this is the first algorithm that has the ability to learn continuous-valued
latent-variable graphical models in a completely nonparametric fashion without making any assumptions
about the function forms of the variables’ distributions.

6

Experiments

We designed three sets of experiments to evaluate the
performance of our proposed algorithm, including both
synthetic data and real data and ranging from both
discrete domain and continuous domain.
6.1

gorithm and average across all possible joint realizations of the variables (g, h, e). We report the results
in Figure 6 below. From the plot, we see that the
average KL divergence between our algorithm’s results and the ground truth posterior quickly decreases
and approaches 0 as the size of the training data increases. This result demonstrates that our algorithm
learns quickly to perform accurate inference over latent graphical models.

Synthetic Dataset

Figure 6: Our algorithm’s learning performance on the
synthetic dataset. We test our algorithm on 6 different
models. Each model has a different set of randomly
initialized parameters.

In this experiment, we test the performance of our algorithm on the task of learning and running inference
on the discrete-valued latent-variable graphical model
depicted in Figure 5 using artifically generated synthetic data and compare it with the EM algorithm.
G

J

H

D

A

Figure 7: Comparison between our algorithm and the
EM algorithm on the synthetic dataset.

?
F

B

C

E

Figure 5: A toy directed latent-variable graphical
model. Green nodes indicate observable variables and
red nodes indicate latent variables.

We also run the EM algorithm to learn the same model
with the same synthetic training dataset, and compare its performance and training time with our algorithm. The results are plotted in Figure 7. Our
algorithm achieves comparable learning performance,
but is much faster to train than EM.
6.2

We randomly set all the local conditional probability
tables in this directed graphical model as the ground
truth model parameters, and then sample a dataset
of 30000 sets of joint observations for the observable
variables. Next we apply our new algorithm to learn
this model, and evaluate its performance on the task of
inferring the posterior distribution of variable D given
the observed values at variables G, H, and E. We
compute the Kullback-Leibler divergence between our
algorithm’s inferred posterior and the ground truth
posterior calculated using the exact Shafer–Shenoy al-

DNA Splice Dataset

Here we consider the computational biology task of
classifying splice junctions in DNA sequences, and test
our learning algorithm using the DNA splice dataset
on the UCI machine learning repository (Asuncion &
Newman, 2007). This dataset contains a total of 3190
length-60 DNA sequences. Each sequence is labeled as
one of the three different categories: Intron/Exon site,
Exon/Intron site or Neither. Our goal is to learn to
classify unseen DNA instances. We adopt a generative
approach and use second-order nonhomogeneous Hid-

Borui Wang, Geoffrey J. Gordon

Figure 8: The performance of our learning algorithm
on the DNA splice dataset.

cal models with loopy structures, since the feature descriptors are continuous-valued and have complex distributions. But as we discussed in section 5, our new
algorithm can smoothly handle the learning problem
of such models through RKHS embeddings. In our
experiment, using Gaussian radial basis function kernels with bandwidth parameter σ = 3, regularization
parameter λ = 0.01 and model length K = 30, we
are able to achieve an overall recognition accuracy of
75.69%. Figure 10 plots the recognition accuracies and
the normalized confusion matrix of our results.

dern Markov Models to model the DNA splice junction
sequences. For each of the three categories, we use our
algorithm to learn a different second-order nonhomogeneous HMM. At test time, we compute the probabilities that a test instance is generated from each of
the three second-order HMMs, choosing the one with
the highest probability as the predicted category. We
are able to achieve an overall classification accuracy of
87.97%, and the detailed results are reported in Figure
8.
6.3

Human Action Recognition in Videos

In this experiment we consider the computer vision
problem of recognizing human actions from videos, using the classic KTH human action dataset (Schuldt
et al., 2004). The KTH dataset contains a total of 2391
video sequences of 6 different categories of actions:
boxing, handclapping, handwaving, jogging, running
and walking. Here we choose to use the second-order
non-homogeneous state space model to model the generative process of human action videos, as depicted
in Figure 9. For each video episode in the dataset,
we take K key frames evenly spaced across the whole
time frame, and then extract a 2800-dimensional histogram of oriented gradients (HOG) feature vector and
a 180-dimensional histogram of optical flow (HOF)
feature vector for each of the key frames. During
training, we concatenate the HOG and HOF feature
vectors into a 2980-dimensional vector that summarizes each key frame, and then use dimensionality
reduction to reduce it to 5 dimensions. These 5dimensional feature vectors are the observable variables Xt in our latent graphical model. For each of
the 6 action categories, we use our algorithm to learn
a different second-order nonhomogeneous state space
model. At test time, we compute the probabilities
P (X1 = x1 , X2 = x2 , ..., XK = xK | class label) that
a test video sequence is generated from each of the
6 second-order state space models, choosing the one
with the highest probability as the predicted action
category. This formulation poses a difficult problem
of learning continuous-valued latent-variable graphi-

Figure 9: Illustration of our second-order nonhomogeneous state space model for human action
recognition from videos.

Figure 10: The performance of our learning algorithm
on the KTH human action recognition dataset.

7

Conclusion

We have developed a new algorithm for learning general latent-variable graphical models using predictive
belief propagation, two-stage instrumental regression
and RKHS embeddings of distributions. We proved
that the algorithm gives a consistent estimator of the
inference results. We evaluate the algorithm’s learning performance on both synthetic and real datasets,
showing that it learns different types of latent graphical model efficiently and achieves good inference results in both discrete and continuous domains. We
believe that our algorithm provides a powerful and
flexible new learning framework for latent graphical
models.

Learning General Latent-Variable Graphical Models with Predictive Belief Propagation and HSE

References
Anandkumar, A., Ge, R., Hsu, D., Kakade, S., and
Telgarsky, M. Tensor decompositions for learning
latent variable models. Journal of Machine Learning
Research, 15:2773–2832, 2014.
Anderson, S., Olesen, K., Jensen, F., and Jensen, F.
HUGIN - a shell for building Bayesian belief universes for expert systems. In Proceedings of the 18th
International Joint Conference on Artificial Intelligence (IJCAI), pp. 1080–1085, 1989.
Asuncion, A. and Newman, D. J. UCI machine
learning repository.
https://archive.ics.
uci.edu/ml/datasets/Molecular+Biology+
(Splice-junction+Gene+Sequences), 2007.
Boots, B., Siddiqi, S., and Gordon, G. Closing the
learning-planning loop with predictive state representations. The International Journal of Robotics
Research, 30(7):954–966, 2011.
Boots, B., Gretton, A., and Gordon, G. Hilbert space
embeddings of predictive state representations. In
Proceedings of the 29th Conference on Uncertainty
in Artificial Intelligence (UAI 2013), pp. 92–101,
2013.
Dempster, A., Laird, N., and Rubin, D. Maximum
likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B, 39(1):1–38, 1977.
Hefny, A., Downey, C., and Gordon, G. Supervised
learning for dynamical system learning. In Advances
in Neural Information Processing Systems 28 (NIPS
2015), 2015.
Hsu, D., Kakade, S., and Zhang, T. A spectral algorithm for learning hidden Markov models. In Proceedings of the 22nd Annual Conference on Learning
Theory, 2009.
Lauritzen, S. and Spiegelhalter, D. Local computations with probabilities on graphical structures and
their application to expert systems. Journal of the
Royal Statistical Society, Series B, 50(2):157–224,
1988.
Littman, M., Sutton, R., and Singh, S. Predictive representations of state. In Advances in Neural Information Processing Systems 14 (NIPS 2001), 2001.
Parikh, A., Song, L., Ishteva, M., Teodoru, G., and
Xing, E. A spectral algorithm for latent junction
trees. In Proceedings of the 28th Conference on Uncertainty in Artificial Intelligence (UAI 2012), pp.
675–684, 2012.

Schuldt, C., Laptev, I., and Caputo, B. Recognizing
human actions: A local svm approach. In Proceedings of the 17th International Conference on Pattern
Recognition (ICPR 2004), 2004.
Shenoy, P. and Shafer, G. Axioms for probability and
belief-function propagation. In Proceedings of the
6th Conference on Uncertainty in Artificial Intelligence (UAI 1990), pp. 169–198, 1990.
Singh, S., James, M., and Rudary, M. Predictive state
representations: A new theory for modeling dynamical systems. In Proceedings of the 20th Conference
on Uncertainty in Artificial Intelligence (UAI 2004),
pp. 512–519, 2004.
Stock, J. H. and Watson, M. W. Introduction to
Econometrics. Addison-Wesley, 3rd edition, 2011.

