Weighted Model Integration with Orthogonal Transformations

1 Introduction

Motivation Weighted model counting (WMC) is a natural problem to which we can reduce many probabilistic inference tasks, for instance, in Bayesian networks and probabilistic programs [Chavira and Darwiche, 2008; Chistikov et al., 2015; Fierens et al., 2015]. In WMC, we are given a propositional formula, where each satisfying assignment has a weight-e.g., a probability of occurrence-and our goal is to compute the total weight of the formula. Recently, a number of works have leveraged state-of-the-art satisfiability (SAT) solvers as oracles for solving WMC problems [Chakraborty et al., 2014; 2013].

To handle richer domains, researchers from the artificial intelligence, as well as the formal methods, community have also started tackling the weighted model integration (WMI) problem [Belle et al., 2015b; 2015a; 2016; Chistikov et al., 2015; 2014], where formulas are over infinite domains, e.g., real arithmetic. The idea is that inference in, for instance, hybrid Markov networks and rich probabilistic programs, can be reduced to WMI problems. Recent works tackling the WMI

problem have primarily targeted the first-order theory of linear real arithmetic (LRA), which is well supported in state-ofthe-art satisfiability modulo theories (SMT) solvers [Barrett et al., 2009].

Problem setting In this paper, we continue the investigation of WMI algorithms, with an emphasis on a versatile class of techniques for formulas in LRA. Given a formula ϕ in real arithmetic with k free variables, we view ϕ as a region in Rn . We assume the value of each free variable x in ϕ is drawn independently from some probability distribution Dx. Our goal is to compute the probability of satisfaction of ϕ.

A number of recent works [Albarghouthi et al., 2017a; 2016; Sankaranarayanan et al., 2013; Chistikov et al., 2015] have attacked this problem by decomposing ϕ into a set of hyperrectangles (rectangles in Rn ). The idea is that computing the weight of a hyperrectangular region is simple, as all dimensions have constant bounds. Of course, there may be infinitely many hyperrectangles in the decomposition of ϕ. By considering more and more hyperrectangles in ϕ, we approach the true weight of ϕ from below, and converge in the limit. Similarly, we can approach the weight of the negation of ϕ (¬ϕ) from below, thus approaching the weight of ϕ from above. A key feature of this class of techniques is that we can terminate the process at any point and retrieve exact lower and upper bounds on ϕ's weight. This is often desirable in verification of programs with uncertainty, where we are interested in proving some upper/lower bound on the probability of a bad event [Sankaranarayanan et al., 2013; Bouissou et al., 2016; Sampson et al., 2014].

Challenges addressed The foremost challenge in hyperrectangular decomposition is that it is hard to fill a slanted region with hyperrectangles; that is, we might require a large number of hyperrectangles, making con-

vergence incredibly slow. (See illustration on the right where gray rectangles are used to fill the region between the two slanted faces.) To tackle this challenge, we apply a linear transformation (e.g., a rotation) to the given formula ϕ in order to make its faces axis-aligned. In doing so, we must (i) understand how transforming ϕ affects its weight, and (ii) ensure that hyperrectangular decomposition can still be applied to the transformed formula. The problem is that, in general, transforming a formula may change its weight; may produce irrational coefficients (preventing the use of SMT solvers); or may produce stochastic dependence between its free variables (undermining a key premise of hyperrectangular decomposition).

Contributions We make the following key contributions:

• We argue that hyperrectangular decomposition on LRA formulas can benefit from linear transformations, but show that this is only true for formulas with Gaussian free variables. Additionally, we show the desirable properties of orthogonal transformations in this setting. (Section 3)
• We show how to construct orthogonal transformation matrices with rational entries, using Givens rotations and Pythagorean triples. We also provide a method to obtain suitable transformations for formulas. (Section 4)
• We implement and evaluate our approach using the Z3 SMT solver; on representative benchmarks, our approach often dramatically improves performance of weighted model integration tasks. (Section 5)

2 Motivation and Illustration

As a motivating example, suppose we have the independent random variables x1 ∼ N(0, 1) and x2 ∼ N(0, 1). Let ϕ be a first-order formula in the theory of linear real arithmetic (LRA) over variables x1 and x2:

Intuitively, the satisfying assignments of ϕ can be viewed as a region in R2 , as illustrated in Figure 1(a), where ϕ is the region inside the parallelogram. (For simplicity, we illustrate our approach on a convex polyhedron; however, the class of formulas we consider, LRA, admits formulas with disjunctions, i.e., unions of polyhedra.)

Hyperrectangular decomposition Our goal is to compute the probability that ϕ is satisfied. Formally, we want to evaluate the integral ∫ϕ p(x1, x2)dx1dx2, where p(x1, x2) is the joint probability density function (PDF) of x1 and x2.

A number of recent techniques use SMT solvers to decompose the region ϕ into an infinite set of non-overlapping rectangles H1,H2, . . ., such that ∨∞i=1 Hi ≡ ϕ. The idea is that computing the weight of a hyperrectangular region is easy, as all dimensions have constant upper and lower bounds. Intuitively, the weight of ϕ is the sum of weights of all Hi: ∑∞i=1(∫Hi p(x1, x2)dx1dx2). Thus, recent techniques [Sankaranarayanan et al., 2013; Albarghouthi et al., 2017a] iteratively compute the weight of more and more hyperrectangles Hi, converging to the weight of ϕ in the limit.

For the purposes of our work here, we are primarily concerned with the quality of the hyperrectangles that the algorithm considers. Intuitively, the larger the hyperrectangles, the faster the convergence. In the case of our formula ϕ, the largest possible hyperrectangle in terms of weight is shown in Figure 1(b). It is easy to see that this rectangle covers a relatively small subset of ϕ. This is due to the two slanted faces of ϕ, which force us to find hyperrectangles that only touch the faces with two of their corners.

Applying linear transformations Ideally, we would decompose ϕ into a series of hyperrectangles where most of the weight is concentrated in a small number of rectangles, thus accelerating convergence to the true weight of ϕ. To do so, we observe that we can apply a linear transformation T to the region ϕ, resulting in a new region Tϕ that admits a better hyperrectangular decomposition. In our example, we might choose T to be a rotation ϕ that aligns the two slanted faces with an axis. The resulting formula Tϕ is shown in Figure 1(c). The red region is the rectangle with the largest weight that fits in Tϕ. Observe how this rectangle covers the densest region-around the origin.

In this illustrative example, we assumed that we can rotate the formula such that some of its faces are axis-aligned. This is not always possible. It may be that the "best" transformation matrix contains irrational entries, leading to irrational constants in the transformed formula. Meanwhile, SMT solvers expect formulas to have rational constants. Later, we will show how to construct transformation matrices that make formulas better-conditioned for hyperrectangular decomposition, without introducing irrational constants.

3 ProblemDefinition

We now formalize our weighted model integration problem. Formulas We consider formulas in the quantified linear real arithmetic (LRA) fragment of first-order logic. LRA is supported by state-of-the-art SMT solvers like Z3 [De Moura and Bjørner, 2008].

A formula ϕ is constructed using the following grammar:

where ci, c ∈ R. We use the column vector x to denote the variables that are free in ϕ (i.e., not bound by a ∃ or ∀ quantifier). Note that LRA is closed under quantifier elimination: for every formula ϕ, there is an equivalent LRA formula ϕ′ that is quantifier free.

We will use ϕ[x ↦→ x′ ] to denote ϕ with occurrences of variables x replaced by respective elements of x′ .

Weighted model integration Given a formula ϕ with k free variables xT = [x1, . . . , xk], we assume x is a vector of independently distributed random variables.

We can view the formula ϕ as a region in Rk . We define the weight of ϕ as follows: WMI(ϕ) = ∫ϕ p(x)dx1 . . . dxk, where p(x) is the joint probability density function of x.

Hyperrectangular decomposition We describe our algorithm for weighted model integration. Its details are not important for our exposition-only the high-level idea. We refer the reader to past works for more information. [Sankaranarayanan et al., 2013; Albarghouthi et al., 2017a; 2016].

We will use H for an LRA formula defining a hyperrectangular region in Rk ; i.e., H is of the form ∧ki=1 li ⩽ xi ⩽ ui, where li, ui ∈ R. A key point is that WMI(H) is easily computed using the cumulative density functions (CDFs) F1, . . . , Fk of x1, . . . , xk, since the upper and lower bounds of each dimension are constant, and the xi are independent: WMI(H) = ∏i (Fi(ui) - Fi(li)).

A hyperrectangular decomposition of a formula ϕ is an infinite sequence of hyperrectangles H1,H2, . . . such that (i) ∨i Hi ≡ ϕ and (ii) Hi ∧ Hj ≡ false for any i ̸= j.

Recent approaches [Albarghouthi et al., 2017a; Sankaranarayanan et al., 2013] use hyperrectangular decomposition to compute the weight of ϕ, by summing the hyperrectangles' individual weights.

At any point during this summation, the sum is a lower bound on WMI(ϕ), and convergence to the actual weight is guaranteed in the limit. 1 In order to additionally compute an upper bound on WMI(ϕ), one could simply invoke this process on ¬ϕ, since we know that WMI(ϕ)+WMI(¬ϕ) = 1. Permissible distributions Note that hyperrectangular decomposition requires the formula's free variables to be independently distributed. The following theorem, credited to Skitovitch and Darmois, characterizes the class of distributions that remain independent under linear transformation [Skitovitch, 1953; Darmois, 1953].

Theorem 1. Let L1 = ∑i αixi and L2 = ∑i βixi, with xi independent real random variables and αi, βi ∈ R. If L1 and L2 are independent, then all xj with nonzero αj, βj must be Gaussian.

In short: nontrivial linear combinations of independent random variables are independent only if those random variables are Gaussian. It follows that only formulas with Gaussian free variables can be subjected to a linear transformation, and then integrated via hyperrectangular decomposition.

Canonical forms Without further loss of generality, we assume all of our weighted model integration problems are over formulas ϕ where x ∼ Nk(0, I), and I is the identity matrix.

Note that we can transform any formula ϕ with variables x ∼ Nk(µ, Σ) into a canonical form ϕ′ with variables x′ ∼ Nk(0, I), such that WMI(ϕ) = WMI(ϕ′ ).

To perform this transformation, we exploit the following well-known property of multivariate Gaussians: There exists a matrix C ∈ Rk×k such that if z = C -1 (x - µ), then z ∼ Nk(0, I). The following theorem characterizes correctness of this transformation. The idea is that we can replace the variables x of ϕ with the new variables z, creating a new formula ϕ′ whose free variables are z.

1 In the case ϕ is decomposable into finitely many hyperrectangles, then the algorithm can converge in finitely many steps.

Henceforth, we shall assume that all weighted model integration problems are over formulas in canonical form.

Example 1. Consider the simple formula ϕ ≡ x > 0, where x ∼ N(10, 20). We can construct the canonical form ϕ′ ≡ ∃x. z = (x - 10)/20 ∧ x > 0, where z ∼ N (0, 1). We have WMI(ϕ) = WMI(ϕ′ ). ■

Orthogonal transformations While theorem 1 shows that only Gaussian distributions can retain independence under linear transformations, it does not characterize the linear transformations which make that possible. We show that orthogonal transformations-rotations and reflectionspreserve independence in our setting.

Suppose we have a formula ϕ in canonical form. Let p(x) denote the joint PDF of its free variables, and let T be a linear transformation. Then we know the joint density p(x) ∝ exp(- xTx 2 ), and therefore the transformed density p(Tx) ∝ exp(- xT T T Tx 2 ). If the matrix T is orthogonal, then T TT = I, and p(Tx) = p(x); hence, the free variables will remain independent.

Orthogonal transformations also have the convenient property of preserving weight. If T is an orthogonal matrix, then | detT| = 1. It follows from integration by substitution of variables that WMI(ϕ) = WMI(Tϕ).

4 WMI with Orthogonal Transformations

As suggested previously, orthogonal transformations can make formulas better-conditioned for WMI. However, this strategy comes with challenges. First, SMT solvers expect formulas with rational coefficients, but orthogonal matrices will generally have irrational entries. A naïve approximation of irrational entries with rationals would generally produce a non-orthogonal transformation that preserves neither independence nor weight. Second, choosing a transformation that maximizes the efficiency of WMI for a formula entails a difficult optimization problem.

In this section, we present (i) a technique for constructing orthogonal matrices with rational entries and (ii) a technique for choosing a favorable matrix.

4.1 Orthogonal Matrices with Rational Entries

Our method for constructing orthogonal matrices with rational entries relies on Givens rotations. A Givens rotation is an (orthogonal) matrix of the form

where c = √a2 + b2 . This matrix rotates vectors in Rk , but only acts on their ith and j th components. Its entries ac and bc can be interpreted as cos θ and sin θ, respectively, where θ is the angle swept by the rotation.

Any orthogonal matrix Q can be constructed as a product of Givens rotations (with an additional trivial reflection in the case that Q is a reflection); in this sense, the Givens rotations are primitives of the orthogonal transformations.

Notice that the entries of Gij are rational whenever a, b, and c are integers; i.e., whenever (a, b, c) form a Pythagorean triple. The Pythagorean triples (a, b, c) generate a set of unit vectors [ ac , bc ]T ∈ Q2 that is dense on the unit circle in R2 [Shiu, 1983]; see Figure 2(c). In consequence, any angle θ can be approximated to arbitrary precision by arccos ac = arcsin bc , where (a, b, c) is a Pythagorean triple. We use Shiu's method (1983) for obtaining (a, b, c), given θ and a precision requirement. In practice, we limit our precision in terms of the lengths of integers a, b, and c-formulas with large integers impose higher computational cost on the SMT solver during hyperrectangular decomposition.

Since a Pythagorean triple (a, b, c) can be generated to approximate any angle θ ≃ arccos( ac ) = arcsin( bc ) to arbitrary precision, it follows that we can approximate any Givens rotation Gij ∈ Rk×k to arbitrary precision with a rational Givens rotation Gq;ij ∈ Qk×k . Since any orthogonal transformation can be constructed from Givens rotations, we can approximate any orthogonal transformation with a rational one. The next section describes our method for choosing a useful transformation and building its approximation from rational Givens rotations.

4.2 Composing a Suitable Transformation

The following example formula will guide our discussion:

with x1, x2 ∼ N(0, 1) iid (see Figure 2(a)). Note that in its current form, ϕex would be poorly suited for hyperrectangular decomposition: most of its weight is near its boundary (around the origin), and only the corners of sampled rectangles would capture it. It may be that an orthogonal transformation would make ϕex better suited for WMI. In this subsection, we give a method for choosing that transformation.

Aligning one of the faces of the boundary of ϕex with one of the axes would allow the weight near that face to be efficiently captured. However, note that the faces have differing amounts of weight in their proximity; one could say that the face passing through the origin (x1+x2 ⩽ 0) has more weight than the other (x1 + 2x2 ⩽ 1), which passes through regions of lower density.

From these observations, we derive the following methodology for choosing an orthogonal transformation, which we exemplify below for ϕex:

1. On each face of ϕex, numerically integrate the joint PDF of x1, x2; this surface integral will be called the weight of a face. Let the normal vectors of the faces of ϕex be called a1 = [1, 1]T , a2 = [1, 2]T ; numerically, we find their weights to be w1 ≈ 0.37 and w2 ≈ 0.18.
2. Since w1 > w2, we construct an orthogonal matrix that aligns a1 with one of the axes. The result is T = 1 √2 [ 1 -1 1 1 ].
3. Since T contains irrational entries, we cannot use it to transform ϕex. Instead, we approximate T with rational Tq. Our example is in two dimensions, so it suffices for Tq to be a single rational Givens rotation:

1. Transformϕex with Tq to yield

illustrated in Figure 2(b).

In two dimensions it suffices to choose one face of the formula and rotate it into alignment with an axis; this fully defines the rotation. In higher dimensions, orthogonal transformations have more degrees of freedom and it becomes more challenging to fully specify them.

We generalize steps 1-3 in the above procedure to higher dimensions in the following way:

• Construct matrix A, whose columns a1, . . . , an are normal unit vectors for the faces of the formula, ordered by their faces' weights.
• Compute a rational approximation of A's QRfactorization; A = QR where Q is orthogonal and R is upper-triangular. Specifically, we use a variant of QR-factorization by Givens rotations, a well-known method that composes Q from Givens rotations. Our variant is modified to use only rational Givens rotations, and results in a rational approximation Qq of Q.
• Let Tq = Qq.

We give a brief description of QR-factorization, and some rationale for using it in this fashion. Every matrix A has a QR-factorization; i.e. A = QR where the columns q1, . . . , qn of Q are orthonormal and R is upper-triangular. This factorization is unique, and the columns qi have the following property:

(where normalization is implicit on the RHS).

In effect, Q is orthogonal and transforms faces into alignment with axes, prioritized by the faces' weights. Specifically, QTa1 is aligned perfectly with the first axis, QTa2 is aligned approximately with the second axis, and so on. The rational Qq approximates this property, while preserving the formula's rationality.

The computational cost of QR-factorization by Givens rotations is O(n3 ). In our context, this polynomial cost is dwarfed by the worst-case exponential cost incurred by the SMT solver.

5 Implementation and Evaluation

Implementation We implemented our orthogonal transformations technique atop an implementation of the algorithm proposed in [Albarghouthi et al., 2017a], which uses the Z3 SMT solver [De Moura and Bjørner, 2008] to perform hyperrectangular decomposition. Constructing a single hyperrectangle involves making a Z3 query. We use the approach presented in Section 4 to construct an orthogonal transformation, after using Redlog 2 to eliminate quantifiers from the formula.

Evaluation Our primary goal in evaluation is to understand the effect of applying orthogonal transformations on SMT formulas prior to performing weighted model integration. To this end, we considered two classes of benchmarks: (i) First, we generated random benchmarks of quantifier-free formulas of increasing size, where the task is to compute the weight of the formula, where variables are independently and identically distributed. Our primary finding is that transformations can dramatically speed up convergence, particularly as we increase the number of dimensions. (ii) Second, we adapted random-walk benchmarks from the programming languages literature [Chatterjee et al., 2016]. Given a probabilistic program simulating a random walk in Rn , the goal is to compute the probability of landing in a specific region of Rn after k steps. We encoded these probabilistic programs as SMT formulas, using the methods of Chistikov et al. [2015].

Results on synthetic problems For the synthetic benchmarks, we randomly generated formulas using two parameters: (i) the number of free variables (3, 5, or 7), and (ii) the number of conjuncts in the formula (1, 2, or 3). For each combination of number of free variables and conjuncts, we generated 40 benchmarks. The benchmarks were conjunctions of linear inequalities. The coefficients of each inequality were generated as uniformly distributed unit vectors in 3-, 5-, or 7-dimensional space; the formulas can be understood as conjunctions of randomly oriented halfspaces. We let each formula ϕ run for 100 seconds, where we, in parallel, ran WMI(ϕ) and WMI(¬ϕ). This is with the exception of the 7variable/3-conjunct benchmarks, which tended to be too large for the hyperrectangular decomposition method.

Figure 4 shows our results on the synthetic benchmarks. The x-axes represent the sum of WMI(ϕ) + WMI(¬ϕ) after 100 seconds for the case without orthogonal transformations; the y-axes represent results with orthogonal transformations. Recall that the longer our WMI algorithm runs, the closer it arrives to the actual volume. The closer the value is to 1 the better. We observe that for lower dimensions (3), orthogonal transformations do not provide a significant improvement in computed volume; on benchmarks with 3 free variables and 1 conjunct, an orthogonal transformation gave only a 17% increase to the weight on average. However, as we consider larger dimensions, the effect of orthogonal transformations becomes more pronounced; on average, orthogonal transformations gave a 130% increase in weighted volume for the 5-variable/2-conjunct benchmarks and a 1500% increase for the 7-variable/1-conjunct benchmarks. This dependence on dimension results from the corners of hyperrectangles becoming increasingly prominent in higher dimensions.

2 http://www.redlog.eu/

For a closer picture of the difference in performance with and without orthogonal transformations, consider Figure 5. Here, we pick two representative benchmarks and plot the value of 1-(WMI(ϕ)+WMI(¬ϕ)), as computed by the algorithm over the course of 100 seconds (lower value is better). Our results show that WMI with transformations consistently provides better bounds for the duration of execution. Thus, at any point in execution, we can halt weighted model integration and derive better exact bounds in the presence of orthogonal transformations.

Results on random-walk problems For the random walk problems, we generated a set of benchmarks named nd-ks-t, where n is the space Rn in which the random walk takes place; k is the number of steps considered; and t is the type of walk: simp is a program without conditionals, while b and b2 are more complex walks involving branching.

The goal is to compute the probability of an event ϕ encoding the position of the walk after k steps. Using the technique of Chistikov et al., we encoded programs as quantified LRA formulas. The sizes of the generated formulas are shown in the table in Figure 3. The largest formula we considered has 90 atomic predicates of the form ∑i cixi ⩽ c, 15 free variables, and 67 existentially quantified variables.

The bar chart in Figure 3 shows the lower bound on the probability after 100s of execution. We notice that orthogonal transformations yielded consistently better results, with the exception of 2 benchmarks. Our results here display a similar pattern to the synthetic benchmarks: the effects of orthogonal transformations are more pronounced for larger dimensions.

6 RelatedWork

We now compare our work with existing WMI techniques.

Exact weighted model integration Our work is closely related to exact model integration for formulas in SMT theories. To our knowledge, the first such technique is due to [Ma et al., 2009], where the goal was to compute the volume of a formula ϕ in LRA. There, the authors iteratively made calls to an off-the-shelf tool-LattE [De Loera et al., 2012]-for computing the volume of a polytope. This technique decomposes the formula into the set of polytopes (effectively, its DNF form). Recent work by Belle et al. [Belle et al., 2015a; 2016] also utilizes LattE as a backend tool, but generalizes the problem to piecewise-polynomial weight specifications. In their setting, weights defined by Gaussian distributions are approximated using piecewise-polynomial functions.

In comparison with the aforementioned works, our approach is not restricted to volume computation or piecewisepolynomial weights. Hyperrectangular decomposition can integrate arbitrary weight functions, although the orthogonal transformation technique described in this paper must be restricted to a formula's Gaussian variables.

Recent work introduced probabilistic inference modulo theories [Braz et al., 2016]. The instantiations of their DPLLlike algorithm have only been over discrete bounded domains. In contrast, our work targets weighted model integration over the unbounded, dense domain of real arithmetic.

Approximate weighted model integration Weighted model counting is a #P-complete problem. Recently, there has been interest in providing (ϵ, δ)-guarantees in model counting using polynomially many calls to an NP oraclea SAT or SMT solver. A number of works have studied the propositional counting setting, e.g., [Chakraborty et al., 2014; 2013; Ermon et al., 2014; 2013].

In the weighted model integration setting, recent work has shown how to use SMT solvers to approximate volumes of LRA formulas [Belle et al., 2015b; Chistikov et al., 2015]. In comparison to these works, our goal is to provide an exact (non-probabilistic) guarantee. Specifically, we provide exact upper and lower bounds on the desired probability by decomposing an LRA formula into hyperrectangles. Hyperrectangular decomposition is also used by Chistikov et al. in the approximate setting.

7 Conclusion

We addressed the class of WMI algorithms that decompose a LRA formula into a set of disjoint hyperrectangular regions. We argued that the hyperrectangular decomposition approach benefits from transformations on formulas' Gaussian variables, and showed how to to construct orthogonal transformations over rationals, preserving the independence of Gaussian variables without introducing irrational coefficients. We implemented and evaluated our approach, demonstrating the substantial improvements available through orthogonal transformations. As a result, downstream applications such as repair and synthesis [Albarghouthi et al., 2017b] are improved.

Acknowledgements

We thank Jin-Yi Cai for suggesting exploring rotations to improve performance. This material is based upon work supported by the National Science Foundation under Grant numbers 1566015 and 1652140.