SAT-Metropolis: Combining Markov Chain Monte Carlo with SAT/SMT Sampling

Probabilistic inference with Markov Chain Monte Carlo (MCMC) is at the core of statistical analysis with a myriad of applications in domains such as biology, medicine and physics. Probabilistic inference for problems with soft constraints has seen notable progress [14], but inference in the presence of hard constraints remains a difficult task [12]. Yet hard constraints appear in a wide range of real-world applications such as road routing, genetics or computer security; just to mention a few of them.

Probabilistic inference via MCMC is a technique to sample data or compute expected values from distributions whose explicit density functions are unknown, but we can evaluate them up to a multiplicative factor. For instance, in probabilistic Bayesian inference $P(x|y)=P(y|x)P(x)/P(y)$, the term $P(y|x)P(x)$ (the likelihood and the prior) can be evaluated at any point of the state space of the target distribution, but computing the normalizing factor $1/P(y)$ is infeasible in practice. MCMC algorithms tackle this problem by following a random walk on the state space guided by the term $P(y|x)P(x)$ [4].

Unfortunately, MCMC algorithms often fail in the presence of hard constraints that make the state space sparse or even partition it into disjoint parts of positive probability separated by areas of zero probability. Since MCMC algorithms reject zero probability samples, they struggle to explore such distributions completely, even though they are correct in the limit. Classic algorithms like Metropolis-Hastings [19, 11] and Gibbs sampling [9] cannot effectively reach all positive probability states under these conditions [12]. They fail to converge or must be run for a prohibitively long time.

Existing work addressing this problem focuses on the design of sampling algorithms for specific types of probabilistic inference problems [12, 13, 22]. For instance, Hazelton and coauthors propose a method to sample from a $\mathbb{Z}$-polytope defined by inverse linear counting problems – i.e., problems of the form $𝐲=\mathrm{𝐀𝐱}$ where $𝐲$ is a vector of counts, $𝐀$ is a contingency matrix and $𝐱$ the vector of variables of interest [12]. Although efficient, their method is limited to problems of this particular form. Poon and Domingos pioneered using SAT samplers in probabilistic inference. They propose an inference method for problems with propositional constraints that encode relational properties [22]. They create a slice sampler using SAT sampling to perform probabilistic inference over Markov Logic Networks [23]. Unfortunately, this requires encoding probabilistic models as Markov Logic Networks, a model not used commonly in statistics and limited when continuous variables are needed.

In this work, we extend a classic MCMC method with SAT/SMT sampling to support probabilistic inference under mixed soft and hard constraints. In particular, we use SAT/SMT samplers as proposal distributions for the Metropolis algorithm [19]. Our method works on standard probabilistic models widely used in statistics. The theory of discrete integer linear arithmetic with multiplication allows us to handle multivariate polynomial hard constraints. We evaluate it using multiple modern SAT/SMT samplers that provide an excellent basis to “jump” between the set of states with positive probability and avoid probability-zero regions of the state space by design. Finally, we discuss what properties of SAT/SMT samplers help to effectively address probabilistic inference. In summary, our contributions are:

• $\mathrm{■}$

  A MCMC algorithm for probabilistic inference problems under mixed soft and hard constraints, supporting multivariate polynomial hard constraints,

• $\mathrm{■}$

  An implementation of the algorithm using three different SAT/SMT samplers,

• $\mathrm{■}$

  An evaluation of convergence of the inference for different encodings and samplers,

• $\mathrm{■}$

  An evaluation of the scalability of the algorithm in realistic case studies.

Let $T$ be a theory over a set of variables with some concrete domain, logical symbols, predicate symbols, function symbols and a fixed interpretation. Let $\phi$ denote a formula in $T$ and $𝐱$ be a vector with the variables appearing in $\phi$. A satisfying assignment is an assignment of domain values to the variables in $𝐱$ such that $\phi$ evaluates to true. We use $\{𝐱\mid 𝐱\vDash \phi \}$ to denote the set of the satisfying assignments of $\phi$.

A Satisfiability Modulo Theory (SMT) sampler is an algorithm that generates samples from the set of solutions of a formula $\phi$, denoted ${𝐱}_i\stackrel{\$}{\leftarrow }𝒟(\{𝐱\mid 𝐱\vDash \phi \})$, according to a probability distribution $𝒟$. When $\phi$ is a Boolean formula, we call it a SAT sampler. A SAT/SMT sampler is uniform if the probability of obtaining any satisfiable assignment is $1/|\{𝐱\mid 𝐱\vDash \phi \}|$, denoted $𝒰(\{𝐱\mid 𝐱\vDash \phi \})$. A SAT/SMT sampler is almost uniform if the probability of sampling any satisfiable assignments is at least $1-ϵ/|\{𝐱\mid 𝐱\vDash \phi \}|$ and at most $1+ϵ/|\{𝐱\mid 𝐱\vDash \phi \}|$ for some $ϵ>0$.

In this paper, we use the theory $T_{\mathrm{𝑀𝐼𝐴}}$ of integer linear arithmetic extended with multiplication to handle multivariate polynomial constraints. A $T_{\mathrm{𝑀𝐼𝐴}}$ formula includes logical connectives $\{\neg ,\wedge \}$ (and derived), predicates $\{>,=\}$ (and derived) and arithmetic operations $\{+,-,\ast \}$, but no division. The interpretation of the symbols is standard. Variables have finite integer domains. We use bit-blasting to encode formulae in Boolean logics for SAT-samplers.

The Metropolis algorithm samples from a target distribution $P(𝐲)=P^{\ast }(𝐲)/Z$ by using evaluations of $P^{\ast }(𝐲)$, a function different by an unknown constant factor $1/Z$. The vector $𝐲$ denotes the list of parameters for the target distribution. When applied to Bayesian inference, the target distribution is the posterior defined by the Bayes rule $P(𝐱\mid D)=P(D\mid 𝐱)P(𝐱)/Z$. The term $P^{\ast }(𝐲)$ is instantiated as $P(D\mid 𝐱)P(𝐱)$, i.e., the likelihood times the prior. When no ambiguity arises, we write $P(\cdot )$ to denote a probability measure on a measurable space that is compatible with its parameters.

The likelihood function $P(D\mid 𝐱)$ is used to model soft and hard constraints in the underlying probabilistic model. Let $D=\{c_1,\mathrm{\dots },c_n\}$ be a set of events that model constraint parameters, e.g., observations from real-world data or logical formulae. The set $D$ is commonly known as a set of observations. The likelihood is defined as $P(D\mid 𝐱)={\prod }_{i}P(c_i\mid 𝐱)$. A hard constraint is a constraint with 0–1 likelihood. For example, for a model with parameter $x$, we can impose $x>0$, i.e., $c_i=\{x|x>0\}$, by defining $P(x>0\mid 𝐱)={\mathrm{𝟏}}_{x>0}(x)$ where ${\mathrm{𝟏}}_{x>0}(x)$ is the indicator function; its value is $1$ if $x>0$ and $0$ otherwise. Note that the posterior probability of values of $x$ that do not satisfy the constraint is zero, . Soft constraints, on the other hand, allow for specifying noisy properties in the model. For example, consider that value $c_i=v$ has been collected using a recording device with a known measurement error of ${\sigma }_v\in {ℝ}_{>0}$. It is possible to model the measurement error using a soft constraint $P(v\mid x)=𝒩(v;\mu =x,\sigma ={\sigma }_v)$. In this example, we use a Gaussian distribution ($𝒩$) to model that values of $x$ close to $v$ have high likelihood, and other values, though less likely, are also possible; due to the measurement error of the recording device.

To sample values from the posterior, the Metropolis algorithm takes a random walk over the parameter space starting in a random point $𝐱$ with a positive posterior probability, $P(𝐱\mid D)>0$. Given a state ${𝐱}_i$, a proposal distribution, $q({𝐱}^{\ast }|𝐱)$, is used to sample a candidate next state in the random walk. The Metropolis algorithm requires $q(\cdot )$ to be symmetric, i.e., $q({𝐱}^{\ast }|𝐱)=q(𝐱|{𝐱}^{\ast })$. An extension of the algorithm, Metropolis-Hastings, relaxes this constraint [11, 16]. The probabilistic inference presented in Section 5 does not require asymmetric proposal distributions. Given a state ${𝐱}_i$ and proposal ${𝐱}^{\ast }$, the next state ${𝐱}_{i+1}$ is selected according to the following acceptance ratio:

In particular, ${𝐱}_{i+1}$ is selected as

That is, the proposed state is accepted with probability $\min (1,r(𝐱,{𝐱}^{\ast }))$. Otherwise, the proposal is rejected and the random walk stays in the current state.

If the proposal distribution is positive, i.e., $q({𝐱}^{\ast }\mid 𝐱)>0$ for all $𝐱,{𝐱}^{\ast }$, the probability distribution induced by a collection of $n$ samples ${\{{𝐱}_i\}}_{i=1..n}$ is known to converge to $P(𝐱|D)$ as $n\to \mathrm{\infty }$ [11, 16]. The theoretical guarantees aside, the choice of proposal distribution drastically affects the speed of convergence in practice. Much of the work on MCMC addresses designing efficient (in terms of convergence speed) proposal distributions such as Gibbs sampling or Hamiltonian Monte Carlo (HMC) [4], predominantly for continuous variables. Despite these advances, it remains a challenge to sample from posterior distributions with hard discrete constraints.

The popular MCMC algorithms, such as Metropolis, Metropolis-Hastings, Gibbs sampling and HMC, are not suitable for sampling from distributions with hard constraints. The problem is that the proposal distribution is unaware of the hard constraints and, consequently, proposes states violating them. Such proposals are promptly rejected – for any ${𝐱}^{\ast }$ violating hard constraints, we have $P({𝐱}^{\ast }\mid D)=0$, the acceptance ratio $r(𝐱,{𝐱}^{\ast })=0/P(D|𝐱)P(𝐱)=0$, and the probability of accepting ${𝐱}^{\ast }$ is $\min (1,0)=0$. Thus, the random walk gets stuck in its current state ($𝐱$).

This problem gets much worse for sparse constraints when the state space increases. It is intractable to randomly sample a satisfying solution to a hard constraint. As we show in Section 6, problems do not have to get very complex for the sampler to get stuck completely.

The key idea to address the above problems is to use SAT/SMT sampling as the proposal distribution. We present, sat-metropolis, our modification of the Metropolis sampler as Algorithm 1 below. The algorithm takes four arguments: (i) a possibly not normalized probability density $P$ to approximate, (ii) a set of hard constraints $\mathrm{\Phi }$, (iii) an initial state for the random walk, and (iv) the number of samples $n$ to generate. As common for the sampling algorithms, the initial state must have positive probability $P({𝐱}_{\mathrm{init}})>0$, which in this case means that it must also satisfy the hard constraints. The output is a set of samples ${\{x_i\}}_{i=1}^n$. When used for Bayesian inference, $P$ is instantiated as likelihood times prior, i.e., the unnormalized posterior distribution (cf. Section 3).

We begin in the provided initial state (line 2). The main novelty is found in line 4. The proposal distribution, $𝒰(\{𝐱|𝐱\vDash \mathrm{\Phi }\})$ is a uniform distribution over all states that satisfy the hard constraints. A uniform proposal distribution over non-zero probability states reduces autocorrelation among samples, which helps reducing bias and improving speed of convergence [17]. Furthermore, it is compatible with the inner workings of existing SAT/SMT samplers, which work by producing batches of uniform samples. However, the proposal distribution in Metropolis can be non-uniform and conditional on a state. This type of proposal distributions can also achieve low autocorrelation, if they are properly configured. Since SAT/SMT sampling from conditional and non-uniform is not supported by existing tools, we do not consider it in this version of sat-metropolis. The remaining steps in Algorithm 1 proceed as usual in Metropolis (cf. Section 3). Lines 5 and 6 compute the acceptance probability. Lines 7 to 11 select the new proposal ${𝐱}^{\ast }$ with probability $\alpha$. Otherwise we remain in the current state. The procedure is repeated $n$ times.

The proposal distribution in sat-metropolis is symmetric and positive; as it is independent of the current state and it assigns non-zero probability to any pair of states. Therefore, it trivially follows that sat-metropolis converges in the limit (cf. Section 3). However, the practical use of this algorithm will depend on how efficiently and uniformly we can sample from the proposal distribution.

We implemented sat-metropolis in Python, using three engines for sampling solutions to hard constraints: SPUR [2], CMSGen [10], and MegaSampler [21]. We selected these engines as they vary in: the uniformity of generated samples, speed of sampling and amount of pre-processing. These features are the main factors that affect the speed of convergence and scalability of our algorithm (cf. section 6). Furthermore, these engines are at the forefront of the research in SAT/SMT sampling. The implementation is publicly available at https://github.com/itu-square/sat-metropolis.

Figure 1 gives a high level overview of the implementation. The implementation takes the target distribution as a Python function from the domain of the parameter space ($𝐱\in {𝒳}^m$) to reals $ℝ$. This distribution specifies the soft constraints. To specify hard constraints, the user should provide an SMT problem in abstract syntax constructed using the Z3 API [6].

The implementation starts by requesting $n+1$ proposal samples from an SMT/SAT sampling engine. For simplicity, we use the first sample as the initial state ${𝐱}_{\mathrm{init}}$. This assumes that the target distribution has positive density for all states satisfying the hard constraints, which is the case for the problems we consider. Then, it proceeds as shown in Algorithm 1 with the only difference that the sampling step on line 4 retrieves a sample from the received proposals. This is equivalent to the behavior specified by Algorithm 1, as proposals are sampled independently of the current state. Depending on the type of engine used for generating samples, there are requirements on the input Z3 problem and pre-processing steps. We discuss them below. However all the engines used allow us to handle problems with polynomial hard constraints that was not possible for previous methods for probabilistic inference.

We evaluate sat-metropolis by answering the following research questions.

1. RQ1

   How does the SAT/SMT sampling method affect the convergence of sat-metropolis?

2. RQ2

   How big problems can be handled by sat-metropolis with the contemporary SAT/SMT samplers?

Table 1 lists the probabilistic inference problems used in the evaluation. To answer RQ1, we use two database reconstruction problems, DB CACM and NZ Stats. DB CACM was previously introduced as an example of database reconstruction using SMT solvers [8]. NZ Stats uses data published by Statistics New Zealand. We use the data and the information on the anonymization algorithms used in these cases, to reason about beliefs for the possible reconstructions. These problems are simple enough to solve analytically, and we use the analytical solution as the ground truth for evaluation of convergence. Note that evaluating the convergence of the estimation on large problems without ground truth is not possible.

To answer RQ2, we use the last two problems in the table (Books and Roads). Books and Roads are examples of problems with realistic sizes regarding the number of variables and constraints and domain sizes. They had been used to evaluate specialized methods to solve probabilistic inference under linear constraints [12] (Section 7). Determining how sat-metropolis performs on them is a good indicator on whether modern SAT/SMT-sampling can handle realistic tasks. We consider multiple versions of these two problems by systematically adjusting their size; with the objective of finding the limits of scalability.

In general, we consider problems that range from 10 to 72 variables and 2 to 141 constraints. The formulae include polynomial and linear constraints (cf. Table 1). The structure of the formulae is problem dependent and we describe it the sections below. The SAT version of these problems ranges from 148 to 10252 variables and from 499 to 47299 clauses. Clauses contain from 2 to 4 literals and the clause-to-var ratio ranges from 3.06 to 5.08. Table 4 in Appendix A contains the complete list of problems that we consider in this paper.

The experiments were run on Ubuntu 24.04 with an Intel Core i7-1355U CPU and 16GB RAM. The experiment package is available at https://github.com/itu-square/sat-metropolis.

Recent work on MCMC based probabilistic inference with hard constraints targets discrete inverse linear problems [12, 18, 7, 3]. These are problems of the form $𝐲=\mathrm{𝐀𝐱}$ where $𝐲$ is an observed counts vector and $𝐀$ is a contingency matrix that defines the constraints on the variables in the vector $𝐱$ w.r.t. the observed counts. A major advance in tackling this type of problems was the use of Markov (sub-)basis (see [3] for an overview). Let ${ℱ}_{𝐲}=\{𝐱|𝐲=\mathrm{𝐀𝐱}\}$ be the set of count vectors satisfying the hard constraints – this set is known as the $𝐲$-fiber. Markov basis methods produce a Markov chain such that all points within the $𝐲$-fiber are reachable from any state – recall that this is a requirement for convergence of MCMC. However, computing the Markov basis for a given problem can be computationally expensive. This has motivated works suppressing the need for computing the full Markov basis. Dobra [7] designed a method to generate sampling directions dynamically while ensuring that the union of all directions is within the $𝐲$-fiber. More recently, the work by Hazelton et al. [12] introduces a method based on dynamic lattice basis which allows for sampling from the $𝐲$-fiber without computing the full Markov basis, which notably improves sampling performance. Our method does not achieve the same level of performance as these specialized methods. However, sat-metropolis can tackle a larger class of problems. Our method can handle problems with polynomial hard constraints whereas discrete inverse linear problems only feature linear constraints. Nevertheless, despite being able to tackle a more general class of problems, our evaluation showed that sat-metropolis with CMSGen can solve the book ratings problem of Hazelton.

We are not the first to propose the use of SAT samplers for probabilistic inference. Poon and Domingos [22] define a slice sampler to tackle probabilistic inference problems specified as Markov Logic Networks [23]. Markov Logic Networks are useful to express statistical problems with relational constraints, but they require uncommon encodings for other problems such as linear or polynomial constraints. Instead of designing a sampler for a specific type of probabilistic models, namely Markov Logic Networks, we adapt Metropolis, which a standard and widely used algorithm in statistics. sat-metropolis takes probabilistic inference problems expressed as an unnormalized density function (typically the numerator in Bayes rule: prior $\times$ likelihood), which is common place in the statistics literature. The slice sampler uses SampleSAT [24] to sample solutions for the relational constraints (expressed as propositional formulae). SampleSAT does not have uniformity guarantees and works on SAT problems. As a consequence, it cannot be used to analyze the effect of uniformity or bit-blasting, as we did in our work; with the backends SPUR (perfect uniformity) or MegaSampler (which works directly on SMT problems and can be compared to SPUR/CMSGen that require bit-blasting). Furthermore, the most scalable backend of sat-metropolis (cf. Section 6.2), CMSGen, has been recently shown to be one of the most scalable SAT sampling tools [10]. Thus, our work provides a more detailed and updated perspective on the use of SAT/SMT sampling methods for probabilistic inference than Poon and Domingos [22].

We have demonstrated that SAT/SMT sampling technology can effectively be used to tackle probabilistic inference problems with multivariate polynomial hard constraints. SAT/SMT sampling presents itself as a promising technology for tackling the problem of probabilistic inference via MCMC, where hard and soft constraints are mixed. To this end, we have adapted the Metropolis algorithm to use SAT/SMT samplers as proposal distributions. This prevents proposing zero-probability states, which lead the standard Metropolis algorithm to fail to explore the target distribution well. We have implemented sat-metropolis with three SAT/SMT samplers; SPUR, CMSGen and MegaSampler, and have evaluated convergence and scalability for these three variants. The results indicate that uniformity impacts convergence of sat-metropolis positively; as SPUR is the backend that exhibits the best convergence. At the same time, CMSGen is the most scalable sampler in this application. Even though sat-metropolis is generic with respect to the constraint class supported, it can handle problems used to benchmark specialized methods.