Approximate Probabilistic Inference with Bounded Error for Hybrid Probabilistic Logic Programming

1 Introduction

Probabilistic logics are gaining popularity as modelling and reasoning tools, since they combine the power of logic to represent knowledge with the ability of probability theory to deal with uncertainty. Especially, probabilistic languages based on logic programming (LP) have become popular, as they enable knowledge representation by employing similar structure as the popular Bayesian Networks (BNs) [Pearl, 1988], but lift this from a propositional to a first-order level, thus dramatically increasing expressive power. Examples of such languages are ICL [Poole, 2008] and ProbLog [Fierens et al., 2015]. The need for such expressive languages emerged from the fact that in many areas more and more data become available, which does not only imply uncertainty, but often provides rich structure in terms of relations between entities. Probabilistic logic methods have been applied to a wide range of problem domains. Examples include link and node prediction in metabolic networks [Kimmig and Costa, 2012], dealing with a potentially unknown number of relations between multiple objects by a robot [Moldovan et al., 2012] and information fusion in the safety and security domain [Michels et al., 2013].

Probabilistic inference is an inherently computationally hard problem. However, it is well known that by exploiting logical structure as present in many distributions, exact inference is rendered feasible for many realistic problems [Chavira and Darwiche, 2008]. This is achieved by viewing the problem as a weighted model counting (WMC) problem, utilising previous work on SAT solving and model counting. These inference methods are especially suitable for probabilistic LP, as here models always have a rich logical structure (e.g. [Fierens et al., 2015]). Unfortunately, the above-mentioned structure-sensitive inference methods are usually restricted to discrete, finite distributions. Many real-world problems, however, also require continuous variables, i.e. they are best modelled using a hybrid distribution. Inference for hybrid distributions is however either made possible by significantly limiting the expressivity of the languages, e.g. [Gutmann et al., 2010; Chistikov et al., 2015], or by resorting to approximate inference methods, e.g. [Gutmann et al., 2011].

Sampling algorithms - the most widely used methods for approximate inference - are however known to converge slowly when dealing with (near) deterministic probabilities and cannot handle observed rare events well. Also inference methods based on the approximate representation of the actual continuous distribution, e.g. [Sanner and Abbasnejad, 2012; Belle et al., 2015a], suffer from observed rare events. All such methods furthermore share the drawback that they only provide weak guarantees about the quality of the estimates. In domains such as diagnosis and crisis management, where wrong decisions may have a huge impact, guarantees about the quality of approximations are highly desirable. The need to deal with rare events is also characteristic for such domains.

In this paper, we extend probabilistic LP with continuous distributions in a way that poses little restriction on the definition and use of continuous variables. Although the inference task is not computable, we show that we can provide an approximation with bounded error, which essentially provides more information about the quality of the result than the approximation methods mentioned above. The algorithm iteratively computes approximations with decreasing error, adopting a coarse-to-fine approach, making use of information from previous runs and also exploiting the logical structure of the problem to effectively decrease the error of an approximation. A comparison is made with sampling methods, showing that in spite of the much stronger guarantees on the quality of the result, our method is competitive and even superior for hybrid problems that possess a rich logical structure and concern observed rare events.

2 Hybrid Probabilistic Logic Programming

We here introduce the language used, which is similar to Hybrid ProbLog [Gutmann et al., 2010], but allows for more expressive constraints on continuous variables.

Example 1 We consider a diagostic problem. Suppose there are three reasons why a component i can fail: the component breaks with some probability, the temperature got too high, or a sub-component fails. Additionally, there is a cooler that cools down the temperature of all components by an amount that is uncertain as well. This cooler can also fail, indicated by the variable NoC. We can express this with the following rules:

For i = 0, the component does not depend on other components, so it can only fail because of the first two reasons. We can use arbitrary distributions for the continuous variables, for example:

The distributions for Break(i) and Limit(i) are assumed to be the same for all components i. Assuming that a given component fails, a natural query is what is the probability that the first component also fails, e.g. P(fails(0)♣fails(9)).

Formally, we consider a finite number of random variables X↪Y↪Z↪ ▷ ▷ ▷, which have a support, i.e. the range of values on which the distribution is defined. These ranges are in principle arbitrary; in this paper we consider binary and realvalued random variables as typical representatives of discrete and continuous variables. For the continuous variables, we assume that their cumulative distributions functions (CDFs) can be computed exactly, thereby ignoring the potential error of typical implementations of such functions as such errors are typically negligible and predictable. We assume all random variables to be independent initially, and express dependencies by the logical structure, as common in probabilistic LP.

We use LP rules with constraints on the random variables as building blocks for the bodies, with the only restriction that satisfiability of those constraints remains decidable. Concretely, for real-numbered random variables we allow inequalities on combinations of variables. These combinations are not necessarily linear as for instance required for Conditional Linear Gaussians [Lauritzen, 1992], because constraints that consist of multiplications with real numbers are decidable. Note that this is a much weaker requirement than e.g. made for Hybrid ProbLog, which only allows a single random variable in a constraint, such that not only the constraint is decidable, but the probability of the constraint is computable as well.

As usual for probabilistic LP, the inference task we consider is the computation of the probability that a grounded predicate holds, given a conjunction of other grounded predicates.

3 Iterative Hybrid Probabilistic Model Counting

We here present a novel algorithm, the iterative hybrid probabilistic model counting (IHPMC) algorithm; it is an anytime algorithm that provides approximations with bounded error for the language outlined above.

3.1 Hybrid Probability Trees

Queries are converted to propositional formulas as is often done for probabilistic LP inference. Details of this transformation are not discussed in this paper as it is essentially the same as for ProbLog [Fierens et al., 2015]. The only difference is that the propositional formulas also include constraints on random variables. For a fixed i 6= 0, the rules in Example 1 can for instance be represented as: failsi $ Breaki = true _Temp-Cooling > Limiti _(Temp > Limiti ^ NoC = true) _ fails i-1.

To perform inference on such formulas we introduce the concept of Hybrid Probability Trees (HPTs). They are similar to the principles used in solvers for binary weighted model counting. Whereas in previous work the concept is applied to discrete, finite ranges, we extend the concept to continuous random variables. The basic idea is to split the variable's range into two parts at each edge and condition the formula on the split made. To support continuous variables, we allow further splits of the same variable at deeper levels. Then, for each choice taken at a certain child, the probability is assigned that the random variable takes a value in the chosen partition, conditioned on the range chosen for the random variable earlier on the path. Such probabilities can be computed for continuous random variables by their CDFs, because we assume independence of the random variables. The splitting is repeated until the formula simplifies to ? or >, which we for now assume to happen in a finite number of steps. By attaching a corresponding probability to each branch, this can be used to compute the probability of the event represented by the formula.

Definition 1 (Hybrid Probability Trees (HPT)) An HPT is a binary tree with at each node n a propositional formula 'n and for all random variables X a range denoted by range(n↪X). If r is the root node, it holds that range(r↪X) equals the support of X. Now let n be some node with children ¶c1↪ c2♦. To each edge n ! ci, i 2 ¶1↪ 2♦, a particular range ⌧ni is associated to a fixed random variable Y such that ⌧n1 [ ⌧n2 = range(n↪Y) and ⌧n1 \ ⌧n2 = ?. It holds that range(ci↪Y) = ⌧ni and range(ci↪Z) = range(n↪Z) if Y 6= Z, with i 2 ¶1↪ 2♦. Furthermore, the formula 'i associated to ci equals 'n simplified by the restrictions on the

?

Figure 1: Example HPT

range on Y at ci, i.e. by observing that some primitive constraints can be replaced by > or ?. Furthermore, to each edge n ! ci, with i 2 ¶1↪ 2♦, a probability pni is assigned, such that pni = P(Y 2 ⌧ni ♣ Y 2 range(n↪Y)). Finally, if l is a leaf node, then it holds that 'l = > or 'l = ?.

Given a particular HPT, it is straightforward to compute the probability of the event represented by the formula at the root node. First, the probability of a leaf is by definition 0 or 1. Furthermore, due to the properties of the tree it is easy to show that for each non-leaf node n the probability of 'n given the ranges of random variables can be computed by pn = pn1 · p1 + pn2 · p2, where p1 and p2 are the probabilities associated to n's children.

Example 2 Consider a simplified diagnosis problem:

Here NoC represents that the cooling fails, implying a lower allowed temperature for the component. The temperature is modelled by T ⇠ N(20▷0↪ 5▷0). The event that the component fails is then represented by the formula (NoC = true ^ T > 20▷0) _ T > 30▷0.

A possible HPT is given in Figure 1. At each node we simplify the formula as much as possible given the choices made. Note that the probability assigned to the edge with choice T 2 (-1↪ 20] is not P (T 2 (-1↪ 20]) = 0▷5, as the range of T is already restricted on the path before. Therefore, it is P (T 2 (-1↪ 20] ♣ T 2 (-1↪ 30]) ⇡ 0▷5116. From this tree, the probability of component failure can be computed by 0▷9772 · 0▷01 · 0▷4884 + 0▷0228 ⇡ 0▷0276.

3.2 Partially Evaluated Hybrid Probability Trees

The concept of HPT is restricted to problems for which exact inference is possible. For continuous events, which cannot be represented by hyperrectangles, we can never simplify all leaves to ? or >. This problem is tackled by using partially evaluated HPTs (PHPTs). For PHPTs we drop the requirement that all leaves must contain ? or > and allow arbitrary formulas. One can therefore further extend a PHPT either ad infinitum or until one finds an HPT. In the next section we show how PHPTs can be used to compute approximations with known maximal error.

>

Example 3 Suppose the limit of the temperature the component can bear is uncertain as well, modelled by L ⇠ N(30▷0↪ 5▷0). A PHPT for the event that the temperature is above the limit is depicted in Figure 2. It is impossible to find a finite HPT for this problem, as the event does not have the shape of a hyperrectangle. Note that the formula T > L cannot be simplified directly, in case the range of only one of the variables is split. This is in contrast to binary versions of model counting algorithms, for which the choices of variable values are completely reflected in the simplified formula.

3.3 Approximating Probabilities by Partially Evaluated Hybrid Probability Trees

Bounds on Event Probabilities

Each leaf of a PHPT with a formula which is not ? or > corresponds to an area of the sample space which is only partially part of the event associated to the query. While this does not provide an exact probability of this event, it does provide a bound on this probability by assuming that the probability of such areas is 0▷0 or 1▷0 respectively. In this way, PHPTs can be used to compute probability bounds, denoted by P and P .

Example 4 The PHPT in Figure 2 contains only one leaf with >, whereas all other leaves contribute at most partially to the event's probability. The lower bound of the event's probability is therefore P(T > L) = 0▷5 · 0▷9772 = 0▷4886. Analogously, as there is only one path which certainly does not contribute to the probability, the probability's upper bound is P (T > L) = 1 - 0▷5 · 0▷0228 = 0▷9886.

Proposition 1 For any event e and PHPT for e: P(e)  P (e)  P (e).

All proofs are given in [Michels, 2016, Chapter 6, Appendix A]. Furthermore, we can always achieve arbitrary precision, by evaluating the PHPT sufficiently deep.

Proposition 2 For every event e and maximal error ✏, there is a PHPT, such that: P (e) - P (e)  ✏ ^ P (e) - P (e)  ✏.

Bounds on Conditional Event Probabilities

Conditional probabilities are defined as P(e ♣ o) = P(e ^ o)◁P(o). As this definition requires the non-conditional

Input: Events e and o, maximal error ✏ Result: Approximation of P(e ♣ o) with max. error ✏ 1 hpt p = PHPT with single root node e ^ o 2 hpt n = PHPT with single root node ≻e ^ o 3 p = 0▷0; p = 1▷0 4 while (p - p)◁2 > ✏ do 5 eval hpt = choose hpt(hpt p↪ hpt n) 6 eval node = choose node(eval hpt) 7 branch rv = choose rv(eval node) 8 (rv l ↪ rv r) = choose part(branch rv↪ eval node) 9 add children to eval node according to (rv l ↪ rv r) 10 p = lower(hpt p)◁(lower(hpt p) +upper(hpt n)) 11 p = upper(hpt p)◁(upper(hpt p) + lower(hpt n)) 12 end 13 return (p+p)◁2

probabilities of two events, we cannot compute bounds on this probability using a single PHPT, but have to use a second one. We can furthermore compute tighter bounds by making use of the fact that the bounds of e ^ o and e are not independent: an area cannot at the same time be outside of e, but within e ^ o. This can be exploited by expressing the bounds of conditional probabilities in terms of bounds on two nonconditional probabilities [Michels et al., 2015]:

Proposition 3

Also for conditional probabilities we can therefore achieve arbitrary precisions.

Proposition 4 For all events e and o and every maximal error ✏, there are PHPTs, such that: P(e ♣ o) - P(e ♣ o)  ✏ ^P(e ♣ o) -P(e ♣ o)  ✏.

3.4 Anytime Inference by Iterative Hybrid Probability Tree Evaluation

The anytime approximation IHPMC algorithm based on PHPTs is given in Algorithm 1. We start with two PHPTs with initial root nodes e ^ o and ≻e ^ o, which we use to compute in each iteration an approximation according to Proposition 3. In case no evidence is present, we can use a single PHPT with formula e. In the general case, at each iteration we first choose one of the two PHPTs and then, in this PHPT, we choose a non-evaluated node to evaluate further. From the formula of this node we choose a random variable that will be used for splitting and finally we choose a partitioning for this random variable.

Note that in the algorithm we could also compute components of the formula independently if they do not share variables, as is usually done in WMC solvers. For brevity this is left out here. While the basic algorithm is straightforward, there are four non-deterministic choices at each iteration which largely determines the effectiveness of the algorithm. In the following, we describe general heuristics which seem to be very effective for most problems.

As the first choice, we pick the PHPT which has the highest difference between lower and upper bound, i.e. , the largest error. Second, we choose the non-evaluated node which represents the sub-space with the highest probability mass, as these nodes have the most potential for reducing the error. A more sophisticated heuristic, that tries to predict how much effort is required to reduce the error of the node's formula, offered no improvement in experiments. The reason is that an accurate prediction takes as much time as just trying to reduce the error and continuing with another branch if the error for the firstly chosen one cannot be reduced quickly.

In contrast to the choice of the node, the heuristic for choosing a random variable must be a good predictor for how much this simplifies the formula, as this choice influences all descendants of the node in the tree. We use a heuristics preferring variables occurring more frequently, as also common for binary model counting [Sang et al., 2005, Section 3.2]. For continuous variables we however have to adapt the heuristic to make sure that continuous variables are not selected repeatedly without an opportunity to eliminate a primitive constraint, and thereby simplifying the formula. For instance, in Example 3, simplification of the constraint can only be achieved if both variables are selected in an alternating order.

For the binary variables there is obviously only a single partitioning possible (¶true♦ and ¶false♦). The partitioning of a continuous random variable amounts to picking a point, on which the range is split further. Initially, we find points that can lead to a simplified sub-formula. For instance, in Figure 2 the choice L 2 (20↪ 1) in the left branch can be used to simplify the formula to ?, given the choice for T above in the tree. In case there are multiple of those points, we count how many sub-formulas can be simplified by choosing this point. There are also cases in which there is no such point, as at the top or the rightmost evaluated node in Figure 2. In this case we use a heuristic trying to maximise the probability of the path eliminating the constraint.

Proposition 5 IHPMC (Algorithm 1) terminates, i.e. it always finds an approximation with the desired error bound in finite time, given the heuristics as described above.

4 Experiments

4.1 Set-up

The algorithm was implemented in the functional programming language Haskell 1 . A comparison was made with the sampler implementations of BLOG 0.8 2 [Milch et al., 2005] and Distributional clauses 3 [Gutmann et al., 2011] (DC); both are optimised for problems with logical structure. We compare to the Likelihood weighted (BLOG LW) and the naive rejection sampler (BLOG RJ) of BLOG; the MCMC sampler fails completely on the problem. To allow for a comparison, we ignore the bounded error guarantee of IHPMC and just consider the error of approximations. We measure the squared error and for the sampling approaches the mean squared error over 50 runs versus inference time. All experiments were run on a Laptop with a Intel Core i3 2.4 GHz processor and 4GB of RAM. We use the following diagnostic problem:

1 http://www.steffen-michels.de/ihpmc

2 https://bayesianlogic.github.io

3 https://github.com/davidenitti/DC

The number of components (n) and probability of Break(i) (p) are varied in the experiments. For the continuous variables we use the following distributions: Temp ⇠ N(20▷0↪ 5▷0) and Limit(i) ⇠ N(θ↪ 5▷0), where θ is a parameter that we also vary during the experiments. Note that since we use a single variable Temp for all components, dependencies are created such that probabilities cannot be computed exactly in a straightforward way.

4.2 Results

Non-Conditional Probabilities With Varying Parameters

The results depicted in Figure 3a are based on a run with parameters, where the continuous cause has a much higher probability than for a realistic diagnostic problem, but is used to illustrate behaviour of the algorithm. In the first milliseconds the error of IHPMC is very unstable. Note that the actual error depicted in the figure is not the guaranteed error bound, which is usually much larger. Actually, it takes about 2▷6 seconds before the error is guaranteed to be below 0▷01. Decreasing the error bound can temporarily result in an increased error, which explains that the error does not decrease monotonically. In contrast the errors of the sampling algorithms decrease more or less monotonically, as the errors were averaged over several runs.

IHPMC can profit from the case in which the discrete cause has more impact, as shown by the experiment in Figure 3b. The error bound drops below 0▷01 already after about 6 milliseconds. The sampling algorithms can also profit from the discrete structure, but are clearly outperformed by IHPMC. Similarly, IHPMC can profit if the event is more rare, which is more realistic for a diagnostic problem (Figure 3c). Here is takes about 28 milliseconds to get an error bound of 0▷01. In this case, IHPMC cannot make use that much of the discrete structure, but in the first iterations a single path with a high probability is found that implies Break0 = false ^Temp  Limit0 ^ · · · ^ Breakn-1 = false ^ Temp  Limitn-1. This disproves failsn-1, which means the upper probability bound jumps to a small value.

Scalability

To show scalability of IHPMC we performed an experiment with a larger set of components (n = 100), where it takes about 2▷3 seconds to get a guaranteed error below 0▷01. As the results in Figure 3d show, in the beginning two of the sampling implementations outperform IHPMC, as sampling algorithms by their nature do not suffer from a high number of dimensions. However, recall that IHPMC is outperformed by algorithms providing much less information, as no guarantees about the error are provided by the sampling algorithms. After about 1▷5 second however the results of IHPMC become competitive.

Computing Probabilities Conditioned on Rare-Event

Finally, in Figure 3e, we show that IHPMC is superior if probabilities conditioned on rare events are computed, which is quite natural for a diagnostic setting. The performance of IHPMC is also affected by conditioning on the rare event; it takes about 250 milliseconds to guarantee an error below 0▷01, while for the probability of the evidence only, this takes about 15 milliseconds. The performance of the samplers however decreases significantly more. Error measures only tell half of the story as IHPMC provides a different kind of result as sampling approaches. This is illustrated in Figure 3f, where the bounds computed by IHPMC have been visualised together with point approximations produced by a sampler. Each point represents a single run of the sampler, which can or cannot be within the bounds computed by IHPMC. However, as a user of a sampling method, one only obtains a single point with, if at all, only a very weak guarantee about its error.

5 RelatedWork

Probabilistic LP has been extended with continuous distributions, in a way that still allows for exact inference, by either restricting the events expressible [Gutmann et al., 2010] or the distributions employable [Islam et al., 2012]. There are also representations of continuous distributions that are powerful enough to approximate arbitrary distributions, but still allow for exact inference, e.g. piecewise polynomials [Sanner and Abbasnejad, 2012]. A recent algorithm for inference on such representation, which exploits the logical structure of distributions was presented in [Belle et al., 2015a]. Such approximations can in principle resemble intended distributions very closely, but there is no guarantee about the approximation error of event probabilities, as a small probability of evidence can significantly increase the error made when computing conditional probabilities.

Distributional clauses [Gutmann et al., 2011] is a probabilistic logic programming language, that supports arbitrary continuous distributions and events, as also supported by various other probabilistic programming languages, such as BLOG [Milch et al., 2005]. Inference for such languages is only possible by sampling approaches, usually Markov chain Monte Carlo (MCMC) methods, which can be very effective if they are hand-tailored for specific problems. Generic MCMC frameworks are however suboptimal and can provide poor results when rare events are observed. A significant difference with our work is that even under perfect circumstances, when detailed knowledge about the approximated distributions is available, guarantees for approximations obtained by MCMC algorithms can only be given in terms of for instance the Monte Carlo standard error, which gives a probabilistic, but no hard guarantee.

Recent hashing-based sampling methods, which exploit logical structure of distributions, have been extended for hybrid distributions as well [Belle et al., 2015b; Chistikov et al., 2015]. Such methods provide a probabilistic guarantee on the quality of approximations, but are restricted to exactly computable representations of continuous distributions (piecewise-polynomials [Belle et al., 2015b] and uniform distributions [Chistikov et al., 2015]). Moreover, this

(e) P(fails(9) ♣ fails(0)), n = 10, p = 0▷0001, θ = 60▷0

(f) P(fails(9) ♣ fails(0)), n = 10, p = 0▷0001, θ = 60▷0

Figure 3: Experiment Results

remains a probabilistic guarantee on non-conditional probabilities, while IHPMC provides a hard guarantee for computed conditional probabilities.

Several approximation methods have been proposed that provide hard guarantees for discrete distributions [Poole, 1993]. Recently, very effective methods based on stateof-the-art knowledge compilation techniques have been proposed [Renkens et al., 2014; 2015; Vlasselaer et al., 2015]. In contrast to IHPMC, these algorithms are restricted to discrete problems and do not handle evidence, which are significant limitations in practice. In [Michels et al., 2015], we provided a theoretical basis for approximate inference with continuous distributions, but did not provide an efficient and practical algorithm for actually computing good approximations.

6 Conclusions

In this paper, we propose a hybrid extension for probabilistic logic programming, which allows for a much wider class of continuous distributions than existing extensions allowing for exact inference. At the same time, our extension allows one to compute approximations with bounded and arbitrarily small error. We propose a novel anytime algorithm for computing such approximations, thereby exploring the logical and continuous structure of the problem by utilising information from previous iterations. We experimentally show that our algorithm is competitive with state-of-the-art sampling algorithms, for typical relational problems, which are characterised by a moderate number of dimensions and a significant amount of structure, in spite of the fact that it provides much stronger guarantees. Our method outperforms sampling methods by far if rare events with deterministic structure are provided as evidence. Moreover, even when our algorithm converges slowly for a certain problem, this is clearly visible to the user, thereby preventing wrong conclusions. MCMC methods do not offer similar safeguards.