Probabilistic software analysis methods extend classic static analysis techniques to consider the effects of probabilistic uncertainty, whether explicitly embedded within the code – as in probabilistic programs – or externalized in a probabilistic input distribution. Analogously to their classic counterparts, these analyses aim at inferring the probability of a target event to occur during execution, e.g., reaching a program state or triggering an exception. Besides helping understand the behaviour of software programs that interact with the uncertain, real-world inputs, they also have a profound impact as verification tools for safety-critical systems. Providing efficient, scalable algorithms that are general to various applications remain one of the main challenges that prevent the wide adoption of these analysis tools. In our recent paper published in the ACM Joint European Software Engineering Conference and Symposium on the Foundations of Software Engineering (ESEC/FSE ‘21), we introduce a new algorithm that combines advances in symbolic execution and probabilistic inference to enable these analysis tools to be applied to more complex and general problems.

Background

Probabilistic Symbolic Execution

Probabilistic symbolic execution (PSE) is a static analysis technique aiming at quantifying the probability of a target event occurring during execution. It uses a symbolic execution engine to extract conditions on the values of inputs or specified random variables that lead to the occurrence of the target event. It then computes the probability of such constraints being satisfied given a probability distribution over the inputs or specified random variables. These constraints are called path conditions because they uniquely identify the execution path induced by an input satisfying them.

To illustrate with an example, consider the following snippet for a probabilistic program modelling the behaviour of a safety controller for a flying vehicle.

// Probabilistic profile p(x)
altitude ::= Gaussian (8000, 100);
obstacle_x, obstacle_y ::= Gaussian (
    [ -2 , -2],
    [[0.2 , 0.1] , [0.1 , 0.2]]) ;
// Controller program
if (altitude <= 9000) { ...
    if (Math.pow(obstacle_x, 2) + Math.pow (obstacle_y, 2) <= 1) {
        callSupervisor ();
        ...
    }
} else { callSupervisor (); }

The purpose of the program is to detect environmental conditions – excessive altitude or collision distance of an obstacle – that may compromise the crew’s safety and call for a supervisor’s intervention. The purpose of the analysis is to estimate the probability of invoking callSupervisor at any point in the code. Safety-critical applications may require this probability to be very small (e.g., $< 10^{-7}$) and to be estimated with high accuracy. The symbolic execution of the snippet, where random variables are marked as symbolic, would return the following two path conditions (PCs), corresponding to the possible invocations of callSupervisor:

• $PC_0$: $\texttt{altitude} > 9000$; and
• $PC_1$: $\texttt{altitude} \leq 9000 \land \texttt{obstacle_x}^2 + \texttt{obstacle_y}^2 \leq 1$.

The probability of satisfying a path condition $PC$ can be computed based as $$ \begin{aligned} p_{PC} &:= \Pr(\mathbf{x} \models PC) = \int_{\mathbf{x}} \mathbb{1}{PC}(\mathbf{x}) p(\mathbf{x}) \, d\mathbf{x} \\ &\approx \frac{1}{N} \sum{i=1}^N \mathbb{1}_{PC}(\mathbf{x}^{(i)}) =: \hat{p}^{DMC}, \; \text{where} \, \mathbf{x}^{(i)} \sim p(\mathbf{x}) \end{aligned} $$

where $\mathbb{1}_{PC}(\mathbf{x})$ denotes the indicator function. For clarity, we will use $\bar{p}(\mathbf{x})$ to denote the truncated distribution satisfying the constraints, i.e., $\bar{p} := \mathbb{1} (\mathbf{x}) p(\mathbf{x})$.

Monte Carlo Methods

Because analytical solutions to the integral are in general intractable or infeasible, Monte Carlo methods are used to approximate $p_{PC}$. A naive approach to this problem would be Direct Monte Carlo. In DMC, we generate independent samples $\mathbf{x}^{(i)}$ from the prior $p(\mathbf{x})$ and obtain an estimate of $p_{PC}$ by counting the number of samples that satisfy the constraints. Unfortunately, this method is very inefficient for rare-event problems. When the target event has a probability of less than $10^{-6}$, obtaining a single instance that satisfies the constraints require millions of samples. It is impractical to obtain accurate estimates of the event probability. Symbolic execution methods, on the other hand, can find feasible solutions without brute force sampling. However, symbolic methods alone are unable to provide probabilistic quantification. In our work, we are interested in designing algorithms that can combine the best of both worlds.

qCoral

Borges et al.¹ proposed qCoral, a compositional Monte Carlo method to estimate the probability of satisfying a path condition over non-linear numerical constraints with arbitrary small estimation variance. qCoral combines insights from program analysis, interval constraint propagation, and stratified sampling to mitigate the complexity of the integration problem and reduce the variance of its estimates. To further reduce the variance of the probability estimates of each independent constraint, qCoral uses interval constraint propagation and branch-and-bound methods to find a disjoint union of n-dimensional boxes that reliably enclose all the solutions of a constraint – where $n$ is the number of variables in the constraint. Regions of the input domain outside the boxes are guaranteed to contain no solutions of the constraint ($\mathbb{1}{\cdot}(\cdot)=0$). A box is classified as either an inner box, which contains only solutions, or an outer box, which may contain both solutions and non-solutions.

Stratified sampling with interval constraint propagation can lead to a significant variance reduction in the aggregated estimate by reducing the uncertainty only to the regions of the domain enclosed within the outer boxes, potentially avoiding sampling from large regions of the domain that can be analytically determined as including only or no solutions. However, qCoral can only produce scalable and accurate estimates for the satisfaction probability for constraints that

• have low dimensionality or can be reduced to low-dimensional subproblems via constraints slicing,
• are amenable to scalable and effective interval constraint propagation, and
• whose input distribution have CDFs in analytical form and allows efficient sampling from their truncated distributions.

SYMPAIS

In our paper, we take an alternative approach based on importance sampling² (IS) to reduce the variance in our estimates of rare event probability. IS has been successfully applied to many applications such as reliability analysis³. An advantage of IS over the stratified sampling approach as used in qCoral is that we are not restricted to only input distributions which tractable CDFs. An importance-sampling based approach can be applied to any problem with tractable PDFs, including distributions with correlated dimensions. With a properly chosen $q(\mathbf{x})$, IS also offers better variance reduction compared to stratified sampling.

While any distribution with $q(\mathbf{x}) > 0$ over the entire domain guarantees the estimate will eventually converge to the correct value, an optimal choice of $q(\mathbf{x})$ determines the convergence rate of the process and its practical efficiency. In our context of estimating the probability of satisfying path conditions PC, the optimal proposal distribution $q^(\mathbf{x})$ is exactly the truncated, normalized distribution $p(\mathbf{x})$ satisfying PC, $$ q^(\mathbf{x}) = \frac{1}{p_{PC}}p(\mathbf{x})\mathbb{1}{PC}(\mathbf{x}). $$ In general, it is infeasible to sample from $q^*(\mathbf{x})$ as it requires the calculation of $p{PC}$, which is exactly our target. For practical applications, specialized proposals are proposed to achieve high sample efficiency in domain-specific applications. In our work, we are interested in designing an automated procedure that iteratively finds a proposal distribution that is tailored to the given probabilistic program.

To do this, we adapt the PI-MAIS algorithm⁴ and tailor it to the problem of probabilistic program analysis. PI-MAIS is a probabilistic inference algorithm that runs a Markov Chain Monte Carlo (MCMC) algorithm (e.g. with the Random-Walk Rosenbluth-Metropolis-Hastings or Hamiltonian Monte Carlo) to sample from an intractable posterior distribution. In this case, the intractable posterior is our intractable $\bar{p}(\mathbf{x})$. The samples from the Markov chains are used to construct an adaptive proposal distribution that can be used for the program analysis above. Concretely, given samples $\mathbf{\mu}i$ from the $N$ parallel chains at the $t$‘th iteration, we use a Mixture-of-Gaussians (MoG) $$ \begin{equation*} q(\cdot) = \frac{1}{N} \sum{i=1}^{N} \mathcal{N}(\cdot, \mathbf{\mu}_i, \Sigma), \end{equation*} $$ as the importance proposal. The proposal distribution is multi-modal and, with enough number of chains, can be used to approximate arbitrary $q^*(\mathbf{x})$. The figure below gives a graphical illustration of SYMPAIS.

SYMPAIS further utilizes information from the symbolic execution of the program to further improve sample efficiency.

First, the MCMC algorithm requires feasible solutions to initialize the chains. This is typically done in the Bayesian inference literature by sampling from the prior $p(\mathbf{x})$. However, generating diverse initial solutions that cover the support of $\bar{p}(\mathbf{x})$ by sampling from prior is prohibitively expensive. In SYMPAIS, we address this by combing an Interval Constraint Propagator (ICP) to generate diverse initial solutions. The diverse initial solutions are further re-sampled by an importance-reweighting procedure that further reduces the warm-up iterations required for the MCMC samples to converge to the typical set. In addition, we further use the contracted interval bounds from ICP to truncate the modes in the MoG proposal distribution.

These optimizations significantly improve SYMPAIS’s sample efficiency, making it a practical algorithm for high-dimensional applications without sacrificing convergence guarantees offered from the original PI-MAIS paper.

Evaluation

We provide an open-source implementation of SYMPAIS. The open-source code includes a general implementation of PI-MAIS as well as code specializing PI-MAIS to perform probabilistic program analysis as done in SYMPAIS. We also provide notebooks for you to interactively explore SYMPAIS.

The reference Python implementation uses Google’s JAX library. JAX provides excellent performance that allows us to implement a parallel-chain MCMC sampler that can take advantage of the SPMD parallelism that enables us to run hundreds of Markov chains easily with a laptop CPU.

We evaluated SYMPAIS on a range of benchmark tasks, including artificial examples such as high-dimensional spheres; tasks processing non-convex constraints such as torus. We also reported performance on problems from the volComp benchmark as well as the performance on analyzing constraints from quantifying the adversarial robustness of ACAS Xu neural networks. Please refer to our paper for a complete discussion of the result and evaluation.

Sphere

The RAE results are illustrated in the figure above. As expected: DMC achieves the worst performance throughout all tests. For low-dimensional problems ($d \leq 4$), qCoral is the most efficient, while its performance deteriorates significantly when the $d$ increases and RealPaver fails to prune out large portions of the domain that contain no solutions. SYMPAIS’s performance is comparable to qCoral in low dimensions, but up to one order of magnitude more accurate when the dimensionality grows ($d \geq 8$).

The figure on the right shows the convergence rate of RAE for different methods over sample size for $d=8$. SYMPAIS achieves the final RAE of DMC with $<5%$ of the sampling budget and the final RAE of qCoral with $<10%$. The improvement in sample efficiency becomes more significant for $d > 8$.

Torus

Summary

We introduce SYMbolic Parallel Adaptive Importance Sampling (SYMPAIS), a general algorithm for probabilistic program analysis that lifts some of the assumptions of prior work and is more scalable for higher-dimensional and non-linear problems. We believe that our work provides an interesting combination of symbolic execution with statistical inference for practical applications. While statistical methods provide better scalability compared to symbolic methods, incorporating information from a symbolic execution method can further improve sample efficiency. We look forward to exploring more interesting hybrid approaches that can combine the best of both worlds.

References