Journal of Artiﬁcial Intelligence Research 15 (2001) 391-454

Submitted 6/18; published 12/01

Parameter Learning of Logic Programs for
Symbolic-statistical Modeling
Taisuke Sato
Yoshitaka Kameya

sato@mi.cs.titech.ac.jp
kame@mi.cs.titech.ac.jp

Dept. of Computer Science, Graduate School of Information
Science and Engineering, Tokyo Institute of Technology
2-12-1 Ookayama Meguro-ku Tokyo Japan 152-8552

Abstract
We propose a logical/mathematical framework for statistical parameter learning of parameterized logic programs, i.e. deﬁnite clause programs containing probabilistic facts witha parameterized distribution. It extends the traditional least Herbrand model semantics inlogic programming to distribution semantics, possible world semantics with a probability
distribution which is unconditionally applicable to arbitrary logic programs including onesfor HMMs, PCFGs and Bayesian networks.
We also propose a new EM algorithm, the graphical EM algorithm, that runs for aclass of parameterized logic programs representing sequential decision processes where eachdecision is exclusive and independent. It runs on a new data structure called support graphs
describing the logical relationship between observations and their explanations, and learns
parameters by computing inside and outside probability generalized for logic programs.
The complexity analysis shows that when combined with OLDT search for all explanations for observations, the graphical EM algorithm, despite its generality, has the sametime complexity as existing EM algorithms, i.e. the Baum-Welch algorithm for HMMs, theInside-Outside algorithm for PCFGs, and the one for singly connected Bayesian networks
that have been developed independently in each research ﬁeld. Learning experiments with
PCFGs using two corpora of moderate size indicate that the graphical EM algorithm cansigniﬁcantly outperform the Inside-Outside algorithm.

1. Introduction
Parameter learning is common in various ﬁelds from neural networks to reinforcement learning to statistics. It is used to tune up systems for their best performance, be they classiﬁers
or statistical models. Unlike these numerical systems described by mathematical formulas however, symbolic systems, typically programs, do not seem amenable to any kind ofparameter learning. Actually there has been little literature on parameter learning of programs.
This paper is an attempt to incorporate parameter learning into computer programs.
The reason is twofold. Theoretically we wish to add the ability of learning to computer
programs, which the authors believe is a necessary step toward building intelligent systems.
Practically it broadens the class of probability distributions, beyond traditionally used numerical ones, which are available for modeling complex phenomena such as gene inheritance,
consumer behavior, natural language processing and so on.

c°2001
AI Access Foundation and Morgan Kaufmann Publishers. All rights reserved.

Sato & Kameya

The type of learning we consider here is statistical parameter learning applied to logic
programs.1 We assume that facts (unit clauses) in a program are probabilistically true andhave a parameterized distribution.2 Other clauses, non-unit deﬁnite clauses, are always
true as they encode laws such as “if one has a pair of blood type genes a and b, one’s
blood type is AB”. We call logic programs of this type a parameterized logic program anduse for statistical modeling in which ground atoms3 provable from the program represent
our observations such as “one’s blood type is AB” and the parameters of the program are
inferred by performing ML (maximum likelihood) estimation on the observed atoms.
The probabilistic ﬁrst-order framework sketched above is termed statistical abduction(Sato & Kameya, 2000) as it is an amalgamation of statistical inference and abduction
where probabilistic facts play the role of abducibles, i.e. primitive hypotheses.4 Statisticalabduction is powerful in that it not only subsumes diverse symbolic-statistical frameworks
such as HMMs (hidden Markov models, Rabiner, 1989), PCFGs (probabilistic context freegrammars, Wetherell, 1980; Manning & Schütze, 1999) and (discrete) Bayesian networks(Pearl, 1988; Castillo, Gutierrez, & Hadi, 1997) but gives us freedom of using arbitrarily
complex logic programs for modeling.5
The semantic basis for statistical abduction is distribution semantics introduced by Sato
(1995). It deﬁnes a parameterized distribution, actually a probability measure, over the setof possible truth assignments to ground atoms and enables us to derive a new EM algorithm6for ML estimation called the graphical EM algorithm (Kameya & Sato, 2000).
Parameter learning in statistical abduction is done in two phases, search and EM learning. Given a parameterized logic program and observations, the ﬁrst phase searches for all
explanations for the observations. Redundancy in the ﬁrst phase is eliminated by tabulating
partial explanations using OLDT search (Tamaki & Sato, 1986; Warren, 1992; Sagonas, T.,& Warren, 1994; Ramakrishnan, Rao, Sagonas, Swift, & Warren, 1995; Shen, Yuan, You, &Zhou, 2001). It returns a support graph which is a compact representation of the discovered
explanations. In the second phase, we run the graphical EM algorithm on the support graph
1. In this paper, logic programs mean deﬁnite clause programs. A deﬁnite clause program is a set of deﬁnite
clauses. A deﬁnite clause is a clause of the form A ← L1 , . . . , Ln (0 ≤ n) where A, L1 , . . . , Ln are atoms.
A is called the head, L1 , . . . , Ln the body. All variables are universally quantiﬁed. It reads if L1 and· · · and Ln hold, then A holds. In case of n = 0, the clause is called a unit clause. A general clause is
one whose body may contain negated atoms. A program including general clauses is sometimes called ageneral program (Lloyd, 1984; Doets, 1994).
2. Throughout this paper, for familiarity and readability, we will somewhat loosely use “distribution” as asynonym for “probability measure”.
3. In logic programming, the adjective “ground” means no variables contained.
4. Abduction means inference to the best explanation for a set of observations. Logically, it is formalized asa search for an explanation E such that E, KB ` G where G is an atom representing our observation, KB
a knowledge base and E a conjunction of atoms chosen from abducibles, i.e. a class of formulas allowed
as primitive hypotheses (Kakas, Kowalski, & Toni, 1992; Flach & Kakas, 2000). E must be consistent
with KB.
5. Existing symbolic-statistical modeling frameworks have restrictions and limitations of various types compared with arbitrary logic programs (see Section 7 for details). For example, Bayesian networks do notallow recursion. HMMs and PCFGs, stochastic grammars, allow recursion but lack variables and data
structures. Recursive logic programs are allowed in Ngo and Haddawy’s (1997) framework but they
assume domains are ﬁnite and function symbols seem prohibited.
6. “EM algorithm” stands for a class of iterative algorithms for ML estimation with incomplete data
(McLachlan & Krishnan, 1997).

392

Parameter Learning of Logic Programs for Symbolic-statistical Modelingand learn the parameters of the distribution associated with the program. Redundancy inthe second phase is removed by the introduction of inside and outside probability for logic
programs computed from the support graph.
The graphical EM algorithm has accomplished, when combined with OLDT search for
all explanations, the same time complexity as the specialized ones, e.g. the Baum-Welchalgorithm for HMMs (Rabiner, 1989) and the Inside-Outside algorithm for PCFGs (Baker,1979), despite its generality. What is surprising is that, when we conducted learning experiments with PCFGs using real corpora, it outperformed the Inside-Outside algorithm byorders of magnitudes in terms of time for one iteration to update parameters. These experimental results enhance the prospect for symbolic-statistical modeling by parameterized
logic programs of even more complex systems than stochastic grammars whose modelinghas been diﬃcult simply because of the lack of an appropriate modeling tool and their sheer
complexities. The contributions of this paper therefore are• distribution semantics for parameterized logic programs which uniﬁes existing symbolicstatistical frameworks,• the graphical EM algorithm (combined with tabulated search), a general yet eﬃcient
EM algorithm that runs on support graphs and• the prospect suggested by the learning experiments for modeling and learning complex
symbolic-statistical phenomena.
The rest of this paper is organized as follows. After preliminaries in Section 2, a probability space for parameterized logic programs is constructed in Section 3 as a mathematicalbasis for the subsequent sections. We then propose a new EM algorithm, the graphicalEM algorithm, for parameterized logic programs in Section 4. Complexity analysis of thegraphical EM algorithm is presented in Section 5 for HMMs, PCFGs, pseudo PCSGs andsc-BNs.7 Section 6 contains experimental results of parameter learning with PCFGs by thegraphical EM algorithm using real corpora that demonstrate the eﬃciency of the graphicalEM algorithm. We state related work in Section 7, followed by conclusion in Section 8. The
reader is assumed to be familiar with the basics of logic programming (Lloyd, 1984; Doets,1994), probability theory (Chow & Teicher, 1997), Bayesian networks (Pearl, 1988; Castillo
et al., 1997) and stochastic grammars (Rabiner, 1989; Manning & Schütze, 1999).

2. Preliminaries
Since our subject intersects logic programming and EM learning which are quite diﬀerent
in nature, we separate preliminaries.
2.1 Logic Programming and OLDT
In logic programming, a program DB is a set of deﬁnite clauses8 and the execution is search
for an SLD refutation of a given goal ← G. The top-down interpreter recursively selects the
7. Pseudo PCSGs (probabilistic context sensitive grammars) are a context-sensitive extension of PCFGs
proposed by Charniak and Carroll (1994). sc-BN is a shorthand for a singly connected Bayesian network(Pearl, 1988).
8. We do not deal with general logic programs in this paper.

393

Sato & Kameya

next goal and unfolds it (Tamaki & Sato, 1984) into subgoals using a nondeterministicallychosen clause. The computed result by the SLD refutation, i.e. a solution, is an answer
substitution (variable binding) θ such that DB ` Gθ.9 Usually there is more than one
refutation for ← G, and the search space for all refutations is described by an SLD tree
which may be inﬁnite depending on the program and the goal (Lloyd, 1984; Doets, 1994).
More often than not, applications require all solutions. In natural language processing
for instance, a parser must be able to ﬁnd all possible parse trees for a given sentence as
every one of them is syntactically correct. Similarly in statistical abduction, we need to
examine all explanations to determine the most likely one. All solutions are obtained bysearching the entire SLD tree, and there is a choice of the search strategy. In Prolog, thestandard logic programming language, backtracking is used to search for all solutions inconjunction with a ﬁxed search order for goals (textually from left-to-right) and clauses(textually top-to-bottom) due to the ease and simplicity of implementation.
The problem with backtracking is that it forgets everything until up to the previous
choice point, and hence it is quite likely to prove the same goal again and again, resulting in
exponential search time. One answer to avoid this problem is to store computed results andreuse them whenever necessary. OLDT is such an instance of memoizing scheme (Tamaki& Sato, 1986; Warren, 1992; Sagonas et al., 1994; Ramakrishnan et al., 1995; Shen et al.,2001). Reuse of proved subgoals in OLDT search often drastically reduces search time
for all solutions, especially when refutations of the top goal include many common subrefutations. Take as an example a logic program coding an HMM. For a given string s,
there exist exponentially many transition paths that output s. OLDT search applied tothe program however only takes time linear in the length of s to ﬁnd all of them unlike
exponential time by Prolog’s backtracking search.
What does OLDT have to do with statistical abduction? From the viewpoint of statistical abduction, reuse of proved subgoals, or equivalently, structure sharing of sub-refutations
for the top-goal G brings about structure sharing of explanations for G, in addition to the
reduction of search time mentioned above, thereby producing a highly compact representation of all explanations for G.
2.2 EM Learning
Parameterized distributions such as the multinomial distribution and the normal distribution provide convenient modeling devices in statistics. Suppose a random sample x1 , . . . , xTof size T on a random variable X drawn from a distribution P (X = x | θ) parameterized
by unknown θ, is observed. The value of θ is determined by ML estimation as the MLE
Q(maximum likelihood estimate) of θ, i.e. as the maximizer of the likelihood 1≤i≤T P (xi | θ).
Things get much more diﬃcult when data are incomplete. Think of a probabilistic
relationship between non-observable cause X and observable eﬀect Y such as one between
diseases and symptoms in medicine and assume that Y does not uniquely determine thecause X. Then Y is incomplete in the sense that Y does not carry enough information tocompletely determine X. Let P (X = x, Y = y | θ) be a parameterized joint distribution
over X and Y . Our task is to perform ML estimation on θ under the condition that X is
9. By a solution we ambiguously mean both the answer substitution θ itself and the proved atom Gθ, as
one gives the other.

394

Parameter Learning of Logic Programs for Symbolic-statistical Modeling

non-observable while Y is observable. Let y1 , . . . , yT be a random sample of size T drawn
Pfrom the marginal distribution P (Y = y | θ) = x P (X = x, Y = y | θ). The MLE of θ isQ
obtained by maximizing the likelihood 1≤i≤T P (yi | θ) as a function of θ.
While mathematical formulation looks alike in both cases, the latter, ML estimation withincomplete data, is far more complicated and direct maximization is practically impossible
in many cases. People therefore looked to indirect approaches to tackle the problem of
ML estimation with incomplete data to which the EM algorithm has been a standard
solution (Dempster, Laird, & Rubin, 1977; McLachlan & Krishnan, 1997). It is an iterativealgorithm applicable to a wide class of parameterized distributions including the multinomial
distribution and the normal distribution such that the MLE computation is replaced by theiteration of two easier, more tractable steps. At n-th iteration, it ﬁrst calculates the valueof Q function introduced below using current parameter value θ(n) (E-step)10 :def X

Q(θ | θ(n) ) =P (x | y, θ(n) ) ln P (x, y | θ).

(1)x

Next, it maximizes Q(θ | θ(n) ) as a function of θ and updates θ(n) (M-step):θ(n+1) = argmaxθ Q(θ | θ(n) ).

(2)

Since the old value θ(n) and the updated value θ(n+1) do not necessarily coincide, the E-steps
and M-steps are iterated until convergence, during which the (log) likelihood is assured toincrease monotonically (McLachlan & Krishnan, 1997).
Although the EM algorithm merely performs local maximization, it is used in a variety
of settings due to its simplicity and relatively good performance. One must notice howeverthat the EM algorithm is just a class name, taking diﬀerent form depending on distributions
and applications. The development of a concrete EM algorithm such as the Baum-Welchalgorithm for HMMs (Rabiner, 1989) and the Inside-Outside algorithm for PCFGs (Baker,1979) requires individual eﬀort for each case.

10. Q function is related to ML estimation as follows. We assume here only one data, y, is observed. From
Jensen’s inequality (Chow & Teicher, 1997) and the concavity of ln function, it follows thatX

P (x | y, θ (n) ) ln P (x | y, θ) −x

XP (x | y, θ(n) ) ln P (x | y, θ (n) ) ≤ 0x

and hence that
Q(θ | θ(n) ) − Q(θ(n) | θ(n) )=XP (x | y, θ(n) ) ln P (x | y, θ) −P (x | y, θ (n) ) ln P (x | y, θ(n) ) + ln P (y | θ) − ln P (y | θ(n) )

x

x

≤

X

ln P (y | θ) − ln P (y | θ(n) ).

Consequently, we haveQ(θ | θ(n) ) ≥ Q(θ(n) | θ(n) ) ⇒ ln p(y | θ) ≥ ln p(y | θ(n) ) ⇒ p(y | θ) ≥ p(y | θ(n) ).

395

Sato & Kameya

3. Distribution Semantics
In this section, we introduce parameterized logic programs and deﬁne their declarative semantics. The basic idea is as follows. We start with a set F of probabilistic facts (atoms)
and a set R of non-unit deﬁnite clauses. Sampling from F determines a set F 0 of true
atoms, and the least Herbrand model of F 0 ∪ R determines the truth value of every atom
in DB = F ∪ R. Hence every atom can be considered as a random variable, taking on 1
(true) or 0 (false). In what follows, we formalize this process and construct the underlying
probability space for the denotation of DB.
3.1 Basic Distribution PF
Let DB = F ∪R be a deﬁnite clause program in a ﬁrst-order language L with countably many
variables, function symbols and predicate symbols where F is a set of unit clauses (facts)
and R a set of non-unit clauses (rules). In the sequel, unless otherwise stated, we consider
for simplicity DB as the set of all ground instances of the clauses in DB, and assume that
F and R consist of countably inﬁnite ground clauses (the ﬁnite case is similarly treated).
We then construct a probability space for DB in two steps. First we introduce a probability
space over the Herbrand interpretations11 of F i.e. the truth assignments to ground atoms
in F . Next we extend it to a probability space over the Herbrand interpretations of all
ground atoms in L by using the least model semantics (Lloyd, 1984; Doets, 1994).
Let A1 , A2 , . . . be a ﬁxed enumeration of atoms in F . We regard an inﬁnite vector ω =hx1 , x2 , . . .i of 0s and 1s as a Herbrand interpretation of F in such a way that for i = 1, 2, . . .
Ai is true (resp. false) if and only if xi = 1 (resp. xi = 0). Under this isomorphism, the setof all possible Herbrand interpretations of F coincides with the Cartesian product:
∞
def Y

ΩF =

{0, 1}i .

i=1

We construct a probability measure PF over the sample space ΩF 12 from a collection of(n)ﬁnite joint distributions PF (A1 = x1 , . . . , An = xn ) (n = 1, 2, . . . , xi ∈ {0, 1}, 1 ≤ i ≤ n)such that








(n)

0 ≤ PF (A1 = x1 , . . . , An = xn ) ≤ 1(n)x1 ,...,xn PF (A1 = x1 , . . . , An = xn ) = 1
P
(n+1)
(n)
(A1 = x1 , . . . , An+1 = xn+1 ) = PF (A1 = x1 , . . . , An = xn ).
xn+1 PFP(3)

The last equation is called the compatibility condition. It can be proved (Chow & Teicher,1997) from the compatibility condition that there exists a probability space (ΩF , F, PF )where PF is a probability measure on F, the minimal σ algebra containing open sets of ΩF ,such that for any n,(n)

PF (A1 = x1 , . . . , An = xn ) = PF (A1 = x1 , . . . , An = xn ).
11. A Herbrand interpretation interprets a function symbol uniquely as a function on ground terms andassigns truth values to ground atoms. Since the interpretation of function symbols is common to all
Herbrand interpretations, given L, they have a one-to-one correspondence with truth assignments toground atoms in L. So we do not distinguish them.
12. We regard ΩF as a topological space with the product topology such that each {0, 1} is equipped withthe discrete topology.

396

Parameter Learning of Logic Programs for Symbolic-statistical Modeling

We call PF a basic distribution.13(n)
The choice of PF is free as long as the compatibility condition is met. If we want all(n)interpretations to be equiprobable, we should set PF (A1 = x1 , . . . , An = xn ) = 1/2n forevery hx1 , . . . , xn i. The resulting PF is a uniform distribution over ΩF just like the one
over the unit interval [0, 1]. If, on the other hand, we stipulate no interpretation except
ω0 = hc1 , c2 , . . .i should be possible, we put, for each n,((n)
PF (A1 = x1 , . . . , An = xn ) =1 if∀i xi = ci (1 ≤ i ≤ n)
0 o.w.

Then PF places all probability mass on ω0 and gives probability 0 to the rest.
Deﬁne a parameterized logic program as a deﬁnite clause program14 DB = F ∪ R where
F is a set of unit clauses, R is a set of non-unit clauses such that no clause head in R is
uniﬁable with a unit clause in F and a parameterized basic distribution PF is associated with
F . A parameterized PF is obtained from a collection of parameterized joint distributions(n)satisfying the compatibility condition. Generally, the more complex PF ’s are, the more
ﬂexible PF is, but at the cost of tractability. The choice of parameterized ﬁnite distributions
made by Sato (1995) was simple:
(2n)

PF

=

(ON 1 = x1 , OFF 2 = x2 , . . . , OFF 2n = x2n | θ1 , . . . , θn )nYPbs (ON 2i−1 = x2i−1 , OFF 2i = x2i | θi )

i=1

where
Pbs (ON 2i−1 = x2i−1 , OFF 2i = x2i | θi )0
if x2i−1 = x2i
θiif x2i−1 = 1, x2i = 0= 1 − θ if xi2i−1 = 0, x2i = 1.

(4)

Pbs (ON 2i−1 = x2i−1 , OFF 2i = x2i | θi ) (1 ≤ i ≤ n) represents a probabilistic binary switch,i.e. a Bernoulli trial, using two exclusive atoms ON 2i−1 and OFF 2i in such a way that either
one of them is true on each trial but never both. θi is a parameter specifying the probability
that the switch i is on. The resulting PF is a probability measure over the inﬁnite product of
independent binary outcomes. It might look too simple but expressive enough for Bayesian
networks, Markov chains and HMMs (Sato, 1995; Sato & Kameya, 1997).
3.2 Extending PF to PDB
In this subsection, we extend PF to a probability measure PDB over the possible world s
for L, i.e. the set of all possible truth assignments to ground atoms in L through the least
13. This naming of PF , despite its being a probability measure, partly reﬂects the observation that it behaves
like an inﬁnite joint distribution PF (A1 = x1 , A2 = x2 , . . .) for an inﬁnite random vector hA1 , A2 , . . .i(n)of which PF (A1 = x1 , . . . , An = xn ) (n = 1, 2, . . .) are marginal distributions. Another reason is
intuitiveness. These considerations apply to PDB deﬁned in the next subsection as well.
14. Here clauses are not necessarily ground.

397

Sato & Kameya

Herbrand model (Lloyd, 1984; Doets, 1994). Before proceeding however, we need a coupleof notations. For an atom A, deﬁne Ax by(

Ax = Aif x = 1
Ax = ¬A if x = 0.

Next take a Herbrand interpretation ν ∈ ΩF of F . It makes some atoms in F true andothers false. Let Fν be the set of atoms made true by ν. Then imagine a deﬁnite clauseprogram DB 0 = R ∪ Fν and its least Herbrand model MDB 0 (Lloyd, 1984; Doets, 1994).
MDB 0 is characterized as the least ﬁxed point of a mapping TDB 0 (·) below¯)¯ there is some A ← B , . . . , B ∈ DB 0 (0 ≤ k)¯1kA¯
¯ such that {B1 , . . . , Bk } ⊆ I(
def

TDB 0 (I) =where I is a set of ground atoms.15 Or equivalently, it is inductively deﬁned by
I0 = ∅
In+1 = TDB 0 (In )
MDB 0

=

[

In .

n

Taking into account the fact that MDB 0 is a function of ν ∈ ΩF , we henceforth employ afunctional notation MDB (ν) to denote MDB 0 .
Turning back, let A1 , A2 , . . . be again an enumeration, but of all ground atoms in L.16
Form ΩDB , similarly to ΩF , as the Cartesian product of denumerably many {0, 1}’s and identify it with the set of all possible Herbrand interpretations of the ground atoms A1 , A2 , . . .
in L, i.e. the possible world s for L. Then extend PF to a probability measure PDB over ΩDB(n)
as follows. Introduce a series of ﬁnite joint distributions PDB (A1 = x1 , . . . , An = xn ) forn = 1, 2, . . . by
[Ax1 1 ∧ · · · ∧ Axnn ]Fdef={ν ∈ ΩF | MDB (ν) |= Ax1 1 ∧ · · · ∧ Axnn }PDB (A1 = x1 , . . . , An = xn )defPF ([Ax1 1 ∧ · · · ∧ Axnn ]F ).

(n)

=

(n)

Note that the set [Ax1 1 ∧ · · · ∧ Axnn ]F is PF -measurable and by deﬁnition, PDB ’s satisfy thecompatibility condition
X

(n+1)

PDB

(n)

(A1 = x1 , . . . , An+1 = xn+1 ) = PDB (A1 = x1 , . . . , An = xn ).

xn+1

Hence there exists a probability measure PDB over ΩDB which is an extension of PF suchthat
PDB (A1 = x1 , . . . , An = xn ) = PF (A1 = x1 , . . . , An = xn )15. I deﬁnes, mutually, a Herbrand interpretation such that a ground atom A is true if and only if A ∈ I.
A Herbrand model of a program is a Herbrand interpretation that makes every ground instance of every
clause in the program true.
16. Note that this enumeration enumerates ground atoms in F as well.

398

Parameter Learning of Logic Programs for Symbolic-statistical Modelingfor any ﬁnite atoms A1 , . . . , An in F and for every binary vector hx1 , . . . , xn i (xi ∈ {0, 1}, 1 ≤i ≤ n). Deﬁne the denotation of the program DB = F ∪ R w.r.t. PF to be PDB . The denotational semantics of parameterized logic programs deﬁned above is called distribution
semantics. As remarked before, we regard PDB as a kind of inﬁnite joint distribution
PDB (A1 = x1 , A2 = x2 , . . .). Mathematical properties of PDB are listed in Appendix A
where our semantics is proved to be an extension of the standard least model semantics inlogic programming to possible world semantics with a probability measure.
3.3 Programs as Distributions
Distribution semantics views parameterized logic programs as expressing distributions. Traditionally distributions have been expressed by using mathematical formulas but the use of
programs as (discrete) distributions gives us far more freedom and ﬂexibility than mathematical formulas in the construction of distributions because they have recursion and arbitrary composition. In particular a program can contain inﬁnitely many random variables
as probabilistic atoms through recursion, and hence can describe stochastic processes that
potentially involve inﬁnitely many random variables such as Markov chains and derivationsin PCFGs (Manning & Schütze, 1999).17
Programs also enable us to procedurally express complicated constraints on distributions
such as “the sum of occurrences of alphabets a or b in an output string of an HMM must bea multiple of three”. This feature, procedural expression of arbitrarily complex (discrete)
distributions, seems quite helpful in symbolic-statistical modeling.
Finally, providing mathematically sound semantics for parameterized logic programs
is one thing, and implementing distribution semantics in a tractable way is another. In
the next section, we investigate conditions on parameterized logic programs which makeprobability computation tractable, thereby making them usable as a means for large scale
symbolic-statistical modeling.

4. Graphical EM Algorithm
According to the preceding section, a parameterized logic program DB = F ∪ R in aﬁrst-order language L with a parameterized basic distribution PF (· | θ) over the Herbrandinterpretations of ground atoms in F speciﬁes a parameterized distribution PDB (· | θ) overthe Herbrand interpretations for L. In this section, we develop, step by step, an eﬃcient EM
algorithm for the parameter learning of parameterized logic programs by interpreting PDB
as a distribution over the observable and non-observable events. The new EM algorithm istermed the graphical EM algorithm. It is applicable to arbitrary logic programs satisfying
certain conditions described later provided the basic distribution is a direct product ofmulti-ary random switches, which is a slight complication of the binary ones introduced in
Section 3.1.
From this section on, we assume that DB consists of usual deﬁnite clauses containing(universally quantiﬁed) variables. Deﬁnitions and changes relating to this assumption are
17. An inﬁnite derivation can occur in PCFGs. Take a simple PCFG {p : S → a, q : S → SS} where S is a
start symbol, a a terminal symbol, p + q = 1 and p, q > 0. In this PCFG, S is rewritten either to a withprobability p or to SS with probability q. The probability of the occurrence of an inﬁnite derivation is
calculated as max {0, 1 − (p/q)} which is non-zero when q > p (Chi & Geman, 1998).

399

Sato & Kameya

listed below. For a predicate p, we introduce iﬀ (p), the iﬀ deﬁnition of p byiﬀ (p)def=∀x (p(x) ↔ ∃y1 (x = t1 ∧ W1 ) ∨ · · · ∨ ∃yn (x = tn ∧ Wn )) .

Here x is a vector of new variables of length equal to the arity of p, p(ti ) ← Wi (1 ≤ i ≤n, 0 ≤ n), an enumeration of clauses about p in DB, and yi , a vector of variables occurringin p(ti ) ← Wi . Then deﬁne comp(R) as follows.
head(R)def={B | B is a ground instance of a clause head appearing in R}

iﬀ (R)

def={iﬀ (p) | p appears in a clause head in R}

Eq

def

{f (x) = f (y) → x = y | f is a function symbol}=∪ {f (x) 6= g(y) | f and g are diﬀerent function symbols}∪ {t 6= x | t is a term properly containing x}comp(R)def=

iﬀ (R) ∪ Eq

Eq , Clark’s equational theory (Clark, 1978), deductively simulates uniﬁcation. Likewisecomp(R) is a ﬁrst-order theory which deductively simulates SLD refutation with the helpof Eq by replacing a clause head atom with the clause body (Lloyd, 1984; Doets, 1994).
We here introduce some deﬁnitions which will be frequently used. Let B be an atom.
An explanation for B w.r.t. DB = F ∪ R is a conjunction S such that S, R ` B, and as a
set comprised of its conjuncts, S ⊂ F holds and no proper subset of S satisﬁes this. The
set of all explanations for B is called the support set for B and designated by ψDB (B).18
4.1 Motivating Example
First of all, we review distribution semantics by a concrete example. Consider the following
program DB b = Fb ∪ Rb in Figure 1 modeling how one’s blood type is determined by blood
type genes probabilistically inherited from the parents.19
The ﬁrst four clauses in Rb state a blood type is determined by a genotype, i.e. a pair ofblood type genes a, b and o. For instance, btype(’A’):- (gtype(a,a) ; gtype(a,o) ;gtype(o,a)) says that one’s blood type is A if his (her) genotype is ha, ai, ha, oi or ho, ai.
These are propositional rules.
Succeeding clauses state general rules in terms of logical variables. The ﬁfth clausesays that regardless of the values of X and Y, event gtype(X,Y) (one’s having genotype
hX, Yi) is caused by two events, gene(father,X) (inheriting gene X from the father) andgene(mother,Y) (inheriting gene Y from the mother). gene(P,G):- msw(gene,P,G) is aclause connecting rules in Rb with probabilistic facts in Fb . It tells us that the gene G
is inherited from a parent P if a choice represented by msw(gene,P,G)20 is made. The
18. This deﬁnition of a support set diﬀers from the one used by Sato (1995) and Kameya and Sato (2000).
19. When we implicitly emphasize the procedural reading of logic programs, Prolog conventions are employed(Sterling & Shapiro, 1986). Thus, ; stands for “or”, , “and” :- “implied by” respectively. Strings
beginning with a capital letter are (universally quantiﬁed) variables, but quoted ones such as ’A’ areconstants. The underscore is an anonymous variable.
20. msw is an abbreviation of “multi-ary random switch” and msw(·, ·, ·) expresses a probabilistic choice fromﬁnite alternatives. In the framework of statistical abduction, msw atoms are abducibles from which
explanations are constructed as a conjunction.

400

Parameter Learning of Logic Programs for Symbolic-statistical Modeling











Rb

btype(’A’)
btype(’B’)
btype(’O’)=
btype(’AB’)
gtype(X,Y)
gene(P,G)::::::-(gtype(a,a) ; gtype(a,o) ; gtype(o,a)).
(gtype(b,b) ; gtype(b,o) ; gtype(o,b)).
gtype(o,o).
(gtype(a,b) ; gtype(b,a)).
gene(father,X), gene(mother,Y).
msw(gene,P,G).

Fb= {msw(gene,father,a), msw(gene,father,b), msw(gene,father,o),msw(gene,mother,a), msw(gene,mother,b), msw(gene,mother,o)}Figure 1: ABO blood type program DB b

genetic knowledge that the choice of G is by chance and made from {a, b, o} is expressed by
specifying a joint distribution Fb as follows.
def

PFb (msw(gene,t,a) = x, msw(gene,t,b) = y, msw(gene,t,o) = z | θa , θb , θo ) = θax θby θoz
where x, y, z ∈ {0, 1}, x + y + z = 1, θa , θb , θo ∈ [0, 1], θa + θb + θo = 1 and t is either
father or mother. Thus θa is the probability of inheriting gene a from a parent. Statisticalindependence of the choice of gene, once from father and once from mother, is expressedby putting
PFb ( msw(gene,father,a) = x, msw(gene,father,b) = y, msw(gene,father,o) = z,
msw(gene,mother,a) = x0 , msw(gene,mother,b) = y 0 , msw(gene,mother,o) = z 0| θa , θ b , θ o )= PFb (x, y, z | θa , θb , θo )PFb (x0 , y 0 , z 0 | θa , θb , θo ).
In this setting, atoms representing our observation are obs(DB b ) = {btype(’A’), btype(’B’),
btype(’O’), btype(’AB’)}. We observe one of them, say btype(’A’), and infer a possible
explanation S, i.e. a minimal conjunction of abducibles msw(gene,·,·) such that
S, Rb ` btype(’A’).
S is obtained by applying a special SLD refutation procedure to the goal ← btype(’A’)
which preserves msw atoms resolved upon in the refutation. Three explanations are found.
S1 = msw(gene,father,a) ∧ msw(gene,mother,a)S2 = msw(gene,father,a) ∧ msw(gene,mother,o)S3 = msw(gene,father,o) ∧ msw(gene,mother,a)So ψDB b (btype(a)), the support set for btype(a), is {S1 , S2 , S3 }. The probability of eachexplanation is respectively computed as PFb (S1 ) = θa2 and PFb (S2 ) = PFb (S3 ) = θa θo . From
Proposition A.2 in Appendix A, it follows that PDB b (btype(’A’)) = PDB b (S1 ∨ S2 ∨ S3 ) =PFb (S1 ∨ S2 ∨ S3 ) and that
PDB b (btype(’A’) | θa , θb , θo ) = PFb (S1 ) + PFb (S2 ) + PFb (S3 )= θa2 + 2θa θo .
401

Sato & Kameya

Here we used the fact that S1 , S2 and S3 are mutually exclusive as the choice of gene is
exclusive. Parameters, i.e. θa , θb and θo are determined by ML estimation performed on arandom sample such as {btype(’A’), btype(’O’), btype(’AB’)} of btype as follows.
hθa , θb , θo i = argmaxhθa ,θb ,θo i PDB b (btype(’A’))PDB b (btype(’O’))PDB b (btype(’AB’))= argmaxhθa ,θb ,θo i (θa2 + 2θa θo )θo2 θa θb
This program contains neither function symbol nor recursion though our semantics
allows for them. Later we see an example containing both, a program for an HMM (Rabiner& Juang, 1993).
4.2 Four Simplifying Conditions
DB b in Figure 1 is simple and probability computation is easy. This is not generally thecase. Since our primary interest is learning, especially eﬃcient parameter learning of parameterized logic programs, we hereafter concentrate on identifying what property of a program
makes probability computation easy like DB b , thereby makes eﬃcient parameter learning
possible.
To answer this question precisely, let us formulate the whole modeling process. Suppose
there exist symbolic-statistical phenomena such as gene inheritance for which we hope
to construct a probabilistic computational model. We ﬁrst specify a target predicate p
whose ground atom p(s) represents our observation of the phenomena. Then to explain
the empirical distribution of p, we write down a parameterized logic program DB = F ∪ R
having a basic distribution PF with parameter θ that can reproduce all observable patternsof p(s). Finally, observing a random sample p(s1 ), . . . , p(sT ) of ground atoms of p, weQadjust θ by ML estimation, i.e. by maximizing the likelihood L(θ) = Tt=1 PDB (p(st ) | θ) sothat PDB (p(·) | θ) approximates as closely to the empirically observed distribution of p as
possible.
At ﬁrst sight, this formulation looks right, but in reality it is not. Suppose two events
p(s) and p(s0 ) (s 6= s0 ) are observed. We put L(θ) = PDB (p(s) | θ)PDB (p(s0 ) | θ). But thiscannot be a likelihood at all simply because in distribution semantics, p(s) and p(s0 ) are
two diﬀerent random variables, not two realizations of the same random variable.
A quick remedy is to note that in the case of blood type program DB b where obs(DB b ) ={btype(’A’), btype(’B’), btype(’O’), btype(’AB’)} are observable atoms, only one ofthem is true for each observation, and if some atom is true, others must be false. In otherwords, these atoms collectively behave as a single random variable having the distribution
PDB b whose values are obs(DB b ).
Keeping this in mind, we introduce the following condition. Let obs(DB) (⊂ head(R))
be a set of ground atoms which represent observable events. We call them observable atoms.
Uniqueness condition:P
PDB (G ∧ G0 ) = 0 for any G 6= G0 ∈ obs(DB), andG∈obs(DB) PDB (G) = 1.

402

Parameter Learning of Logic Programs for Symbolic-statistical Modeling

The uniqueness condition enables us to introduce a new random variable Yo representing
our observation. Fix an enumeration G1 , G2 , . . . of observable atoms in obs(DB) and deﬁne
Yo by21
Yo (ω) = k iﬀ ω |= Gk for ω ∈ ΩDB (k ≥ 1).

(5)
QT

Let Gk1 , Gk2 , . . . , GkT ∈ obs(DB) be a random sample of size T . Then L(θ) = t=1 PDB (Gkt |Qθ) = Tt=1 PDB (Yo = kt | θ) qualiﬁes for the likelihood function w.r.t. Yo .
The second condition concerns the reduction of probability computation to addition.
Take again the blood type exmaple. The computation of PDB b (btype(’A’)) is decomposed
into a summation because explanations in the support set are mutualy exclusive. So we
introduce
Exclusiveness condition:
For every G ∈ obs(DB) and the support set ψDB (G), PDB (S ∧ S 0 ) = 0 for any S 6=
S 0 ∈ ψDB (G).
Using the exclusiveness condition (and Proposition A.2 in Appendix A), we haveX

PDB (G) =PF (S).

S∈ψDB (G)From a modeling point of view, it means that while a single event, or a single observation,
G, may have several (or even inﬁnite) explanations ψDB (G), only one of ψDB (G) is allowed
to be true for each observation.
Now introduce ΨDB , i.e. the set of all explanations relevant to obs(DB) by
[

def

ΨDB =

ψDB (G)G∈obs(DB)and ﬁx an enumeration S1 , S2 , . . . of explanations in ΨDB . It follows from Proposition A.2,
the uniqueness condition and the exclusiveness condition that
PDB (Si ∧ Sj ) = 0 for i 6= j
X

X

PDB (S) =

S∈ΨDB

X

and
PDB (S)G∈obs(DB) S∈ψDB (G)X=

PDB (G)G∈obs(DB)= 1.
So we are able to introduce under the uniqueness condition and the exclusiveness condition
yet another random variable Xe , representing an explanation for G, deﬁned byXe (ω) = k iﬀ ω |= Sk for ω ∈ ΩDB .

(6)

The third condition concerns termination.

21.

P

P (G) = 1 only guarantees that the measure of { ω | ω |= Gk for some k (≥ 1)} is one, soG∈obs(DB) DB
there can be some ω satisfying no Gk ’s. In such case, we put Yo (ω) = 0. But values on a set of measurezero do not aﬀect any part of the discussion that follows. This also applies to the deﬁnition of Xe in (6).

403

Sato & Kameya

Finite support condition:
For every G ∈ obs(DB) ψDB (G) is ﬁnite.
PDB (G) is then computed from the support set ψDB (G) = {S1 , . . . , Sm } (0 ≤ m), withP
the help of the exclusiveness condition, as a ﬁnite summation m
i=1 PF (Si ). This condition
prevents an inﬁnite summation that is hardly computable.
The fourth condition simpliﬁes the probability computation to multiplication. Recall
that an explanation S for G ∈ obs(DB) is a conjunction a1 ∧ · · · ∧ am of some abducibles{a1 , . . . , am } ⊂ F (1 ≤ m). In order to reduce the computation of PF (S) = PF (a1 ∧· · ·∧am )to the multiplication PF (a1 ) · · · PF (am ), we assume
Distribution condition:
F is a set Fmsw of ground atoms with a parameterized distribution Pmsw speciﬁed below.
Here atom msw(i,n,v) is intended to simulate a multi-ary random switch whose name is i
and whose outcome is v on trial n. It is a generalization of primitive probabilistic events
such as coin tossing and dice rolling.
1. Fmsw consists of probabilistic atoms msw(i,n,v). The arguments i, n and v are groundterms called switch name, trial-id and a value (of the switch i), respectively. We
assume that a ﬁnite set Vi of ground terms called the value set of i is associated witheach i, and v ∈ Vi holds.
2. Write Vi as {v1 , v2 , . . . , vm } (m = |Vi |). Then, one of the ground atoms { msw(i,n,v1 ),msw(i,n,v2 ), . . . , msw(i,n,vm )} becomes exclusively true (takes on value 1) on eachP
trial. With each i, a parameter θi,v ∈ [0, 1] such that v∈Vi θi,v = 1 is associated. θi,v
is the probability of msw(i,·,v) being true (v ∈ Vi ).
3. For each ground terms i, i0 , n, n0 , v ∈ Vi and v 0 ∈ Vi0 , random variable msw(i,n,v) is
independent of msw(i0 ,n0 ,v 0 ) if n 6= n0 or i 6= i0 .
In other words, we introduce a family of parameterized ﬁnite distributions P(i,n) such that
P(i,n) (msw(i,n,v1 ) = x1 , . . . , msw(i,n,vm ) = xm | θi,v1 , . . . , θi,vm )
(

def

=

x1
xm
θi,v
· · · θi,v
m
1
0

if
o.w.

Pm

k=1 xk = 1(7)where m = |Vi |, xk ∈ {0, 1} (1 ≤ k ≤ m), and deﬁne Pmsw as their inﬁnite productdef Y

Pmsw =

P(i,n) .

i,nUnder this condition, we can compute Pmsw (S), the probability of an explanation S, as the
product of parameters. Suppose msw(ij ,n,v) and msw(ij 0 ,n0 ,v 0 ) are diﬀerent conjuncts inan explanation S = msw(i1 ,n1 ,v1 ) ∧ · · · ∧ msw(ik ,nk ,vk ). If either j 6= j 0 or n 6= n0 holds,
they are independent by construction. Else if j = j 0 and n = n0 but v 6= v 0 , they are notindependent but Pmsw (S) = 0 by construction. As a result, whichever condition may hold,Pmsw (S) is computed from the parameters.
404

Parameter Learning of Logic Programs for Symbolic-statistical Modeling

4.3 Modeling Principle
Up to this point, we have introduced four conditions, the uniqueness condition, the exclusiveness condition, the ﬁnite support condition and the distribution condition, to simplify
probability computation. The last one is easy to satisfy. We just adopt Fmsw together with
Pmsw . So, from here on, we always assume that Fmsw has a parameterized distribution Pmsw
introduced in the previous subsection. Unfortunately the rest are not satisﬁed automatically. According to our modeling experiences however, it is only mildly diﬃcult to satisfy
the uniqueness condition and the exclusiveness condition as long as we obey the following
modeling principle.
Modeling principle: DB = Fmsw ∪R describes a sequential decision process(modulo auxiliary computations) that uniquely produces an observable atomG ∈ obs(DB) where each decision is expressed by some msw atom.22
Translated into programming level, it says that we must take care when writing a program so that for any sample F 0 from Pmsw , there must uniquely exist goal ← G (G ∈ obs(DB))
which has a successful refutation from DB 0 = F 0 ∪ R. We can conﬁrm the principle by theblood type program DB b = Fb ∪ Rb . It describes a process of gene inheritance, and foran arbitrary sample Fb0 from Pmsw , say Fb0 = {msw(gene,father,a), msw(gene,mother,o)},there exists a unique goal, ← btype(’A’) in this case, that has a successful SLD refutationfrom Fb0 ∪ Rb .
The idea behind this principle is that a decision process always produces some result (anobservable atom), and diﬀerent decision processes must diﬀer at some msw thereby entailingmutually exclusive observable atoms. So the uniqueness condition and the exclusivenesscondition will be automatically satisﬁed.
Satisfying the ﬁnite support condition is more diﬃcult as it is virtually equivalent towriting a program DB for which all solution search for ← G (G ∈ obs(DB)) always terminates. Apparently we have no general solution to this problem, but as far as speciﬁc models
such as HMMs, PCFGs and Bayesian networks are concerned, it can be met. All programsfor these models satisfy the ﬁnite support condition (and other conditions as well).
4.4 Four Conditions Revisited
In this subsection, we discuss how to relax the four simplifying conditions introduced in Subsection 4.2 for the purpose of ﬂexible modeling. We ﬁrst examine the uniqueness condition
considering its crucial role in the adaptation of the EM algorithm to our semantics.
The uniqueness condition guarantees that there exists a (many-to-one) mapping fromexplanations to observations so that the EM algorithm is applicable (Dempster et al., 1977).
It is possible, however, to relax the uniqueness condition while justifying the applicationof the EM algorithm. We assume the MAR (missing at random) condition introduced by
Rubin (1976) which is a statistical condition on how a complete data (explanation) becomes an incomplete data (observation), and is customarily assumed implicitly or explicitly
in statistics (see Appendix B). By assuming the MAR condition, we can apply our EM

22. Decisions made in the process are a ﬁnite subset of Fmsw .

405

Sato & KameyaPalgorithm to non-exclusive observations O such that O P (O) ≥ 1 where the uniquenesscondition is seemingly destroyed.
Let us see the MAR condition in action with a simple example. Imagine we walk alonga road in front of a lawn. We occasionally observe their state such as “the road is dry but
the lawn is wet”. Assume that the lawn is watered by a sprinkler running probabilistically.
The program DB rl = Rrl ∪ Frl in Figure 2 describes a sequential process which outputs
an observation observed(road(x),lawn(y)) (“the road is x and the lawn is y”) wherex, y ∈ {wet, dry}.
Rrl

=

Frl

=

{ observed(road(X),lawn(Y)):msw(rain,once,A),( A = yes, X = wet, Y = wet; A = no, msw(sprinkler,once,B),( B = on, X = dry, Y = wet; B = off, X = dry, Y = dry ) ). }{ msw(rain,once,yes), msw(rain,once,no),msw(sprinkler,once,on), msw(sprinkler,once,off) }Figure 2: DB rl

The basic distribution over Frl is speciﬁed like PFb (·) in Subsection 4.1, so we omit it.
msw(rain,once,A) in the program determines whether it rains (A = yes) or not (A = no),whereas msw(sprinkler,once,B) determines whether the sprinkler works ﬁne (B = on)or not (B = off). Since for each sampled values of A = a (a ∈ {yes, no}) and B = b(b ∈ {on, off}), there uniquely exists an observation observed(road(x),lawn(y)) (x, y ∈{wet, dry}), there is a many-to-one mapping χ : χ(a, b) = hx, yi. In other words, we
can apply the EM algorithm to the observations observed(road(x),lawn(y)) (x, y ∈{wet, dry}). What would happen if we observe exclusively either a state of the road orthat of the lawn? Logically, this means we observe ∃y observed(road(x),lawn(y)) or∃x observed(road(x),lawn(y)). Apparently the uniqueness condition is not met, because∃y observed(road(wet),lawn(y)) and ∃x observed(road(x),lawn(wet)) are compatible(they are true when it rains). Despite the non-exclusiveness of the observations, we can still
apply the EM algorithm to them under the MAR condition, which in this case translates
into that we observe either the lawn or the road randomly regardless of their state.
We now brieﬂy check other conditions. Basically they can be relaxed at the cost ofincreased computation. Without the exclusiveness condition for instance, we would need an
additional process of transforming the support set ψDB (G) for a goal G into a set of exclusive
explanations. For instance, if G has explanations {msw(a,n,v), msw(b,m,w)}, we have totransform it into {msw(a,n,v), ¬msw(a,n,v) ∧ msw(b,m,w)} and so on.23 Clearly, thistransformation is exponential in the number of msw atoms and eﬃciency concern leads toassuming the exclusiveness condition.
The ﬁnite support condition is in practice equivalent to the condition that the SLD treefor ← G is ﬁnite. So relaxing this condition might induce inﬁnite computation.
23. ¬msw(a,n,v) is further transformed to a disjunction of exclusive msw atoms like406

W
v 0 6=v,v 0 ∈Va

msw(a,n,v 0 ).

Parameter Learning of Logic Programs for Symbolic-statistical Modeling

Relaxing the distribution condition and accepting probability distributions other thanPmsw serve to expand the horizon of the applicability of parameterized logic programs. In
particular the introduction of parameterized joint distributions P (v1 , . . . , vk ) like Boltzmann distributions over switches msw1 , . . . , mswk where v1 , . . . , vk are values of the switches,
makes them correlated. Such distributions facilitate writing parameterized logic programsfor complicated decision processes in which decisions are not independent but interdependent. Obviously, on the other hand, they increase learning time, and whether the added
ﬂexibility of distributions deserves the increased learning time or not is yet to be seen.
4.5 Naive Approach to EM Learning
In this subsection, we derive a concrete EM algorithm for parameterized logic programs
DB = Fmsw ∪ R assuming that they satisfy the uniqueness condition, the exclusivenesscondition and the ﬁnite support condition.
To start, we introduce Yo , a random variable representing our observations according
to (5) based on a ﬁxed enumeration of observable atoms in obs(DB). We also introduce
another random variable Xe representing their explanations according to (6) based on someﬁxed enumeration of explanations in ΨDB . Our understanding is that Xe is non-observable
while Yo is observable, and they have a joint distribution PDB (Xe = x, Yo = y | θ) whereθ denotes relevant parameters. It is then immediate, following (1) and (2) in Section 2, todef P
derive a concrete EM algorithm from the Q function deﬁned by Q(θ | θ0 ) = x PDB (x |y, θ0 ) ln PDB (x, y | θ) whose input is a random sample of observable atoms and whose output
is the MLE of θ.
In the following, for the sake of readability, we substitute an observable atom G (G ∈obs(DB)) for Yo = y and write PDB (G | θ) instead of PDB (Yo = y | θ). Likewise we
substitute an explanation S (S ∈ ΨDB ) for Xe = x and write PDB (S, G | θ) instead of
PDB (Xe = x, Yo = y | θ). Then it follows from the uniqueness condition that(

PDB (S, G | θ) =0
if S 6∈ ψDB (G)Pmsw (S | θ) if S ∈ ψDB (G).

We need yet another notation here. For an explanation S, deﬁne the count of msw(i,n,v)
in S bydef
σi,v (S) = |{ n | msw(i,n,v) ∈ S}| .
We have done all preparations now. Suppose we make some observations G = G1 , . . . , GT
where Gt ∈ obs(DB) (1 ≤ t ≤ T ). Put
I

def

=

{i | msw(i,n,v) ∈ S ∈ ψDB (Gt ), 1 ≤ t ≤ T }θdef{θi,v | msw(i,n,v) ∈ S ∈ ψDB (Gt ), 1 ≤ t ≤ T }.

=I is a set of switch names that appear in some explanation for one of the Gt ’s and θ denotes
parameters associated with these switches. θ is ﬁnite due to the ﬁnite support condition.

407

Sato & Kameya

Various probabilities and the Q function are computed by using Proposition A.2 inAppendix A together with our assumptions as follows.
PDB (Gt | θ)Pmsw (S | θ)

=

PDB

³_

Y

=

¯ ´
¯

ψDB (Gt ) ¯ θ=XPmsw (S | θ)

(8)S∈ψDB (Gt )σ (S)θi,vi,vi∈I,v∈Vi

Q(θ | θ 0 )

T
X
X

def

=

PDB (S | Gt , θ 0 ) ln PDB (S, Gt | θ)

t=1 S∈ΨDBX=

0η(i, v, θ ) ln θi,v ≤i∈I,v∈Vi

X
i∈I,v∈ViÃη(i, v, θ 0 )η(i, v, θ ) ln P
0 0v 0 ∈Vi η(i, v , θ )
0

!

(9)

where
η(i, v, θ)X
1
Pmsw (S | θ)σi,v (S)P (Gt | θ) S∈ψ (G )t=1 DB

T
X

def

=

t

DB

P

Here we used Jensen’s inequality to obtain (9). Note that PDB (Gt | θ)−1S∈ψDB (Gt )Pmsw (S | θ)σi,v (S) is the expected count of msw(i,·,v) in an SLD refutation of Gt . Speaking
Qof the likelihood function L(θ) = Tt=1 PDB (Gt | θ), it is already shown in Subsection 2.2(footnote) that Q(θ | θ 0 ) ≥ Q(θ 0 | θ 0 ) implies L(θ) ≥ L(θ 0 ). Hence from (9), we reachthe procedure learn-naive(DB,G) below that ﬁnds the MLE of the parameters. The arrayvariable η[i, v] stores η(i, v, θ) under the current θ.
1: procedure learn-naive(DB, G) begin
2: Initialize θ with appropriate values and ε with a small positive number ;PT
3: λ(0) :=% Compute the log-likelihood.
t=1 ln PDB (Gt | θ);
4: repeat5:
foreach i ∈ I, v ∈ Vi do
6:
7:
8:

η[i, v] :=X
1
Pmsw (S | θ)σi,v (S);P (Gt | θ) S∈ψ (G )t=1 DBT
X

foreach i ∈ I, v ∈ Vi do

θi,v := Pη[i, v];
η[i, v 0 ]
i

DB

t

% Update the parameters.

v 0 ∈V9:m := m +1;
P
10:
λ(m) := Tt=1 ln PDB (Gt | θ)
11: until λ(m) − λ(m−1) < ε
12: end% Compute the log-likelihood again.
% Terminate if converged.

This EM algorithm is simple and correctly calculates the MLE of θ, but the calculation of PDB (Gt | θ) and η[i, v](Line 3, 6 and 10) may suﬀer a combinatorial explosion of
explanations. That is, |ψDB (Gt )| often grows exponentially in the complexity of the model.
For instance, |ψDB (Gt )| for an HMM with N states is O(N L ), exponential in the length Lof an input/output string. Nonetheless, suppressing the explosion to realize eﬃcient computation in a polynomial order is possible, under suitable conditions, by avoiding multiple
computations of the same subgoal as we see next.
408

Parameter Learning of Logic Programs for Symbolic-statistical Modeling

4.6 Inside Probability and Outside Probability for Logic Programs
In this subsection, we generalize the notion of inside probability and outside probability(Baker, 1979; Lari & Young, 1990) to logic programs. Major computations in learn-naive(DB,G)P
are those of two terms in Line 6, PDB (Gt | θ) and S∈ψDB (Gt ) Pmsw (S | θ)σi,v (S). Computational redundancy lurks in the naive computation of both terms. We show it by an example.
Suppose there is a propositional program DB p = Fp ∪ Rp where Fp = {a, b, c, d, m} and

f←a∧g





 f←b∧g

Rp =

g←c




g←d∧h




(10)

h ← m.

Here f is an observable atom. We assume that a, b, c, d and m are independent and alsothat {a, b} and {c, d} are pair-wise exclusive. Then the support set for f is calculated asψDB p (f) = {a ∧ c, a ∧ d ∧ m, b ∧ c, b ∧ d ∧ m }.
Hence, in light of (8), we may compute PDB p (f) as
PDB p (f) = PFp (a ∧ c) + PFp (a ∧ d ∧ m) + PFp (b ∧ c) + PFp (b ∧ d ∧ m).

(11)

This computation requires 6 multiplications (because PFp (a ∧ c) = PFp (a)PFp (c) etc.) and3 additions. On the other hand, it is possible to compute PDB p (f) much more eﬃciently byfactoring out common computations. Let A be a ground atom. Deﬁne the inside probability
β(A) of A asdefβ(A) = PDB (A | θ).24(12)
Then by applying Theorem A.1 in Appendix A tocomp(Rp ) ` f ↔ (a ∧ g) ∨ (b ∧ g), g ↔ c ∨ (d ∧ h), h ↔ m(13)which unconditionally holds in our semantics, and by using the independent and the exclusiveness assumption made on Fp , the following equations about inside probability arederived.
 β(f) = β(a)β(g) + β(b)β(g)β(g) = β(c) + β(d)β(h)(14)
 β(h) = β(m)PDB p (f)(= β(f)) is obtained by solving (14) about β(f), for which only 3 multiplications
and 2 additions are required.
It is quite straightforward to generalize (14) but before proceeding, look at a program
DB q = {m} ∪ {g:-m ∧ m, g:-m} where g is an observable atom and m the only msw atom.
We have g ↔ (m ∧ m) ∨ m in our semantics, but to compute P (g) = P (m)P (m) + P (m) is
clearly wrong as it ignores the fact that clause bodies for g, i.e. m∧m and m are not mutually
exclusive, and atoms in the clause body m∧m are not independent (here P (·) = PDB q (·)).
Similarly, if we set a = b = c = d = m, the equation (14) will be totally incorrect.
24. Note that if A is a fact in F , β(A) = Pmsw (A | θ).

409

Sato & Kameya

We therefore add, temporarily in this subsection, two assumptions on top of the exclusiveness condition and the ﬁnite support condition so that equations like (14) become
mathematically correct. The ﬁrst assumption is that “clause” bodies are mutually exclusive i.e. if there are two clauses B ← W and B ← W 0 , PDB (W ∧ W 0 | θ) = 0, and thesecond assumption is that body atoms are independent, i.e. if A ← B1 ∧ · · · ∧ Bk is a rule,
PDB (B1 ∧ · · · ∧ Bk | θ) = PDB (B1 | θ) · · · PDB (Bk | θ) holds.
Please note that “clause” used in this subsection has a special meaning. It is intended tomean G ← τ where G is a goal and τ is a tabled explanation for G obtained by OLDT searchboth of which will be explained in the next subsection.25 In other words, these additional
conditions are not imposed on a source program but on the result of OLDT search. So
clauses for auxiliary computations do not need to satisfy them.
Now suppose clauses about A occur in DB likeA ← B1,1 ∧ · · · ∧ B1,i1···
A ← BL,1 ∧ · · · ∧ BL,iL
where Bh,j (1 ≤ h ≤ L, 1 ≤ j ≤ ih ) is an atom. Theorem A.1 in Appendix A and the above
assumptions ensureβ(A) =i1Yβ(B1,j ) + · · · +

j=1

iL
Y

β(BL,j ).

(15)

j=1

(15) suggests that β(Gt ) can be considered as a function of β(A) if these equations about
inside probabilities are hierarchically organized in such a way that β(Gt ) belongs to the top
layer and any β(A) appearing on the left hand side only refers to β(B)’s which belong to thelower layers. We refer to this condition as the acyclic support condition. Under the acyclic
support condition, equations of the form (15) have a unique solution, and the computationof PDB (G | θ) via inside probabilities allows us to take advantage of reusing intermediate
results stored as β(A), thereby contributing to faster computation of PDB (Gt | θ).
P

Next we tackle a more intricate problem, the computation of S∈ψDB (Gt ) Pmsw (S |P Pθ)σi,v (S). Since the sum equals n msw(i,n,v)∈S∈ψDB (Gt ) Pmsw (S | θ), we concentrateon the computation ofX
defµ(Gt , m) =Pmsw (S | θ)m∈S∈ψDB (Gt )where m = msw(i,n,v). First we note that if an explanation S contains m like S = a1 ∧ · · · ∧ah ∧ m, then we have β(S) = β(a1 ) · · · β(ah )β(m). So µ(Gt , m) is expressed asµ(Gt , m) = α(Gt , m)β(m)

(16)

t ,m)
where α(Gt , m) = ∂µ(G
∂β(m) and α(Gt , m) does not depend on β(m). Generalizing this observation to arbitrary ground atoms, we introduce the outside probability of ground atom A
w.r.t. Gt by
def ∂β(Gt )α(Gt , A) =∂β(A)25. The logical relationship (13) corresponds to (20) where f, g and h are table atoms.

410

Parameter Learning of Logic Programs for Symbolic-statistical Modelingassuming the same conditions as inside probability. In view of (16), the problem of computing µ(Gt , m) is now reduced to that of computing α(Gt , m), which is recursively computable
as follows. Suppose A occurs in the ground program DB likeB1 ← A ∧ W1,1 , · · · , B1 ← A ∧ W1,i1···
BK ← A ∧ WK,1 , · · · , BK ← A ∧ WK,iK .
As β(Gt ) is a function of β(B1 ), . . . , β(BK ) by our assumption, the chain rule of derivatives
leads toµα(Gt , A) =∂β(Gt )∂β(B1 )¶µ∂β(A ∧ W1,1 )
∂β(A)

¶

µ

+ ··· +

∂β(Gt )
∂β(BK )¶µ∂β(A ∧ WK,iK )
∂β(A)

¶

and hence to26α(Gt , Gt ) = 1(17)α(Gt , A) = α(Gt , B1 )i1
Xβ(W1,j ) + · · · + α(Gt , BK )

j=1

iK
X

β(WK,j ).

(18)

j=1

Therefore if all inside probabilities have already been computed, outside probabilities arerecursively computed from the top (17) using (18) downward along the program layers. In
the case of DB p with f and m being chosen atoms, we computeα(f, f) = 1 α(f, g) = β(a) + β(b)α(f, h) = α(f, g)β(d)




(19)

α(f, m) = α(f, h).

From (19), the desired sum µ(f, m) is calculated asµ(f, m) = α(f, m)β(m) = (β(a) + β(b))β(d)β(m)which requires only two multiplications and one addition compared to four multiplications
and one addition in the naive computation.
Gains obtained by computing inside and outside probability may be small for this case,but as the problem size grows, they become enormous, and compensate enough for additional restrictions imposed on the result of OLDT search.
4.7 OLDT Search
To compute inside and outside probability recursively like (15) or (17) and (18), we needat programming level a tabulation mechanism for structure-sharing of partial explanations
26. Because of the independence assumption on body atoms, Wh,j (1 ≤ h ≤ K, 1 ≤ j ≤ ih ) and A are
independent. Therefore∂β(A)β(Wh,j )∂β(A ∧ Wh,j )== β(Wh,j ).
∂β(A)
∂β(A)

411

Sato & Kameya

between subgoals. We henceforth deal with programs DB in which a set table(DB) of tablepredicates are declared in advance. A ground atom containing a table predicate is calleda table atom. The purpose of table atoms is to store their support sets and eliminate theneed of recomputation, and by doing so, to construct hierarchically organized explanations
made up of the table atoms and the msw atoms.
Let DB = Fmsw ∪ R be a parameterized logic program which satisﬁes the ﬁnite supportcondition and the uniqueness condition. Also let G1 , G2 , . . . , GT be a random sample ofobservable atoms in obs(DB). We make the following additional assumptions.
Assumptions:t } of table atoms associatedFor each t (1 ≤ t ≤ T ), there exists a ﬁnite set {τ1t , . . . , τKtt (0 ≤ k ≤ K , 1 ≤ j ≤ m ) such thatwith conjunctions Sek,jtkcomp(R) `³t ∨ ··· ∨ S
etGt ↔ Se0,1
0,m0

³

´
´

³

t
t ∨ ··· ∨ Setetet∧ τ1t ↔ Se1,1
1,m1 ∧ · · · ∧ τKt ↔ SKt ,1 ∨ · · · ∨ SKt ,mK

´

(20)

t

where
t
• each Sek,j(0 ≤ k ≤ Kt , 1 ≤ j ≤ mk ) is, as a set, a subset of Fmsw ∪ {τk+1 , . . . , τKt }(acyclic support condition). As a convention, we put τ0 = Gt and call respectivelyt deft } the set of table atoms for G and Set (k ≥ 0) a t-explanationτDB= {τ0 , τ1t , . . . , τKtk,jtfor τkt .27 The set of all t-explanations for τk is denoted by ψeDB (τkt ) and we consider
ψeDB (·) as a function of table atoms.
t ∧S
et 0 ) =• t-explanations are mutually exclusive, i.e. for each k (0 ≤ k ≤ Kt ), PDB (Sek,jk,j0 (1 ≤ j 6= j 0 ≤ mk ) (t−exclusiveness condition).
t (0 ≤ k ≤ K , 1 ≤ j ≤ m ) is a conjunction of independent atoms (independent• Sek,jtkcondition).28

These assumptions are aimed at eﬃcient probability computation. Namely, the acyclic
support condition makes dynamic programming possible, the t-exclusiveness condition reduces PDB (A ∨ B) to PDB (A) + PDB (B) and the independent condition reduces PDB (A ∧ B)
to PDB (A)PDB (B). There is one more point concerning eﬃciency however. Note that the
t 29 imcomputation in dynamic programming proceeds following the partial order on τDB
posed by the acyclic support condition and access to the table atoms will be much simpliﬁedt
if they are linearly ordered. We therefore topologically sort τDBrespecting the said partialtorder and call the linearized τDB satisfying the three assumptions (the acyclic support condition, the t-exclusiveness condition and the independent condition) a hierarchical systemt = hτ t , τ t , . . . , τ t i (τ = G ) assuming ψeDB (·) isof t-explanations for Gt . We write it as τDB
0
t
0 1
Kt
30
implicitly given. Once a hierarchical system of t-explanations for Gt is successfully built
27. Preﬁx “t-” is an abbreviation of “tabled-”.
28. The independence mentioned here only concerns positive propositions. For B1 , B2 ∈ head(DB), we say
B1 and B2 are independent if PDB (B1 ∧ B2 | θ) = PDB (B1 | θ)PDB (B2 | θ) for any θ.
29. τi precedes τj if and only if the top-down execution of τi w.r.t. DB invokes τj directly or indirectly.
30. So now it holds that if τi precedes τj then i < j.

412

Parameter Learning of Logic Programs for Symbolic-statistical Modelingfrom the source program, equations on inside probability and outside probability such as(14) and (19) are automatically derived and solved in time proportional to the size of theequations. It plays a central role in our approach to eﬃcient EM learning.
One way to obtain such t-explanations is to use OLDT search (Tamaki & Sato, 1986;
Warren, 1992), a complete refutation method for logic programs. In OLDT search, whena goal G is called for the ﬁrst time, we set up an entry for G in a solution table and store
its answer substitutions Gθ there. When a call to an instance G0 of G occurs later, we stopsolving G0 and instead try to retrieve an answer substitution Gθ stored in the solution table
by unifying G0 with Gθ. To record the remaining answer substitutions of G, we prepare alookup table for G0 and hold a pointer to them.
For self-containedness, we look at details of OLDT search using a sample program
DB h = Fh ∪Rh in Figure 431 which depicts an HMM32 in Figure 3. This HMM has two states{s0, s1}. At a state transition, it probabilistically chooses the next destination from {s0, s1}

a,b
s1

s0

a,b

a,b

a,b
Figure 3: Two state HMM

Fh f1: values(init, [s0,s1]).
f2: values(out(_),[a,b]).
=
f3: values(tr(_), [s0,s1]).
h1: hmm(Cs):% To generate a string (chars) Cs,





msw(init,once,Si),%
Set initial state to Si, and thenhmm(1,Si,Cs).
% Enter the loop with clock = 1.





 h2: hmm(T,S,[C|Cs]):- T=<3, % Loop:

Rh =

















msw(out(S),T,C),
msw(tr(S),T,NextS),
T1 is T+1,
hmm(T1,NextS,Cs).
h3: hmm(T,_,[]):- T>3.

% Output C in state S.
% Transit from S to NextS.
% Put the clock ahead.
% Continue the loop (recursion).
% Finish the loop if clock > 3.

Figure 4: Two state HMM program DB h
31. f1, f2, f3, h1, h2 and h3 are temporary marks, not part of the program.
32. An HMM deﬁnes a probability distribution over strings in the given set of alphabets, and works as astochastic string generator (Rabiner & Juang, 1993) such that an output string is a sample from thedeﬁned distribution.

413

Sato & Kameya

and also an alphabet from {a, b} to emit. Note that to specify a fact set Fh and the associated distribution compactly, we introduce here a new notation values(i,[v1 ,...,vm ]). It
declares that Fh contains msw atoms of the form msw(i,n,v) (v ∈ {v1 , . . . , vm }) whose distribution is P(i,n) given by (7) in Subsection 4.2. For example, (f3), values(tr( ),[s0,s1])
introduces msw(tr(t),n,v) atoms into the program such that t can be any ground term,v ∈ {s0, s1} and for a ground term n, they have a distributionyxP(tr(t),n) (msw(tr(t),n,s0) = x, msw(tr(t),n,s1) = y | θi,s0 , θi,s1 ) = θi,s0
θi,s1

where i = tr(t), x, y ∈ {0, 1} and x + y = 1.
This program runs like a Prolog program. For a non-ground top-goal ← hmm(S), it functions as a stochastic string generator returning a list of alphabets such as [a,b,a] in thevariable S as follows. The top-goal calls clause (h1) and (h1) selects the initial state by executing subgoal msw(init,once,Si)33 which returns in Si an initial state probabilisticallychosen from {s0, s1}. The second clause (h2) is called from (h1) with ground S and ground
T. It makes a probabilistic choice of an output alphabet C by asking msw(out(S),T,C) andthen determines NextS, the next state, by asking msw(tr(S),T,NextS). (h3) is there tostop the transition. For simplicity, the length of output strings is ﬁxed to three. This way
of execution is termed as sampling execution because it corresponds to a random samplingfrom PDB h . If the top-goal is ground like ← hmm([a,b,a]), it works as an acceptor, i.e.
returning success (yes) or failure (no).
If all explanations for hmm([a,b,a]) are sought for, we keep all msw atoms resolved upon
during the refutation as a conjunction (explanation), and repeat this process by backtracking until no more refutation is found. If we need t-explanations however, backtracking must
be abandoned because sharing of partial explanations through t-explanations, the purposeof t-explanations itself, becomes impossible. We therefore instead use OLDT search for all
t1: top_hmm(Cs,Ans):- tab_hmm(Cs,Ans,[]).
t2: tab_hmm(Cs,[hmm(Cs)|X],X):- hmm(Cs,_,[]).
t3: tab_hmm(T,S,Cs,[hmm(T,S,Cs)|X],X):- hmm(T,S,Cs,_,[]).
t4: e_msw(init,T,s0,[msw(init,T,s0)|X],X).
t4’: e_msw(init,T,s1,[msw(init,T,s1)|X],X).
:t7: hmm(Cs,X0,X1):- e_msw(init,once,Si,X0,X2), tab_hmm(1,Si,Cs,X2,X1).
t8: hmm(T,S,[C|Cs],X0,X1):T=<3, e_msw(out(S),T,C,X0,X2), e_msw(tr(S),T,NextS,X2,X3),
T1 is T+1, tab_hmm(T1,NextS,Cs,X3,X1).
t9: hmm(T,S,[],X,X):- T>3.
Figure 5: Translated program of DB h
33. If msw(i,n,V) is called with ground i and ground n, V, a logical variable, behaves like a random variable.
It is instantiated to some term v with probability θi,v selected from the value set Vi declared by a values
atom. If, on the other hand, V is a ground term v when called, the procedural semantics of msw(i,n,v)
is equal to that of msw(i,n,V) ∧ V = v.

414

Parameter Learning of Logic Programs for Symbolic-statistical Modelingt-explanation search. In the case of the HMM program for example, to build a hierarchicalsystem of t-explanations for hmm([a,b,a]) by OLDT search, we ﬁrst declare hmm/1 andhmm/3 as table predicate.34 So a t-explanation will be a conjunction of hmm/1 atoms, hmm/3
atoms and msw atoms. We then translate the program into another logic program, analogously to the translation of deﬁnite clause grammars (DCGs) in Prolog (Sterling & Shapiro,1986). We add two arguments (which forms a D-list) to each predicate for the purpose ofaccumulating msw atoms and table atoms as conjuncts in a t-explanation. The translation
applied to DB h yields the program in Figure 5.
In the translated program, clause (t1) corresponds to the top-goal ← hmm(l) with an
input string l, and a t-explanation for the table atom hmm(l) will be returned in Ans. (t2)and (t3) are auxiliary clauses to add to the callee’s D-list a table atom of the form hmm(l)and hmm(t,s,l) respectively (t: time step, s: state). In general, if p/n is a table predicate
in the original program, p/(n + 2) becomes a table predicate in the translated program andan auxiliary predicate tab p/(n + 2) is inserted to signal the OLDT interpreter to check thesolution table for p/n, i.e. to check if there already exist t-explanations for p/n. Likewiseclauses (t4) and (t4’) are a pair corresponding to (f1) which insert msw(init,T,·) to thecallee’s D-list with T = once. Clauses (t7), (t8) and (t9) respectively correspond to (h1),(h2) and (h3).
hmm([a,b,a]):[hmm([a,b,a])][ [msw(init,once,s0), hmm(1,s0,[a,b,a])],[msw(init,once,s1), hmm(1,s1,[a,b,a])] ]hmm(1,s0,[a,b,a]):[hmm(1,s0,[a,b,a])][ [msw(out(s0),1,a), msw(tr(s0),1,q0), hmm(2,s0,[b,a])],[msw(out(s0),1,a), msw(tr(s0),1,s1), hmm(2,s1,[b,a])] ]hmm(1,s1,[a,b,a]):[hmm(1,s1,[a,b,a])][ [msw(out(s1),1,a), msw(tr(s1),1,s0), hmm(2,s0,[b,a])],[msw(out(s1),1,a), msw(tr(s1),1,s1), hmm(2,s1,[b,a])] ]hmm(2,s0,[b,a]):[hmm(2,s0,[b,a])][ [msw(out(s0),2,b), msw(tr(s0),2,s0), hmm(3,s0,[a])],[msw(out(s0),2,b), msw(tr(s0),2,s1), hmm(3,s1,[a])] ]hmm(2,s1,[b,a]):[hmm(2,s1,[b,a])][ [msw(out(s1),2,b), msw(tr(s1),2,s0), hmm(3,s0,[a])],[msw(out(s1),2,b), msw(tr(s1),2,s1), hmm(3,s1,[a])] ]hmm(3,s0,[a]):[hmm(3,s0,[a])][ [msw(out(s0),3,a), msw(tr(s0),3,s0), hmm(4,s0,[])],[msw(out(s0),3,a), msw(tr(s0),3,s1), hmm(4,s1,[])] ]hmm(3,s1,[a]):[hmm(3,s1,[a])][ [msw(out(s1),3,a), msw(tr(s1),3,s0), hmm(4,s0,[])],[msw(out(s1),3,a), msw(tr(s1),3,s1), hmm(4,s1,[])] ]hmm(4,s0,[]):[hmm(4,s0,[])][[]]
hmm(4,s1,[]):[hmm(4,s1,[])][[]]

Figure 6: Solution table
34. In general, p/n means a predicate p with arity n. So although hmm/1 and hmm/3 share the predicate namehmm, they are diﬀerent predicates.

415

Sato & Kameya

Then after translation, we apply OLDT search to ← top hmm([a,b,a],Ans) while noting (i) the added D-list does not inﬂuence the OLDT procedure, and (ii) we associate witheach solution of a table atom in the solution table a list of t-explanations. The resulting
solution table is shown in Figure 6. The ﬁrst row reads that a call to hmm([a,b,a]) occurred and entered the solution table and its solution, hmm([a,b,a]) (no variable binding generated), has two t-explanations, msw(init,once,s0) ∧ hmm(1,s0,[a,b,a]) andmsw(init,once,s1) ∧ hmm(1,s1,[a,b,a]). The remaining task is the topological sorting of the table atoms stored in the solution table respecting the acyclic support condition.
This can be done by using depth-ﬁrst search (trace) of t-explanations from the top-goal for
example. Thus we obtain a hierarchical system of t-explanations for hmm([a,b,a]).
4.8 Support Graphs
Looking back, all we need to compute inside and outside probability is a hierarchical systemof t-explanations, which essentially is a boolean combination of primitive events (msw atoms)
and compound events (table atoms) and as such can be more intuitively representable as agraph. For this reason, and to help visualizing our learning algorithm, we introduce a new
data-structure termed support graphs, though the new EM algorithm in the next subsection
itself is described solely by the hierarchical system of t-explanations.
As illustrated in Figure 7 (a), the support graph for Gt is a graphical representation oftt i (τ t = G ) for G in (20).
the hierarchical system of t-explanations τDB= hτ0t , τ1t , . . . , τKtt
0t
It consists of totally ordered disconnected subgraphs, each of which is labeled with thetcorresponding table atom τkt in τDB(0 ≤ k ≤ Kt ). A subgraph labeled τkt comprises two
special nodes (the start node and the end node) and explanation graphs, each correspondingt in ψeDB (τ t ) (1 ≤ j ≤ m ).
to a t-explanation Sek,jkkte
An explanation graph of Sk,j is a linear graph in which a node is labeled either with at . They are called a table node and a switch
table atom τ or with a switch msw(·,·,·) in Sek,jnode respectively. Figure 7 (b) is the support graph for hmm([a,b,a]) obtained from thesolution table in Figure 6. Each table node labeled τ refers to the subgraph labeled τ , sodata-sharing is achieved through the distinct table nodes referring to the same subgraph.
4.9 Graphical EM Algorithm
We describe here an eﬃcient EM learning algorithm termed the graphical EM algorithm(Figure 8) introduced by Kameya and Sato (2000), that runs on support graphs. Suppose
we have a random sample G = G1 , . . . , GT of observable atoms. Also suppose support
graphs for Gt (1 ≤ t ≤ T ), i.e. hierarchical systems of t-explanations satisfying the acyclic
support condition, the t-exclusiveness condition and the independent condition, have beensuccessfully constructed from a parameterized logic program DB satisfying the uniquenesscondition and the ﬁnite support condition.
The graphical EM algorithm reﬁnes learn-naive(DB,G) by introducing two subroutines,
get-inside-probs(DB, G) to compute inside probabilities and get-expectations(DB, G) to compute outside probabilities. They are called from the main routine learn-gEM(DB,G). Whenlearning, we prepare four arrays for each support graph for Gt in G:• P[t, τ ] for the inside probability of τ , i.e. β(τ ) = PDB (τ | θ) (see (12))
416

Parameter Learning of Logic Programs for Symbolic-statistical Modeling

τk

(a)

explanation graph

msw

Gt:
start

τk

msw

msw

msw

end

msw

:
τk :
start

end
msw

msw

:

(b)
msw(init,once,s0)

hmm(1,s0,[a,b,a])hmm([a,b,a]):

start

end

msw(init,once,s1)

msw(out(s0),1,a)hmm(1,s1,[a,b,a])msw(tr(s0),1,s0)

hmm(2,s0,[b,a])hmm(1,s0,[a,b,a]):

end

start

msw(out(s0),1,a)msw(tr(s0),1,s1)

hmm(2,s1,[b,a])msw(out(s1),1,a)msw(tr(s1),1,s0)

hmm(2,s0,[b,a])hmm(1,s1,[a,b,a]):

end

start

msw(out(s1),1,a)msw(tr(s1),1,s1)

hmm(2,s1,[b,a])

Figure 7: A support graph (a) in general form, (b) for Gt = hmm([a,b,a]) in the HMM
program DB h . A double-circled node refers to a table node.

• Q[t, τ ] for the outside probability of τ w.r.t. Gt , i.e. α(Gt , τ ) (see (17) and (18))
e for the explanation probability of Se (∈ ψeDB (τ t )), i.e. PDB (Se | θ)
• R[t, τ, S]k417

Sato & Kameya

1: procedure learn-gEM (DB, G)
2: begin
3:
Select some θ as initial
4:
5:
6:
7:
8:
9:

parameters;
get-inside-probs(DB,G);P
λ(0) := Tt=1 ln P[t, Gt ];1: procedure get-inside-probs(DB, G)
2: begin3:
for t := 1 to T do begin
4:
Let τ0t = Gt ;5:
for k := Kt downto 0 do begin
6:
P[t, τkt ] := 0;

repeat7:

get-expectations(DB, G);8:

foreach i ∈ I, v ∈ Vi do

9:η[i, v]P:=T
t=1 η[t, i, v]/P[t, Gt ];foreach i ∈ I, v ∈ Vi doP
θi,v := η[i, v]/ v0 ∈Vi η[i, v 0 ];

10:
11:
12:

10:
11:
12:
get-inside-probs(DB, G);13:m := m +1;
P
14:
λ(m) := Tt=1 ln P[t, Gt ]15:
until λ(m) − λ(m−1) < ε
16: end.

13:
14:

foreach Se ∈ ψeDB (τkt ) do begin
Let Se = {A1 , A2 , . . . , A|Se| };e := 1;R[t, τkt , S]e dofor l := 1 to |S|

if Al = msw(i,·,v) thene ∗ = θi,v
R[t, τkt , S]e ∗ = P[t, Al ];else R[t, τkt , S]te
P[t, τk ] += R[t, τkt , S]
end /* foreach Se */15:
16:
end /* for k */17:
end /* for t */18: end.

1: procedure get-expectations(DB, G) begin
2:
for t := 1 to T do begin
3:
foreach i ∈ I, v ∈ Vi do η[t, i, v] := 0;
4:
Let τ0t = Gt ; Q[t, τ0t ] := 1.0;
5:
for k := 1 to Kt do Q[t, τkt ] := 0;
6:
for k := 0 to Kt do
7:
foreach Se ∈ ψeDB (τkt ) do begin
8:

Let Se = {A1 , A2 , . . . , A|Se| };e dofor l := 1 to |S|e
10:
if Al = msw(i,·,v) then η[t, i, v] += Q[t, τkt ] · R[t, τkt , S]
t
t
e
11:
else Q[t, Al ] += Q[t, τk ] · R[t, τk , S]/P[t, Al ]12:
end /* foreach Se */13:
end /* for t */14: end.
9:

Figure 8: graphical EM algorithm.
• η[t, i, v] for the expected count of msw(i,·,v), i.e.

P

S∈ψDB (Gt ) Pmsw (S | θ)σi,v (S)and call the procedure learn-gEM(DB,G) in Figure 8. The main routine learn-gEM(DB,G) initially computes all inside probabilities (Line 4) and enters a loop in which get-expectations(DB,G)
is called ﬁrst to compute the expected count η[t, i, v] of msw(i,·,v) and parameters are updated (Line 11). Inside probabilities are renewed by using the updated parameters before
entering the next loop (Line 12).
418

Parameter Learning of Logic Programs for Symbolic-statistical Modeling

The subroutine get-inside-probs(DB,G) computes the inside probability β(τ ) = PDB (τ | θ)(and stores it in P[t, τ ]) of a table atom τ from the bottom layer up to the topmost layer τ0 =Gt (Line 4) of the hierarchical system of t-explanations for Gt (see (20) in Subsection 4.6).
It takes a t-explanation Se in ψeDB (τkt ) one by one (Line 7), decomposes Se into conjuncts andmultiplies their inside probabilities which are either known (Line 12) or already computed(Line 13).
The other subroutine get-expectations(DB,G) computes outside probabilities following therecursive deﬁnitions (17) and (18) in Subsection 4.6 and stores the outside probabilityα(Gt , τ ) of a table atom τ in Q[t, τ ]. It ﬁrst sets the outside probability of the top-goalτ0 = Gt to 1.0 (Line 4) and computes the rest of outside probabilities (Line 6) going down
the layers of the t-explanation for Gt described by (20) in Subsection 4.6. (Line 10) adds
e = α(Gt , τ t ) · β(S)e to η[t, i, v], the expected count of msw(i,·,v), as aQ[t, τkt ] · R[t, τkt , S]kcontribution of msw(i,·,v) in Se through τkt to η[t, i, v]. (Line 11) increments the outside
probability Q[t, Al ] = α(Gt , Al ) of Al according to the equation (18). Notice that Q[t, τkt ]e
has already been computed and R[t, τkt , S]/P[t,Al ] = β(W ) for Se = Al ∧ W . As shown in
Subsection 4.5, learn-naive(DB,G) is the MLE procedure, hence the following theorem holds.
Theorem 4.1 Let DB be a parameterized logic program, and G = G1 , . . . , GT a random sample of observable atoms. Suppose the ﬁve conditions (uniqueness, ﬁnite support(Subsection 4.2), acyclic support, t-exclusiveness and independence (Subsection 4.7)) are
met. Then learn-gEM (DB, G) ﬁnds the MLE θ ∗ which (locally) maximizes the likelihood
Q
L(G | θ) = Tt=1 PDB (Gt | θ).
(Proof) Sketch.35 Since the main routine learn-gEM(DB,G) is the same as learn-naive(DB,G)P
P
except the computation of η[i, v] = Tt=1 η[t, i, v], we show that η[t, i, v] = S∈ψDB (Gt ) Pmsw (S |P Pθ)σi,v (S) (= n msw(i,n,v)∈S∈ψDB (Gt ) Pmsw (S | θ)). However,η[t, i, v] =

X

X

X

e
α(Gt , τkt )β(S)0≤k≤Kt n msw(i,n,v)∈Se∈ψeDB (τkt )(see (Line 10) in get-expectations(DB, G))=X

α(Gt , msw(i,n,v))β(msw(i,n,v))

n

=

X
n

=

X

µ(Gt , msw(i,n,v)) (see the equation (16))
X

Pmsw (S | θ).

Q.E.D.

n msw(i,n,v)∈S∈ψDB (Gt )e =Here we used the fact that if Se contains msw(i,n,v) like Se = Se0 ∧ msw(i,n,v), β(S)0β(Se )β(msw(i,n,v)) holds, and hencee = α(Gt , τ t )β(Se0 )β(msw(i,n,v))α(Gt , τkt )β(S)k= (contribution of msw(i,n,v) in Se through τkt to α(Gt , msw(i,n,v)))β(msw(i,n,v)).

35. A formal proof is given by Kameya (2000). It is proved there that under the common parameters θ, η[i, v]
in learn-naive(DB,G) coincides with η[i, v] in learn-gEM(DB,G) . So, the parameters are updated tothe same values. Hence, starting with the same initial values, the parameters converge to the samevalues.

419

Sato & Kameya

The ﬁve conditions on the applicability of the graphical EM algorithm may look hardto satisfy at once. Fortunately, the modeling principle in Section 4.3 still stands, and withdue care in modeling, it is likely to lead us to a program that meets all of them. Actually,
we will see in the next section, programs for standard symbolic-statistical frameworks suchas Bayesian networks, HMMs and PCFGs all satisfy the ﬁve conditions.

5. Complexity
In this section, we analyze the time complexity of the graphical EM algorithm applied
to various symbolic-statistical frameworks including HMMs, PCFGs, pseudo PCSGs andBayesian networks. The results show that the graphical EM algorithm is competitive withthese specialized EM algorithms developed independently in each research ﬁeld.
5.1 Basic Property
Since the EM algorithm is an iterative algorithm and since we are unable to predict whenit converges, we measure time complexity by the time taken for one iteration. We therefore
estimate time per iteration on the repeat loop of learn-gEM (DB, G) (G = G1 , . . . , GT ). We
observe that in one iteration, each support graph for Gt (1 ≤ t ≤ T ) is scanned twice, onceby get-inside-probs(DB, G) and once by get-expectations(DB, G). In the scan, addition is
performed on the t-explanations, and multiplication (possibly with division) is performedon the msw atoms and table atoms once for each. So time spent for Gt per iteration by thegraphical EM algorithm is linear in the size of the support graph, i.e. the number of nodes
in the support graph for Gt . Put
et
∆
DB

def

=

[

ψeDB (τ )tτ ∈τDBξnum

defξmaxsize

def

=
=

et |
max |∆
DB

1≤t≤T

max

e∈∆
e tDB
1≤t≤T,Se|S|.

t is the set of table atoms for G , and hence ∆e t is the set of all t-explanations
Recall that τDBtDB
appearing in the right hand side of (20) in Subsection 4.7. So ξnum is the maximum numberof t-explanations in a support graph for the Gt ’s and ξmaxsize the maximum size of a texplanation for the Gt ’s respectively. The following is obvious.

Proposition 5.1 The time complexity of the graphical EM algorithm per iteration is linear
in the total size of support graphs, O(ξnum ξmaxsize T ) in notation, which coincides with thespace complexity because the graphical EM algorithm runs on support graphs.
This is a rather general result, but when we compare the graphical EM algorithm withother EM algorithms, we must remember that the input to the graphical EM algorithm issupport graphs (one for each observed atom) and our actual total learning time is
OLDT time + (the number of iterations) ×O(ξnum ξmaxsize T )420

Parameter Learning of Logic Programs for Symbolic-statistical Modelingwhere “OLDT time” denotes time to construct all support graphs for G. It is the sum oftime for OLDT search and time for the topological sorting of the table atoms, but because
the latter is part of the former order-wise,36 we represent “OLDT time” by time for OLDT
search. Also observe that the total size of support graphs does not exceed time for OLDT
search for G order-wise.
To evaluate OLDT time for a speciﬁc class of models such as HMMs, we need to know
time for table operations. Observe that our OLDT search in this paper is special in thesense that table atoms are always ground when called and there is no resolution with solved
goals. Accordingly a solution table is used only
• to check if a goal G already has an entry in the solution table, i.e. if it was called
before, and• to add a new searched t-explanation for G to the list of discovered t-explanations
under G’s entry.
The time complexity of these operations is equal to that of table access which depends
both on the program and on the implementation of the solution table.37 We ﬁrst suppose
programs are carefully written in such a way that the arguments of table atoms used as indecies for table access are integers. Actually all programs used in the subsequent complexity
analysis (DB h in Subsection 4.7, DB g and DB g0 in Subsection 5.3, DB Gτ in Subsection 5.5)
satisfy or can satisfy this condition by replacing non-integer terms with appropriate integers. We also suppose that the solution table is implemented using an array so that table
access can be done in O(1) time.38
In what follows, we present a detailed analysis of the time complexity of the graphicalEM algorithm applied to HMMs, PCFGs, pseudo PCSGs and Bayesian networks, assumingO(1) time access to the solution table. We remark by the way that their space complexity
is just the total size of solution tables (support graphs).
5.2 HMMs
The standard EM algorithm for HMMs is the Baum-Welch algorithm (Rabiner, 1989; Rabiner & Juang, 1993). An example of HMM is shown in Figure 3 in Subsection 4.7.39 Given
T observations w1 , . . . , wT of output string of length L, it computes in O(N 2 LT ) time int (q) = P (ot ot · · · ot
each iteration the forward probability αm
1 2
m−1 , q | θ) and the backwardtttt
probability βm (q) = P (om om+1 · · · oL | q, θ) for each state q ∈ Q, time step m (1 ≤ m ≤ L)
and a string wt = ot1 ot2 · · · otL (1 ≤ t ≤ T ), where Q is the set of states and N the number ofstates. The factor N 2 comes from the fact that every state has N possible destinations and
36. Think of OLDT search for a top-goal Gt . It searches for msw atoms and table atoms to create a solution table, while doing some auxiliary computations. Therefore its time complexity is never lessthan O(|the number of msw atoms and table atoms in the support graph for Gt |), which coincides withthe time we need to topologically sort table atoms in the solution table by depth-ﬁrst search from τ0 = Gt .
37. Sagonas et al. (1994) and Ramakrishnan et al. (1995) discuss about the implementation of OLDT.
38. If arrays are not available, we may be able to use balanced trees, giving O(log n) access time where nis the number data in the solution table, or we may be able to use hashing, giving average O(1) time
access under a certain condition (Cormen, Leiserson, & Rivest, 1990).
39. We treat here only “state-emission HMMs” which emit a symbol depending on the state. Another type,“arc-emission HMMs” in which the emitted symbol depends on the transition arc, is treated similarly.

421

Sato & Kameya

we have to compute the forward and backward probability for every destination and everyt (q)’s and β t (q)’s, parameters are updated. So, the total
state. After computing all αmmcomputation time in each iteration of the Baum-Welch algorithm is estimated as O(N 2 LT )(Rabiner & Juang, 1993; Manning & Schütze, 1999).
To compare this result with the graphical EM algorithm, we use the HMM program
DB h in Figure 4 with appropriate modiﬁcations to L, the length of a string, Q, thestate set, and declarations in Fh for the output alphabets. For a string w = o1 o2 · · · oL ,hmm(n,q,[om , om+1 , . . . , oL ]) in DB h reads that the HMM is in state q ∈ Q at time n andhas to output [om ,om+1 ,...,oL ] until it reaches the ﬁnal state. After declaring hmm/1 andhmm/3 as table predicate and translation (see Figure 5), we apply OLDT search to the goal
← top hmm([o1 ,...,oL ],Ans) w.r.t. the translated program to obtain all t-explanations
for hmm([o1 ,...,oL ]). For a complexity argument however, the translated program andDB h are the same, so we talk in terms of DB h for the sake of simplicity. In the search,
we ﬁx the search strategy to multi-stage depth-ﬁrst strategy (Tamaki & Sato, 1986). We
assume that the solution table is accessible in O(1) time.40 Since the length of the list in
the third argument of hmm/3 decreases by one on each recursion, and there are only ﬁnitely
many choices of the state transition and the output alphabet, the search terminates, leavingﬁnitely many t-explanations in the solution table like Figure 6 that satisfy the acyclic support condition respectively. Also the sampling execution of ← hmm(L) w.r.t. DB h is nothing
but a sequential decision process such that decisions made by msw atoms are exclusive,
independent and generate a unique string, which means DB h satisﬁes the t-exclusivenesscondition, the independence condition and the uniqueness condition respectively. So, thegraphical EM algorithm is applicable to the set of hierarchical systems of t-explanations forhmm(wt ) (1 ≤ t ≤ T ) produced by OLDT search for T observations w1 , . . . , wT of output
string. Put wt = ot1 ot2 · · · otL . It follows fromt
= {hmm(m,q,[otm ,...,otL ]) | 1 ≤ m ≤ L + 1, q ∈ Q} ∪ {hmm([ot1 ,...,otL ])}
τDB
h

(

ψeDB h (hmm(m,q,[otm ,...,otL ])) =(1 ≤ m ≤ L)¯

msw(out(q),m,om ), msw(tr(q),m,q 0 ), ¯¯ 0
¯q ∈ Q¯
hmm(m + 1,q 0 ,[otm+1 ,...,otL ]))that for a top-goal hmm([ot1 ,...,otL ]), there are at most O(N L) calling patterns of hmm/3
and each call causes at most N calls to hmm/3, implying there occur O(N L · N ) = O(N 2 L)
calls to hmm/3. Since each call is computed once due to the tabling mechanism, we haveξnum = O(N 2 L). Also ξmaxsize = 3. Applying Proposition 5.1, we reach
Proposition 5.2 Suppose we have T strings of length L. Also suppose each table operationin OLDT search is done in O(1) time. OLDT time by DB h is O(N 2 LT ) and the graphicalEM algorithm takes O(N 2 LT ) time per iteration where N is the number of states.
O(N 2 LT ) is the time complexity of the Baum-Welch algorithm. So the graphical EM
algorithm runs as eﬃciently as the Baum-Welch algorithm.41
40. O(1) is possible because in the translated program DB h in Section 4.7, we can identify a goal pattern ofhmm(·,·,·,·,·) by the ﬁrst two arguments which are constants (integers).
41. Besides, the Baum-Welch algorithm and the graphical EM algorithm whose input are support graphs
generated by DB h update parameters to the same value if initial values are the same.

422

Parameter Learning of Logic Programs for Symbolic-statistical Modeling

By the way, the Viterbi algorithm (Rabiner, 1989; Rabiner & Juang, 1993) provides for
HMMs an eﬃcient way of ﬁnding the most likely transition path for a given input/output
string. A similar algorithm for parameterized logic programs that determines the most
likely explanation for a given goal can be derived. It runs in time linear in the size of the
support graph, thereby O(N 2 L) in the case of HMMs, the same complexity as the Viterbi
algorithm (Sato & Kameya, 2000).
5.3 PCFGs
We now compare the graphical EM algorithm with the Inside-Outside algorithm (Baker,
1979; Lari & Young, 1990). The Inside-Outside algorithm is a well-known EM algorithm
for PCFGs (Wetherell, 1980; Manning & Schütze, 1999).42 It takes a grammar in Chomsky
normal form. Given N nonterminals, a production rule in the grammar takes the formi → j, k (1 ≤ i, j, k ≤ N ) (nonterminals are named by numbers from 1 to N and 1 is a
starting symbol) or the form i → w where 1 ≤ i ≤ N and w is a terminal. In each iteration,
it computes the inside probability and the outside probability of every partial parse treeof the given sentence to update parameters for these production rules. Time complexity ismeasured by time per iteration, and is described by N , the number of nonterminals, andL, the number of terminals in a sentence. It is O(N 3 L3 T ) for T observed sentences (Lari& Young, 1990).
To compare the graphical EM algorithm with the Inside-Outside algorithm, we startfrom a propositional program DB g = Fg ∪ Rg below representing the largest grammar
containing all possible rules i → j, k in N nonterminals where nonterminal 1 is a starting
symbol, i.e. sentence.

Fg

Rg

=

{msw(i,[d,d0 ],[j,k]) | 1 ≤ i, j, k ≤ N, d, d0 are numbers}∪ {msw(i,d,w) | 1 ≤ i ≤ N, d is a number, w is a terminal}


¯

¯


 q(i,d0 ,d2 ) :- msw(i,[d0 ,d2 ],[j,k]), ¯ 1 ≤ i, j, k ≤ N,¯
q(j,d0 ,d1 ),=¯¯ 0 ≤ d0 < d1 < d2 ≤ L ¯q(k,d1 ,d2 ).
S n

¯
¯

o

q(i,d,d+1) :- msw(i,d,wd+1 ). ¯ 1 ≤ i ≤ N, 0 ≤ d ≤ L − 1
Figure 9: PCFG program DB g

DB g is an artiﬁcial parsing program whose sole purpose is to measure the size of an
OLDT tree43 created by the OLDT interpreter when it parses a sentence w1 w2 · · · wL . So
42. A PCFG (probabilistic context free grammar) is a backbone CFG with probabilities (parameters) assigned to each production rule. For a nonterminal A having
Pn n production rules {A → αi | 1 ≤ i ≤ n}, aprobability pi is assigned to A → αi (1 ≤ i ≤ n) wherep = 1. The probability of a sentence s is
i=1 i
the sum of probabilities of each (leftmost) derivation of s. The latter is the product of probabilities ofrules used in the derivation.
43. To be more precise, an OLDT structure, but in this case, it is a tree because DB g contains only constants(Datalog program) and there never occurs the need of creating a new root node.

423

Sato & Kameya

(1)

Td

q(1,d,L)
2 j N

2q(1,d,d+1),q(1,d+1,L)q(1,d+1,L)
(1)

Td+1
[Note] q
1

i

k

N

2

k

N

q(1,d,d+1),q(k,d+1,L)q(j,d,d+1),q(1,d+1,L)q(j,d,d+1),q(k,d+1,L)q(k,d+1,L)q(1,d+1,L)q(k,d+1,L)

(k)

Td+1
q(i,d’,d’’) already appears
d+1 d’ d’’ L, 1 id’’-d’ L-d-2qNd+2 e L-1
1 j N

2 k
q(j,d,e),q(1,e,L)N
q(j,d,e),q(k,e,L)d+2 e’ e
1 i’ N,
1 j’ N,q(i’,d,e’),q(j’,e’,e),q(1,e,L)Td’q(k,e,L)N
q(j’,e’,e),q(1,e,L)q(1,e,L)...
p(i)p(1) p(2) p(N)Figure 10: OLDT tree for the query ← q(1,d,L)
the input sentence w1 w2 · · · wL is embedded in the program separately as msw(i,d,wd+1 )(0 ≤ d ≤ L−1) in the second clauses of Rg (this treatment does not aﬀect the complexity argument). q(i,d0 ,d1 ) reads that the i-th nonterminal spans from position d0 to position d1 ,i.e. the substring wd0 +1 · · · wd1 . The ﬁrst clauses q(i,d0 ,d2 ) :- msw(·,·,·), q(j,d0 ,d1 ),q(k,d1 ,d2 ) are supposed to be textually ordered according to the lexicographic order fortuples hi, j, k, d0 , d2 , d1 i. As a parser, the top-goal is set to ← q(1,0,L).44 It asks theparser to parse the whole sentence w1 w2 · · · wL as the syntactic category “1” (sentence).
We make an exhaustive search for this query by OLDT search.45 As before, the multistage depth-ﬁrst search strategy and O(1) time access to the solution table are assumed.
Then the time complexity of OLDT search is measured by the number of nodes in the(k)(1)
OLDT tree. Let Td be the OLDT tree for ← q(k,d,L). Figure 10 illustrates Td for d(0 ≤ d ≤ L − 3) where msw atoms are omitted. As can be seen, the tree has many similar
subtrees, so we put them together (see Note in Figure 10). Due to the depth-ﬁrst strategy,(1)(1)
Td has a recursive structure and contains Td+1 as a subtree. Nodes whose leftmost atom
is not underlined are solution nodes, i.e. they solve their leftmost atoms for the ﬁrst time inthe entire refutation process. The underlined atoms are already computed in the subtrees
to their left.46 They only check the solution table if there are their entries (= already
44. L here is not a Prolog variable but a constant denoting the sentence length.
45. q is a table predicate.
(1)
46. It can be inductively proved that Td+1 contains every computed q(i,d0 ,d00 ) (0 ≤ d ≤ L − 3, d + 1 ≤ d0 <00000d ≤ L, 1 ≤ i ≤ N, d − d ≤ L − d − 2).

424Parameter Learning of Logic Programs for Symbolic-statistical Modelingcomputed) in O(1) time. Since all clauses are ground, such execution only generates asingle child node.
(1)
(1)
(k)
We enumerate hd , the number of nodes in Td but not in Td+1 (1 ≤ k ≤ N ). From
(1)

(k)

Figure 10, we see hd = O(N 3 (L − d)2 ).47 Let hd (2 ≤ k ≤ N ) be the number of nodes in(k)(1)
Td+1 not contained in Td+1 . It is estimated as O(N 2 (L − d − 2)). Consequently, the number
(1)

(1)

P

(k)

of nodes that are newly created in Td is hd + N
h = O(N 3 (L − d)2 ). As a result,PL−3 k=2 d 3 3 48
total time for OLDT search is computed as d=0 hd = O(N L ) which is also the size ofthe support graph.
We now consider a non-propositional parsing program DB g0 = Fg0 ∪ Rg0 in Figure 11
whose ground instances constitute the propositional program DB g . DB g0 is a probabilistic
variant of DCG program (Pereira & Warren, 1980) in which q’/1, q’/6 and between/3 are
declared as table predicate. Semantically DB g0 speciﬁes a probability distribution over theatoms of the form {q’(l) | l is a list of terminals}.

Fg0{msw(si ,t,[sj ,sk ]) | 1 ≤ i, j, k ≤ N, t is a number}∪ {msw(si ,t,w) | 1 ≤ i ≤ N, t is a number, w is a terminal}=
q’(S) :- length(S,D), q’(s1 ,0,D,0, ,S-[]).









q’(I,D0,D2,C0,C2,L0-L2) :- between(D0,D1,D2),




Rg 0

=















msw(I,C0,[J,K]),q’(J,D0,D1,s(C0),C1,L0-L1),q’(K,D1,D2,s(C1),C2,L1-L2).

q’(I,D,s(D),C0,s(C0),[W|X]-X) :- msw(I,C0,W).
Figure 11: Probabilistic DCG like program DB g0

The top-goal to parse a sentence S = [w1 , . . . , wL ] is ← q’([w1 , . . . , wL ]). It invokesq’(s1 ,0,D,0, ,[w1 ,. . .,wL ]-[]) after measuring the length D of an input sentence S bycalling length/2. 49 50 In general, q’(i,d0 ,d2 ,c0 ,c2 ,l0 -l2 ) works identically to q(i,d0 ,d2 )but three arguments, c0 , c2 and l0 -l2 , are added. c0 supplies a unique trial-id for msws to be
used in the body, c2 the latest trial-id in the current computation, and l0 -l2 a D-list holding
a substring between d0 and d2 . Since the added arguments do not aﬀect the shape of the¯¯

47. We here focus on the subtree Td0 . j, i0 and j 0 range from 1 to N , and ¯{(e, e0 ) | d + 2 ≤ e0 < e ≤ L − 1 }¯ =O((L − d)2 ). Hence, the number of nodes in Td0 is O(N 3 (L − d)2 ). The number of nodes in Td but(1)(1)
neither in Td+1 nor in Td0 is negligible, therefore hd = O(N 3 (L − d)2 ).
(1)(1)
48. The number of nodes in TL−1 and TL−2 is negligible.
49. To make the program as simple as possible, we assume that an integer n is represented by a ground term
(1)

(n)

def

z }| {

sn = s(· · ·s (0)· · ·). We also assume that when D0 and D2 are ground, the goal ← between(D0, D1, D2)
returns an integer D1 between them in time proportional to |D1 − D0|.
50. We omit an obvious program for length(l,sn ) which computes the length sn of a list l in O(|l|) time.

425

Sato & Kameya

search tree in Figure 10 and the extra computation caused by length/2 is O(L) and theone by the insertion of between(D0,D1,D2) is O(N L3 ) respectively,51 OLDT time remains
O(N 3 L3 ), and hence so is the size of the support graph.
To apply the graphical EM algorithm correctly, we need to conﬁrm the ﬁve conditions
on its applicability. It is rather apparent however that the OLDT refutation of any topgoal of the form ← q’([w1 ,. . .,wL ]) w.r.t. DB g0 terminates, and leaves a support graphsatisfying the ﬁnite support condition and the acyclic support condition. The t-exclusivenesscondition and the independent condition also hold because the refutation process faithfully
simulates the leftmost stochastic derivation of w1 · · · wL in which the choice of a production
rule made by msw(si ,sc ,[sj ,sk ]) is exclusive and independent (trial-ids are diﬀerent on
diﬀerent choices).
What remains is the uniqueness condition. To conﬁrm it, let us consider another program DB g00 , a modiﬁcation of DB g0 such that the ﬁrst goal length(S,D) in the body of theﬁrst clause and the ﬁrst goal between(D0,D1,D2) in the second clause of Rg0 are moved to
the last position in their bodies respectively. DB g00 and DB g0 are logically equivalent, andsemantically equivalent as well from the viewpoint of distribution semantics. Then think of
the sampling execution by the OLDT interpreter of a top-goal ← q’(S) w.r.t. DB g00 where
S is a variable, using the multi-stage depth-ﬁrst search strategy. It is easy to see ﬁrst thatthe execution never fails, and second that when the OLDT refutation terminates, a sentence[w1 , . . . , wL ] is returned in S, and third that conversely, the set of msw atoms resolved uponin the refutation uniquely determines the output sentence [w1 , . . . , wL ].52 Hence, if thesampling execution is guaranteed to always terminate, every sampling from PFg00 (= PFg0 )uniquely generates a sentence, an observable atom, so the uniqueness condition is satisﬁedby DB g00 , and hence by DB g0 .
Then when is the sampling execution guaranteed to always terminate? In other words,
when does the grammar only generate ﬁnite sentences? Giving a general answer seems
diﬃcult, but it is known that if the parameter values in a PCFG are obtained by learningfrom ﬁnite sentences, the stochastic derivation by the PCFG terminates with probability
one (Chi & Geman, 1998). In summary, assuming appropriate parameter values, we cansay that the parameterized logic program DB g0 for the largest PCFG with N nonterminal
symbols satisﬁes all applicability conditions, and the OLDT time for a sentence of length
L is O(N 3 L3 )53 and this is also the size of the support graph. From Proposition 5.1, weconcludeProposition 5.3 Let DB be a parameterized logic program representing a PCFG with N
nonterminals in the form of DB 0g in Figure 11, and G = G1 , G2 , . . . , GT be the sampled atoms
representing sentences of length L. We suppose each table operation in OLDT search is donein O(1) time. Then OLDT search for G and one iteration in learn-gEM are respectively
done in O(N 3 L3 T ) time.
(1)

PL−3

51. between(D0,D1,D2) is called O(N (L − d)2 ) times in Td . So it is called d=0 O(N (L − d)2 ) = O(N L3 )(1)
times in T0 .
52. Because the trial-ids used in the refutation record which rule is used at what step in the derivation of
w1 · · · wL .
53. In DB g0 , we represent integers by ground terms made out of 0 and s(·) to keep the program short. If
we use integers instead of ground terms however, the ﬁrst three arguments of q’(·,·,·,·,·,·) are enoughto check whether the goal is previously called or not, and this check can be done in O(1) time.

426

Parameter Learning of Logic Programs for Symbolic-statistical ModelingO(N 3 L3 T ) is also the time complexity of the Inside-Outside algorithm per iteration(Lari & Young, 1990), hence our algorithm is as eﬃcient as the Inside-Outside algorithm.
5.4 Pseudo PCSGs
PCFGs can be improved by making choices context-sensitive, and one of such attempts is
pseudo PCSGs (pseudo probabilistic context sensitive grammars) in which a rule is chosen probabilistically depending on both the nonterminal to be expanded and its parent
nonterminal (Charniak & Carroll, 1994).
A pseudo PCSG is easily programmed. We add one extra-argument, N, representing
the parent node, to the predicate q’(I,D0,D2,C0,C2,L0-L2) in Figure 11 and replace
msw(I,C0,[J,K]) with msw([N,I],C0,[J,K]). Since the (leftmost) derivation of a sentencefrom a pseudo PCSG is still a sequential decision process described by the modiﬁed program,the graphical EM algorithm applied to the support graphs generated from the modiﬁed
program and observed sentences correctly performs the ML estimation of parameters in thepseudo PCSG.
A pseudo PCSG is thought to be a PCFG with rules of the form [n, i] → [i, j][i, k](1 ≤ n, i, j, k ≤ N ) where n is the parent nonterminal of i, so the arguments in the previous
subsection are carried over with minor changes. We therefore have (details omitted)
Proposition 5.4 Let DB be a parameterized logic program for a pseudo PCSG with N
nonterminals as shown above, and G = G1 , G2 , . . . , GT the observed atoms representing
sampled sentences of length L. Suppose each table operation in OLDT search can be donein O(1) time. Then OLDT search for G and each iteration in learn-gEM is completed inO(N 4 L3 T ) time.
5.5 Bayesian Networks
A relationship between cause C and its eﬀect E is often probabilistic such as the one between diseases and symptoms, and as such it is mathematically captured as the conditional
probability P (E = e | C = c) of eﬀect e given the cause c. What we wish to know however isthe inverse, i.e. the probability of a candidate cause c given evidence e, i.e. P (C = c | E = e)Pwhich is calculated by Bayes’ theorem as P (E = e | C = c)P (C = c)/ c0 P (E = e | C =c0 )P (C = c0 ). Bayesian networks are a representational/computational framework that ﬁts
best this type of probabilistic inference (Pearl, 1988; Castillo et al., 1997).
A Bayesian network is a graphical representation of a joint distribution P (X1 = x1 , . . . ,XN = xN ) of ﬁnitely many random variables X1 , . . . , XN . The graph is a dag (directed
acyclic graph) such as ones in Figure 12, and each node is a random variable.54
In the graph, a conditional probability table (CPT) representing P (Xi = xi | Πi = ui )(1 ≤ i ≤ N ) is associated with each node Xi where Πi represents Xi ’s parent nodes and ui
their values. When Xi has no parent, i.e. a topmost node in the graph, the table is just amarginal distribution P (Xi = xi ). The whole joint distribution is deﬁned as the product of

54. We only deal with discrete cases.

427

Sato & Kameya

A

B

A

D

C

B

C

E

F

D
E

( G 1 ) Singly-connectedF( G 2 ) Multiply-connected

Figure 12: Bayesian networks
these conditional distributions:
P (X1 = x1 , . . . , XN = xN )55 =NY

P (Xi = xi | Πi = ui ).

(21)

i=1

Thus the graph G1 in Figure 12 deﬁnes
PG1 (a, b, c, d, e, f ) = PG1 (a)PG1 (b)PG1 (c | a)PG1 (d | a, b)PG1 (e | d)PG1 (f | d)where a, b, c, d, e and f are values of corresponding random variables A, B, C, D, E
and F , respectively.56 As mentioned before, one of the basic tasks of Bayesian networks
is to compute marginal probabilities. For example, the marginal distribution PG1 (c, d) is
computed either by (22) or (23) below.
PG1 (c, d) =X

PG1 (a)PG1 (b)PG1 (c | a)PG1 (d | a, b)PG1 (e | d)P (f | d)a,b,e,f



= 

X



PG1 (a)PG1 (b)PG1 (c | a)PG1 (d | a, b) 

a,b

X

(22)


PG1 (e | d)PG1 (f | d) (23)e,f(23) is clearly more eﬃcient than (22). Observe that if the graph were like G2 in
Figure 12, there would be no way to factorize computations like (23) but to use (22) requiringexponentially many operations. The problem is that computing marginal probabilities is
NP-hard in general, and factorization such as (23) is assured only when the graph is singlyconnected like G1 , i.e. has no loop when viewed as an undirected graph. In such case, thecomputation is possible in O(|V |) time where V is the set of vertices in the graph (Pearl,1988). Otherwise, the graph is called multiply-connected, and might need exponential time
to compute marginal probabilities. In the sequel, we show the following.
• For any discrete Bayesian network G deﬁning a distribution PG (x1 , . . . , xN ), there is aparameterized logic program DB G for a predicate bn(·) such that PDB G (bn(x1 ,. . .,xN ))= PG (x1 , . . . , xN ).
55. Thanks to the acyclicity of the graph, without losing generality, we may assume that if Xi is an ancestor
node of Xj , then i < j holds.
56. For notational simplicity, we shall omit random variables when no confusion arises.

428

Parameter Learning of Logic Programs for Symbolic-statistical Modeling• For arbitrary factorizations and their order to compute a marginal distribution, there
exists a tabled program that accomplishes the same computation in the speciﬁed way.
• When the graph is singly connected and evidence e is given, there exists a tabled
program DB Gτ such that OLDT time for ← bn(e) is O(|V |), and hence the timecomplexity per iteration of the graphical EM algorithm is O(|V |) as well.
Let G be a Bayesian network deﬁning a joint distribution PG (x1 , . . . , xN ) and {PG (Xi =xi | Πi = ui ) | 1 ≤ i ≤ N, xi ∈ val(Xi ), ui ∈ val(Πi )} the conditional probabilities
associated with G where val(Xi ) is the set of Xi ’s possible values and val(Πi ) denotes
the set of possible values of the parent nodes Πi as a random vector, respectively. We
construct a parameterized logic program that deﬁnes the same distribution PG (x1 , . . . , xN ).
Our program DB G = FG ∪ RG is shown in Figure 13.
FG= { msw(par(i,ui ),once,xi ) | 1 ≤ i ≤ N, ui ∈ val(Πi ), xi ∈ val(Xi ) }RG = { bn(X1 ,. . .,XN ):-VN

i=1 msw(par(i,Πi ),once,Xi ). }Figure 13: Bayesian network program DB G
FG is comprised of msw atoms of the form msw(par(i,ui ),once,xi ) whose probability isexactly the conditional probability PG (Xi = xi | Πi = ui ). When Xi has no parents, ui isthe empty list []. RG is a singleton, containing only one clause whose body is a conjunctionof msw atoms which corresponds to the product of conditional probabilities. Note that we
intentionally identify random variables X1 , . . . , XN with logical variables X1 , . . . , XN for
convenience.
Proposition 5.5 DB G denotes the same distributions as G.
(Proof) Let hx1 , . . . , xN i be a realization of the random vector hX1 , . . . , XN i. It holds byconstruction thatPDB G (bn(x1 ,. . .,xN )) =
=

N
Y
h=1
N
Y

Pmsw (msw(par(i,ui ),once,xi ))PG (Xi = xi | Πi = ui )h=1

= PG (x1 , . . . , xN ).

Q.E.D.

In the case of G1 in Figure 12, the program becomes57bn(A,B,C,D,E,F):-

msw(par(’A’,[]),once,A),msw(par(’C’,[A]),once,C),
msw(par(’E’,[D]),once,E),

57. 0 A0 , 0 B0 , . . . are Prolog constants used in place of integers.

429

msw(par(’B’,[]),once,B),
msw(par(’D’,[A,B]),once,D),msw(par(’F’,[D]),once,F).

Sato & Kameya

and the left-to-right sampling execution gives a sample realization of the random vectorh A, B, C, D, E, F i. A marginal distribution is computed from bn(x1 ,. . .,xN ) by adding a newclause to DB G . For example, to compute PG1 (c, d), we add bn(C,D):- bn(A,B,C,D,E,F)
to DB G1 (let the result be DB 0G1 ) and then compute PDB 0G (bn(c,d)) which is equal to1
PG1 (c, d) because
PDB 0G (bn(c,d)) = PDB G1 (∃ a, b, e, f bn(a,b,c,d,e,f ))1=

X

PDB G1 (bn(a,b,c,d,e,f ))a,b,e,f= PG1 (c, d).
Regrettably this computation corresponds to (22), not to the factorization (23). Eﬃcient
probability computation using factorization is made possible by carrying out summations
in a proper order.
We next sketch by an example how to carry out speciﬁed summations in a speciﬁedorder by introducing new clauses. Suppose we have a joint distribution P (x, y, z, w) =φ1 (x, y)φ2 (y, z, w)φ3 (x, z, w) such that φ1 (x, y), φ2 (y, z, w) and φ3 (x, z, w) are respectively
computed by atoms p1 (X,Y), p2 (Y,Z,W) and p3 (X,Z,W). Suppose also that we hope tocompute the sumÃ!
P (x) =X

φ1 (x, y)yX

φ2 (y, z, w)φ3 (x, z, w)

z,w

in which we ﬁrst eliminate z, w and then y. Corresponding to each elimination, we introduceP
two new predicates, q(X,Y) to compute φ4 (x, y) = z,w φ2 (y, z, w)φ3 (x, z, w) and p(X) toP
compute P (x) = y φ1 (x, y)φ4 (x, y) as follows.
p(X)
q(X,Y)::-

p1 (X,Y), q(X,Y).
p2 (Y,Z,W), p3 (X,Z,W).

Note that the clause body of q/2 contains Z and W as (existentially quantiﬁed) local variables
and the clause head q(X,Y) contains variables shared with other atoms. In view of thePcorrespondence between
and ∃, it is easy to conﬁrm that this program realizes therequired computation. It is also easy to see by generalizing this example, though we do
not prove here, that there exists a parameterized logic program that carries out the givensummations in the given order for an arbitrary Bayesian network, in particular we areable to simulate VE (variable elimination, Zhang & Poole, 1996; D’Ambrosio, 1999) in our
approach.
Eﬃcient computation of marginal distributions is not always possible but there is awell-known class of Bayesian networks, singly connected Bayesian networks, for which there
exists an eﬃcient algorithm to compute marginal distributions by message passing (Pearl,
1988; Castillo et al., 1997). We here show that when the graph is singly connected, we canconstruct an eﬃcient tabled Bayesian network program DB Gτ assigning a table predicate
to each node. To avoid complications, we explain the construction procedure informallyand concentrate on the case where we have only one interested variable. Let G be a singly
430

Parameter Learning of Logic Programs for Symbolic-statistical Modelingconnected graph. First we pick up a node U whose probability PG (u) is what we seek. We
construct a tree Gτ with the root node U from G, by letting other nodes dangling from U .
Figure 14 shows how G1 is transformed to a tree when we select node B as the root node.
B

D

A

E

F

C
τ

Transformed graph G1

Figure 14: Transforming G1 to a tree
Then we examine each node in Gτ one by one. We add for each node X in the graph acorresponding clause to DB Gτ whose purpose is to visit all nodes connected to X except theone that calls X. Suppose we started from the root node U1 in Figure 15 where evidence
u is given, and have generated clause (24). Now we proceed to an inner node X (U1 calls
X). In the original graph G, X has parent nodes {U1 , U2 , U3 } and child nodes {V1 , V2 }. U3
is a topmost node in G.
U1
X
U2

V2
U3

V1

Tree G τ

Figure 15: General situation
For node X in Figure 15, we add clause (25). When it is called from the parent nodeU1 with U1 being ground, we ﬁrst generate possible values of U2 by calling val U2 (U2 ),and then call call X U2 (U2 ) to visit all nodes connected to X through U2 . U3 is similary
treated. After visiting all nodes in G connecting to X through the parent nodes U2 andU3 (nodes connected to U1 have already been visited), the value of random variable X isdetermined by sampling the msw atom jointly indexed by ’X’ and the values of U1 , U2 and431

Sato & Kameya

U3 . Then we visit X’s children, V1 and V2 . For a topmost node U3 in the original graph,
we add clause (26).

tbn(U1 ) :- msw(par(’U1 ’,[]),once,U1 ), call U1 X(U1 ).

(24)call U1 X(U1 ) :- val U2 (U2 ), call X U2 (U2 ),val U3 (U3 ), call X U3 (U3 ),msw(par(’X’,[U1 ,U2 ,U3 ]),once,X),
call X V1 (X), call X V2 (X).
call X U3 (U3 ) :- msw(par(’U3 ’,[]),once,U3 ).

(25)
(26)

Let DB Gτ be the ﬁnal program containing clauses like (24), (25) and (26). Apparently
DB Gτ can be constructed in time linear in the number of nodes in the network. Also
note that successive unfolding (Tamaki & Sato, 1984) of atoms of the form call ...(·)
in the clause bodies that starts from (24) yields a program DB 0G similar to the one in
Figure 13 which contains msw atoms but no call ...(·)’s. As DB Gτ and DB 0G deﬁne thesame distribution,58 it can be proved from Proposition 5.5 that PG (u) = PDB 0G (bn(u)) =PDB Gτ (tbn(u)) holds (details omitted). By the way, in Figure 15 we assume the construction
starts from the topmost node U1 where the evidence u is given, but this is not necessary.
Suppose we change to start from the inner node X. In that case, we replace clause (24)with call X U1 (U1 ) :- msw(par(’U1 ’,[]),once,U1 ) just like (26). At the same time wereplace the head of clause (25) with tbn() and add a goal call X U1 (u) to the bodyand so on. For the changed program DB 00 Gτ , it is rather straightforward to prove that
PDB 00 Gτ (tbn()) = PG (u) holds. It is true that the construction of the tabled program
DB Gτ shown here is very crude and there is a lot of room for optimization, but it suﬃces
to show that a parameterized logic program for a singly connected Bayesian network runsin O(|V |) time where V is the set of nodes.
To estimate time complexity of OLDT search w.r.t. DB Gτ , we declare tbn and every
predicate of the form call ...(·) as table predicate and verify the ﬁve conditions on theapplicability of the graphical EM algorithm (details omitted). We now estimate the timecomplexity of OLDT search for the goal ← tbn(u) w.r.t. DB Gτ .59 We notice that calls occur
according to the pre-order scan (parents – the node – children) of the tree Gτ , and calls
to call Y X(·) occur val(Y ) times. Each call to call Y X(·) invokes calls to the rest ofnodes, X’s parents and X’s children in the graph Gτ except the caller node, with diﬀrent
set of variable instantiations, but from the second call on, every call only refers to solutions
stored in the solution table in O(1) time. Thus, the number of added computation steps in
58. Since distribution semantics is based on the least model semantics, and because unfold/fold transformation (Tamaki & Sato, 1984) preserves the least Herbrand model of the transformed program, unfold/foldtransformation applied to parameterized logic programs preserves the denotation of the transformed
program.
59. DB Gτ is further transformed for the OLDT interpreter to collect msw atoms like the case of the HMM
program.

432

Parameter Learning of Logic Programs for Symbolic-statistical Modeling

OLDT search by X is bounded from above, by constant O(val(U 1)val(U 2)val(U 3)val(X))
in the case of Figure 15. As a result OLDT time is proportional to the number of nodes
in the original graph G (and so is the size of the support graph) provided that the numberof edges connecting to a node, and that of values of a random variable are bounded fromabove. So we have
Proposition 5.6 Let G be a singly connected Bayesian network deﬁning distribution PG ,V the set of nodes, and DB Gτ the tabled program derived as above. Suppose the number ofedges connecting to a node, and that of values of a random variable are bounded from above
by some constant. Also suppose table access can be done in O(1) time. Then, OLDT time
for computing PG (u) for an observed value u of a random variable U by means of DB Gτ isO(|V |) and so is time per iteration required by the graphical EM algorithm. If there are T
observations, time complexity is O(|V |T ).
O(|V |) is the time complexity required to compute a marignal distribution for a singlyconnected Bayesian network by a standard algorithm (Pearl, 1988; Castillo et al., 1997),
and also is that of the EM algorithm using it. We therefore conclude that the graphicalEM algorithm is as eﬃcient as a specialzed EM algorithm for singly connected Bayesian
networks.60 We must also quickly add that the graphical EM algorithm is applicable toarbitrary Bayesian networks,61 and what Proposition 5.6 says is that an explosion of the
support graph can be avoided by appropriate programming in the case of singly connected
Bayesian networks.
To summarize, the graphical EM algorithm, a single generic EM algorithm, is proved tohave the same time complexity as specialized EM algorithms, i.e. the Baum-Welch algorithmfor HMMs, the Inside-Outside algorithm for PCFGs, and the one for singly connected
Bayesian networks that have been developed independently in each research ﬁeld.
Table 1 summarizes the time complexity of EM learning using OLDT search and thegraphical EM algorithm in the case of one observation. In the ﬁrst column, “sc-BNs”
represents singly connected Bayesian networks. The second column shows a program to use.
DB h is an HMM proram in Subsection 4.7, DB g0 a PCFG program in Subsection 5.3 andDB Gτ a transformed Bayesian network program in Subsection 5.5, respectively. OLDT timein the third column is time for OLDT search to complete the search of all t-explanations.
gEM in the fourth column is time in one iteration taken by the graphical EM algorithm
to update parameters. We use N , M , L and V respectively for the number of states in
an HMM, the number of nonterminals in a PCFG, the length of an input string and thenumber of nodes in a Bayesian network. The last column is a standard (specialized) EMalgorithm for each model.
60. When a marginal distribution of PG for more than one variable is required, we can construct a similartabled program that computes marginal probabilities still in O(|V |) time by adding extra-argumentsthat convey other evidence or by embedding other evidnece in the program.
61. We check the ﬁve conditions with DB G in Figure 13. The uniqueness condition is obvious as samplingalways uniquely generates a sampled value for each random variable. The ﬁnite support condition issatisﬁed because there are only a ﬁnite number of random variables and their values. The acyclic supportcondition is immediate because of the acyclicity of Bayesian networks. The t-exclusiveness condition andthe independent condition are easy to verify.

433

Sato & Kameya

Model
HMMs
PCFGs
sc-BNsuser model

Program

OLDT timegEM

Specialized EM

DB h
DB g0
DB Gτ

O(N 2 L)O(N 2 L)O(M 3 L3 )O(M 3 L3 )O(|V |)
O(|OLDT tree|)O(|V |)
O(|support graph|)

Baum-Welch
Inside-Outside(Castillo et al., 1997)

Table 1: Time complexity of EM learning by OLDT search and the graphical EM algorithm

5.6 Modeling Language PRISM
We have been developing a symbolic-statistical modeling laguage PRISM since 1995 (URL
= http://mi.cs.titech.ac.jp/prism/) as an implementation of distribution semantics(Sato, 1995; Sato & Kameya, 1997; Sato, 1998). The language is intented for modelingcomplex symbolic-statistical phenomena such as discourse interpretation in natural language
processing and gene inheritance interacting with social rules. As a programming language,
it looks like an extension of Prolog with new built-in predicates including the msw predicate
and other special predicates for manipulating msw atoms and their parameters.
A PRISM program is comprised of three parts, one for directives, one for modeling andone for utilities. The directive part contains declarations such as values, telling the system
what msw atoms will be used in the execution. The modeling part is a set of non-unit deﬁnite
clauses that deﬁne the distribution (denotation) of the program by using msw atoms. The
last part, the utility part, is an arbitary Prolog program which refers to predicates deﬁned
in the modeling part. We can use in the utility part learn built-in predicate to carry out
EM learning from observed atoms.
PRISM provides three modes of execution. The sampling execution correponds to arandom sampling drawn from the distribution deﬁned by the modeling part. The second
one computes the probability of a given atom. The third one returns the support set for agiven goal. These execution modes are available through built-in predicates.
We must report however that while the implementation of the graphical EM algorithmwith a simpiﬁed OLDT search mechanism has been under way, it is not completed yet. So
currently, only Prolog search and learn-naive(DB, G) in Section 4 are available for EM learning though we realized, partially, structrure sharing of explanations in the implementionof learn-naive(DB, G). Putting computational eﬃciecy aside however, there is no problem
in expressing and learning HMMs, PCFGs, pseudo PCSGs, Bayesian networks and otherprobailistic models by the current version. The learning experiments in the next sectionused a parser as a substitute for the OLDT interpreter, and the independently implemented
graphical EM algorithm.

6. Learning Experiments
After complexity analysis of the graphical EM algorithm for popular symbolic-probabilistic
models in the previous section, we look at an actual behavior of the graphical EM algorithmwith real data in this section. We conducted learning experiments with PCFGs using two
434Parameter Learning of Logic Programs for Symbolic-statistical Modelingcorpora which have contrasting characters, and compared the performance of the graphicalEM algorithm against that of the Inside-Outside algorithm in terms of time per iteration(= time for updating parameters). The results indicate that the graphical EM algorithm
can outperform the Inside-Outside algorithm by orders of magnigude. Detalis are reported
by Sato, Kameya, Abe, and Shirai (2001). Before proceeding, we review the Inside-Outsidealgorithm for completeness.
6.1 The Inside-Outside Algorithm
The Inside-Outside algorithm was proposed by Baker (1979) as a generalization of the
Baum-Welch algorithm to PCFGs. The algorithm is designed to estimate parameters for a
CFG grammar in Chomsky normal form containing rules expressed by numbers like i → j, k
(1 ≤ i, j, k ≤ N for N nonterminals, where 1 is a starting symbol). Suppose an input
sentence w1 , . . . , wL is given. In each iteration, it ﬁrst computes in a bottom up manner∗inside probabilities e(s, t, i) = P (i ⇒ ws , . . . , wt ) and then computes outside probabilities∗
f (s, t, i) = P (S ⇒ w1 , . . . , ws−1 i wt+1 , . . . , wL ) in a top-down manner for every s, t andi (1 ≤ s ≤ t ≤ L, 1 ≤ i ≤ N ). After computing both probabilities, parameters are
updated by using them, and this process iterates until some predetermined criterion suchas a convergence of the likelihood of the input sentence is achieved. Although Baker did
not give any analysis of the Inside-Outside algorithm, Lari and Young (1990) showed that it
takes O(N 3 L3 ) time in one iteration and Laﬀerty (1993) proved that it is the EM algorithm.
While it is true that the Inside-Outside algorithm has been recognized as a standard EM
algortihm for training PCFGs, it is notoriously slow. Although there is not much literature
explicitly stating time required by the Inside-Outside algorithm (Carroll & Rooth, 1998;
Beil, Carroll, Prescher, Riezler, & Rooth, 1999), Beil et al. (1999) reported for example
that when they trained a PCFG with 5,508 rules for a corpus of 450,526 German subordinate clauses whose average ambiguity is 9,202 trees/clause using four machines (167MHzSun UltraSPARC×2 and 296MHz Sun UltraSPARC-II×2), it took 2.5 hours to complete
one iteration. We discuss later why the Inside-Outside algorithm is slow.
6.2 Learning Experiments Using Two Corpora
We report here parameter learning of existing PCFGs using two corpora of moderate sizeand compare the graphical EM algorithm against the Inside-Outside algorithm in termsof time per iteration. As mentioned before, support graphs, input to the garphical EM
algorithm, were generated by a parser, i.e. MSLR parser.62 All measurements were madeon a 296MHz Sun UltraSPARC-II with 2GB memory under Solaris 2.6 and the thresholdfor an increase of the log likelihood of input sentences was set to 10−6 as a stopping criterion
for the EM algorithms.
In the experiments, we used ATR corpus and EDR corpus (each converted to a POS(part of speech)-tagged corpus). They are similar in size (about 10,000) but contrasting intheir characters, sentence length and ambiguity of their grammars. The ﬁrst experimentemployed ATR corpus which is a Japanese-English corpus (we used only the Japanese part)developed by ATR (Uratani, Takezawa, Matsuo, & Morita, 1994). It contains 10,995 short
62. MSLR parser is a Tomita (Generalized LR) parser developed by Tanaka-Tokunaga Laboratory in Tokyo
Institute of Technology (Tanaka, Takezawa, & Etoh, 1997).

435

Sato & Kameya

conversational sentences, whose minimum length, average length and maximum length arerespectively 2, 9.97 and 49. As a skeleton of PCFG, we employed a context free grammar
Gatr comprising 860 rules (172 nonterminals and 441 terminals) manually developed for
ATR corpus (Tanaka et al., 1997) which yields 958 parses/sentence.
Because the Inside-Outside algorithm only accepts a CFG in Chomsky normal form, weconverted Gatr into Chomsky normal form G∗atr . G∗atr contains 2,105 rules (196 nonterminals and 441 terminals). We then divided the corpus into subgroups of similar length
like (L = 1, 2), (L = 3, 4), . . . , (L = 25, 26), each containing randomly chosen 100 sentences.
After these preparations, we compare at each length the graphical EM algorithm applied to
Gatr and G∗atr against the Inside-Outside algorithm applied to G∗atr in terms of time periteration by running them until convergence.
(sec)

(sec)

(sec)

60
I-O
50

0.7

0.04

0.6

0.035

0.5

40

0.4

0.03

I-O
gEM (original)gEM (Chomsky NF)

0.025
0.02

30
0.3
20

0.015

0.2

10

0.01

0.1

0.005

5

10

15

20

25

L

L

L
0

gEM (original)gEM (Chomsky NF)

0

5

10

15

20

25

0

5

10 15 20 25

Figure 16: Time per iteration : I-O vs. gEM (ATR)
Curves in Figure 16 show the learning results where an x-axis is the length L of an input
sentence and a y-axis is average time taken by the EM algorithm in one iteration to update
all parameters contained in the support graphs generated from the chosen 100 sentences(other parameters in the grammar do not change). In the left graph, the Inside-Outsidealgorithm plots a cubic curve labeled “I-O”. We omitted a curve drawn by the graphicalEM algorithm as it drew the x-axis. The middle graph magniﬁes the left graph. The curvelabeled “gEM (original)” is plotted by the graphical EM algorithm applied to the original
grammar Gatr whereas the one labeled “gEM (Chomsky NF)” used G∗atr . At length 10, theaverage sentence length, it is measured that whichever grammar is employed, the graphicalEM algorithm runs several hundreds times faster (845 times faster in the case of Gatr
and 720 times faster in the case of G∗atr ) than the Inside-Outside algorithm per iteration.
The right graph shows (almost) linear dependency of updating time by the graphical EM
algorithm within the measuared sentence length.
Although some diﬀerence is anticipated in their learning speed, the speed gap between
the Inside-Outside algorithm and the graphical EM algorithm is unexpectedly large. The
most conceivable reason is that ATR corpus only contains short sentences and Gatr is not436

Parameter Learning of Logic Programs for Symbolic-statistical Modelingmuch ambiguous so that parse trees are sparse and generated support graphs are small,which aﬀects favorably the perforamnce of the graphical EM algorithm.
We therefore conducted the same experiment with another corpus which contains much
longer sentences using a more ambiguous grammar that generates dense parse trees. We
used EDR Japanese corpus (Japan EDR, 1995) containing 220,000 Japanese news article
sentences. It is however under the process of re-annotation, and only part of it (randomly
sampled 9,900 sentences) has recently been made available as a labeled corpus. Comparedwith ATR corpus, sentences are much longer (the average length of 9,900 sentences is 20,
the minimum length 5, the maximum length 63) and a CFG grammar Gedr (2,687 rules,
converted to Chomsky normal form grammar G∗edr containing 12,798 rules) developed for
it is very ambiguous (to keep a coverage rate), having 3.0 × 108 parses/sentence at length
20 and 6.7 × 1019 parses/sentence at length 38.
(sec)

5000

(sec)

(sec)
3

10

6000
I-O

8

I-O
gEM (original)

2.5

gEM (original)

2

4000

6
1.5

3000

4
1

2000

2

1000
0

L
5 10 15 20 25 30 35 40

0

0.5

L
5 10 15 20 25 30 35 40

0L
5 10 15 20 25 30 35 40

Figure 17: Time per iteration : I-O vs. gEM (EDR)
Figure 17 shows the obtained curves from the experiments with EDR corpus (the graphical EM algorithm applied to Gedr vs. the Inside-Outside algorithm applied to G∗edr ) under
the same condition as ATR corpus, i.e. plotting average time per iteration to process 100
sentences of the designated length, except that the plotted time for the Inside-Outside algorithm is the average of 20 iterations whereas that for the graphical EM algorithm is theaverage of 100 iterations. As is clear from the middle graph, this time again, the graphicalEM algorithm runs orders of magnitude faster than the Inside-Outside algorithm. At average sentence length 20, the former takes 0.255 second whereas the latter takes 339 seconds,
giving a speed ratio of 1,300 to 1. At sentence length 38, the former takes 2.541 seconds but
the latter takes 4,774 seconds, giving a speed ratio of 1,878 to 1. Thus the speed ratio evenwidens compared to ATR corpus. This can be explained by the mixed eﬀects of O(L3 ),time complexity of the Inside-Outside algorithm, and a moderate increase in the total sizeof support graphs w.r.t. L. Notice that the right graph shows how the total size of support
graphs grows with sentence length L as time per iteration by the graphical EM algorithm
is linear in the total size of support graphs.

437

Sato & Kameya

Since we implemented the Inside-Outside algorithm faithfully to Baker (1979), Lari andYoung (1990), there is much room for improvement. Actually Kita gave a reﬁned InsideOutside algorithm (Kita, 1999). There is also an implementation by Mark Johnson of the
Inside-Outside algorithm down-loadable from http://www.cog.brown.edu/%7Emj/. The
use of such implementations may lead to diﬀerent conclusions. We therefore conducted
learning experiments with the entire ATR corpus using these two implementations andmeasured updating time per iteration (Sato et al., 2001). It turned out that both implementations run twice as fast as our naive implementation and take about 630 seconds periteration while the graphical EM algorithm takes 0.661 second per iteration, which is stillorders of magnitude faster than the former two. Regrettably a similar comparison using
the entire EDR corpus available at the moment was abandoned due to memory overﬂow
during parsing for the construction of support graphs.
Learning experiments so far only compared time per iteration which ignore extra time
for search (parsing) required by the graphical EM algorithm. So a question naturally arises
w.r.t. comparison in terms of total learning time. Assuming 100 iterations for learning
ATR corpus however, it is estimated that even considering parsing time, the graphicalEM algorithm combined with MSLR parser runs orders of magnitude faster than the three
implementations (ours, Kita’s and Johnson’s) of the Inside-Outside algorithm (Sato et al.,2001). Of course this estimation does not directly apply to the graphical EM algorithm
combined with OLDT search, as the OLDT interpreter will take more time than a parser
and how much more time is needed depends on the implementaiton of OLDT search.63
Conversely, however, we may be able to take it as a rough indication of how far our approach,
the graphical EM algorithm combined with OLDT search via support graphs, can go in thedomain of EM learning of PCFGs.
6.3 Examing the Performance Gap
In the previous subsection, we compared the performance of the graphical EM algorithm
against the Inside-Outside algorithm when PCFGs are given, using two corpora and three
implementations of the Inside-Outside algorithm. In all experiments, the graphical EM
algorithm considerably outperformed the Inside-Outside algorithm despite the fact thatboth have the same time complexity. Now we look into what causes such a performance
gap.
Simply put, the Inside-Outside algorithm is slow (primarily) because it lacks parsing.
Even when a backbone CFG grammar is explicitly given, it does not take any advantage ofthe constraints imposed by the grammar. To see it, it might help to review how the insideprobability e(s, t, A), i.e. P(nonterminal A spans from s-th word to t-th word) (s ≤ t), is
calculated by the Inside-Outside algorithm for the given grammar.
e(s, t, A) =

X

r=t−1
X

P(A → BC)e(s, r, B)e(r + 1, t, C)

B,C s.t. A→BC in the grammar r=s

Here P(A → BC) is a probability associated with a production rule A → BC. Note that fora ﬁxed triplet (s, t, A), it is usual that the term P(A → BC)e(s, r, B)e(r+1, t, C) is non-zero
63. We cannnot answer this question right now as the implementation of OLDT search is not completed.

438

Parameter Learning of Logic Programs for Symbolic-statistical Modeling

only for a relatively small number of (B, C, r)’s determined from successful parses and therest of combinations always give 0 to the term. Nonetheless the Inside-Outside algorithm
attempts to compute the term in every iteration for all possible combinations of B, C andr and this is repeated for every possible (s, t, A), resulting in a lot of redundancy. The samekind of redundancy occurs in the computation of outside probability by the Inside-Outsidealgorithm.
The graphical EM algorithm is free of such redundancy because it runs on parse trees (aparse forest) represented by the support graph.64 It must be added, on the other hand, thatsuperiority in learning speed of the graphical EM algorithm is realized at the cost of spacecomplexity because while the Inside-Outside algorithm merely requires O(N L2 ) space for
its array to store probabilities, the graphical EM algorithm needs O(N 3 L3 ) space to store
the support graph where N is the number of nonterminals and L is the sentence length.
This trade-oﬀ is understandable if one notices that the graphical EM algorithm applied
to a PCFG can be considered as partial evaluation of the Inside-Outside algorithm by thegrammar (and the introduction of appropriate data structure for the output).
Finally we remark that the use of parsing as a preprocess for EM learning of PCFGs is
not unique to the graphical EM algorithm (Fujisaki, Jelinek, Cocke, Black, & Nishino, 1989;
Stolcke, 1995). These approaches however still seem to contain redundancies compared withthe graphical EM algorithm. For instance Stolcke (1995) uses an Earley chart to computeinside and outside probability, but parses are implicitly reconstructed in each iteration
dynamically by combining completed items.

7. Related Work and Discussion
7.1 Related Work
The work presented in this paper is at the crossroads of logic programming and probability
theory, and considering an enormous body of work done in these ﬁelds, incompleteness isunavoidable when reviewing related work. Having said that, we look at various attempts
made to integrate probability with computational logic or logic programming.65 In reviewing, one can immediately notice there are two types of usage of probability. One type,
constraint approach, emphasizes the role of probability as constraints and does not necessarily seek for a unique probability distribution over logical formulas. The other type,
distribution approach, explicitly deﬁnes a unique distribution by model theoretical meansor proof theoretical means, to compute various probabilities of propositions.
A typical constraint approach is seen in the early work of probabilistic logic by Nilsson(1986). His central problem, “probabilistic entailment problem”, is to compute the upperand lower bound of probability P(φ) of a target sentence φ in such a way that the bounds
are compatible with a given knowledge base containing logical sentences (not necessarily
logic programs) annotated with a probability. These probabilities work as constraints on
64. We emphasize that the diﬀerence between the Inside-Outside algorithm and the graphical EM algorithm
is solely computational eﬃciency, and they converge to the same parameter values when starting from thesame initial values. Linguistic evaluations of the estimated parameters by the graphical EM algorithm
are also reported by Sato et al. (2001).
65. We omit literature leaning strongly toward logic. For logic(s) concerning uncertainty, see an overview
by Kyburg (1994).

439

Sato & Kameya

the possible range of P(φ). He used the linear programming technique to solve this problem
that inevitably delimits the applicability of his approach to ﬁnite domains.
Later Lukasiewicz (1999) investigated the computational complexity of the probabilistic
entailment problem in a slightly diﬀerent setting. His knowledge base comprises statementsof the form (H | G)[u1 , u2 ] representing u1 ≤ P(H | G) ≤ u2 . He showed that inferring“tight” u1 , u2 is NP-hard in general, and proposed a tractable class of knowledge base called
conditional constraint trees.
After the inﬂuential work of Nilsson, Frish and Haddawy (1994) introduced a deductive system for probabilistic logic that remedies “drawbacks” of Nilsson’s approach, thatof computational intractability and the lack of a proof system. Their system deduces aprobability range of a proposition by rules of probabilistic inferences about unconditionaland conditional probabilities. For instance, one of the rules infers P (α | δ) ∈ [0 y] fromP (α ∨ β | δ) ∈ [x y] where α,β and δ are propositional variables and [x y] (x ≤ y) designatesa probability range.
Turning to logic programming, probabilistic logic programming formalized by Ng andSubrahmanian (1992) and Dekhtyar and Subrahmanian (1997) was also a constraint approach. Their program is a set of annotated clauses of the form A : µ ← F1 : µ1 , . . . , Fn : µnwhere A is an atom, Fi (1 ≤ i ≤ n) a basic formula, i.e. a conjunction or a disjunction of
atoms, and µj (0 ≤ j ≤ n) a sub-interval of [0, 1] indicating a probability range. A query← ∃ (F1 : µ1 , . . . , Fn : µn ) is answered by an extension of SLD refutation. On formalization,
it is assumed that their language contains only a ﬁnite number of constant and predicate
symbols, and no function symbol is allowed.
A similar framework was proposed by Lakshmanan and Sadri (1994) under the same syntactic restrictions (ﬁnitely many constant and predicate symbols but no function symbols)cin a diﬀerent uncertainty setting. They used annotated clauses of the form A ← B1 , . . . , Bn
where A and Bi (1 ≤ i ≤ n) are atoms and c = h[α, β], [γ, δ]i, a conﬁdence level, represents
a belief interval [α, β] (0 ≤ α ≤ β ≤ 1) and a doubt interval [γ, δ] (0 ≤ γ ≤ δ ≤ 1), which
an expert has in the clause.
As seen above, deﬁning a unique probability distribution is of secondary or no concern
to the constraint approach. This is in sharp contrast with Bayesian networks as the wholediscipline rests on the ability of the networks to deﬁne a unique probability distribution(Pearl, 1988; Castillo et al., 1997). Researchers in Bayesian networks have been seeking fora way of mixing Bayesian networks with a logical representation to increase their inherently
propositional expressive power.
Breese (1992) used logic programs to automatically build a Bayesian network from aquery. In Breese’s approach, a program is the union of a deﬁnite clause program and a setof conditional dependencies of the form P(P | Q1 ∧ · · · ∧ Qn ) where P and Qi s are atoms.
Given a query, a Bayesian network is constructed dynamically that connects the query andrelevant atoms in the program, which in turn deﬁnes a local distribution for the connected
atoms. Logical variables can appear in atoms but no function symbol is allowed.
Ngo and Haddawy (1997) extended Breese’s approach by incorporating a mechanism
reﬂecting context. They used a clause of the form P(A0 | A1 , . . . , An ) = α ← L1 , . . . , Lk ,where Ai ’s are called p-atoms (probabilistic atoms) whereas Lj ’s are context atoms disjoint
from p-atoms, and computed by another general logic program (satisfying certain restric440

Parameter Learning of Logic Programs for Symbolic-statistical Modelingtions). Given a query, a set of evidence and context atoms, relevant ground p-atoms are
identiﬁed by resolving context atoms away by SLDNF resolution, and a local Bayesian network is built to calculate the probability of the query. They proved the soundness andcompleteness of their query evaluation procedure under the condition that programs areacyclic66 and domains are ﬁnite.
Instead of deﬁning a local distribution for each query, Poole (1993) deﬁned a global distribution in his “probabilistic Horn abduction”. His program consists of deﬁnite clauses anddisjoint declarations of the form disjoint([h1 :p1 ,...,hn :pn ]) which speciﬁes a probability distribution over the hypotheses (abducibles) {h1 , . . . , hn }. He assigned probabilities toall ground atoms with the help of the theory of logic programming, and furthermore provedthat Bayesian networks are representable in his framework. Unlike previous approaches, his
language contains function symbols, but the acyclicity condition imposed on the programsfor his semantics to be deﬁnable seems to be a severe restriction. Also, probabilities are notdeﬁned for quantiﬁed formulas.
Bacchus et al. (1996) used a much more powerful ﬁrst-order probabilistic language thanclauses annotated with probabilities. Their language allows a statistically quantiﬁed termsuch as k φ(x)|θ(x) kx to denote the ratio of individuals in a ﬁnite domain satisfying φ(x) ∧θ(x) to those satisfying θ(x). Assuming that every world (interpretation for their language)is equally likely, they deﬁne theof a¶sentence ϕ under the given knowledgeµ probabilityτ#worldsN (ϕ∧KB)
where #worldsτN (χ) is the number ofbase KB as the limit limN →∞τ#worldsN (KB)
possible worlds containing N individuals satisfying χ, and τ parameters used in judging
approximations. Although the limit does not necessarily exist and the domain must be ﬁnite,they showed that their method can cope with diﬃculties arising from “direct inference” anddefault reasoning.
In a more linguistic vein, Muggleton (1996, and others) formulated SLPs (stochastic
logic programs) procedurally, as an extension of PCFGs to probabilistic logic programs.
So, a clause C, which must be range-restricted,67 is annotated with a probability p likep : C. The probability of a goal G is the product of such ps appearing in its refutation but
with a modiﬁcation such that if a subgoal g can invoke n clauses, pi : Ci (1 ≤ i ≤ n) atPsome refutation step, the probability of choosing k-th clause is normalized to pk / ni=1 pi .
More recently, Cussens (1999, 2001) enriched SLPs by introducing a special class of
log-linear models for SLD refutations w.r.t. a given goal. He for example considers all
possible SLD refutations for the most general goal ← s(X) and deﬁnes probability P(R)P
of a refutation R as P(R) = Z −1 exp ( i λi ν(R, i)). Here λi is a number associated witha clause Ci and ν(R, i) is a feature, i.e. the number of occurrences of Ci in R. Z is thenormalizing constant. Then, the probability assigned to s(a) is the sum of probabilities ofrefutation for ← s(a).

66. The condition says that every ground atom A must be assigned a unique integer n(A) such that n(A) >n(B1 ), . . . , n(Bn ) holds for any ground instance of a clause of the form A ← B1 , . . . , Bn . Under thiscondition, when a program includes p(X) ← q(X, Y ), we cannot write recursive clauses about q such asq(X, [H|Y ]) ← q(X, Y ).
67. A syntactic property that variables appearing in the head also appear in the body of a clause. A unit
clause must be ground.

441

Sato & Kameya

7.2 Limitations and Potential Problems
Approaches described so far have more or less similar limitations and potential problems.
Descriptive power conﬁned to ﬁnite domains is the most common limitation, which is dueto the use of the linear programming technique (Nilsson, 1986), or due to the syntactic
restrictions not allowing for inﬁnitely many constant, function or predicate symbols (Ng& Subrahmanian, 1992; Lakshmanan & Sadri, 1994). Bayesian networks have the same
limitation as well (only a ﬁnite number of random variables are representable).68 Also there
are various semantic/syntactic restrictions on logic programs. For instance the acyclicity
condition imposed by Poole (1993) and Ngo and Haddawy (1997) prevents the unconditionaluse of clauses with local variables, and the range-restrictedness imposed by Muggleton(1996) and Cussens (1999) excludes programs such as the usual membership Prolog program.
There is another type of problem, the possibility of assigning conﬂicting probabilities
to logically equivalent formulas. In SLPs, P(A) and P(A ∧ A) do not necessarily coincidebecause A and A∧A may have diﬀerent refutations (Muggleton, 1996; Cussens, 1999, 2001).
Consequently in SLPs, we would be in trouble if we naively interpret P(A) as the probabilityof A’s being true. Also assigning probabilities to arbitrary quantiﬁed formulas seems out of
scope of both approaches to SLPs.
Last but not least, there is a big problem common to any approach using probabilities:
where do the numbers come from? Generally speaking, if we use n binary random variables ina model, we have to determine 2n probabilities to completely specify their joint distribution,and fulﬁlling this requirement with reliable numbers quickly becomes impossible as n grows.
The situation is even worse when there are unobservable variables in the model such as
possible causes of a disease. Apparently parameter learning from observed data is a natural
solution to this problem, but parameter learning of logic programs has not been well studied.
Distribution semantics proposed by Sato (1995) was an attempt to solve these problems
along the line of the global distribution approach. It deﬁnes a distribution (probability
measure) over the possible interpretations of ground atoms for an arbitrary logic programin any ﬁrst order language and assigns consistent probabilities to all closed formulas. Also
distribution semantics enabled us to derive an EM algorithm for the parameter learning oflogic programs for the ﬁrst time. As it was a naive algorithm however, dealing with large
problems was diﬃcult when there are exponentially many explanations for an observation
like HMMs. We believe that the eﬃciency problem is solved to a large extent by thegraphical EM algorithm presented in this paper.
7.3 EM Learning
Since EM learning is one of the central issues in this paper, we separately mention work
related to EM learning for symbolic frameworks. Koller and Pfeﬀer (1997) used in theirapproach to KBMC (knowledge-based model construction) EM learning to estimate parameters labeling clauses. They express probabilistic dependencies among events by deﬁnite clauses annotated with probabilities, similarly to Ngo and Haddawy’s (1997) approach,and locally build a Bayesian network relevant to the context and evidence as well as the
68. However, RPMs (recursive probability models) proposed by Pfeﬀer and Koller (2000) as an extensionof Bayesian networks allow for inﬁnitely many random variables. They are organized as attributes ofclasses and a probability measure over attribute values is introduced.

442

Parameter Learning of Logic Programs for Symbolic-statistical Modelingquery. Parameters are learned by applying to the constructed network the specialized EM
algorithm for Bayesian networks (Castillo et al., 1997).
Dealing with a PCFG by a statically constructed Bayesian network was proposed Pynadath and Wellman (1998), and it is possible to combine the EM algorithm with their method
to estimate parameters in the PCFG. Unfortunately, the constructed network is not singlyconnected, and time complexity of probability computation is potentially exponential in thelength of an input sentence.
Closely related to our EM learning is parameter learning of log-linear models. Riezler (1998) proposed the IM algorithm in his approach to probabilistic constraint programming. The IM algorithm is a general parameter estimation algorithm from incomplete data for log-linear models whose probability function P(x) takes the form P(x) =P
Z −1 exp ( ni=1 λi νi (x)) p0 (x) where (λ1 , . . . , λn ) are parameters to be estimated, νi (x) thei-th feature of an observed object x and Z the normalizing constant. Since a feature canbe any function of x, the log-linear model is highly ﬂexible and includes our distribution
Pmsw as a special case of Z = 1. There is a price to pay however; the computational costof Z. It requires a summation over exponentially many terms. To avoid the cost of exact
computation, approximate computation by a Monte Carlo method is possible. Whichever
one may choose however, learning time increases compared to the EM algorithm for Z = 1.
The FAM (failure-adjusted maximization) algorithm proposed by Cussens (2001) is an
EM algorithm applicable to pure normalized SLPs that may fail. It deals with a specialclass of log-linear models but is more eﬃcient than the IM algorithm. Because the statisticalframework of the FAM is rather diﬀerent from distribution semantics, comparison with thegraphical EM algorithm seems diﬃcult.
Being slightly tangential to EM learning, Koller et al. (1997) developed a functionalmodeling language deﬁning a probability distribution over symbolic structures in which
they showed “cashing” of computed results leads to eﬃcient probability computation of
singly connected Bayesian networks and PCFGs. Their cashing corresponds to the computation of inside probability in the Inside-Outside algorithm and the computation of outside
probability is untouched.
7.4 Future Directions
Parameterized logic programs are expected to be a useful modeling tool for complex symbolicstatistical phenomena. We have tried various types of modeling, besides stochastic grammars and Bayesian networks, such as the modeling of gene inheritance in the Kariera tribe(White, 1963) where the rules of bi-lateral cross-cousin marriage for four clans interact withthe rules of genetic inheritance (Sato, 1998). The model was quite interdisciplinary, but
the ﬂexibility of combining msw atoms by means of deﬁnite clauses greatly facilitated themodeling process.
Although satisfying the ﬁve conditions in Section 4• the uniqueness condition (roughly, one cause yields one eﬀect)• the ﬁnite support condition (there are a ﬁnite number of explanations for one observation)• the acyclic support condition (explanations must not be cyclic)
443

Sato & Kameya• the t-exclusiveness condition (explanations must be mutually exclusive)• the independence condition (events in an explanation must be independent)
for the applicability of the graphical EM algorithm seems daunting, our modeling experiences so far tell us that the modeling principle in Section 4 eﬀectively guides us to successfulmodeling. In return, we can obtain a declarative model described compactly by a high level
language whose parameters are eﬃciently learnable by the graphical EM algorithm as shownin the preceding section.
One of the future directions is however to relax some of the applicability conditions,especially the uniqueness condition that prohibits a generative model from failure or fromgenerating multiple observable events. Although we pointed out in Section 4.4 that the MAR
condition in Appendix B adapted to our semantics can replace the uniqueness condition andvalidates the use of the graphical EM algorithm even when a complete data does not uniquelydetermine the observed data just like the case of “partially bracketed corpora” (Pereira &Schabes, 1992), we feel the need to do more research on this topic. Also investigating
the role of the acyclicity condition seems theoretically interesting as the acyclicity is often
related to the learning of logic programs (Arimura, 1997; Reddy & Tadepalli, 1998).
In this paper we only scratched the surface of individual research ﬁelds such as HMMs,
PCFGs and Bayesian networks. Therefore, there remains much to be done about clarifying
how experiences in each research ﬁeld are reﬂected in the framework of parameterized logic
programs. For example, we need to clarify the relationship between symbolic approaches
to Bayesian networks such as SPI (Li, Z. & D’Ambrosio, B., 1994) and our approach.
Also it is unclear how a compiled approach using the junction tree algorithm for Bayesian
networks can be incorporated into our approach. Aside from exact methods, approximate
methods of probability computation specialized for parameterized logic programs must also
be developed.
There is also a direction of improving learning ability by introducing priors instead of ML
estimation to cope with data sparseness. The introduction of basic distributions that makeprobabilistic switches correlated seems worth trying in the near future. It is also importantto take advantage of the logical nature of our approach to handle uncertainty. For example,
it is already shown by Sato (2001) that we can learn parameters from negative examples
such as “the grass is not wet” but the treatment of negative examples in parameterized
logic programs is still in its infancy.
Concerning developing complex statistical models based on the “programs as distributions” scheme, stochastic natural language processing which exploits semantic information
seems promising. For instance, uniﬁcation-based grammars such as HPSGs (Abney, 1997)
may be a good target beyond PCFGs because they use feature structures logically describable, and the ambiguity of feature values seems to be expressible by a probability
distribution.
Also building a mathematical basis for logic programs with continuous random variables
is a challenging research topic.

444

Parameter Learning of Logic Programs for Symbolic-statistical Modeling

8. Conclusion
We have proposed a logical/mathematical framework for statistical parameter learning ofparameterized logic programs, i.e. deﬁnite clause programs containing probabilistic facts
with a parameterized probability distribution. It extends the traditional least Herbrandmodel semantics in logic programming to distribution semantics, possible world semantics
with a probability distribution over possible worlds (Herbrand interpretations) which is
unconditionally applicable to arbitrary logic programs including ones for HMMs, PCFGs
and Bayesian networks.
We also have presented a new EM algorithm, the graphical EM algorithm in Section 4,
which learns statistical parameters from observations for a class of parameterized logic programs representing a sequential decision process in which each decision is exclusive andindependent. It works on support graphs, a new data structure specifying a logical relationship between an observed goal and its explanations, and estimates parameters by computing
inside and outside probability generalized for logic programs.
The complexity analysis in Section 5 showed that when OLDT search, a complete tabledrefutation method for logic programs, is employed for the support graph construction andtable access is done in O(1) time, the graphical EM algorithm, despite its generality, has
the same time complexity as existing EM algorithms, i.e. the Baum-Welch algorithm for
HMMs, the Inside-Outside algorithm for PCFGs and the one for singly connected Bayesian
networks that have been developed independently in each research ﬁeld. In addition, forpseudo probabilistic context sensitive grammars with N nonterminals, we showed that thegraphical EM algorithm runs in time O(N 4 L3 ) for a sentence of length L.
To compare actual performance of the graphical EM algorithm against the InsideOutside algorithm, we conducted learning experiments with PCFGs in Section 6 using two
real corpora with contrasting characters. One is ATR corpus containing short sentences forwhich the grammar is not much ambiguous (958 parses/sentence), and the other is EDR
corpus containing long sentences for which the grammar is rather ambiguous (3.0 × 108
at average sentence length 20). In both cases, the graphical EM algorithm outperformed
the Inside-Outside algorithm by orders of magnitude in terms of time per iteration, whichsuggests the eﬀectiveness of our approach to EM learning by the graphical EM algorithm.
Since our semantics is not limited to ﬁnite domains or ﬁnitely many random variables
but applicable to any logic programs of arbitrary complexity, the graphical EM algorithm is
expected to give a general yet eﬃcient method of parameter learning for models of complex
symbolic-statistical phenomena governed by rules and probabilities.

Acknowledgments
The authors wish to thank three anonymous referees for their comments and suggestions.
Special thanks go to Takashi Mori and Shigeru Abe for stimulating discussions and learning
experiments, and also to Tanaka-Tokunaga Laboratory for kindly allowing them to use
MSLR parser and the linguistic data.

445

Sato & Kameya

Appendix A. Properties of PDB
In this appendix, we list some properties of PDB deﬁned by a parameterized logic program
DB = F ∪ R in a countable ﬁrst-order language L.69 First of all, PDB assigns consistentprobabilities70 to every closed formula φ in L bydef

PDB (φ) = PDB ({ω ∈ ΩDB | ω |= φ})while guaranteeing continuity in the sense thatlimn→∞ PDB (φ(t1 ) ∧ · · · ∧ φ(tn )) = PDB (∀xφ(x))
limn→∞ PDB (φ(t1 ) ∨ · · · ∨ φ(tn )) = PDB (∃xφ(x))
where t1 , t2 , . . . is an enumeration of ground terms in L.
The next proposition, Proposition A.1, relates PDB to the Herbrand model. To proveit, we need some terminology. A factor is a closed formula in prenex disjunctive normal
form Q1 · · · Qn M where Qi (1 ≤ i ≤ n) is either an existential quantiﬁcation or a universal
quantiﬁcation and M a matrix. The length of quantiﬁcations n is called the rank of thefactor. Deﬁne Φ as a set of formulas made out of factors, conjunctions and disjunctions.
Associate with each formula φ in Φ a multi-set r(φ) of ranks by

r(φ) =∅
if φ is a factor with no quantiﬁcation{n}
if φ is a factor with rank n r(φ ) ] r(φ ) if φ = φ ∨ φ or φ = φ ∧ φ .
1
2
1
2
1
2

Here ] stands for the union of two multi-sets. For instance {1, 2, 3}]{2, 3, 4} = {1, 2, 2, 3, 3, 4}.
We use the multi-set ordering in the proof of Proposition A.1 because the usual induction
on the complexity of formulas does not work.
Lemma A.1 Let φ be a boolean formula made out of ground atoms in L. PDB (φ) =PF ({ν ∈ ΩF | MDB (ν) |= φ}).
(Proof) We have only to prove the lemma about a conjunction of atoms of the form D1x1 ∧· · · ∧ Dnxn (xi ∈ {0, 1}, 1 ≤ i ≤ n).
PDB (D1x1 ∧ · · · ∧ Dnxn ) = PDB ({ω ∈ ΩDB | ω |= D1x1 ∧ · · · ∧ Dnxn })= PDB (D1 = x1 , . . . , Dn = xn )= PF ({ν ∈ ΩF | MDB (ν) |= D1x1 ∧ · · · ∧ Dnxn })Q.E.D.

Proposition A.1 Let φ be a closed formula in L. PDB (φ) = PF ({ν ∈ ΩF | MDB (ν) |= φ}).

69. For deﬁnitions of ΩF , PF , MDB (ν), ΩDB , PDB and others used below, see Section 3.
70. By consistent, we mean probabilities assigned to logical formulas respect the laws of probability such as
0 ≤ P (A) ≤ 1, P (¬A) = 1 − P (A) and P (A ∨ B) = P (A) + P (B) − P (A ∧ B).

446

Parameter Learning of Logic Programs for Symbolic-statistical Modeling(Proof) Recall that a closed formula has an equivalent prenex disjunctive normal form
that belongs to Φ. We prove the proposition for formulas in Φ by using induction on themulti-set ordering over {r(φ) | φ ∈ Φ}. If r(φ) = ∅, φ has no quantiﬁcation. So the
proposition is correct by Lemma A.1. Suppose otherwise. Write φ = G[Q1 Q2 · · · Qn F ]where Q1 Q2 · · · Qn F indicates a single occurrence of a factor in G.71 We assume Q1 = ∃x
(Q1 = ∀x is similarly treated). We also assume that bound variables are renamed to avoid
name clash. Then G[∃xQ2 · · · Qn F ] is equivalent to ∃xG[Q2 · · · Qn F ] in light of the validityof (∃xA) ∧ B = ∃x(A ∧ B) and (∃xA) ∨ B = ∃x(A ∨ B) when B contains no free x.
PDB (φ) = PDB (G[Q1 Q2 · · · Qn F ])= PDB (∃xG[ Q2 · · · Qn F [x] ])===lim PDB (G[ Q2 · · · Qn F [t1 ]] ∨ · · · ∨ G[Q2 · · · Qn F [tk ] ])k→∞lim PDB (G[ Q2 · · · Qn F [t1 ] ∨ · · · ∨ Q2 · · · Qn F [tk ] ])k→∞lim PF ({ν ∈ ΩF | MDB (ν) |= G[ Q2 · · · Qn F [t1 ] ∨ · · · ∨ Q2 · · · Qn F [tk ] ]})k→∞(by induction hypothesis)= PF ({ν ∈ ΩF | MDB (ν) |= ∃xG[Q2 · · · Qn F [x]]})= PF ({ν ∈ ΩF | MDB (ν) |= φ})

Q.E.D.

We next prove a theorem on the iﬀ deﬁnition introduced in Section 4. Distribution
semantics considers the program DB = F ∪ R as a set of inﬁnitely many ground deﬁnite
clauses such that F is a set of facts (with a probability measure PF ) and R a set of rules,
and no clause head in R appears in F . Put
defhead(R) = {B | B appears in R as a clause head}.
For B ∈ head(R), let B ← Wi (i = 1, 2, . . .) be an enumeration of clauses about B in R.
Deﬁne iﬀ (B), the iﬀ (if-and-only-if ) form of rules about B in DB 72 bydef

iﬀ (B) = B ↔ W1 ∨ W2 ∨ · · ·Since MDB (ν) is a least Herbrand model, the following is obvious.
Lemma A.2 For B in head(R) and ν ∈ ΩF , MDB (ν) |= iﬀ (B).
Theorem A.1 below is about iﬀ (B). It states that at general level, both sides of the iﬀ
deﬁnition p(x) ↔ ∃y1 (x = t1 ∧ W1 ) ∨ · · · ∨ ∃yn (x = tn ∧ Wn ) of p(·) coincide as randomvariables whenever x is instantiated to a ground term.
Theorem A.1 Let iﬀ (B) = B ↔ W1 ∨ W2 ∨ · · · be the iﬀ form of rules about B ∈ head(R).
PDB (iﬀ (B)) = 1 and PDB (B) = PDB (W1 ∨ W2 ∨ · · ·).
71. For an expression E, E[γ] means that γ may occur in the speciﬁed positions of E. If γ1 ∨ γ2 in E[γ1 ∨ γ2 ]indicates a single occurrence of γ1 ∨ γ2 in a positive boolean formula E, E[γ1 ∨ γ2 ] = E[γ1 ] ∨ E[γ2 ] holds.
72. This deﬁnition is diﬀerent from the usual one (Lloyd, 1984; Doets, 1994) as we are here talking at groundlevel. W1 ∨ W2 ∨ · · · is true if and only if one of the disjuncts is true.

447

Sato & Kameya(Proof)
PDB (iﬀ (B)) = PDB ({ω ∈ ΩDB | ω |= B ∧ (W1 ∨ W2 ∨ · · ·)})+PDB ({ω ∈ ΩDB | ω |= ¬B ∧ ¬(W1 ∨ W2 ∨ · · ·)})=

lim PDB ({ω ∈ ΩDB | ω |= B ∧

k→∞

k
_

Wi })

i=1

+ lim PDB ({ω ∈ ΩDB | ω |= ¬B ∧ ¬
k→∞

=

k
_

Wi })

i=1

lim PF ({ν ∈ ΩF | MDB (ν) |= B ∧

k→∞

k
_

Wi })

i=1

+ lim PF ({ν ∈ ΩF | MDB (ν) |= ¬B ∧ ¬
k→∞

k
_

Wi })

i=1

(Lemma A.1)
= PF ({ν ∈ ΩF | MDB (ν) |= iﬀ (B)})= PF (ΩF ) (Lemma A.2)
= 1
It follows from PDB (iﬀ (B)) = 1 that
PDB (B) = PDB (B ∧ iﬀ (B)) = PDB (W1 ∨ W2 ∨ · · ·).

Q.E.D.

We then prove a proposition useful in probability computation. Let ψDB (B) be thesupport set for an atom B introduced in Section 4 (it is the set of all explanations for B). InW
the sequel, B is a ground atom. Write ψDB (B) = {S1 , S2 , . . .} and ψDB (B) = S1 ∨S2 ∨· · ·73
Deﬁne a set ΛB by_def
ΛB = {ω ∈ ΩDB | ω |= B ↔ ψDB (B)}.
W

Proposition A.2 For every B ∈ head(R), PDB (ΛB ) = 1 and PDB (B) = PDB ( ψDB (B)).
(Proof) We ﬁrst prove PDB (ΛB ) = 1 but the proof exactly parallels that of Theorem A.1
except that W1 ∨ W2 ∨ · · · is replaced by S1 ∨ S2 ∨ · · · using the fact that B ↔ S1 ∨ S2 ∨ · · ·is true in every least Herbrand model of the form MDB (ν). Then from PDB (ΛB ) = 1, we
have
PDB (B) = PDB (B ∧ (B ↔
= PDB (__

ψDB (B)).

ψDB (B)))
Q.E.D.

Finally, we show that distribution semantics is a probabilistic extension of the traditional
least Herbrand model semantics in logic programming by proving Theorem A.2. It says thatthe probability mass is distributed exclusively over possible least Herbrand models.
Deﬁne Λ as the set of least Herbrand models generated by ﬁxing R and varying a subset
of F in the program DB = F ∪ R. In symbols,
73. For a set K = {E1 , E2 , . . .} of formulas,W

K denotes a (-n inﬁnite) disjunction E1 ∨ E2 ∨ · · ·448

Parameter Learning of Logic Programs for Symbolic-statistical ModelingdefΛ = {ω ∈ ΩDB | ω = MDB (ν) for some ν ∈ ΩF }.
Note that as Λ is merely a subset of ΩDB , we cannot conclude PDB (Λ) = 1 a priori, but thenext theorem, Theorem A.2, states PDB (Λ) = 1, i.e. distribution semantics distributes theprobability mass exclusively over Λ, i.e. possible least Herbrand models.
To prove the theorem, we need some preparations. Recalling that atoms outside head(R)∪F have no chance of being proved from DB, we introduce
Λ0 = {ω ∈ ΩDB | ω |= ¬D for every ground atom D 6∈ head(R) ∪ F }.
def

For a Herbrand interpretation ω ∈ ΩDB , ω|F (∈ ΩF ) is the restriction of ω to those atoms
in F .
Lemma A.3 Let ω ∈ ΩDB be a Herbrand interpretation.
Wω = MDB (ν) for some ν ∈ ΩF iﬀ ω ∈ Λ0 and ω |= B ↔ ψDB (B) for every B ∈ head(R).
(Proof) Only-if part is immediate from the property of the least Herbrand model. For
if-part, suppose ω satisﬁes the right hand side. We show that ω = MDB (ω|F ). As ω and
MDB (ω|F ) coincide w.r.t. atoms not in head(R), it is enough to prove that they also giveW
the same truth values to atoms in head(R). Take B ∈ head(R) and write ψDB (B) =S1 ∨ S2 ∨ · · · Suppose ω |= B ↔ S1 ∨ S2 ∨ · · · Then if ω |= B, we have ω |= Sj for some j,thereby ω|F |= Sj , and hence MDB (ω|F ) |= Sj , which implies MDB (ω|F ) |= B. Otherwise
ω |= ¬B. So ω |= ¬Sj for every j. It follows that MDB (ω|F ) |= ¬B. Since B is arbitrary,
we conclude that ω and MDB (ω|F ) agree on the truth values assigned to atoms in head(R)as well.
Q.E.D.
Theorem A.2 PDB (Λ) = 1.
(Proof) From Lemma A.3, we have
Λ = {ω ∈ ΩDB | ω = MDB (ν) for some ν ∈ ΩF }= Λ0 ∩\

ΛB .

B∈head(R)PDB (ΛB ) = 1 by Proposition A.2. To prove PDB (Λ0 ) = 1, let D1 , D2 , . . . be an enumerationof atoms not belonging to head(R) ∪ F . They are not provable from DB = F ∪ R, andhence false in every least Herbrand model MDB (ν) (ν ∈ ΩF ). So
PDB (Λ0 ) ==lim PDB ({ω ∈ ΩDB | ω |= ¬D1 ∧ · · · ∧ ¬Dm })m→∞lim PF ({ν ∈ ΩF | MDB (ν) |= ¬D1 ∧ · · · ∧ ¬Dm })m→∞= PF (ΩF ) = 1.
Since a countable conjunction of measurable sets of probability measure one has also
probability measure one, it follows from PDB (ΛB ) = 1 for every B ∈ head(R) and PDB (Λ0 ) =1 that PDB (Λ) = 1.
Q.E.D.

449

Sato & Kameya

Appendix B. The MAR (missing at random) Condition
In the original formulation of the EM algorithm by Dempster et al. (1977), it is assumed
that there exists a many-to-one mapping y = χ(x) from a complete data x to an incomplete(observed) data y. In the case of parsing, x is a parse tree and y is the input sentence and x
uniquely determines y. In this paper, the uniqueness condition ensures the existence of sucha many-to-one mapping from explanations to observations. We however sometimes face asituation where there is no such many-to-one mapping from complete data to incomplete
data but nonetheless we wish to apply the EM algorithm.
This dilemma can be solved by the introduction of a missing-data mechanism whichmakes a complete data incomplete. The missing-data mechanism, m, has a distribution
gφ (m | x) parameterized by φ and y, the observed data, is described as y = χm (x). It says
x becomes incomplete y by m. The correspondence between x and y, i.e. {hx, yi | ∃m(y =χm (x))} naturally becomes many-to-many.
Rubin (1976) derived two conditions on gφ (data are missing at random and data areobserved at random) collectively called the MAR (missing at random) condition, and showedthat if we assume a missing-data mechanism behind our observations that satisﬁes the MAR
condition, we may estimate parameters of the distribution over x by simply applying the
EM algorithm to y, the observed data.
We adapt the MAR condition to parameterized logic programs as follows. We keep agenerative model satisfying the uniqueness condition that outputs goals G such as parsetrees. We further extend the model by additionally inserting a missing-data mechanism
m between G and our observation O like O = χm (G) and assume m satisﬁes the MAR
condition. Then the extended model has a many-to-many correspondence between explanations and observations, and generates non-exclusive observations such that P (O ∧O0 ) > 0
PP(O 6= O0 ), which causes O P (O) ≥ 1 where P (O) = G:∃m O=χm (G) PDB (G). Thanks to
the MAR condition however, we are still allowed to apply the EM algorithm to such nonexclusive observations. Put it diﬀerently, even if the uniqueness condition is seemingly
destroyed, the EM algorithm is applicable just by (imaginarily) assuming a missing-datamechanism satisfying the MAR condition.

References
Abney, S. (1997). Stochastic attribute-value grammars. Computational Linguistics, 23 (4),597–618.
Arimura, H. (1997). Learning acyclic ﬁrst-order horn sentences from entailment. In
Proceedings of the Eighth International Workshop on Algorithmic Learning Theory.
Ohmsha/Springer-Verlag.
Bacchus, F., Grove, A., Halpern, J., & Koller, D. (1996). From statistical knowledge bases
to degrees of belief. Artiﬁcial Intelligence, 87, 75–143.
Baker, J. K. (1979). Trainable grammars for speech recognition. In Proceedings of Spring
Conference of the Acoustical Society of America, pp. 547–550.

450

Parameter Learning of Logic Programs for Symbolic-statistical Modeling

Beil, F., Carroll, G., Prescher, D., Riezler, S., & Rooth, M. (1999). Inside-Outside estimationof a lexicalized PCFG for German. In Proceedings of the 37th Annual Meeting of the
Association for Computational Linguistics (ACL’99), pp. 269–276.
Breese, J. S. (1992). Construction of belief and decision networks. Computational Intelligence, 8 (4), 624–647.
Carroll, G., & Rooth, M. (1998). Valence induction with a head-lexicalized PCFG. In Proceedings of the 3rd Conference on Empirical Methods in Natural Language Processing(EMNLP 3).
Castillo, E., Gutierrez, J. M., & Hadi, A. S. (1997). Expert Systems and Probabilistic
Network Models. Springer-Verlag.
Charniak, E., & Carroll, G. (1994). Context-sensitive statistics for improved grammatical language models. In Proceedings of the 12th National Conference on Artiﬁcial
Intelligence (AAAI’94), pp. 728–733.
Chi, Z., & Geman, S. (1998). Estimation of probabilistic context-free grammars. Computational Linguistics, 24 (2), 299–305.
Chow, Y., & Teicher, H. (1997). Probability Theory (3rd ed.). Springer.
Clark, K. (1978). Negation as failure. In Gallaire, H., & Minker, J. (Eds.), Logic andDatabases, pp. 293–322. Plenum Press.
Cormen, T., Leiserson, C., & Rivest, R. (1990). Introduction to Algorithms. The MIT Press.
Cussens, J. (1999). Loglinear models for ﬁrst-order probabilistic reasoning. In Proceedings ofthe 15th Conference on Uncertainty in Artiﬁcial Intelligence (UAI’99), pp. 126–133.
Cussens, J. (2001). Parameter estimation in stochastic logic programs. Machine Learning,
44 (3), 245–271.
D’Ambrosio, B. (1999). Inference in Bayesian networks. AI Magazine, summer, 21–36.
Dekhtyar, A., & Subrahmanian, V. S. (1997). Hybrid probabilistic programs. In Proceedingsof the 14th International Conference on Logic Programming (ICLP’97), pp. 391–405.
Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from incomplete
data via the EM algorithm. Royal Statistical Society, B39 (1), 1–38.
Doets, K. (1994). From Logic to Logic Programming. The MIT Press.
Flach, P., & Kakas, A. (Eds.). (2000). Abduction and Induction – Essays on Their Relation
and Integration. Kluwer Academic Publishers.
Frish, A., & Haddawy, P. (1994). Anytime deduction for probabilistic logic. Journal of
Artiﬁcial Intelligence, 69, 93–122.
Fujisaki, T., Jelinek, F., Cocke, J., Black, E., & Nishino, T. (1989). A probabilistic parsing
method for sentence disambiguation. In Proceedings of the 1st International Workshop
on Parsing Technologies, pp. 85–94.
Japan EDR, L. (1995). EDR electronic dictionary technical guide (2nd edition). Technicalreport, Japan Electronic Dictionary Research Institute, Ltd.

451

Sato & Kameya

Kakas, A. C., Kowalski, R. A., & Toni, F. (1992). Abductive logic programming. Journalof Logic and Computation, 2 (6), 719–770.
Kameya, Y. (2000). Learning and Representation of Symbolic-Statistical Knowledge (inJapanese). Ph. D. dissertation, Tokyo Institute of Technology.
Kameya, Y., & Sato, T. (2000). Eﬃcient EM learning for parameterized logic programs.
In Proceedings of the 1st Conference on Computational Logic (CL2000), Vol. 1861 of
Lecture Notes in Artiﬁcial Intelligence, pp. 269–294. Springer.
Kita, K. (1999). Probabilistic Language Models (in Japanese). Tokyo Daigaku Syuppan-kai.
Koller, D., McAllester, D., & Pfeﬀer, A. (1997). Eﬀective Bayesian inference for stochastic programs. In Proceedings of 15th National Conference on Artiﬁcial Intelligence(AAAI’97), pp. 740–747.
Koller, D., & Pfeﬀer, A. (1997). Learning probabilities for noisy ﬁrst-order rules. In Proceedings of the 15th International Joint Conference on Artiﬁcial Intelligence (IJCAI’97),
pp. 1316–1321.
Kyburg, H. (1994). Uncertainty logics. In Gabbay, D., Hogger, C., & Robinson, J. (Eds.),
Handbook of Logics in Artiﬁcial Intelligence and Logic Programming, pp. 397–438.
Oxford Science Publications.
Laﬀerty, J. (1993). A derivation of the Inside-Outside Algorithm from the EM algorithm.
Technical report, IBM T.J.Watson Research Center.
Lakshmanan, L. V. S., & Sadri, F. (1994). Probabilistic deductive databases. In Proceedingsof the 1994 International Symposium on Logic Programming (ILPS’94), pp. 254–268.
Lari, K., & Young, S. J. (1990). The estimation of stochastic context-free grammars using
the Inside-Outside algorithm. Computer Speech and Language, 4, 35–56.
Li, Z., & D’Ambrosio, B. (1994). Eﬃcient inference in Bayes networks as a combinatorial
optimization problem. International Journal of Approximate Reasoning, 11, 55–81.
Lloyd, J. W. (1984). Foundations of Logic Programming. Springer-Verlag.
Lukasiewicz, T. (1999). Probabilistic deduction with conditional constraints over basic
events. Journal of Artiﬁcial Intelligence Research, 10, 199–241.
Manning, C. D., & Schütze, H. (1999). Foundations of Statistical Natural Language Processing. The MIT Press.
McLachlan, G. J., & Krishnan, T. (1997). The EM Algorithm and Extensions. Wiley
Interscience.
Muggleton, S. (1996). Stochastic logic programs. In de Raedt, L. (Ed.), Advances inInductive Logic Programming, pp. 254–264. IOS Press.
Ng, R., & Subrahmanian, V. S. (1992). Probabilistic logic programming. Information andComputation, 101, 150–201.
Ngo, L., & Haddawy, P. (1997). Answering queries from context-sensitive probabilisticknowledge bases. Theoretical Computer Science, 171, 147–177.
Nilsson, N. J. (1986). Probabilistic logic. Artiﬁcial Intelligence, 28, 71–87.
452

Parameter Learning of Logic Programs for Symbolic-statistical Modeling

Pearl, J. (1988). Probabilistic Reasoning in Intelligent Systems. Morgan Kaufmann.
Pereira, F. C. N., & Schabes, Y. (1992). Inside-Outside reestimation from partially bracketed
corpora. In Proceedings of the 30th Annual Meeting of the Association for Computational Linguistics (ACL’92), pp. 128–135.
Pereira, F. C. N., & Warren, D. H. D. (1980). Deﬁnite clause grammars for language analysis— a survey of the formalism and a comparison with augmented transition networks.
Artiﬁcial Intelligence, 13, 231–278.
Pfeﬀer, A., & Koller, D. (2000). Semantics and inference for recursive probability models. In
Proceedings of the Seventh National Conference on Artiﬁcial Intelligence (AAAI’00),
pp. 538–544.
Poole, D. (1993). Probabilistic Horn abduction and Bayesian networks. Artiﬁcial Intelligence, 64 (1), 81–129.
Pynadath, D. V., & Wellman, M. P. (1998). Generalized queries on probabilistic context-freegrammars. IEEE Transaction on Pattern Analysis and Machine Intelligence, 20 (1),65–77.
Rabiner, L. R. (1989). A tutorial on hidden markov models and selected applications inspeech recognition. Proceedings of the IEEE, 77 (2), 257–286.
Rabiner, L. R., & Juang, B. (1993). Foundations of Speech Recognition. Prentice-Hall.
Ramakrishnan, I., Rao, P., Sagonas, K., Swift, T., & Warren, D. (1995). Eﬃcient tablingmechanisms for logic programs. In Proceedings of the 12th International Conference
on Logic Programming (ICLP’95), pp. 687–711. The MIT Press.
Reddy, C., & Tadepalli, P. (1998). Learning ﬁrst-order acyclic horn programs from entailment. In Proceedings of the 15th International Conference on Machine Learning;(and Proceedings of the 8th International Conference on Inductive Logic Programming). Morgan Kaufmann.
Riezler, S. (1998). Probabilistic Constraint Logic Programming. Ph.D. thesis, UniversitätTübingen.
Rubin, D. (1976). Inference and missing data. Biometrika, 63 (3), 581–592.
Sagonas, K., T., S., & Warren, D. (1994). XSB as an eﬃcient deductive database engine.
In Proceedings of the 1994 ACM SIGMOD International Conference on Managementof Data, pp. 442–453.
Sato, T. (1995). A statistical learning method for logic programs with distribution semantics.
In Proceedings of the 12th International Conference on Logic Programming (ICLP’95),
pp. 715–729.
Sato, T. (1998). Modeling scientiﬁc theories as PRISM programs. In Proceedings of ECAI’98
Workshop on Machine Discovery, pp. 37–45.
Sato, T. (2001). Minimum likelihood estimation from negative examples in statistical abduction. In Proceedings of IJCAI-01 workshop on Abductive Reasoning, pp. 41–47.
Sato, T., & Kameya, Y. (1997). PRISM: a language for symbolic-statistical modeling.
In Proceedings of the 15th International Joint Conference on Artiﬁcial Intelligence(IJCAI’97), pp. 1330–1335.
453

Sato & Kameya

Sato, T., & Kameya, Y. (2000). A Viterbi-like algorithm and EM learning for statisticalabduction. In Proceedings of UAI2000 Workshop on Fusion of Domain Knowledge
with Data for Decision Support.
Sato, T., Kameya, Y., Abe, S., & Shirai, K. (2001). Fast EM learning of a family of PCFGs.
Titech technical report (Dept. of CS) TR01-0006, Tokyo Institute of Technology.
Shen, Y., Yuan, L., You, J., & Zhou, N. (2001). Linear tabulated resolution based on Prolog
control strategy. Theory and Practice of Logic Programming, 1 (1), 71–103.
Sterling, L., & Shapiro, E. (1986). The Art of Prolog. The MIT Press.
Stolcke, A. (1995). An eﬃcient probabilistic context-free parsing algorithm that computes
preﬁx probabilities. Computational Linguistics, 21 (2), 165–201.
Tamaki, H., & Sato, T. (1984). Unfold/fold transformation of logic programs. In Proceedingsof the 2nd International Conference on Logic Programming (ICLP’84), Lecture Notes
in Computer Science, pp. 127–138. Springer.
Tamaki, H., & Sato, T. (1986). OLD resolution with tabulation. In Proceedings of the 3rdInternational Conference on Logic Programming (ICLP’86), Vol. 225 of Lecture Notes
in Computer Science, pp. 84–98. Springer.
Tanaka, H., Takezawa, T., & Etoh, J. (1997). Japanese grammar for speech recognition
considering the MSLR method. In Proceedings of the meeting of SIG-SLP (Spoken
Language Processing), 97-SLP-15-25, pp. 145–150. Information Processing Society ofJapan. In Japanese.
Uratani, N., Takezawa, T., Matsuo, H., & Morita, C. (1994). ATR integrated speech andlanguage database. Technical report TR-IT-0056, ATR Interpreting Telecommunications Research Laboratories. In Japanese.
Warren, D. S. (1992). Memoing for logic programs. Communications of the ACM, 35 (3),
93–111.
Wetherell, C. S. (1980). Probabilistic languages: a review and some open questions. Computing Surveys, 12 (4), 361–379.
White, H. C. (1963). An Anatomy of Kinship. Prentice-Hall.
Zhang, N., & Poole, D. (1996). Exploiting causal independence in Bayesian network inference. Journal of Artiﬁcial Intelligence Research, 5, 301–328.

454