Gilles Barthe

Benjamin Grégoire

Justin Hsu

Pierre-Yves Strub

Coupling proofs are probabilistic product programs

Couplings are a powerful mathematical tool for reasoning about pairs of probabilistic processes. Recent developments in formal verification identify a close connection between couplings and pRHL, a relational program logic motivated by applications to provable security, enabling formal construction of couplings from the probability theory literature. However, existing work using pRHL merely shows existence of a coupling and does not give a way to prove quantitative properties about the coupling, which are need to reason about mixing and convergence of probabilistic processes. Furthermore, pRHL is inherently incomplete, and is not able to capture some advanced forms of couplings such as shift couplings. We address both problems as follows.

First, we define an extension of pRHL, called ×pRHL, which explicitly constructs the coupling in a pRHL derivation in the form of a probabilistic product program that simulates two correlated runs of the original program. Existing verification tools for probabilistic programs can then be directly applied to the probabilistic product to prove quantitative properties of the coupling. Second, we equip ×pRHL with a new rule for while loops, where reasoning can freely mix synchronized and unsynchronized loop iterations. Our proof rule can capture examples of shift couplings, and the logic is relatively complete for deterministic programs.

We show soundness of ×pRHL and use it to analyze two classes of examples. First, we verify rapid mixing using different tools from coupling: standard coupling, shift coupling, and path coupling, a compositional principle for combining local couplings into a global coupling. Second, we verify (approximate) equivalence between a source and an optimized program for several instances of loop optimizations from the literature.

Introduction

The coupling method [START_REF] Levin | Markov chains and mixing times[END_REF][START_REF] Lindvall | Lectures on the coupling method[END_REF][START_REF] Thorisson | Coupling, Stationarity, and Regeneration[END_REF][START_REF] Villani | Optimal transport: old and new[END_REF] is an elegant mathematical tool for analyzing the relationship between probabilistic processes. Informally, couplings correlate outputs of two probabilistic processes by specifying how corresponding sampling statements are correlated; reasoning about the correlated processes can then imply interesting properties of the original processes.

A classical application of couplings is showing that two probabilistic processes converge in distribution. Consider, for instance, a symmetric simple random walk over Z: starting from some initial position p, the process repeatedly samples a value s uniformly in [Copyright notice will appear here once 'preprint' option is removed.] {-1, 1} and updates its position to p + s. As the process evolves, the distribution on position spreads out from its initial position, and converges to a limit distribution that is the same for all initial positions. This property can be proved by constructing a coupling where the probability that the two coupled walks end in the same position approaches 1 as we run process.

Beyond merely showing convergence, typically we are interested in how quickly the processes converge. The rate of convergence, sometimes measured by the mixing time of a process, is highly important [START_REF] Levin | Markov chains and mixing times[END_REF]. A common use of probabilistic processes is to efficiently sample from a distribution that approximates a complicated distribution (e.g., Metropolis et al. [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF]). For instance, the Glauber dynamics [START_REF] Jerrum | A very simple algorithm for estimating the number of k-colorings of a low-degree graph[END_REF] approximately samples from the uniform distribution on proper colorings of a graph-a hard distribution to compute-by maintaining a coloring and randomly re-coloring a single vertex at a time. More generally, the Glauber dynamics is an example of Markov Chain Monte-Carlo (MCMC), a family of techniques that underlie many computational simulations in science, and machine learning algorithms for performing probabilistic inference (see Brooks et al. [START_REF] Brooks | Handbook of Markov Chain Monte Carlo[END_REF] for a survey). The mixing time determines how many steps we need to run the process before we are approximately sampling from the target distribution.

Recently, Barthe et al. [START_REF] Barthe | Relational reasoning via probabilistic coupling[END_REF] noticed a close connection between couplings and probabilistic relational Hoare logic (pRHL), originally designed for reasoning about the computational security of cryptographic constructions [START_REF] Barthe | Formal certification of code-based cryptographic proofs[END_REF]. Namely, every valid derivation in pRHL implies the existence of a coupling of the output distributions of the two programs. Using this observation, Barthe et al. [START_REF] Barthe | Relational reasoning via probabilistic coupling[END_REF] verify pRHL judgments that imply properties of random walks, the ballsin-bins process, and the birth-death process. Barthe et al. [START_REF] Barthe | Proving differential privacy via probabilistic couplings[END_REF] extend the connection to approximate couplings in approximate probabilistic relational Hoare logic (apRHL), introduced for reasoning about differentially private computations [START_REF] Barthe | Probabilistic relational reasoning for differential privacy[END_REF], and exploit this connection to prove differential privacy of several examples that had previously been beyond the reach of formal verification.

While pRHL is a useful tool for constructing couplings, it has two limitations. First, it cannot directly reason about the two coupled processes. This poses a problem for proving mixing and convergence properties. For instance, Barthe et al. [START_REF] Barthe | Relational reasoning via probabilistic coupling[END_REF] prove that if a certain property P holds on the coupled samples, then two random walks meet under the coupling. By a theorem about random walks, this means that the (total variation) distance between the two distributions is at most the probability of P under the coupling. However, we do not know what this probability is or how it grows as we run the random walk for more iterations, since pRHL cannot reason about the coupled process. We run into similar difficulties if we try to prove convergence using path coupling [START_REF] Bubley | Path coupling: A technique for proving rapid mixing in Markov chains[END_REF], a general construction that shows fast convergence by upper-bounding the expected distance when we make a transition from two coupled states. While we can express the transition function in code, pRHL cannot reason about expectations.

Second, pRHL cannot express some natural classes of couplings because the rule for while loops requires that both loops execute in lockstep. For instance, shift couplings [START_REF] Aldous | Shift-Coupling[END_REF] allow the two processes to meet at a random time shift, e.g., the first process at time t1 could track the second process at some other time t2 = t1 + δ. This incompleteness of pRHL is also a limitation for more standard applications of pRHL to program analysis, like validation of loop transformations, and more generally for any optimization that alters the control flow of the programs.

We address these problems as follows. First, we deepen the connection between pRHL and couplings with an observation that is reminiscent of proofs-as-programs: not only do pRHL judgements correspond to couplings, but pRHL proofs encode a probabilistic product program that constructs the distribution witnessing the coupling. This program is similar to existing product program constructions [START_REF] Barthe | Relational verification using product programs[END_REF][START_REF] Zaks | CoVaC: Compiler Validation by Program Analysis of the Cross-Product[END_REF] in that it simulates two program executions with a single program. However, the probabilistic product also coordinates the samplings in the two executions, as specified by the coupling encoded in the original pRHL derivation. Second, we propose a general rule for while loops; our rule subsumes several existing rules and is sufficiently expressive to capture several examples of shift coupling. The resulting system, which we call product pRHL or ×pRHL for short, has several benefits.

First, we obtain a simple, algorithmic procedure to construct the probabilistic product given a pRHL derivation. The product directly simulates the coupled processes, so we can prove quantitative properties about probabilities or expected values of this coupling using existing probabilistic verification techniques. Moreover, intermediate assertions in a pRHL derivation to be soundly transported to the probabilistic product in ×pRHL derivations. For instance, we may prove synchronized control flow in pRHL, and directly assume this fact in the probabilistic product. Many facts are often practically easier to manipulate in pRHL, since pRHL works purely on non-probabilistic assertions.

Second, we obtain a powerful logic that can reason about many examples from the coupling literature, especially shift couplings, and from the translation validation literature, especially loop optimizations. On the foundational side, we prove that our logic is relatively complete for deterministic programs.

Summary of contributions

We make the following contributions.

• We present a proof-relevant program logic ×pRHL that extracts a probabilistic product program from a valid derivation, and prove (using the Coq proof assistant) that the logic is sound.

• We propose new rules for loops and random sampling in ×pRHL.

We also prove relative completeness for deterministic programs.

• We demonstrate several applications of ×pRHL to showing convergence of probabilistic processes; several of these examples use path coupling, a compositional tool for constructing couplings which bounds convergence in terms of expected properties of the coupling, and shift coupling, a generalization of coupling where the two processes are allowed to meet at different times. We also show how to validate some common loop transformations.

Preliminaries

We begin by giving a bird's eye view on probabilistic couplings, which take output distributions from two probabilistic processes and place them in the same probabilistic space.

In the following, we will work with sub-distributions over discrete (finite or countable) sets. Definition 1. A discrete sub-distribution over a set A is completely defined by its mass function µ : A → [0, 1] such that a∈A µ(a) is defined and bound by 1. The quantity a∈A µ(a) ∈ [0, 1] is called the weight of µ and denoted by |µ|. The support supp(µ) of µ is defined as {x ∈ A | µ(x) = 0} and is discrete (i.e., countable) by construction. A distribution is a sub-distribution with weight 1. The probability of an event P w.r.t. µ, written Prµ[P ] (or Pr[P ] when µ is clear from the context), is defined as {x∈A|P (x)} µ(x).

One can equip distributions with a monadic structure using the Dirac distributions 1 x for the unit1 and distribution expectation Ex∼µ[M (x)] for the bind, where

Ex∼µ[M (x)] : x → a µ(a) • M (a)(x).
When working with sub-distributions over tuples, the probabilistic versions of the usual projections on tuples are called marginals. For distributions over pairs, we define the first and second marginals π1(µ) and π2(µ) of a distribution µ over A × B by π1(µ)(a) = b∈B µ(a, b) and π2(µ)(b) = a∈A µ(a, b). For a distribution µ over memories and a set of variables X ⊆ Vars, we define the X-marginal distribution of µ as

πX (µ)(x) = v∈Vars\X µ(x, v).
We can also construct a sub-distribution over tuples from two distributions: the product sub-distribution of sub-distributions µ1 and µ2 over A and B respectively is defined by

(µ1 × µ2)(a, b) = µ1(a) • µ2(b).
We are now ready to formally define coupling. Definition 2. Two sub-distributions µ1, µ2 resp. over A and B are said to be coupled with witness µ ∈ Distr(A × B), written

µ µ1 & µ2 , iff π1(µ) = µ1 and π2(µ) = µ2. For Ψ ⊆ A × B, we write µ Ψ µ1 & µ2 if µ
µ1 & µ2 and moreover supp(µ) ⊆ Ψ. We will often abuse notation and call µ a coupling of µ1 and µ2.

As an example, suppose that A = B and µ1 = µ2 are the uniform distributions. Then, any bijection f : A → A gives a coupling of µ1 and µ2; we call the resulting coupling D f (A) so that D f (A)

µ1 & µ2 . The coupling D f (A) assigns positive probability only to pairs (v, f v) with v ∈ A.

For another example, suppose again that A = B and µ1 = µ2, but the distributions are not necessarily uniform. Then, the identity map id : A → A always gives a coupling of µ1 and µ2, correlating samples from both distributions to be the same. We will write D = (µ1) for the resulting coupling, so that D = (µ1)

µ1 & µ2 . If µ1 is the uniform distribution we will sometime write D = (A) instead of D = (µ1). Note that D = (µ1) assigns positive probability only to pairs (v, v) with µ1(v), µ2(v) = 0.

To reason about convergence of probabilistic processes, we will use the total variation distance on distributions (also known as statistical distance). Definition 3. Let µ1 and µ2 be sub-distributions over a countable set A. The total variation (TV) distance between µ1 and µ2 is defined by:

TV(µ1, µ2) = 1 2 a∈A |µ1(a) -µ2(a)|.
To bound this distance, it is enough to find a coupling and bound the probability that the two coupled variables differ; this is the fundamental idea underlying the coupling method.

Theorem 4 (Fundamental theorem of couplings (e.g., [START_REF] Lindvall | Lectures on the coupling method[END_REF])). Let µ1 and µ2 be distributions over a countable set. Then for any coupling µ of µ1, µ2, we have

TV(µ1, µ2) ≤ Pr (x,x )∼µ [x = x ].

Product programs

Language

We will work with a core, probabilistic imperative language with a command for random sampling from primitive distributions. The set of commands is defined as follows:

C ::= skip noop | abort abort | X ← E deterministic assignment | X $ ← D(E) probabilistic assignment | C; C sequencing | if E then C else C conditional | while E do C while loop
Here, X is a set of variables, E is a set of (deterministic) expressions, and D is a set of distribution expressions. Variables and expressions are typed, ranging over booleans, integers, lists, etc. The expression grammar is entirely standard, and we omit it.

We will use several shorthands for commands:if e then c for if e then c else skip, abort for the looping command while do skip, and c e,k for the k-fold composition of c restricted to e. Formally, c e,k ı = 0; while (ı < k) ∧ e do c; ı++ where ı is a fresh variable. We allow k to be an arbitrary expression; when k is a constant, c e,k ≡ k times if e then c; . . . ; if e then c for the simple notion of program equivalence introduced in Figure 2.

The denotational semantics of programs is adapted from the seminal work of Kozen [START_REF] Kozen | Semantics of probabilistic programs[END_REF]. We first interpret every ground type T as a set T ; other constructors C are interpreted as functions C that respect their arities. The semantics of expressions and distribution expressions is then parameterized by a state m, and is defined in the usual way. Last, commands are interpreted as a function from memories to sub-distributions over memories, where memories are finite maps from variables to values. Formally, we let StateX denote the set of memories over the finite set of variables X. Moreover, we use to denote the standard disjoint union on finite maps, so : StateX 1 × StateX 2 → StateX 1 ∪X 2 for disjoint finite sets X1 and X2. The interpretation of c, written c , is a function from StateX to Distr(StateX ), where Vars(c) = X. The definition of • is given in Figure 1.

Last, for any predicate over memories Φ, we say that a command is Φ-lossless iff for any memory m s.t. Φ(m), the weight of c m is equal to 1. A command is lossless if it is -lossless, i.e. if | c m | = 1 for every memory m.

Proof system

Our proof system manipulates judgments of the form:

Φ c1 c2 Ψ c
where c1, c2 and c are statements over disjoint variables, and Φ and Ψ are assertions over the variables of c1 and c2. We will not fix a specific syntax for assertions, but one natural choice is first-order formulas over the program variables of both programs.

Informally, the judgment is valid if c is a probabilistic product program for c1 and c2 under the pre-condition Φ, i.e. for every initial memory m = m1 m2 satisfying the pre-condition Φ, the sub-distribution c m is a coupling for c1 m 1 and c2 m 2 , and moreover supp( c m ) only contains states that satisfy Ψ.

Definition 5 (Valid judgment).

• Two commands c1 and c2 are separable iff their variables are disjoint, i.e. Vars(c1) ∩ Vars(c2) = ∅.

• A judgment Φ c1 c2 Ψ c is valid iff the commands c1 and c2 are separable and for every memory m = m1 m2 such that m Φ, we have

c m Ψ ( c1 m 1 ) & ( c2 m 2 ) .
For comparison, judgments in pRHL are of the form

Φ c1 c2 Ψ
and assert that for every initial memory m = m1 m2 such that m Φ, there exists a coupling µ such that

µ Ψ ( c1 m 1 ) & ( c2 m 2 ) .
In contrast, our notion of judgment is proof-relevant, since the derivation of the judgment guides the product construction. We briefly comment on several rules in our system, presented in Fig. 3. The [CONSEQ] rule is similar to the rule of consequence, and reflects that validity is preserved by weakening the post-condition and strengthening the pre-condition; none of the programs is modified in this case. The [FALSE] rule corresponds to the fact that under a false precondition nothing needs to be proved.

The [CASE] rule allows proving a judgment by case analysis; specifically, the validity of a judgment with pre-condition Φ can be established from the validity of two judgments, one where the pre-condition is strengthened with e and the other where the precondition is strengthened with ¬e. The product program is of the form if e then c else c , where c and c correspond to the product programs built from the proofs of the first and second premises respectively.

The [STRUCT] rule allows replacing programs by provably equivalent programs. The rules for proving program equivalence are given in Figure 2, and manipulate judgments of the form Φ c ≡ c , where Φ is a relational assertion. We only introduce equations that are needed for recovering derived rules, striving to keep the notion of structural equivalence as simple as possible.

The [ASSG] rule corresponds to the pRHL rule for assignments; in this case, the product program is simply the sequential composition of the two assignments.

The [RAND] rule informally takes a coupling between the two distributions used for sampling in the left and right program, and produces a product program that draws the pair of correlated samples from the coupling. Since our language supports sampling from distribution expressions and not only distributions, the rule asks for the existence of a coupling for each interpretation of the distribution expressions under a valuation satisfying the pre-condition Φ of the judgment. Furthermore, note that the rule requires that every element in the support of the coupling validates the post-condition. This is similar to pRHL, and is a natural requirement given the notion of valid judgment. However, the rule is strictly more general than the corresponding pRHL rule, which is restricted to the case where Ψ is the graph of some bijection.

The [WHILE] rule for while loops constructs a product program that interleaves synchronous and asynchronous executions of the loop bodies. The first premise establishes that k1 and k2 are strictly positive when the loop invariant Ψ holds. Then, we specify an expression e-which may mention variables from both sides-that holds true exactly when at least one of the guards is true.

skip m = 1 m abort m = 0 x ← e m = 1 m[x:= e m ] x $ ← g m = E v∼ g m [1 m[x:=v] ] c1; c2 m = E ξ∼ c 1 m [ c2 m ] if e then s1 else s2 m = if e m then c1
Next, the notation ⊕{p0, p1, p2} indicates that exactly one of the tests p0, p1, and p2 must hold. These predicates guide the construction of the product program. If p0 holds, then both guards should be equal and we can execute the two sides k1 and k2 iterations respectively, preserving the loop invariant Ψ. If p1 holds and the right loop has not terminated yet, then the left loop also has not terminated yet and we may execute the left loop one iteration. If p2 holds and the left loop has not terminated yet, then the right loop also has not terminated yet and we may execute the right loop one iteration.

The last pair of premises deal with termination. Note that some condition on termination is needed for soundness of the logic: if the left loop terminates with probability 1 while the right loop terminates with probability 0 (i.e., never), it is impossible to construct a valid product program since there is no distribution on pairs that has first marginal with weight 1 and second marginal with weight 0. So, we require that the first and second loops are lossless assuming p1 and p2 respectively. This ensures that with probability 1, there are only finitely many steps where we execute the left or right loop separately. Note that the product program may still fail to terminate with positive probability since there may be infinite sequences of iterations where p0 holds so both loops advance, but both programs will yield sub-distributions with the same weight so we can still find a coupling.

With these premises, the construction of the product program is straightforward. At the top level, the product continues as long as e is true, i.e., while at least one of the two loops can make progress. Each iteration, it performs a case analysis on the three predicates. If p0 holds, then we execute the product from running the two loops k1 and k2 iterations respectively. If p1 holds, then we execute the product from running the left loop one iteration; if p2 then we execute the product from running the right loop one iteration.

Derived rules

Presentations of relational Hoare logic often include one-sided rules, which are based on the analysis of a single program (rather than two). By reasoning about program equivalence, we can derive all one-sided rules and the two-sided rules of conditionals and loops from pRHL, presented in Fig. 4, within our system.2 Proposition 6. All the rules in Figure 4 (and their symmetric version) are derived rules.

We briefly comment on some of the derived rules. The [ASSG-L] rule is the one-sided rule for assignment. It can be derived from its two-sided counterpart, the [ASSG] rule, using the [STRUCT] rule. It can be derived from its two-sided counterpart, the [RAND] rule, using the [STRUCT] rule and the fact that the product distribution µ1 × µ2 is a valid coupling of proper distributions µ1 and µ2.

Φ c ≡ c Φ c1 ≡ c2 Φ c2 ≡ c1 Φ x $ ← 1 x ≡ skip Φ (x1, x2) $ ← µ1 × µ2 ≡ x1 $ ← µ1; x2 $ ← µ2 Φ =⇒ x = e Φ x ← e ≡ skip Φ c; skip ≡ c Φ skip; c ≡ c Φ c1 ≡ c 1 Φ c1; c2 ≡ c 1 ; c2 c2 ≡ c 2 Φ c1; c2 ≡ c1; c 2 Φ =⇒ e Φ if e then c else c ≡ c Φ =⇒ ¬e Φ if e then c else c ≡ c Φ ∧ e c1 ≡ c2 Φ ∧ ¬e c 1 ≡ c 2 Φ if e then c1 else c 1 ≡ if e then c2 else
We have one-sided and two-sided rules for conditionals. The [COND-L] rule is the one-sided version; it can be derived from the [CASE] and [STRUCT] rules. The [COND-S] rule is the two-sided version; the rule assumes that the two guards of the conditional statements are synchronized, so that one must only need to reason about the cases where both statements enter the true branch, and the case where both statements enter the false branch. It can be derived from the [CASE], [STRUCT], and [FALSE] rules.

We also have one-sided and two-sided rules for loops. The [WHILE-L] rule corresponds to the pRHL one-sided rule for while loops; it can be derived from the [WHILE] rule by setting p1 = , p0 = p2 = ⊥, and k1 = k2 = 1, e = e1, e2 = ⊥, c2 = skip and using the [STRUCT] and [FALSE] rules. The rule [WHILE-S] corresponds to the pRHL two-sided rule for while loops. This rule assumes that the two loops are synchronized, i.e., the guards of the two loops are equal assuming the loop invariant. This rule can be derived from the general [WHILE] rule by setting p0 = ,

Φ c1 c2 Ψ c Φ =⇒ Φ Ψ =⇒ Ψ Φ c1 c2 Ψ c FALSE ⊥ c1 c2 Ψ skip STRUCT Φ c1 c2 Ψ c Φ c1 ≡ c 1 Φ c2 ≡ c 2 Φ c ≡ c Φ c 1 c 2 Ψ c CASE Φ ∧ e c1 c2 Ψ c Φ ∧ ¬e c1 c2 Ψ c Φ c1 c2 Ψ if e then c else c ASSG Ψ[e1/x1, e2/x2] x1 ← e1 x2 ← e2 Ψ x1 ← e1; x2 ← e2 RAND ∀m.m Φ =⇒ µ m {(v 1 ,v 2 )|m[x 1 :=v 1 ,x 2 :=v 2 ] Ψ} µ1 m & µ2 m Φ x1 $ ← µ1 x2 $ ← µ2 Ψ (x1, x2) $ ← µ SEQ Φ c1 c2 Ξ c Ξ c 1 c 2 Ψ c Φ c1; c 1 c2; c 2 Ψ c; c WHILE Ψ =⇒ k1 > 0 ∧ k2 > 0 Ψ =⇒ (e1 ∨ e2) = e Ψ ∧ e =⇒ ⊕{p0, p1, p2} Ψ ∧ p0 ∧ e =⇒ e1 = e2 Ψ ∧ p1 ∧ e =⇒ e1 Ψ ∧ p2 ∧ e =⇒ e2 while e1 ∧ p1 do c1 Ψ-lossless while e2 ∧ p2 do c2 Ψ-lossless Ψ ∧ p0 c e 1 ,k 1 1 c e 2 ,k 2 2 Ψ c0 Ψ ∧ e1 ∧ p1 c1 skip Ψ c1 Ψ ∧ e2 ∧ p2 skip c2 Ψ c2
Ψ while e1 do c1 while e2 do c2 Ψ ∧ ¬e1 ∧ ¬e2 while e do if p0 then c0 else if p1 then c1 else c2

Soundness and (relative) completeness

We have formally verified the soundness theorem below in the Coq proof assistant, with Ssreflect.

Theorem 7 (Soundness). Every derivable judgment is valid.

Proof sketch. By induction on the proof derivation. The most interesting case is the loop rule, which we briefly sketch here. We will use the same notations from Fig. 3. Proving that the product program satisfies the post-condition is immediate: by construction, the condition e does not hold after the execution of the loop. Moreover, the premises ensure that the loop body of the product preserves the invariant Ψ. So, we just need to prove that the semantics of the projections of the product program correspond to respective semantics of the original programs. We detail the first projection, i.e. that for any initial memory m s.t. m |= Ψ, we have

w1 π 1 m = π1( w m ) (1)
where w while e do c and wi while ei do ci; the other projection is similar. By definition, we have for any e and c w m = lim n→∞ ce,n m Since ce,n m is monotonic in n, proving Eq. ( 1) amounts to showing

c1 e 1 ,n π 1 m ≤ π1( w m ) (2) π1( ce,n m ) ≤ w1 π 1 m (3)
for all n. We proceed by induction on n. Let mi πi(m) for i ∈ {1, 2}. The base case is immediate and we consider the n + 1 case. We will prove (2) by case analysis on e1 m and e m . If e1 m does not hold, then c1 e 1 ,n m 1 = 1 m 1 . Now if e2 m is false, then e m is also false and (2) holds. Otherwise if e2 m is true, p0 m and p1 m cannot hold, since that would imply that e1 m holds. So in this case, we can conclude (2) since w m = 1 m 1 × w2 m 2 and w2 is Ψ-lossless.

Otherwise, e1 m holds. We proceed by case analysis on p0 m , p1 m and p1 m , as they are mutually exclusive by the premise. ; c1 e 1 ,n m 1 which allows us to conclude using the induction hypothesis, since w m = c0; w m in this case and c e 1 ,k 1

5 2016/9/21 ASSG-L Ψ[e1/x1] x1 ← e1 skip Ψ x1 ← e1 RAND-L µ1 lossless ∀v1 ∈ supp(µ1), Ψ[v1/x1] x1 $ ← µ1 skip Ψ x1 $ ← µ1 COND-L Φ ∧ e1 c1 c2 Ψ c Φ ∧ ¬e1 c 1 c2 Ψ c Φ if e1 then c1 else c 1 c2 Ψ if e1 then c else c COND-S Φ =⇒ e1 = e2 Φ ∧ e1 c1 c2 Ψ c Φ ∧ ¬e1 c 1 c 2 Ψ c Φ if e1 then c1 else c 1 if e2 then c2 else c 2 Ψ if e1 then c else c WHILE-L Ψ ∧ e1 c1 skip Ψ c1 while e1 do c1 Ψ-lossless Ψ while e1 do c1 skip Ψ ∧ ¬e1 while e1 do c1 WHILE-S Ψ =⇒ e1 = e2 Ψ ∧ e1 c1 c2 Ψ c Ψ while e1 do c1 while e2 do c2 Ψ ∧ ¬e1 while e1 do c
1

m 1 = π1( c0 m ) • If p1 m holds, we have c1 e 1 ,n+1 m 1 = c1; c1 e 1 ,n m 1 w m = c1; w m
which allows to conclude using the induction hypothesis.

• If p2 m holds, we have

c1 e 1 ,n+1 π 1 m = E m 1 ∼π 1 while e 2 ∧p 2 do c 2 ;w m [ c1 e 1 ,n+1 m 1 ].
Note that after the executing while e2 ∧ p2 do c2-which is lossless-Ψ and e1 still hold while e2 does not. By the premises of the rule, p2 must be false so p0 or p1 holds on every memory m 1 in the support. By the first two cases with initial memory m 1 , we have:

c1 e 1 ,n+1 m ≤ π1 w m .
We can then conclude:

c1 e 1 ,n+1 π 1 m ≤ E m 1 ∼π 1 while e 2 ∧p 2 do c 2 ;w m [π1 w m ] (linearity) = π1 (E m ∼ while e 2 ∧p 2 do c 2 ;w m [ w m ]) (separable) = π1 w m . (definition)
Where the second to last step is because while e2 ∧ p2 do c2 does not modify π1 m since c1, c2 are separable.

The proof of Eq. (3) follows the same pattern.

Although it is not a primary objective of our work, we briefly comment on completeness of the logic. First, the coupling method is not complete for proving rapid mixing of Markov chains. 3 Second, it is not clear that our proof system is complete with respect to hoisting random assignments out of loops.

However, we note that the deterministic fragment of our logic achieves completeness for programs that satisfy a sufficiently strong termination property; the key is that the new rule for while subsumes self-composition for while loops, provided they terminate on all initial memories satisfying the invariant. More precisely, we can prove the following completeness theorem. Theorem 8. Let c1 and c2 be separable deterministic programs. If {Φ} c1; c2 {Ψ} is derivable using the rules from Figure 9, then Φ c1 c2 Ψ _ is derivable. Therefore, ×pRHL is relatively complete for deterministic programs.

Proof sketch. It suffices to prove that if {Φ} c1 {Ψ} is derivable using Hoare logic, then

Φ c1 skip Ψ _
is derivable. The proof proceeds by induction on the derivation.

Convergence from couplings

The fundamental theorem of couplings (Theorem 4) gives a powerful method to prove convergence of random processes. First, we recast it in terms of ×pRHL.

Proposition 9. Let c1 and c2 be separable programs and assume that the following judgment is valid:

Φ c1 c2 Ψ =⇒ x1 = x2 c
Then for every memory m such that m |= Φ, we have

TV(µ1, µ2) ≤ Pr x← c m [¬Ψ]
where µ1, µ2 are the distributions obtained by sampling x1 from c1 m and x2 from c2 m respectively.

This result is a direct consequence of the soundness of the logic, and allows proving convergence in two parts. First, we use ×pRHL to establish a valid post-condition of the form Ψ =⇒ x1 = x2. Second, we prove that for every memory m satisfying some precondition, the product program c built from the derivation satisfies:

Pr x← c m [¬Ψ] ≤ β
There are multiple approaches for proving properties of this formreasoning directly about the semantics of programs; existing formalisms for bounding probabilities and reasoning about expectations (e.g., [START_REF] Kozen | A probabilistic PDL[END_REF][START_REF] Morgan | Probabilistic predicate transformers[END_REF]); program logics for probabilistic programs (e.g., [START_REF] Barthe | A program logic for union bounds[END_REF][START_REF] Hartog | Probabilistic extensions of semantical models[END_REF][START_REF] Ramshaw | Formalizing the Analysis of Algorithms[END_REF]). We will check the property on pen and paper; mechanizing the proofs is left for future work.

Application: convergence of Markov chains

We now turn to our first group of examples: proving convergence of probabilistic processes. Suppose we have a probabilistic process on a set Ω of possible states. At each time step, the process selects the next state probabilistically. Consider two runs of the same probabilistic process started from two different states in Ω. We would like to know how many steps we need to run before the two distributions on states converges to a common distribution. We consider several classic examples.

Notation. Throughout this section, we consider two copies of the same program. To ensure that the two copies are separable, we tag all the variables of the first copy with 1, and all the variables of the second copy with 2.

Simple, symmetric random walk

x ← s; ı ← 0; l ← []; while ı < T do r $ ← {-1, 1}; x ← x + r; l ← x :: l; ı ← ı + 1; return x Figure 5. Random walk
Our first example is a simple random walk on the integers. Let the state space Ω be Z. At each step, the process chooses uniformly to move left (decreasing the position by 1) or right (increasing the position by 1). The program rwalk in Figure 5 implements the process, executed for T steps. The variable l is a ghost variable. While it does not influence the process, it keeps track of the list of sampled values, and will be used to state assertions when we construct the coupling. Now, consider the random walk started from starting positions s1, s2. We want to show that the TV-distance between the two distributions on positions decreases as we run for more steps T ; roughly, the random walk forgets its initial position. It is not hard to see that if s1 -s2 is an odd integer, then we will not have convergence: at any timestep t, the support of one distribution will be on even integers while the support of the other distribution will be on odd integers, so the TV-distance remains 1.

When s1 -s2 = 2k is even, we can construct a coupling to show convergence. Barthe et al. [START_REF] Barthe | Relational reasoning via probabilistic coupling[END_REF] used pRHL to couple these random walks by mirroring; informally, the coupled walks make mirror moves until they meet, when they make identical moves to stay equal. Specifically, they show that

Φ rwalk1 rwalk2 k ∈ psum(rev(l1)) =⇒ x1 = x2
where Φ s1 + s2 = 2k, rev(l) reverses the list, and psum(l) is the list of partial sums of l (sums over its initial segments). We can lift the judgment to ×pRHL

Φ rwalk1 rwalk2 k ∈ psum(rev(l1)) =⇒ x1 = x2 rwalk0
where rwalk0 is the following product program:

x1 ← s1; x2 ← s2; ı1 ← 0; ı2 ← 0; l1 ← []; l2 ← []; while ı1 < T do if x1 = x2 then (r1, r2) $ ← D = ({-1, 1}); else (r1, r2) $ ← Dopp({-1, 1}); x1 ← x1 + r1; x2 ← x2 + r2; l1 ← x1 :: l1; l2 ← x2 :: l2; ı1 ← ı1 + 1; ı2 ← ı2 + 1; return (x1, x2) where opp x -x.
We briefly sketch the derivation of the product. We start by an application of the [WHILE-S] rule with the invariant

Ψ k ∈ psum(rev(l1)) =⇒ x1 = x2 ∧k / ∈ psum(rev(l1)) =⇒ x1 -x2 = 2k -sum(l1),
where sum(l) is the sum of the list. Then, we apply a [SEQ] rule (consuming the first random sampling on each side) with intermediate assertion

Ξ Ψ ∧ (x1 = x2 =⇒ r1 = r2) ∧ (x1 = x2 =⇒ r1 = -r2)
The sub-proof obligation on tails is straightforward, the interesting one is for the random sampling. We start by using the [CASE] rule with e x1 = x2, which introduce the conditional in the product program. If the equality holds, then the two random values are synchronized using the D = ({-1, 1}); if not, they are mirrored using Dopp({-1, 1}). In our mirror coupling, k ∈ psum(l1) implies that both the have already met, and continue to have the same position. However, the derivation by itself does not tell us how far the two distributions are, as a function of T . To get this information, we will use the probabilistic product construction and the following classical result from the theory of random walks.

Theorem 10 (e.g., [START_REF] Avena | Random walks[END_REF]). Let X0, X1, . . . be a symmetric random walk on the integers with initial position X0 = 0. Then, for any position k ∈ Z, the probability that the walk does not reach k within t steps is at most

Pr[X0, . . . , Xt = k] ≤ ke √ 2 π √ t .
Now, we can analyze how quickly the two walks mix.

Theorem 11. If we perform a simple random walk for T steps from two positions that are 2k apart, then the resulting distributions µ1 and µ2 on final positions satisfy TV(µ1, µ2) ≤ ke √ 2 π √ T . Formally, for every two memories m1 and m2 such that m1(x) -m2(x) = 2k,

7 2016/9/21
we have

TV( rwalk m 1 , rwalk m 2 ) ≤ ke √ 2 π √ T ,
Proof. Conceptually, we can think of the difference x1 -x2 as a random walk which increases by 2 with half probability and decreases by 2 with half probability. By applying Theorem 10 to this random walk, we find that in the product program c0,

Pr[k / ∈ psum(l1)] ≤ ke √ 2 π √ T .
Then, we can conclude by Proposition 9.

The Dynkin process

Our second example models a process called the Dynkin process.

There is a sequence of N concealed cards, each with a number drawn uniformly at random from {1, . . . , 10}. A player starts at some position in {1, . . . , 10}. Repeatedly, the player looks at the number at their current position, and moves forward that many steps. For instance, if the player reveals 2 at their current location, then she moves forward two spaces. The player stops when she passes the last card of the sequence. We want to show fast convergence of this process if we start from any two initial positions. In code, the Dynkin process is captured by the program dynkin defined in Figure 6. Here, l stores the history of positions of the player; this ghost variable will be useful both for writing assertions about the coupling, and for assertions in the product program. Just like for random walks, we can consider the mixing rate of this process, starting from two positions s1, s2.

x ← s; l ← [x]; while x < N do r $ ← [1, 10]; x ← x + r; l ← x :: l; return x Figure 6. Dynkin process We will couple two runs of the Dynkin process from two starting positions asynchronously: we will move whichever process is behind, holding the other process temporarily fixed. If both processes are at the same position, then they move together.

Formally, we prove the following ×pRHL judgment:

dynkin 1 dynkin 2 (l1 ∩ l2 = ∅) =⇒ x1 = x2 dynkin 0
where dynkin 0 is the following program:

x1 ← s1; x2 ← s2; l1 ← [x1]; l2 ← [x2]; while x1 < N ∨ x2 < N do if x1 = x2 then (r1, r2) $ ← D = ([1, 10]); x1 ← x1 + r1; x2 ← x2 + r2; l1 ← x1 :: l1; l2 ← x2 :: l2; else if x1 < x2 then r1 $ ← [1, 10]; x1 ← x1 + r1; l1 ← x1 :: l1; else r2 $ ← [1, 10]; x2 ← x2 + r2; l2 ← x2 :: l2; return (x1, x2)
To couple the two programs, we use the rule [WHILE] with k1 = k2 = 1. To control which process will advance, we define:

p0 (x1 = x2) p1 (x1 < x2) p2 (x1 > x2).
The lossless conditions are satisfied since the distance between x1 and x2 strictly decreases at each iteration of the loops. For the samplings, when p1 or p2 hold we use the one-sided rules for random sampling ([RAND-L] and corresponding [RAND-R]); when p0 holds, we use the identity coupling.

Theorem 12. Let m1 and m2 two memories such that m1(s), m2(s) ∈ [1, 10], and suppose N > 10. Then:

TV( dynkin m 1 , dynkin m 2 ) ≤ ( 9 /10) N/5-2 .
The two distributions converge exponentially fast as N grows.

Proof. In the product program dynkin 0 , we want to bound the probability that l1 and l2 are disjoint; i.e., the probability that the two processes never meet. We proceed in two steps. First, it is not hard to show that |x1 -x2| < 10 is an invariant. Thus, at each iteration, there is a 1 /10 chance that the lagging process hits the leading process. Second, each process moves at most 10 positions each iteration and we finish when both processes reach the end, so there are at least 2( N /10 -1) = N /5 -2 chances to hit. Therefore, in the product we can show

Pr[(l1 ∩ l2) = ∅] ≤ 9 /10) N/5-2 . ( 4
)
We can then conclude by Proposition 9.

To highlight the quantitative information verified by our approach, we note that the corresponding theorem for random walks (Theorem 11) shows that the total variation distance between two random walks decreases as O(1/ √ T ). In contrast, Theorem 12 shows that total variation distance between two Dynkin processes converges as O(0.9 N ), giving much faster mixing (exponentially fast instead of polynomially fast).

Remark. The Dynkin process is inspired by the following twoplayer game. Both players pick a starting position. There is one sequence of random cards that is shared by both players, and players look at the card at their current position and move forward that many spaces. The random cards are shared, so a player samples the card only if the other player has not yet visited the position. If the other player has already landed on the position, the later player looks as the revealed card and moves forward.

While the Dynkin process samples every card that the player lands on, it is not hard to see that the product program exactly models the two-player game. More specifically, the product interleaves the players so that at each turn, the player that is lagging behind makes the next move. By this scheduling, as long as one player is lagging behind, the players have not landed on the same position and so each player lands on unseen cards and draws random samples to reveal. Once the players meet, the product program makes the same move for both players.

For this game, Eq. ( 4) bounds the probability that the players do not land on the same final position from any two initial positions. This result is the basis of the magic trick called Dynkin's card trick or the Kruskal count. If one player is the magician and the other player is a spectator, if the spectator starts at a secret position and runs the process mentally, the magician can guess the correct final position with high probability by starting at any position and counting along.

Path coupling

So far, we have seen how to prove convergence of probabilistic processes by constructing a coupling, reasoning about the probability of the processes not meeting, and applying the coupling theorem. While this reasoning is quite powerful, for more complicated processes it may be difficult to directly construct a coupling that shows fast mixing. For example, it can be difficult to find and reason about a coupling on the distributions from two states s, s if there are many transitions apart in the Markov chain.

To address this problem, Bubley and Dyer [START_REF] Bubley | Path coupling: A technique for proving rapid mixing in Markov chains[END_REF] proposed the path coupling technique, which allows us to consider just pairs of adjacent states, that is states where s can transition in one step to s . Roughly speaking, if we can give a good coupling on the distributions from two adjacent states for every pair of adjacent states, then path coupling shows that the state distributions started from two arbitrary states converge quickly.

As the name suggests, path coupling considers paths of states in a probabilistic process. For this to make sense, we need to equip the state space with additional structure. For the basic setup, let Ω be a finite set of states and suppose that we have a metric d : Ω × Ω → N. We require that d is a path metric: if d(s, s ) > 1, then there exists s = s, s such that d(s, s ) = d(s, s ) + d(s , s ). Two states are said to be adjacent if d(s, s ) = 1. We will assume that the diameter of the state space, i.e. the distance between any two states, is finite: ∆ < ∞. The Markov chain is then defined by iterating a transition function f : Ω → Distr(Ω).

The main idea behind path couplings is that if we can couple the distributions from any two adjacent states, then there exists a coupling for the distributions from two states at an arbitrary distance, constructed by piecing together the couplings along the path. Furthermore, if the expected distance between states contracts under the coupling on adjacent states, i.e., the resulting expected distance is strictly less than 1, then the same holds for the coupling on two states at any distance. More formally, we have the following main theorem of path coupling.

Lemma 13 (Bubley and Dyer [START_REF] Bubley | Path coupling: A technique for proving rapid mixing in Markov chains[END_REF]). Consider a Markov chain with transition function f over a set Ω with diameter at most ∆. Suppose that for any two states s and s such that d(s, s ) = 1, there exists a coupling µ of f (s), f (s ) such that

E (r,r )∼µ [d(r, r )] ≤ β.
Then, starting from any two states s and s and running t steps of the chain, there is a coupling µ of f t (s), f t (s ) such that

TV(f t (s), f t (s )) ≤ Pr (r,r )∼µ [r = r ] ≤ β t ∆.
This lemma applies for all β, but is most interesting for β < 1 when it implies that the Markov chain mixes quickly. With the main theorem in hand, we will show how to verify the conditions for path coupling on two examples from Bubley and Dyer [START_REF] Bubley | Path coupling: A technique for proving rapid mixing in Markov chains[END_REF].

Graph coloring: the Glauber dynamics.

Our first example is a Markov chain to provide approximately uniform samples from the set of colorings of a finite graph; it was first analyzed by Jerrum [START_REF] Jerrum | A very simple algorithm for estimating the number of k-colorings of a low-degree graph[END_REF], our analysis follows Bubley and Dyer [START_REF] Bubley | Path coupling: A technique for proving rapid mixing in Markov chains[END_REF] Recall that a finite graph G is defined by a finite set of vertices V , and a symmetric relation E relating pairs vertices that are connected by an edge; we will let N G(v) ⊆ V denote the set of neighbors of v, i.e the set vertices that have an edge to v in G. Let C be the set of colors; throughout, we assume that C is finite. A coloring w of G is a map from V to C. A coloring is valid (sometimes called proper) if all neighboring vertices have different colors: for all v ∈ N G(v) we have w(v) = w(v ). We write VG(w) if w is a valid coloring.

The following program glauber(T ) models T steps of the Glauber dynamics in statistical physics:

ı ← 0; while ı < T do; v $ ← V ; c $ ← C; if VG(w[v → c]) then w ← w[v → c] ı ← ı + 1
Informally, the algorithm starts from a valid coloring w and iteratively modifies it by sampling uniformly a vertex v and a color c, recoloring the v with c if this continues to be a valid coloring.

We want to measure the statistical distance between two executions of the process starting from two initial colorings w1 and w2. There are two natural approaches. The first option is to couple the two copies of glauber directly, analyze the product program and apply Proposition 9. The problem is that when the two colorings are far apart, it is hard to reason about how the processes might meet under a coupling; Jerrum [START_REF] Jerrum | A very simple algorithm for estimating the number of k-colorings of a low-degree graph[END_REF] adopted this strategy, but the resulting proof is dense and complex.

The second, far simpler option is to apply path coupling. Here, we build a product for just one iteration of the loop, and it suffices to consider cases where the two initial states are adjacent. This drastically simplifies the coupling and analysis of the product program, so we adopt this approach here. For the sake of clarity, we adapt the transition function so that its output is stored in a fresh variable w , and call the resulting program glauber † :

v $ ← V ; c $ ← C; if VG(w[v → c]) then w ← w[v → c] else w ← w
Note that glauber † ; w ← w ; ı ← ı + 1 is semantically equivalent to the loop body of glauber.

To apply the path coupling theorem Lemma 13, we need to define a path metric on Ω and construct a coupling for the process started from two adjacent states. For the path metric, we define the distance d(w1, w2) to be the Hamming distance: the number of vertices where w1 and w2 provide different colors. We say two states are adjacent if d(w1, w2) = 1; these states differ in the color of exactly one vertex. In order to apply path coupling, we need to find a coupling of the transition function on adjacent states such that the expected distance shrinks. We first build the coupling using ×pRHL. Specifically, we prove

d(w1, w2) = 1 glauber † 1 glauber † 2 d(w 1 , w 2 ) ≤ 2 c0 (5)
where c0 is the following program:

v1, v2 $ ← D = (V ); if v1 ∈ N G(v0) then c1, c2 $ ← D π ab (C) else c1, c2 $ ← D = (C) if VG(w1[v1 → c1]) then w 1 ← w1[v1 → c1] else w 1 ← w1 if VG(w2[v2 → c2]) then w 2 ← w2[v2 → c2] else w 2 ← w2
We briefly sketch how to prove the judgment. Note that the two states must agree at all vertices, except at a single vertex v0. Let w1(v0) = a and w2(v0) = b. First, we couple the vertex sampling with the rule [RAND] using the identity coupling, ensuring v1 = v2. Then, we use the rule [CASE] to perform a case analysis on the sampled vertex, call it v. If v is a neighbor of the differing vertex v0, we use the rule [RAND] and the transposition bijection π ab : C → C defined by the clauses:

π ab (a) = b π ab (b) = a π ab (x) = x otherwise
to ensure that c1 = π ab (c2). Otherwise, we use the rule [RAND] and the identity coupling to ensure c1 = c2. By applying the one-sided rules for conditionals to the left and the right sides ([COND-L] and [COND-R]), we can conclude the derivation. Next, we must reason about the expected value of the distance between w 1 and w 2 after executing the product program. Lemma 14. Let n = |V | and k = |C|, and suppose that the graph G has degree bounded by D. That is, for any v ∈ V , there are at most D vertices v such that E(v, v ). If k ≥ D , then there is a coupling µ the distributions after running glauber † on adjacent states such that

E (w 1 ,w 2 )∼µ [d(w 1 , w 2 )] ≤ 1 -1/n + 2D/kn. 9 2016/9/21
Proof. Let w1 and w2 be adjacent states. We must bound the expected distance between the states w 1 , w 2 in the product program. Let d = d(w 1 , w 2 ), we have:

E (w 1 ,w 2 )∼µ [d ] = 0 • Pr[d = 0] + 1 • Pr[d = 1] + 2 • Pr[d = 2] = 1 -Pr[d = 0] + Pr[d = 2] ≤ 1 -Pr[v1 = v0 ∧ VG(w )] + Pr[v1 ∈ N G(v0) ∧ c1 = b] ≤ 1 - 1 n 1 - D k + D nk = 1 - 1 n + 2D nk .
where w = w1[v0 → c1]. The first equality holds because the distance between the two resulting coloring will be at most 2 by judgment Eq. ( 5). The second equality holds since

1 = Pr[d = 0] + Pr[d = 1] + Pr[d = 2].
The second to last step follows since each vertex has at most D neighbors, so there are at least k -D valid colors at any vertex.

Applying the path coupling lemma (Lemma 13), noting that the diameter is n since there are n vertices, proves that the Glauber dynamics mixes quickly if there are sufficiently many colors k.

Theorem 15. Consider the Glauber dynamics on k colors with a graph G with n vertices and degree at most D, and suppose k ≥ 2D + 1. Then, for some constant β < 1, TV( glauber(T ) m 1 , glauber(T ) m 2 ) ≤ β T n for any two initial memories m1 and m2 containing valid colorings.

This theorem recovers the result by Bubley and Dyer [START_REF] Bubley | Path coupling: A technique for proving rapid mixing in Markov chains[END_REF]; this is the key step to showing that running Glauber dynamics for a small number of steps and taking a sample is almost equivalent to drawing a uniformly random sample from all proper colorings of the graph.

Independent sets: the conserved hard-core model

Our second example is from graph theory and statistical physics, modeling the evolution of a physical system in the conserved hardcore lattice gas (CHLG) model [START_REF] Bubley | Path coupling: A technique for proving rapid mixing in Markov chains[END_REF]. Suppose that we have a set P of particles, and we have a graph G. A placement is a map w : P → V , placing each particle at a vertex of the graph. We wish to place the particles so that each vertex has at most one particle, and no two particles are located at adjacent vertices; we call such a placement safe and denote it by SG(w). For a specific graph, there could be multiple safe placements.

If we want to sample a uniformly random safe placement, we can use a simple Markov chain. We take the state space Ω = P → V to be the set of placements. Again, we take G and P to be finite. We start using a safe initial placement. Each step, we sample a particle p from P and a vertex v from V uniformly at random and try to move p to v. If w[p → v] is safe, then we make this the new placement; otherwise, we keep the same placement. We can model T steps of this dynamics with the following program chlg(T ):

ı ← 0; while ı < T do; p $ ← P ; v $ ← V ; if SG(w[p → v]) then w ← w[p → v] ı ← ı + 1
As in the previous example, we adapt the loop body to form chlg † :

p $ ← P ; v $ ← V ; if SG(w[p → v]) then w ← w[p → v] else w ← w
Like the graph coloring sampler, we take the path metric on placements to be Hamming distance and try to find a coupling on the distributions from adjacent initial placements.

E (w 1 ,w 2 )∼µ [d(w 1 , w 2 )] ≤ 1 - 1 s 3(D + 1) n .
Proof. Let chlg † 1 , chlg † 2 be two copies of the transition function, with variables tagged. Consider two adjacent placements w1 and w2. We will sketch how to couple the transitions.

We use rule [RAND] twice to couple the particle and vertex samplings with the identity coupling, ensuring p1 = p2 and v1 = v2. Then, we can apply the one-sided rules for conditionals to the left and the right sides ([COND-L] and [COND-R]) to conclude the following judgment:

d(w1, w2) = 1 chlg † 1 chlg † 2 d(w1, w2) ≤ 2 c0,
where c0 is the following product program:

(p1, p2) $ ← D = (P ); (v1, v2) $ ← D = (V ); if SG(w1[p1 → v1]) then w 1 ← w1[p1 → v1] else w 1 ← w1 if SG(w2[p2 → v2]) then w 2 ← w2[p2 → v2] else w 2 ← w2
Then, we can bound the expected distance between w 1 , w 2 in the product program. Let d = d(w 1 , w 2 ), we have:

E (w 1 ,w 2 )∼µ [d ] = 1 -Pr[d = 0] + Pr[d = 2] = 1 -Pr[p = p0 ∧ SG(w1[p → v])] + Pr[p = p0 ∧ ¬(SG(w1[p → v]) ⇔ SG(w2[p → v]))] ≤ 1 -Pr[p = p0 ∧ SG(w1[p → v])] + Pr[p = p0 ∧ ¬(SG(w1[p → v]) ∧ SG(w2[p → v]))]
We can bound the two probability terms. For the first term, we know that the probability of selecting p = p0 is 1/s, and the probability that p is safe at v if it avoids all other points (at most s-1) and all the neighbors of the other points (at most (s -1)D); this probability is the same for both placements w1 and w2, since the two placements are identical on points besides p0.

For the second term, we know that the probability of selecting p = p0 is 1 -1/s, and p is not safe at v in placement w1 or in w2 if we select the position a, b, or one of their neighbors. Putting everything together, we can conclude:

E (w 1 ,w 2 )∼µ [d(w 1 , w 2 )] ≤ 1 - 1 s 1 - (s -1)(D + 1) n + (s -1) | N G(a)| + 1 + | N G(b)| + 1 sn ≤ 1 - 1 s 1 - (s -1)(D + 1) n + (s -1) 2(D + 1) sn = 1 - 1 s 3(D + 1) n .
Applying the path coupling lemma (Lemma 13) shows that if we iterate the transition function on two initial placements, the resulting distributions on placements converge quickly. Remark. This theorem is slightly weaker than the corresponding result by Bubley and Dyer [START_REF] Bubley | Path coupling: A technique for proving rapid mixing in Markov chains[END_REF], who prove rapid mixing under the weaker condition s ≤ n/2(D + 1) + 1. Roughly, they use the maximal coupling on the two transition distributions, giving a tighter analysis and better bound. It is also possible use the maximal coupling in ×pRHL, but the corresponding specification of the coupling distribution would be proved as part of soundness of the logic, rather than as a property of a probabilistic program.

Application: loop optimizations

Program equivalence is one of the original motivation for relational program logics [START_REF] Benton | Simple relational correctness proofs for static analyses and program transformations[END_REF]. In this section, we demonstrate the effectiveness of our logic using several examples of exact and approximate program equivalence. Our first example is a loop transformation which originates from the literature on parallelizing compilers but also has applications in computer-aided cryptography. Our second example is drawn from the recent literature on approximate computing, and is an instance of loop perforation.

Loop strip-mining

Loop strip-mining (or loop sectioning) is a transformation that turns a loop into a nested loop. While the transformation originates from literature on parallelizing compilers and is primarily used to take advantage of vectorized instructions, it is also useful for formally proving the computational security of certain cryptographic constructions. The following example, Fig. 7, is inspired from a proof of indifferentiability of the SHA3 hash function [START_REF] Canteaut | A unified indifferentiability proof for permutation-or block cipher-based hash functions[END_REF] 4 Using the rule for while loops, we can prove the following:

x1 = x2 c1 c2 x1 = x2 c
The crux of the proof is applying the [WHILE] rule with k1 = 1 and k2 = M , and e = ı1 < N , and p0 = , and p1 = p2 = ⊥ and an invariant Ψ which strengthens the assertion x1 = x2 mainly by adding l2 = ı1 * M . Side conditions using p1 and p2 are trivial to prove (using the FALSE rule) since they have ⊥ in hypothesis. It remains to check the premise for for p0, but we now have two synchronized loops; we can use the [STRUCT] rule to remove the conditional on p0 which is always true in this case.

Loop perforation

Loop perforation [START_REF] Misailovic | Probabilistically accurate program transformations[END_REF][START_REF] Sidiroglou-Douskos | Managing performance vs. accuracy trade-offs with loop perforation[END_REF] is a program transformation that delivers good trade-offs between performance and accuracy, and is practical in many applications, including image and audio processing, simulations and machine-learning. Informally, loop perforation transforms a loop that performs n iterations of its body into a loop that performs m < n iterations of its body, followed by a simple post-processing statement. Figure 8 shows an example of loop perforation inspired from a financial analysis application, called swaptions. In this example, every other loop iteration is skipped, and the post-processing statement simply multiplies by 2 the value s computed by the optimized loop. As for the previous example, we can prove the following judgment:

c1 c2 s1 = s2 c
The product program can be built using the [WHILE] rule. We use e = ı1 < 2 * n, p0 = , p1 = p2 = ⊥ and k1 = 2 and k2 = 1, and the invariant is ı1 = 2 * ı2. The invariant allows to show, using the Nested loop: Finally, in some applications the number of iterations performed by the perforated loop is probabilistic; for instance, the program

ı1 ← 0; while ı1 < N do 1 ← 0; while 1 < M do l1 ← ı1 * M + ; r1 $ ← µ; x1 ← f (l1, x1, r1); 1 ← 1 + 1; ı1 ← ı1 + 1; l1 ← N * M ; Single loop: l2 ← 0; while l2 < N * M do ı2 ← l2 / M ; 2 ← l2 % M ; r2 $ ← µ; x2 ← f (l2, x2, r2); l2 ← l2 + 1; ı2 ← N ; 2 ← M ; Product program: ı1 ← 0; l2 ← 0; while ı1 < N do 1 ← 0; while 1 < M do l1 ← ı1 * M + 1; ı2 ← l2 / M ; 2 ← l2 % M ; (r1, r2) $ ← D = (µ); x1 ← f (l1, x1, r1); x2 ← f (l2, x2, r2); 1 ← 1 + 1; l2 ← l2 + 1; ı1 ← ı1 + 1; l1 ← N * M ; ı2 ← N ; 2 ← M ; Figure 7. Loop strip-mining Original program: s1 ← 0; for (ı1 ← 0, ı1 < 2 * n, ı1 ← ı1 + 1) do x1 $ ← µ; s1 ← s1 + x1; Perforated program: s2 ← 0; for (ı2 ← 0, ı2 < n, ı2 ← ı2 + 1) do x2 $ ← µ; s2 ← s2 + x2; s2 ← 2 * s2; Product program: s1 ← 0; s2 ← 0; ı1 ← 0; ı2 ← 0; while ı1 < 2 * n do x1, x2 $ ← D = (µ); s1 ← s1 + x1; s2 ← s2 + x2; ı1 ← ı1 + 1; x1 $ ← µ; s1 ← s1 + x1; ı1 ← ı1 + 1; ı2 ← ı2 + 1; s2 ← 2 * s2;
k $ ← factors(n); s ← 0; for (ı ← 0, ı < n, i ← i + k) do x $ ← µ; s ← s + x; s ← k * s;
selects uniformly at random a factor k of the original number n of iterations, and performs n/k iterations. It is possible to relate the original and the perforated loop, using the [WHILE] rule as before.

11

2016/9/21

Other optimizations and program transformations

Barthe et al. [START_REF] Barthe | Relational verification using product programs[END_REF] define an inductive method for building valid product programs, and use their method for validating instances of loop optimizations. Their method combines a rule for each program construction and a rule akin to our [STRUCT] rule. Despite this similarity, the two methods are fundamentally different: their treatment of while loops is restricted to synchronized executions. As a consequence, their structural rule is based on a more advanced refinement relation between programs. Nevertheless, we can reproduce all their examples in our formalism, taking advantage of our more powerful rule for loops.

Related work

Relational logics can be seen as a proof-theoretical counterpart of semantics-based relational methods such as logical relations. Under this view, our logic bears strong similarities with proof-relevant logical relations [START_REF] Benton | Proof-relevant logical relations for name generation[END_REF]. As for proof-relevant logical relations, we expect that manipulating explicit witnesses rather than existentially quantified can help developing meta-theoretical studies of our logic.

Much of the recent work on product programs and relational logics has been motivated by applications to security and compiler correctness. For instance, Barthe et al. [START_REF] Barthe | Secure Information Flow by Self-Composition[END_REF] explore self-composition for a variety of programming languages and show that it induces a sound and complete reduction of an information flow policy to a safety property. Independently, Darvas et al. [START_REF] Darvas | A theorem proving approach to analysis of secure information flow[END_REF] consider self-composition and deductive verification based on dynamic logic, also for verifying information flow policies. Later, Terauchi and Aiken [START_REF] Terauchi | Secure information flow as a safety problem[END_REF] introduce the class of 2-safety properties and show a reduction from 2-safety to safety of the self-composed program. Their reduction is more efficient than self-composition as it selectively applies self-composition or a synchronous product construction akin to cross-products (described below). Further improvements appear in Kovács et al. [START_REF] Kovács | Relational abstract interpretation for the verification of 2-hypersafety properties[END_REF], Müller et al. [START_REF] Müller | An analysis of universal information flow based on self-composition[END_REF]. In a related work, Beringer and Hofmann [START_REF] Beringer | Secure information flow and program logics[END_REF] observe that one can encode 2-safety properties in standard Hoare logic, provided that assertions are sufficiently expressive to model the denotational semantics of programs. Beringer [START_REF] Beringer | Relational decomposition[END_REF] further refines this approach, by introducing the notion of relational decomposition.

Zaks and Pnueli [START_REF] Zaks | CoVaC: Compiler Validation by Program Analysis of the Cross-Product[END_REF] define a cross-product construction, which is well-suited for reasoning about programs with identical controlflow. Barthe et al. [START_REF] Barthe | Relational verification using product programs[END_REF] generalize the notion of cross-product by proposing a more general notion of product program which subsumes self-composition and cross-products, and show how it enables validation of common loop optimizations. Specifically, they define an inductive relation for proving that c is a valid product for c1 and c2; informally, their rules closely follow those of our system (for deterministic constructs), except for the general rule for while loops; instead, they use a rule that is closer to the pRHL rule, and a rule akin to our [STRUCT] rule, with a much stronger relationship between programs in order to compensate for the lack of generality of their rule. Barthe et al. [START_REF] Barthe | Product programs and relational program logics[END_REF] carry a more precise study of the relative expressiveness of product program constructions and relational program logics. In a different thread of work, Barthe et al. [START_REF] Barthe | Beyond 2-safety: Asymmetric product programs for relational program verification[END_REF] generalize the notion of product program so that it supports verification of refinement properties (modeled by universal quantification over runs of the first program and existential quantification over runs of the second program), as well as the 2-safety properties (modeled by universal quantification over runs of the first and the second programs). These constructions are focused on non-probabilistic programs. Motivated by applications to differential privacy, Barthe et al. [START_REF] Barthe | Proving differential privacy in Hoare logic[END_REF] define a specialized product construction from probabilistic programs to deterministic programs, so that the original program is differentially private, provided its deterministic product program satisfies some Hoare specification. To the best of our knowledge, this is the sole product construction that goes from a probabilistic language to a deterministic one.

Benton [START_REF] Benton | Simple relational correctness proofs for static analyses and program transformations[END_REF] and Yang [START_REF] Yang | Relational separation logic[END_REF] were among the first to consider relational program logics that support direct reasoning about two programs. Benton [START_REF] Benton | Simple relational correctness proofs for static analyses and program transformations[END_REF] introduces Relational Hoare Logic, proves correctness of several program transformations, and soundly embeds a type system for information flow security into his logic. Yang [START_REF] Yang | Relational separation logic[END_REF] defines Relational Separation Logic and proves the equivalence between Depth-First Search and the Schorr-Waite algorithm. Barthe et al. [START_REF] Barthe | Formal certification of code-based cryptographic proofs[END_REF] develop probabilistic Relational Hoare Logic, and use it for proving computational security of cryptographic constructions. In a follow-up work, Barthe et al. [START_REF] Barthe | Probabilistic relational reasoning for differential privacy[END_REF] develop an approximate variant of probabilistic Relational Hoare Logic, and verify differential privacy of several algorithms. More recently, Sousa and Dillig [START_REF] Sousa | Cartesian Hoare logic for verifying k-safety properties[END_REF] propose Cartesian Hoare Logic, an extension of relational Hoare logic to an arbitrary finite number of executions.

Several authors have considered relational logics for higher-order programs. Nanevski et al. [START_REF] Nanevski | Verification of information flow and access control policies with dependent types[END_REF] develop a relational logic to reason about information flow properties of a higher-order language with mutable state. Ghani et al. [START_REF] Ghani | Comprehensive parametric polymorphism: Categorical models and type theory[END_REF] introduce a relational type theory, and a supporting categorical model, for reasoning about parametricity. Barthe et al. [START_REF] Barthe | Probabilistic relational verification for cryptographic implementations[END_REF] propose a relational extension of a subset of the F * ; in a follow-up work, Barthe et al. [START_REF] Barthe | Higher-order approximate relational refinement types for mechanism design and differential privacy[END_REF] combine a relational refinement type system with a graded monad which they use for modeling differentially private computations.

There are several works that develop more specialized program logics for analyzing relational properties of programs. For instance, Amtoft et al. [START_REF] Amtoft | A logic for information flow in object-oriented programs[END_REF] introduce independence assertions and a supporting program logic for proving information flow security. In a similar way, Chaudhuri et al. [START_REF] Chaudhuri | Continuity analysis of programs[END_REF] propose a logical approach for proving continuity properties of programs, and Carbin et al. [START_REF] Carbin | Proving acceptability properties of relaxed nondeterministic approximate programs[END_REF] develop a logical approach for reasoning about the reliability of approximate computation.

Further afield, there has been a significant amount of work on semantical methods for probabilistic programs and processes initiated by Kozen [START_REF] Kozen | Semantics of probabilistic programs[END_REF], see e.g. [START_REF] Bizjak | Step-indexed logical relations for probability[END_REF][START_REF] Kozen | Kolmogorov extension, martingale convergence, and compositionality of processes[END_REF][START_REF] Lago | On coinductive equivalences for higher-order probabilistic functional programs[END_REF][START_REF] Sangiorgi | Environmental bisimulations for probabilistic higher-order languages[END_REF] for some recent developments. In conjunction with these semantics, research in deductive verification methods for non-relational properties of probabilistic programs is an active area of research; examples include [START_REF] Hartog | Probabilistic extensions of semantical models[END_REF][START_REF] Kozen | A probabilistic PDL[END_REF][START_REF] Morgan | Probabilistic predicate transformers[END_REF][START_REF] Ramshaw | Formalizing the Analysis of Algorithms[END_REF] to cite only a few systems.

Conclusion and future directions

We have introduced ×pRHL, a new program logic that deepens the connection between probabilistic couplings and relational verification of probabilistic programs in two different ways. First, ×pRHL broadens the class of couplings supported by relational verification. Second, ×pRHL derivations explicitly build a probabilistic product program, which can be used to analyze mixing times. We have shown the flexibility of our approach on several examples.

There is ample room for future work. On the theoretical side, it would be interesting to extend ×pRHL to handle continuous distributions as was recently done by Sato [START_REF] Sato | Approximate relational hoare logic for continuous random samplings[END_REF] for apRHL. Also, we believe that we are just scratching the surface of probabilistic product programs; there should be many further applications, notably in relationship with path couplings, in domains such as Brownian motion [START_REF] Lindvall | Lectures on the coupling method[END_REF], molecular evolution [START_REF] Dixit | A finite population model of molecular evolution: Theory and computation[END_REF], and anonymity [START_REF] Gomulkiewicz | Rapid mixing and security of chaum's visual electronic voting[END_REF]. On the more practical side, it would be natural to integrate ×pRHL in EasyCrypt [START_REF] Barthe | Easycrypt: A tutorial[END_REF], a proof-assistant used for reasoning about computational security of cryptographic constructions, and currently the main implementation of pRHL. We expect that several proofs of cryptographic constructions can be simplified using the new loop rule, and we are also planning to use the loop rule in an ongoing formalization of indifferentiability of the SHA3 standard for hash functions.

Hoare Logic for Deterministic Programs

We only consider deterministic programs in this section, i.e. programs without random samplings. In that case, the semantic s m of a deterministic command is either 0 if the program does not terminate, or of the form 1 m if it does. From that remark, we define the operator • * as follow:

s * m m if s m = 1 m ⊥ otherwise
A while loop while e do s is said to terminate certainly under condition Φ, written Φ ⇓ while e do s, iff for any m s.t. Φ(m), we have while e do s * m = ⊥. A Hoare triple is a triple of the form {Φ} s {Ψ} where Φ and Ψ are assertions and s is a deterministic command. A Hoare triple is said to be valid, written |= {Φ} s {Ψ}, iff for any memory m s.t. Φ(m), either s * m = ⊥ or Ψ( s * m ). Figure 9 gives a sound inference system for Hoare triples: If {Φ} s {Ψ} is derivable from the rules of Figure 9, then |= {Φ} s {Ψ}. Moreover, this system is complete for programs that only contain certainly terminating loops.

{Φ} skip {Φ}

  m else c2 m while b do c m = lim n→∞ c b,n m where c b,n k times if b then c; . . . ; if b then c; if b then abort

Figure 1 .

 1 Figure 1. Interpretation of commands

c 2 eFigure 2 .

 22 Figure 2. Equivalence rules

Figure 3 .

 3 Figure 3. Proof rules

Figure 4 .

 4 Figure 4. Derived rules

Lemma 16 .

 16 Let s = |P | and n = |V |, and suppose that the graph G has degree bounded by D. Starting from any two adjacent safe placements w1 and w2, there is a coupling µ on the distributions after one step such that

Theorem 17 .

 17 Consider the conserved lattice gas model with s = |P | particles on a graph G with n = |V | vertices and degree at most D, where s ≤ n/3(D + 1) + 1. Then, for a constant β < 1, TV( chlg m 1 , chlg m 2 ) ≤ β T s for any two initial memories m1, m2 containing safe placements.

Figure 8 .

 8 Figure 8. Loop perforation

  /x]} x ← e {Φ} [ASSIGN] {Φ ∧ e} c1 {Ψ} {Φ ∧ ¬e} c2 {Ψ} {Φ} if e then c1 else c2 {Ψ} [COND] {Φ ∧ e} c {Φ} Φ ⇓ while e do c {Φ} while e do c {Φ ∧ ¬e} [WHILE] {Φ} c1 {Ξ} {Ξ} c2 {Ψ} {Φ} c1; c2 {Ψ} [SEQ] {Φ } c {Ψ } |= Φ =⇒ Φ |= Ψ =⇒ Ψ {Φ} c {Ψ} [CONSEQ]

Figure 9 .

 9 Figure 9. Hoare Logic Rules for Deterministic Programs

1 x is the distribution where x has probability 1 and all other elements have probability 0.

2016/9/21

2016/9/21

In fact, the one-sided and two-sided rules are inter-derivable for all constructions except random assignments and loops.

2016/9/21

Kumar and Ramesh [39] show that the class of causal couplings-which contains all couplings in our logic-are unable to prove rapid mixing for some rapidly-mixing Markov chains[START_REF] Diaconis | Geometric Bounds for Eigenvalues of Markov Chains[END_REF][START_REF] Jerrum | Approximating the permanent[END_REF] 6 2016/9/21

2016/9/21

2016/9/21

for simplicity, the programs use an operator f which takes randomness as an argument (note that the value r is sampled immediately before the assignment using f ), although in the proof of the SHA3 hash function f is a procedure call whose body performs random samplings.

2016/9/21

2016/9/21


Amal Houidhek

Vincent Colotte

Zied Mnasri
email: zied.mnasri@enit.rnu.tn

Denis Jouvet
email: denis.jouvet@loria.fr

Imene Zangar

Statistical modelling of speech units in HMM-based speech synthesis for Arabic

Keywords: speech synthesis, statistical modeling, Arabic language, speech unit modeling

This paper investigates statistical parametric speech synthesis of Modern Standard Arabic (MSA). Hidden Markov Models (HMM)based speech synthesis system relies on a description of speech segments corresponding to phonemes, with a large set of features that represent phonetic, phonologic, linguistic and contextual aspects. When applied to MSA two specific phenomena have to be taken in account, the vowel lengthening and the consonant gemination. This paper studies thoroughly the modeling of these phenomena through various approaches: as for example, the use of different units for modeling short vs. long vowels and the use of different units for modeling simple vs. geminated consonants. These approaches are compared to another one which merges short and long variants of a vowel into a single unit and, simple and geminated variants of a consonant into a single unit (these characteristics being handled through the features associated to the sound). Results of subjective evaluation show that there is no significant difference between using the same unit for simple and geminated consonant (as well as for short and long vowels) and using different units for simple vs. geminated consonants (as well for short vs. long vowels).

Introduction

During the last decade, statistical parametric speech synthesis (SPSS) systems, such as HMM-based system called HTS [START_REF] Black | Statistical parametric speech synthesis[END_REF], have attracted a lot of interest compared to those based on unit selection [START_REF] Hunt | Unit selection in a concatenative speech synthesis system using a large speech database[END_REF]. Actually, HTS system is able to produce a smooth and good quality speech. Besides, it presents the advantages of being trainable, having a small footprint and changing voice characteristics. This approach has been already applied to many languages, such as Japanese [START_REF] Zen | Statistical parametric speech synthesis[END_REF], English [START_REF] Tokuda | An HMMbased speech synthesis system applied to English[END_REF] and French (Le [START_REF] Maguer | Evaluation of contextual descriptors for HMM-based speech synthesis in French[END_REF]. HTS system is split into two main parts corresponding to training and synthesis [START_REF] Zen | Statistical parametric speech synthesis[END_REF]. The training part consists of modeling spectrum features (e.g., Mel-cepstral coefficients and their dynamic features) and excitation features (including log (F0) and its dynamic features) from a speech database. In this stage, both linguistic and prosodic contexts are taken into account. Then, speech features are modeled by context-dependent HMMs. The synthesis part is a multi-step process: first an HMM representing the sentence is built by concatenating the context dependent phone HMMs, then, the HMMs' state durations are calculated through maximizing their probability, and speech level features (e.g., Mel-Cepstral coefficients, log(F0)…etc) are predicted. Finally, these features are processed by the synthesis filter, i.e., MLSA (Mel-Log Spectrum Approximation) to generate a waveform. In HTS, a speech unit is described with a set of contextual features at different levels: phoneme, syllable, word, phrase and utterance. These features are related to different aspects of a speech unit, whether linguistic, phonetic, phonologic or prosodic. These contextual features, which describe the characteristics of the speech segments associated to the phone units, are used during the HTS training process to build decision trees for parameter sharing. The relevance of the contextual features have a considerable effect on the HMM models and thus on synthesis quality. Part of the contextual features is language dependent and thus some descriptors may be added or neglected according to the language. In [START_REF] Tokuda | An HMMbased speech synthesis system applied to English[END_REF] a set of specific descriptors has been presented for English speech synthesis using HTS, whereas in (Le [START_REF] Maguer | Evaluation of contextual descriptors for HMM-based speech synthesis in French[END_REF] a different set of descriptors was introduced for French. However, previous studies focused only on evaluating the impact of a chosen subset of descriptors on acoustic parameter modeling, like in (Le [START_REF] Maguer | Evaluation of contextual descriptors for HMM-based speech synthesis in French[END_REF], and on modeling the duration of consonants and vowels (Silen et al., 2010).

The scope of this paper is to adapt the HTS system to the Arabic language through the introduction and the study of the relevance of some specific features. The originality of this work consists in evaluating the effects of many specific features combinations on the quality of HTS-synthesized Arabic speech. To extract relevant descriptors, a comprehensive review of Arabic language specifications (Al-Ani and Salman, 1970) was referred to.

In fact, MSA presents different specificities like stress (lexical stress), gemination and vowel lengthening. Long vowels are almost twice as long as short vowels and geminated consonants are twice as long as simple consonants [START_REF] Khouja | Durée des consonnes géminées en parole arabe : mesures et comparaison[END_REF]. Whereas stress was introduced in the original labels list of HTS [START_REF] Tokuda | An HMMbased speech synthesis system applied to English[END_REF] specific descriptors for gemination and vowel lengthening have not been so far considered in the standard set of descriptors for HTS. Thus, this paper expands and compares several choices of units that differentiate, or not, long vowels from short vowels, and/or geminated consonants from simple consonants. Using these different models of speech unit, Arabic utterances were produced for each combinations. Then, they were evaluated based on perceptive tests and on an objective evaluation of the duration of the phones.

The paper is organized as follows; section 2 presents different aspects of speech unit modelling in Arabic speech synthesis, particularly in SPSS, and then details the four proposed systems, developed for evaluating unit modelling in Arabic speech synthesis using HTS. Section 3 describes the evaluation protocol and data used in the experiments, then presents and discusses the results of the evaluation.

Speech-unit modeling for Arabic speech synthesis

Speech synthesis has been in continuous progress since the very early modern speech synthesizers [START_REF] Klatt | Software for a cascade/parallel formant synthesizer[END_REF] and (Taylor et al., 1995) until latest techniques, like SPSS. Depending on approaches and languages, different speech units are used. A speech unit corresponds to a phonetic or a phonological segment, which could be a phone, a diphone [START_REF] Moulines | A real-time French text-to-speech system generating high-quality synthetic speech[END_REF], a triphone or a syllable [START_REF] Kishore | Unit size in unit selection speech synthesis[END_REF]. The phoneme is considered as the smallest speech entity, any modification of a phoneme in a word may change its meaning.

According to [START_REF] Taylor | Text-to-speech synthesis[END_REF] a phone is defined as a realization of a phoneme in a certain context.

Arabic speech synthesis

Arabic speech has always been catching up with each novel TTS (Text To Speech) technique, starting from articulatory speech synthesis, formant-based, concatenative and more recently statistical parametric speech synthesis. In [START_REF] Rajouani | Synthesis-by-rule of Arabic language[END_REF], synthesis by rules is used to produce Arabic speech. Based on formants, i.e., maxima of the spectrogram, acoustic speech is generated through formants characteristics (i.e., bandwidth and amplitude) in addition to the rules of evolution of formants between phonemes. Concatenationbased methods consist in putting together a dynamically selected set of natural speech units to produce acoustic signals. Diphone units are speech signal segments going from the middle of one speech sound (phoneme) to the middle of the next one; diphone concatenation-based approaches lead to a better speech quality, compared to previous methods. In [START_REF] Baloul | Développement d'un système automatique de synthèse de la parole à partir du texte arabe standard voyellé[END_REF] diphone concatenation is used to produce Arabic speech. Further researches, used larger speech units in concatenationbased methods, such as di-syllable in [START_REF] Cheffour | [END_REF] and triphones in [START_REF] Ahmed | Réalisation d'un système hybride de synthèse de la parole Arabe utilisant un dictionnaire de polyphones[END_REF]. A particular approach of concatenation method, and maybe the most successful one in terms of quality and naturalness, is based on unit selection. It consists in selecting speech units of different sizes, to be concatenated, in order to produce speech. Applied to Arabic language in [START_REF] Abdelmalek | High quality Arabic Text-to-speech synthesis using unit selection[END_REF], it was based on phonemes and syllables.

Statistical parametric speech synthesis (SPSS)

HTS was adapted to many languages, and different set of descriptors have been used corresponding to the characteristics of each targeted language. For French (Le [START_REF] Maguer | Evaluation of contextual descriptors for HMM-based speech synthesis in French[END_REF], English [START_REF] Tokuda | An HMMbased speech synthesis system applied to English[END_REF] and German [START_REF] Krstulovic | An HMM-based speech synthesis system applied to German and its adaptation to a limited set of expressive football announcements[END_REF] the phoneme was the speech unit. However, Japanese [START_REF] Zen | An example of context-dependent label format for HMM-based speech synthesis in English[END_REF] Swedish [START_REF] Lundgren | An HMM-based Text-to-speech system applied to Swedish[END_REF] and Portuguese HTS-based speech synthesis systems were based on the syllable as a speech unit.

Previous works on HTS system to produce Arabic speech used phonemes as speech units. In [START_REF] Khalil | Arabic HMM-based speech synthesis[END_REF], the basic HTS system was adapted to MSA without explicitly considering the gemination or the vowel lengthening phenomena. [START_REF] Abdel-Hamid | Improving Arabic HMM based speech synthesis quality[END_REF] focused on some modification of the features and of the signal generation processing to improve the quality of produced speech.

Modern Standard Arabic phonology

MSA has 28 consonants (each consonant has a simple and a geminated version) and 3 vowels (they exist in both short and long version). Studies of Arabic phonology present specific characteristics of MSA [START_REF] Newman | The phonetics of Arabic[END_REF], in particular gemination and vowel lengthening.

Gemination ("shadda") [START_REF] Newman | The phonetics of Arabic[END_REF]: In spelling, it is represented by adding a diacritic sign ( ّ) above a consonant followed by a short vowel. Acoustically, it corresponds to making the duration of the geminated consonant twice as long as that of the simple consonant. A geminated consonant is phonetically and phonologically different from the same non-geminated one. This feature is considered as a distinctive factor, for example, a geminated /r/ changes the meaning of the word /darasa/ (he studied) to /darrasa/ (he taught).

Vowel lengthening ("madd") [START_REF] Selouani | Arabic phonetic features recognition using modular connectionist architectures[END_REF]: It consists of lengthening the duration of the vowels. Whereas short vowels are not usually written, long vowels are always indicated by the following graphemes و /uu/, ي /ii/, /اaa/. In addition, vowel lengthening changes the meaning of words, for example, if the initial vowel /a/ in the word /hatafa/ (he shouted) is lengthened, the resulting word /haatafa/ means (he telephoned).

Proposed unit modelling systems

To investigate the modeling of speech units in Arabic speech synthesis (normal vs. geminated consonant, short vowel vs. long vowel), a study of different choices of speech units is conducted. The aim is to investigate whether it is better to consider geminated consonants (resp. long vowels) as fully fledged phonemes (as implied by the Arabic phonology) or if both gemination and vowel lengthening can only be considered as specific linguistic features among the others? Below, four sets of units corresponding to four possible combinations are presented:

1) Differentiating simple vs. geminated consonants and short vs. long vowels (C*2__V*2): In this first system, a simple consonant (e.g., /d/) and its geminated counterpart (e.g., /dd/) are modelled with different units. Similarly, for vowels, a short vowel (e.g., /a/) and a long one (e.g., /aa/) are represented by two distinct units. This leads to the most expanded system of speech unit modelling relying on two units per consonant and two units per vowel.

2) Differentiating simple vs. geminated consonants and merging long and short vowels (C*2__V*1):

In this approach, a long vowel and its short counterpart are modelled with the same unit (for example the same model for /aa/ and /a/). On the other hand, the model of a geminated consonant is kept different from the model of the corresponding simple consonant, thus leading to two units per consonant and one unit per vowel.

3) Merging simple and geminated consonants and differentiating short vs. long vowels (C*1__V*2): The third system separates long vowels from short ones; so for example, the long vowel /aa/ and the short vowel /a/ are modelled by two different units. However, simple consonants and their geminated counterparts are modelled by the same units. This leads to one unit per consonant and two units per vowel.

4) Merging simple and geminated consonant and merging short and long vowels (C*1__V*1):

It is the most compact system. Here a simple consonant (e.g., /d/) and the corresponding geminated consonant (/dd/) are represented by the same unit. The same for vowels, a long vowel (e.g., /aa/) and the corresponding short one (/a/) are modelled by the same unit. This leads to one unit per consonant and one unit per vowel.

In addition, linguistic features are computed over the phoneme sequence of each stimuli in the training corpus. The set of contextual features for defining the phone labels is inspired from the one defined for English [START_REF] Tokuda | An HMMbased speech synthesis system applied to English[END_REF]. Each segment is described with a set of features corresponding to different levels (phoneme, syllable, word, phrase, utterance). However, features related to the accent information and to TOBI [START_REF] Silverman | Tobi: A standard for labeling English prosody[END_REF] are not employed in the proposed set. Indeed, the novelty is adding two additional features to take into account specificities of the Arabic language. Whereas the first feature is used to indicate the possible geminated characteristic (possible values are simple consonant, geminated consonant, or not a consonant), the second one refers to vowel lengthening (possible values are short vowel, long vowel, or not a vowel). Finally, it should be noted that the four systems use the same set of descriptors.

Evaluation

Data

For the development and the evaluation of the systems described above, a Modern Standard Arabic corpus is used [START_REF] Halabi | Phonetic inventory for an Arabic speech corpus[END_REF]. It consists of 1595 separate utterances recorded from a male-speaker reading news bulletin, in a neutral style. The signal is sampled at 48 kHz. For speech synthesis, a speaker-dependent modelling is used, and one HTS model is developed (trained) for each of the four choices of units. For each system, 1565 Arabic utterances from the corpus are used for training the model parameters and 30 utterances are kept for evaluation purpose (test set). Signals are produced for each system using HTS with STRAIGHT vocoder [START_REF] Kawahara | Restructuring speech representations using a pitch-adaptive time frequency smoothing and an instantaneous-frequency-based F0 extraction: Possible role of a repetitive structure in sounds 1[END_REF].

For evaluating and comparing the different approaches, the 30 utterances of the test set were synthesized with each system, thus providing 30 speech stimuli for each system.

Results and discussions

In this experiment, seven native Arabic listeners, who neither are speech specialists nor involved in speech synthesis, participated in the tests.

Global quality and naturalness evaluation

Listeners were asked to assess two features. The first one is the overall quality which refers to the general quality of produced signals and if it was easy to listen to. The second one is the naturalness of the synthesized speech by focusing on the evaluation of the intonation (whether the pitch's evolution is natural) and of the rhythm, (whether the length of phonemes sounds natural too). Signals produced with the four proposed systems are assessed by giving a score on a scale from 1 to 5 ranging from very bad to excellent. Fig. 1 and Fig. 2 display the mean values of MOS listening tests with 95% confidence intervals, with respect to quality and naturalness. According to the MOS scores, the four models lead to rather similar quality of speech signals. The similarity between quality and naturalness results can be explained by the fact that listeners, who neither are specialists in phonetics nor accustomed to speech evaluation, had not seen the difference between the global quality and naturalness questions even though each question was preceded with an introductory explanation. Each pair consists of the same utterance produced with two different systems. During the test, the order of presenting the speech signals is randomly chosen for each trial. Listeners were asked to compare the second signal to the first one and to give a score from 1 to 7 ranging from much worse to much better.

Comparison test

The one-to-one comparison of the four systems in Fig. 3 shows that there is no clear preference for a particular system.

Evaluation of duration

For each system, the average, over the vowels, of the ratios between the mean duration of long vowels (LV) and the mean duration of corresponding short vowels (SV) is calculated as well as the average ratio for geminated consonants (GC) vs. simple consonants (SC). These average ratios are then compared to those obtained for natural speech. Only phonemes with a number of occurrences greater or equal to ten are considered. 1 shows that ratios of long vowel durations to short vowel durations, as well as ratios of geminated consonant durations to simple consonant durations, computed for the different synthesis systems are similar to the ratios calculated on natural speech.

Number of occurrences

To investigate the reasons of the similarity of ratios values of the four models, the duration decision trees were analysed for the four models. Note that there is only a single tree for all phonemes in each model. Fig. 4, which displays the top of the decision tree, shows that questions about the length of the current segment (i.e., geminated consonant or not: "C-Cur-Gem-Cons" and long vowel or not: "C-Seg ==vl") are placed on the top of the tree; this is observed for the four models.

Note that, "C-"refers to the current phoneme, "L-"to the left phoneme and "sil" to silence. Fig. 5 displays the values of RMSE for four phoneme classes (short vowels, long vowels, simple consonants and geminated consonants). Results show that for each class, there is no important difference between RMSE values of the four models. Besides, when considering the mean duration values of each class, NRMSE (normalized root mean square error) values are similar, around 25% for geminated consonant and 35% for the 3 other classes.

Conclusion

In this paper, HMM based speech synthesis (HTS) was applied to Modern Standard Arabic. In Arabic language, two specific phenomena need to be considered: the consonant gemination and the vowel lengthening. Thus, two new features were included in the HTS descriptors to take into account these specificities.

Moreover, several modelling approaches have been investigated with respect to the choice of modelling units, as for example, the use of different units for modelling long vs. short vowels, and/or the use of different units for modelling simple vs. geminated consonants. These approaches were compared to another one, which merges short and long variants of a vowel into a single unit, and simple and geminated variants of a consonant into a single unit. Subjective listening evaluation tests (MOS and preference tests) showed that there is no significant difference between the various modelling systems.

When using HMM-based speech synthesis system for MSA, an identification of geminated consonants and long vowels as fully-fledged phonemes is not necessary, as long as this information exists in the set of the specific phonetic feature.

As DNN are now being used for speech synthesis, it will be interesting to study if, unlike for HMM speech synthesis, neural networks would benefit from the explicit differentiation of geminated vs. simple consonants and long vs. simple vowels.

Acknowledgements

This research work was conducted under PHC-Utique Program in the framework of CMCU (Comité Mixte de Coopération Universitaire) grant N°15G1405.

Fig. 2 .

 2 Fig.1 . Global quality results with 95% confidence intervals.

  Speech signals produced with the four approaches of speech unit modeling A (C*2__V*2), B (C*1__V*1), C (C*1__V*2) and D (C*2__V*1) are compared to each other.

Fig. 3 .

 3 Fig. 3. Results of comparisons tests (percentage of preference)

Fig. 4 .Fig. 5 .

 45 Fig. 4. Duration decision tree of (C*2__V*2) model (top part).