Automated Safety Verification of Posterior Distributions of Probabilistic Programs

1 Introduction

Uncertainty is ubiquitous in real-world systems, and probabilistic behaviours are among the most challenging sources of such uncertainty. To model systems or behaviours that involve probability, probabilistic programs have proven effective for decades [van de Meent et al., 2018; Barthe et al., 2020; Gordon et al., 2014]. As surveyed in [van de Meent et al., 2018], probabilistic programming has enabled applications such as constrained procedural graphics [Ritchie et al., 2015] and recursive multi-agent reasoning [Goodman and Stuhlm¨uller, 2014]. A variety of probabilistic programming languages have been developed to address different real-world needs, including Church [Goodman et al., 2008], Anglican [Wood et al., 2014; Tolpin et al., 2016], and Pyro [Bingham et al., 2019].

Automated verification of probabilistic programs has become a central topic in the formal verification of probabilistic systems. In particular, numerous techniques have been proposed for verifying almost-sure termination and quantitative bounds on expected run-time or outcomes (e.g. [Batz et al., 2023; Ngo et al., 2018; Bao et al., 2025; Chakarov and Sankaranarayanan, 2013; Chatterjee et al., 2020; Li et al., 2025]). It is worth noting that verifying almost-sure termination is, in general, undecidable [Kaminski and Katoen, 2015], highlighting the inherent difficulties and challenges in the verification of probabilistic programs.

In this paper, we address a verification problem that has received comparatively little attention in the literature: establishing a safety property over posterior distributions of probabilistic programs. We illustrate a motivating example below.

Motivating Example. We consider two robots with identical architectures that independently perform autonomous exploration. Their behaviours are modelled by the following probabilistic program. Let C denote the following program.

b := true, t := 0 while (b) { }

{ D; b := false } [1/100] {t += 1} where the program D is given by { {teamA := true} [1/10] {skip} } □ {

{teamB := true} [1/10] {skip} }

Each robot continues exploring until it encounters an issue; for simplicity, we assume that such an issue occurs with probability 1/100 in each iteration. Once a robot detects an issue and terminates its exploration, it may require manual support from either Team A or Team B, with probability 1/10, depending on its assigned task. This uncertainty is modelled via nondeterminism, resolved in an adversarial (demonic) manner. Crucially, the two robots share the task but operate independently, regardless of what the task is.

Now suppose that the two robots indeed share the same task and operate independently. Although both robots may request support from the same team, the team cannot serve both simultaneously due to the cost of assistance. Therefore, we would like to ensure that the probability that both robots request help from the same team remains reasonably small. This situation can be captured by the terminating probabilities P σ teamA and P σ teamB, where σ denotes a scheduler that decides the task assignment by resolving the nondeterminism, and P σ teamA (resp. P σ teamB) is the probability that a robot terminates by requesting support from Team A (resp. Team B) under the scheduler σ. Given a threshold θ > 0, the verification task asks whether

Since the program models the behaviour of a single robot, we take the square of each probability to account for two independent robots.

In this example, by a simple calculation we obtain

from the initial state in which both teamA and teamB are assigned false.

Challenge. While the above example can be verified by hand due to its simplicity, it raises several questions when we generalize the underlying safety problem. In particular, given a program and a safe region R ⊆ [0, 1]n , we ask whether, regardless of how nondeterminism is resolved by a scheduler σ, the resulting terminating probabilities P σ 1 , . . . , P σ n (from an initial state) to designated terminal states, which we collectively call the posterior distribution, always lie within R, that is,

First, formulating this problem in a suitable way while keeping it tractable is non-trivial. Allowing an arbitrary R would render the problem unsolvable in practice, so we must identify a reasonable class for R to preserve tractability. Moreover, a na¨ıve approach would require enumerating all posterior distributions induced by schedulers-an infeasible strategy, as the scheduler space is typically infinite. To the best of our knowledge, no existing work provides an automated verification method for this problem in probabilistic programs, including the motivating example above.

Our Approach. In this paper, we first formulate the verification problem under certain assumptions on the safe region R. Specifically, we assume that R is a downward-closed convex set of the form

where θ ≥ 0 is a threshold and f : [0, 1]n → R≥0 is a polynomial with rational coefficients. The motivating example corresponds to the special case f(p1, p2) := p21 + p22 .

We then propose a novel automated verification method for this problem. Specifically, we develop a counterexampleguided abstraction refinement (CEGAR) framework [Ball and Rajamani, 2002; Clarke et al., 2003] to verify it effectively.

Our approach maintains an abstraction U, represented as a convex polytope obtained as a finite intersection of halfspaces, such that

To resolve the safety problem, it suffices to show that U ⊆ R. We prove that, under our choice of U, checking U ⊆ R is feasible due to topological properties of U and R.

A CEGAR loop requires refinement via counterexamples; in our setting, refinement is triggered when U ̸⊆ R. Since U is a convex polytope, we show that counterexamples can be extracted from its vertices. We then refine U to a strictly smaller polytope U ′ ⊊ U by adding a new half-space induced by such a counterexample. This half-space is derived by synthesizing a suitable invariant, which can be obtained symbolically by existing solvers (e.g. [Chatterjee et al., 2025]).

We implement our CEGAR loop using well-maintained existing solvers, including Z3 [de Moura and Bjørner, 2008], PolyQEnt [Chatterjee et al., 2025], and Apron [Jeannet and Min´e, 2009]. We evaluate the effectiveness of our approach through preliminary experiments.

Structure. In §2, we present preliminaries on probabilistic programs and formally define our verification problem. In §3, we illustrate the key ideas and observations underlying our approach. There, adequacy [Batz et al., 2024] with respect to predicate-transformer semantics [McIver and Morgan, 2005] and the Fenchel-Moreau duality are recalled as crucial ingredients. In §4, we describe our algorithm in detail and establish desirable properties such as progress and soundness. In §5, we report preliminary experimental results obtained based on our implementation. Related work is discussed in §6, and we conclude in §7.

Notation. With a slight abuse of notation, for vectors p, q ∈ [0, 1]n we write p · q for their inner product. We write [m] for the set {1, . . . ,m} for a natural number m.

2 Preliminaries and Problem Setting

In this section, we introduce our target language, and formalize the verification problem of interest.

2.1 Target Language

Our target language is the probabilistic guarded command language (pGCL) [McIver and Morgan, 2005], a simple imperative language. Formally, we fix a countable set of variables Vars := {x, y, z, . . . }, and define the countably infinite set of states by St := {s: Vars → N | the support of s is finite}. Note that restricting variables to take non-negative values is a standard assumption (e.g. [Batz et al., 2023]) and it does not reduce expressiveness [Batz et al., 2021]. A probabilistic program C is then inductively given by

where 1) e : St → Q≥0 is a linear arithmetic expression (such as x and x + r with r ∈ R≥0); 2) ϕ ⊆ St is a predicate that is defined by inequalities, negation, and conjunction; and 3) p is a probability p ∈ Q≥0 ∩ [0, 1]. Here, {C}[p]{C} and {C}□{C} denote probabilistic and nondeterministic branching, respectively. See, e.g., [Batz et al., 2023; Batz et al., 2024] for the full definition of pGCL. As a syntax sugar, we write p · C1 + (1 - p) · C2 for {C1}[p]{C2} and □ {C1} {C2} for {C1}□{C2}.

2.2 Operational Semantics

For the semantics of pGCL, we employ the operational semantics based on Markov decision processes (MDPs) [Puterman, 1994], following the formulation of [Gretz et al., 2014] in the style of [Plotkin, 2004]. The key idea is that program computations (or executions) are modeled as probabilistic transitions in an MDP, while nondeterministic choices are represented as nondeterministic choices of actions. Program termination corresponds to reaching designated target states in the MDP. We also refer to [Batz et al., 2024] for a selfcontained reference on this operational semantics.

We recall standard preliminaries on MDPs. An MDP M is a tuple (Q, T, A, P ) consisting of a countable set of positions Q, a subset of terminal positions T ⊆ Q, a finite set of actions A, and a transition probability function P : Q × A × Q →

[0, 1] satisfying ∑q′∈Q P(q, a, q′ ) = 1 for all q ∈ Q and a ∈ A. We further assume that terminal positions are sinks, i.e. P(t, a, t) = 1 for any t ∈ T and a ∈ A.

Nondeterminism in pGCL is modeled by action choices, which are resolved by schedulers. Formally, a scheduler is a function σ: Q+ → A. Given a scheduler σ and a path q1 · · · qm of positions, the path probability is defined as

Finally, given an initial position q ∈ Q and a scheduler σ, the terminating probability on a target position t ∈ T is defined by

where the summation ranges over paths q · · · qm · t such that none of q, . . . , qm are terminal.

Following [Gretz et al., 2014], a program C induces an MDP M C in which each position corresponds to a configuration (s, C ′ ) consisting of a state s and a program C ′ . Each terminal position of M C is in bijection with a state s, and we write ts for the terminal position corresponding to s.

We define the terminating probability of C for reaching a terminal state s′ from an initial state s under a scheduler σ as the terminating probability TPσ(s,C),ts′ on M C . The resulting posterior distribution of C from s under σ is the family (TPσ(s,C),ts′ )s′∈St over states.

This posterior distribution forms a subdistribution, i.e.

and the sum may be strictly less than 1. In particular, we do not assume almost-sure termination.

2.3 Target Specification and Problem

We conclude this section by formally defining our verification problem. Given a safe region R, the goal is to verify whether the program is always safe; that is, whether every posterior distribution (TPσ(s,C),ts′ )s′∈St lies in R for all schedulers σ.

We formally define the safe region R as follows. First, fix a partition (S1, . . . , Sn) of St. For a scheduler σ, we define the posterior distribution over the partition (S1, . . . , Sn) from an initial state s as

The safe region R is a non-empty convex set of the form

where θ ∈ Q≥0 is a threshold and f : [0, 1]n → R≥0 is a monotone polynomial with rational coefficients.

Note that Rθ is downward-closed: if (p1, . . . , pn) ⪯ (p′1 , . . . , p′n ) and (p′1 , . . . , p′n ) ∈ Rθ, then (p1, . . . , pn) ∈ Rθ. This follows from the monotonicity of f .

Our target problem is now formulated as follows:

Algorithm 1 CEGAR Loop

Require: a program C, an initial state sinit, a partition (S1, . . . , Sn) of St, and a non-empty convex downward closed set Rθ := {(p1, . . . , pn) ∈ [0, 1]n | f(p1, . . . , pn) ≤ θ} with a polynomial function f .

• 1: Initialize h by [0, 1]n that stores results in the iteration.
• 2: while True do
• 3: if h implies supσ f(P C,σ 1 , . . . , P C,σ n ) ≤ θ then
• 4: return Yes and terminate the loop.
• 5: else
• 6: Update h based on the found counter example.
• 7: end if
• 8: end while

Problem Statement: Given a program C and an initial state s ∈ St, the verification problem asks whether, for every scheduler σ, the induced posterior distribution (P σ s,S1 , . . . , P σ s,Sn ) belongs to Rθ. Equivalently, the problem reduces to checking whether the following inequality holds:

3 Overview of Our Approach: Exploiting Adequacy and Duality

To present our automated verification method for the target problem, we introduce two key observations that naturally lead to our approach. Alg. 1 provides an overview of the algorithm. There, we maintain an approximation (abstraction) h of C, whose precise details are not needed at this point.

The essence can be summarized as follows: (1) we attempt to prove Eq. (1) using the current approximation h (line 3); and (2) if h is not precise enough to show Eq. (1), we refine h using a counterexample (line 6).

Specifically, the first point is handled by the adequacy of the predicate transformer semantics with respect to the operational semantics [Batz et al., 2024]. The second point is addressed using the Fenchel-Moreau duality from convex optimization (cf. [Boyd and Vandenberghe, 2004]).

3.1 Adequacy

We first address how sufficiently precise abstractions h allow us to conclude that all posterior distributions lie in Rθ. We explain our approach as follows. Assume we have a finite sequence (w1, u1), . . . , (wm, um) of pairs consisting of a weight vector wi ∈ [0, 1]n and a bound ui ∈ [0, 1] such that

We call such a finite sequence a valid sequence. The following property then holds:

Lemma 1. Given a valid sequence (w1, u1), . . . , (wm, um), if the following property holds:

Then Eq. (1) holds.

Proof. By the definition of a valid sequence, every posterior distribution (P σ s,S1 , . . . , P σ s,Sn ) lies in the set {(p1, . . . , pn) ∈ [0, 1]n ∣∣ wi · (p1, . . . , pn) ≤ ui} for each i ∈ [m]. This yields the desired implication.

Conveniently, Eq. (2) can be checked as follows. Since the set

is a convex polytope in [0, 1]n , it has finitely many vertices p1, . . . , pℓ that generate it, i.e. it is the smallest closed convex set containing these vertices (see, e.g., [Boyd and Vandenberghe, 2004]).

As Rθ is also a closed convex set, it suffices to check whether all vertices p1, . . . , pℓ lie in Rθ, which is straightforward since f is a polynomial, assuming all wi and ui are rational. Such vertices can be extracted by standard algorithms (e.g. via polyhedral libraries including [Jeannet and Min´e, 2009]). In fact, in our CEGAR loop the valid sequence does serve as the abstraction h.

Now an important question is how to obtain such a valid sequence. We achieve this by exploiting the adequacy [Batz et al., 2024] of the predicate transformer semantics (or weakest pre-expectation semantics) [Kozen, 1985; McIver and Morgan, 2005] with respect to the operational semantics.

Formally, for a program C, the predicate transformer semantics yields a function ΦC : [0, 1]St → [0, 1]St , called the predicate transformer of C, which serves as a denotational semantics. Adequacy characterizes how ΦC relates to the operational semantics of C.

Lemma 2 (adequacy). Given a program C and a vector w ∈ [0, 1]n (called weight), we have ΦC (w)(s) = supσ w · (P σ s,S1 , . . . , P σ s,Sn ), where w: St → [0, 1] is given by w(s) = w(i) if s ∈ Si; recall that (S1, . . . , Sn) is a partition of St.

Proof. It follows by Theorem 6 in [Batz et al., 2024].

Thanks to well-known properties of predicate transformers-such as compositionality and least-fixed-point characterizations (see, e.g., [Batz et al., 2024])-there exist automated verification methods (e.g., [Chatterjee et al., 2017; Takisaka et al., 2021; Batz et al., 2023]) that compute an upper bound u on ΦC (w), establishing

Given a weight vector w ∈ [0, 1]n , we may be possible to obtain an upper bound u by these existing automated methods, which in turn becomes an element in a valid sequence.

3.2 Fenchel-Moreau Duality

We now explain how to refine the abstraction h, which consists of a valid sequence (w1, u1), . . . , (wm, um), in the case where h does not imply Eq. (2). Recall that the safe region R is assumed to be a non-empty downward-closed convex subset of the unit cube [0, 1]n . Since f is a polynomial, R is also closed by continuity.

From a duality perspective, R admits a representation as the intersection of half-spaces. Formally, we have the following representation theorem:

Proposition 1 (duality). A safety region R ∈ [0, 1]n can be rewritten by

Proof. This follows by the bijective correspondence of Fenchel conjugates. To see this, get the indicator function δR : [0, 1]n → [0,∞] of R defined by δR(p) := 0 if p ∈ R and δR(p) := ∞ otherwise. Since δR is a proper lowersemi continuous convex function, δR has its Fenchel conjugate (δR)∗ , where the mapping satisfies (δR)∗∗ = δR by the Fenchel-Moreau theorem. By calculating (δR)∗∗ , we have the half-space representation. Note that we only need nonnegative weight vectors w ∈ [0, 1]n since R is downward closed in [0, 1]n . See, e.g., [Boyd and Vandenberghe, 2004] for the detail of the Fenchel conjugate.

An important consequence of the duality theorem concerns certificates for counterexamples p /∈ R, where in our algorithm p is one of the vertices p1, . . . , pℓ generating the convex polytope. By Prop. 1, the fact that p /∈ R guarantees the existence of a weight w ∈ [0, 1]n such that

which serves as a certificate for p /∈ R. Within our CEGAR loop, such a weight w is automatically synthesised based on the discovered counterexample. Once w is obtained, we compute an upper bound u with respect to w such that u ≤ supq∈Rw · q < w · p using existing techniques, and refine h by updating

4 Our Algorithm

In this section, we present our algorithm in full detail. We already sketched the overall workflow in Alg. 1; here, we elaborate on the procedure described in lines 3-7 of Alg. 1.

Specifically, Alg. 2 presents the procedure in detail. Recall that h = (wi, ui)i∈I is a valid sequence, and let CP(h) denote the convex polytope induced by h. To verify whether this valid sequence entails Eq. (1) via Lem. 1, we first compute the vertices p1, . . . , pm of CP(h) (line 1). We then initialize variables ϵ ∈ R≥0 and r,w ∈ [0, 1]n , which will encode a counterexample in case h fails to imply Eq. (1).

For each vertex p ∈ {p1, . . . , pm}, we execute the subroutine in lines 4-9. First, we check whether p ∈ Rθ, where

as shown in line 4. Since f is a polynomial, this membership test is straightforward.

If p /∈ Rθ, we estimate how far p is from the safe region Rθ. To this end, we solve the optimisation problem in line

5. Intuitively, we compute the closest point q ∈ Rθ to p, which yields a weight w separating p and q, meaning that w · q < w · p. Note that the minimiser q is unique because the function ∥ · ∥22 is strictly convex. Moreover, the following desirable property holds:

Lemma 3. Assume that q is the minimiser and w = p - q. We have w ∈ [0, 1]n and

Proof. It is easy to see that w ∈ [0, 1]n because Rθ is downward closed. Since q is the minimiser of ∥q-p∥2 in Rθ, we have (p - q) · (x - q) ≤ 0 for all x ∈ Rθ. This implies (p - q) · x ≤ (p - q) · q, for all x ∈ Rθ. Since q ∈ Rθ, we have the desired equality.

Crucially, Lem. 3 implies that the weight w separates Rθ and p, since maxx∈Rθ w · x < w · p. It is therefore natural to choose w := p - q and, by Lem. 2, attempt to verify ΦC (w)(s) ≤ w · q in order to exclude the point p from the current abstraction (valid sequence) h.

Since verifying ΦC (w)(s) ≤ w · q is often the bottleneck, we only store the minimizer q and the weight w of the farthest point p, as handled in lines 6-8.

In line 11, we check whether ϵ = 0, which indicates that all points p1, . . . , pm lie in Rθ. Otherwise, we attempt to exclude the farthest point p /∈ Rθ by finding an upper bound u on ΦC (w)(s) such that u ≤ w · q.

Whenever such an upper bound is found, our algorithm strictly refines the valid sequence h. We formally state this progress property below:

Theorem 1 (progress). In line 14, suppose that we find an upper bound u that satisfies the condition, then the new valid sequence h′ := (w1, u1), . . . , (w|I|, u|I|), (w, u) satisfies CP(h′ ) ⊊ CP(h), where h is the previous valid sequence.

Proof. Take the farthest point p of CP(h) that induces the weight w. Since the upper bound u satisfies u ≤ w·q < w·p, it follows that p /∈ CP(h′ ). Equivalently, the hyperplane {x | w·x = w·q} separates p from q (and even from the set Rθ). By construction, we clearly have CP(h′ ) ⊆ CP(h).

Remark 1. Computing the exact minimizer q may not be possible, as q is not necessarily rational 1 . Nevertheless, progress can be ensured if, instead of the minimizer, we find a point q ∈ Rθ such that q ⪯ p and obtain an upper bound u with respect to the weight w = p-q satisfying u <w · p. In our implementation, this condition is enforced as an additional constraint.

We summarise this section by showing the soundness of our algorithm.

Theorem 2 (soundness). If Alg. 2 terminates, then Eq. (1) holds.

Proof. It follows by Lem. 1 and the fact that R is a closed convex set.

1 We use Z3 [de Moura and Bjørner, 2008] to solve the minimization problem. To the best of our knowledge, it is not clear whether Z3 ensures optimality of the solutions it returns.

Algorithm 2 Refinement by a counter-example

Require: the Requirement in Alg. 1, and h = (wi, ui)i∈I.

• 1: Compute the set p1, . . . , pm of vertices of the bounded convex polytope defined by h in [0, 1]n .
• 2: Set a threshold ϵ := 0, and initialize r ∈ [0, 1]n and w ∈ [0, 1]n .
• 3: for p in p1, . . . , pm do
• 4: if p ̸∈ Rθ then
• 5: Solve the following problem and synthesize the unique solution q:

subject to q ∈ Rθ is substochastic (∑i qi ≤ 1).

• 6:
• 7: Update ϵ ← ∥p - q∥22 , r ← q, and w ← p - q.
• 8: end if
• 9: end if
• 10: end for
• 11: if ϵ = 0 (implies all p1, . . . , pm ∈ Rθ) then
• 12: return Yes (Eq. (1) is verified)
• 13: else
• 14: Compute an upper bound u that satisfies u ≤ w · r w.r.t. the weight w, and refine h with (w, u).
• 15: end if

Remark 2 (limitation). We conclude by noting some limitations of our approach. First, it cannot decide that Eq. (1) does not hold. Naively, showing this would require synthesizing a posterior distribution (P σ s,S1 , . . . , P σ s,Sn ) that violates Eq. (1), which looks challenging and non-trivial; we leave this direction for future work.

A second limitation concerns the termination of upperbound computation in line 14. Our method relies on invariant synthesis via the Knaster-Tarski theorem and, as is common for probabilistic programs over infinite state spaces, termination is not guaranteed. Nevertheless, the procedure is sound: once it terminates, it yields a valid upper bound.

5 Experimental Evaluation

We evaluate the effectiveness of our approach through preliminary experiments conducted with our prototype implemented in MuValQFL2 . Our results demonstrate that the prototype shows promise in addressing the target problem, while also revealing several limitations that present opportunities for future improvement.

5.1 Setting

We implemented our CEGAR loop (Alg. 1), which uses Alg. 2 for refinement, to verify the target problem. To compute the vertices of convex polytopes (line 1 of Alg. 2), we employ Apron [Jeannet and Min´e, 2009], and we invoke Z3 [de Moura and Bjørner, 2008] to solve the optimization problem in line 5 of Alg. 2. In line 14 of Alg. 2, we perform invariant synthesis on top of the existing solver [Chatterjee et al., 2025]. We attempt to synthesise an inductive invariant-i.e., a certificate for an upper bound justified by the Knaster-Tarski theorem-using piecewise polynomial templates. For this, our prototype tries to synthesize a template that admits an inductive invariant by applying Farkas's lemma [Farkas, 1902] and positivity theorems [Handelman, 1988; Putinar, 1993].

2 Our implementation is available in https://github.com/ hiroshi-unno/coar.

For the Case of Linear Constraints. For the safe region Rθ represented by a linear function f(p1, . . . , pn) = a1p1 + · · ·+anpn, we employ a simpler procedure instead of Alg. 2. By Lem. 2, the problem can be directly verified if we obtain an upper bound u satisfying ΦC (w)(s) ≤ u ≤ θ, where w = (a1, . . . , an). Note that all coefficients ai are nonnegative since f must be monotone. Under this condition, our procedure terminates in at most one iteration for the linear constraint.

5.2 Experimental Results

We constructed 17 benchmarks, including the motivating example in §1, to evaluate the performance of our approach. To the best of our knowledge, no existing benchmarks are available for our target problems. We designed our benchmarks using a variety of programs to assess how the size of partitions for posterior distributions-referred to as the dimension-and the degree of the polynomial f affect the performance of our implementation.

Following [Kura et al., 2025], each benchmark is written as a fixed-point equation system translated from programs written in pGCL, following the procedure in [Kura et al., 2025]. Roughly speaking, this translation amounts to providing the characteristic functions [Batz et al., 2023] of while loops. For each benchmark, we assign a safe region Rθ as a polynomial function f, assuming that Rθ is induced by f satisfying the required assumptions.

The results are summarized in Tab. 1. The experiments were conducted on a machine equipped with an Intel(R) Core(TM) Ultra 5 125U processor and 16 GB of RAM. Dim denotes the dimension of posterior distributions (i.e., the natural number n in (P σ 1 , . . . , P σ n )). Deg denotes the degree of the polynomial defining the safe region; for example, the degree of f(x, y) = x2 + y2 is 2. Num of Itr indicates the number of CEGAR iterations until termination, Inv denotes the running time of invariant synthesis in seconds, and OMT denotes the running time to solve the optimisation problem (by Z3). The total time for the entire procedure is almost the same as the sum Inv (sec) +OMT(sec), so we omit it from the table. This is because the time required to generate the generators of the convex polytope is very short due to the low dimensionality and the small number of iterations.

5.3 Evaluation

Overall, our CEGAR loop appears to work effectively in low-dimensional settings, particularly in two- and threedimensional spaces. We now illustrate how our CEGAR loop works on the benchmarks Choice2 and 3Dim2 in detail.

Example 1 (the benchmark Choice2). The program used in the benchmark is given by

where

C1 = 1/2 * (x:=x) + 1/8 * (x := 0) + 1/4 * (x := 1) + 1/8 * (x := 2), C2 = 1/2 * (x:=x) + 1/4 * (x := 0) + 1/8 * (x := 1) + 1/8 * (x := 2).

The safe region R is given by {(x, y) ∈ [0, 1]2 | x2 + y2 ≤ 9/16}, over the two-dimensional posterior distribution (P σ s,[x=1] , P σ s,[x=0] ). Here [x = n] is the state that assigns the value n to the variable x for each n. From the initial state x = 4, our CEGAR loop consructs the following valid sequence h := (w1, u1), (w2, u2), (w3, u3): 1) w1 = (0, 55/144) and u1 = 55/288; 2) w2 = (63/128, 0) and u2 = 63/256; and 3) w3 = (1/2, 9/16) and u3 = 127/256.

Example 2 (the benchmark 3Dim2). The program used in the benchmark is given by

where

C1 = 1/2 * (x:=x) + 1/16 * (x := 0) + 1/4 * (x := 1) + 1/16 * (x := 2) + 1/8 * (x := 3), C2 = 1/2 * (x:=x) + 1/4 * (x := 0) + 1/16 * (x := 1) + 1/16 * (x := 2) + 1/8 * (x := 3).

The safe region R is given by {(x, y, z) ∈ [0, 1]3 | x2 + y2 + z2 ≤ 5/8}, over the three-dimensional posterior distribution (P σ s,[x=0] , P σ s,[x=1] , P σ s,[x=2] ). From the initial state x = 4, our CEGAR loop consructs the following valid sequence h := (w1, u1), . . . , (w4, u4): 1) w1 = (0, 137/448, 0) and u1 = 137/896; 2) w2 = (991/3072, 0, 0) and u2 = 991/6144; 3) w3 = (95/3072, 0, 1/2) and u3 = 351/4096; and 4) w4 = (3/4, 49/64, 1/2) and u4 = 2527/4096.

We observe that tighter specifications tend to require more complex certificates, which are challenging for existing solvers to obtain. In fact, Choice3 and Choice4 share the same underlying program, and because of the tighter specification, our prototype cannot solve Choice4. Lastly, we remark that all benchmarks that require only one iteration have degree 1.

5.4 Limitations

We discuss several limitations of our approach observed through the implementation. The backend solver used in our prototype [Chatterjee et al., 2025; de Moura and Bjørner, 2008] does not support exponential templates, whereas synthesizing inductive invariants to obtain an upper bound u (line 14 of Alg. 2) often requires such templates. Consequently, we only include benchmarks that do not demand highly expressive templates. In our preliminary experiments, our implementation synthesizes only piecewise affine invariants when it successfully computes upper bounds for the benchmarks.

Another limitation concerns scalability with respect to the dimensionality. As shown in Tab. 1, our implementation begins to struggle from dimension 4. One reason is that increasing the number of partitions significantly enlarges the queries generated by [Chatterjee et al., 2025], making satisfiability checking expensive for SMT solvers (we use Z3 [de Moura and Bjørner, 2008]).

Another difficulty arises from the instability of SMT solvers when solving the optimization problem in line 5 of Alg. 2 if the degree is 2. For the benchmark 4Dim2, Z3 fails to solve an optimization problem generated during our CEGAR loop. Notably, by using OptiMathSAT [Sebastiani and Trentin, 2020] instead of Z3, the benchmark 4Dim2 can be solved with two iterations. However, employing OptiMathSAT as the backend SMT solver prevents our implementation from solving most of the other benchmarks; we list the result with OptiMathSAT in Tab. 2.

6 RelatedWork

In this section, we discuss the relationship between our approach and existing work.

6.1 Verifying Posterior Distributions

Posterior distributions of probabilistic programs have been actively studied as a challenging problem. In particular, posterior distributions have been used to reason about equivalence and similarity between programs [Murawski and Ouaknine, 2005]. Such relational reasoning has received considerable attention, and various logical frameworks and sensitivity analyses have been developed in this context [Culpepper and Cobb, 2017; Albarghouthi and Hsu, 2018; Barthe et al., 2016; Huang et al., 2018].

Notably, the recent work [Chatterjee et al., 2024] introduces an automated method for relational reasoning over probabilistic programs with infinite state spaces. Their framework and ours are largely orthogonal, as they address different verification problems. We remark that their framework does not support nondeterminism in probabilistic programs and relies on technical assumptions such as almostsure termination, whereas our approach supports nondeterminism and does not require almost-sure termination.

6.2 Multi-objective Optimisation

For finite-state MDPs, computing or approximating the set of posterior distributions induced by schedulers has been well studied. Etessami et al. [Etessami et al., 2008] propose multi-objective model-checking algorithms based on approximating the Pareto curve of an MDP. For the multi-objective reachability problem, each point on the Pareto curve corresponds to a Pareto-optimal posterior distribution under a (possibly stochastic) scheduler. Their approach employs the approximation algorithm of Papadimitriou and Yannakakis [Papadimitriou and Yannakakis, 2000], and Forejt et al. [Forejt et al., 2012] later propose an iterative approximation algorithm based on the so-called sandwich algorithm [Rennen et al., 2011; Solanki et al., 1993].

These techniques do not apply to probabilistic programs with infinite state spaces, although our CEGAR loop is inspired by the sandwich algorithm. In particular, the idea of selecting a "good" weight vector w to improve the approximation is central to the sandwich algorithm. In this paper, we develop a novel CEGAR loop for probabilistic programs that synthesizes such a weight w based on violations of the current abstraction with respect to a specification (given as a safe region).

6.3 Healthiness and Duality

Our framework is by no means the first to focus on FenchelMoreau duality in the context of probabilistic programs. Indeed, Morgan et al. [Morgan et al., 1996] and McIver and Morgan [McIver and Morgan, 2005] already observed this duality in their study of healthiness conditions for predicate transformers. Such duality aspects in probabilistic programs are well-studied in studies of denotational semantics of higher-order probabilistic programs [Goubault-Larrecq, 2007; Keimel and Plotkin, 2009]

This duality for predicate transformers has further been generalized categorically [Hino et al., 2016], encompassing predicate transformers for finite-state probabilistic programs. Our contribution is to build on this dual viewpoint to develop a novel automated verification method.

7 Conclusion

We formally define a safety problem for posterior distributions of probabilistic programs. To verify such properties automatically, we introduce a novel CEGAR loop that iteratively refines abstractions represented as convex polytopes. We assess the effectiveness of our approach through a preliminary experiment.

As noted in Rem. 2, it would be desirable to automatically synthesize a posterior distribution that serves as a counterexample to Eq. (1). We leave this seemingly challenging problem for future work.

A bottleneck of our approach is the lack of support for exponential templates, which are often needed to compute upper bounds on ΦC (w). The backend solver we rely on [Chatterjee et al., 2025] currently supports only polynomial templates. A backend capable of handling exponential templates would be highly desirable, potentially by leveraging LLMs to synthesize such invariants; we leave this as future work.

While we focus on proving the safety of probabilistic programs, it would also be useful to be able to disprove the safety of probabilistic programs. This seems to be a more challenging and non-trivial problem, but it may be possible to provide an efficient sound procedure by restricting the language, for example by disallowing non-determinism. We leave this as future work.

Lastly, we aim to improve the scalability of our approach to higher-dimensional problems. In addition to the SMT query size issue observed in §5, several challenges arise in higher dimensions. For example, approximating closed convex sets using half-spaces generally becomes more difficult as the dimensionality increases-a phenomenon also reported in the literature (e.g., [Watanabe et al., 2024]) for finite-state MDPs.

Acknowledgements

We thank the anonymous reviewers for their helpful comments. This study was supported by JSPS KAKENHI Grant Number JP25K24739. K.W. was supported by the JST grant No. JPMJPR25KD.