Gilles Barthe

Thomas Espitau

Luis María Ferrer Fioriti

Justin Hsu

Synthesizing Probabilistic Invariants via Doob's Decomposition

When analyzing probabilistic computations, a powerful approach is to first find a martingale-an expression on the program variables whose expectation remains invariant-and then apply the optional stopping theorem in order to infer properties at termination time. One of the main challenges, then, is to systematically find martingales. We propose a novel procedure to synthesize martingale expressions from an arbitrary initial expression. Contrary to state-of-the-art approaches, we do not rely on constraint solving. Instead, we use a symbolic construction based on Doob's decomposition. This procedure can produce very complex martingales, expressed in terms of conditional expectations. We show how to automatically generate and simplify these martingales, as well as how to apply the optional stopping theorem to infer properties at termination time. This last step typically involves some simplification steps, and is usually done manually in current approaches. We implement our techniques in a prototype tool and demonstrate our process on several classical examples. Some of them go beyond the capability of current semi-automatic approaches.

Introduction

Probabilistic computations are a key tool in modern computer science. They are ubiquitous in machine learning, privacy-preserving data mining, cryptography, and many other fields. They are also a common and flexible tool to model a broad range of complex real-world systems. Not surprisingly, probabilistic computations have been extensively studied from a formal verification perspective. However, their verification is particularly challenging.

In order to understand the difficulty, consider the standard way to infer properties about the final state of a non-probabilistic program using a strong invariant (an assertion which is preserved throughout program execution) and a proof of termination. This proof principle is not easily adapted to the probabilistic case. First, probabilistic programs are interpreted as distribution transformers [START_REF] Kozen | Semantics of probabilistic programs[END_REF] rather than state transformers. Accordingly, assertions (including strong invariants) must be interpreted over distributions. Second, the notion of termination is different for probabilistic programs. We are usually not interested in proving that all executions are finite, but merely that the probability of termination is 1, a slightly weaker notion. Under this notion, there may be no finite bound on the number of steps over all possible executions. So, we cannot use induction to transfer local properties to the end of the program-more complex limiting arguments are needed.

We can avoid some of these obstacles by looking at the average behavior of a program. That is, we can analyze numerical expressions (over program variables) whose average value is preserved. These expressions are known as martingales, and have several technical advantages. First, martingales are easy to manipulate symbolically and can be checked locally. Second, the average value of martingale is preserved at termination, even if the control-flow of the program is probabilistic. This fact follows from the optional stopping theorem (OST), a powerful result in martingale theory.

While martingales are quite useful, they can be quite non-obvious. Accordingly, recent investigation has turned to automatically synthesizing martingales. Stateof-the-art frameworks are based on constraint solving, and require the user to provide either a template expression [START_REF] Katoen | Linear-invariant generation for probabilistic programs: Automated support for proof-based methods[END_REF][START_REF] Chakarov | Probabilistic program analysis with martingales[END_REF] or a limit on the search space [START_REF] Chen | Counterexample-guided polynomial loop invariant generation by lagrange interpolation[END_REF][START_REF] Chakarov | Expectation invariants for probabilistic program loops as fixed points[END_REF]. The main advantage of such approaches is that they are generally completethey find all possible martingales in the search space. However, they have their drawbacks: a slightly wrong template can produce no invariant at all, and a lot of search space may be needed to arrive at the martingale.

We propose a framework that complements current approaches-we rely on purely symbolic methods instead of solving constraints or searching. We require the user to provide a "seed" expression, from which we always generate a martingale. Our approach uses Doob's decomposition theorem, which gives a symbolic method to construct a martingale from any sequence of random values. Once we have the martingale, we can apply optional stopping to reason about the martingale at loop termination. While the martingale and final fact may be quite complex, we can use non-probabilistic invariants and symbolic manipulations to automatically simplify them.

We demonstrate our techniques in a prototype implementation, implementing Doob's decomposition and the Optional Stopping Theorem. Although these proof principles have been long known to probability theory, we are the first to incorporate them into an automated program analysis. Given basic invariants and hints, our prototype generates martingales and facts for a selection of examples.

Mathematical preliminaries

We briefly introduce some definitions from probability theory required for our technical development. We lack the space to discuss the definitions in-depth, but we will explain informally what the various concepts mean in our setting. Interested readers can find a more detailed presentation in any standard probability theory textbook (e.g., Williams [START_REF] Williams | Probability with Martingales[END_REF]).

First, we will need some basic concepts from probability theory.

Definition 1 Let Ω be the set of outcomes.

-A sigma algebra is a set F of subsets of Ω, closed under complements and countable unions, and countable intersections.

-A probability measure is a countably additive mapping P : F → [0, 1] such that P(Ω) = 1. -A probability space is a triple (Ω, F, P).

Next, we can formally define stochastic processes. These constructions are technical but completely standard.

Definition 2 Let (Ω, F, P) be a probability space.

-A (real) random variable is a function

X : Ω → R. X is F-measurable if X -1 ((a, b]) ∈ F for every a, b ∈ R. -A filtration is a sequence {F i } i∈N of sigma algebras such that F i ⊆ F and F i-1 ⊆ F i for every i > 0.
When there is a fixed filtration on F, we will often abuse notation and write F for the filtration. -A stochastic process adapted to filtration F is a sequence of random variables

{X i } i∈N such that each X i is F i -measurable.
Intuitively, we can think of Ω as a set where each element represents a possible outcome of the samples. In our setting, grouping samples according to the loop iteration gives a natural choice for the filtration: we can take F i to be the set of events that are defined by samples in iteration i or before. A stochastic process X is adapted to this filtration if X i is defined in terms of samples from iteration i or before. Sampled variables at step i are independent of previous steps, so they are not F i-1 -measurable.

Expectation. To define martingales, we need to introduce expected values and conditional expectations. The expected value of a random variable is defined as

E[X] Ω X • dP
where is the Lebesgue integral [START_REF] Williams | Probability with Martingales[END_REF]. We say that a random variable is integrable if E[|X|] is finite. Given a integrable random variable X and a sigma algebra G, a conditional expectation of X with respect to G is a random variable Y such that Y is G-measurable, and

E[X • 1 {A} ] = E[Y • 1 {A} ]
for all events A ∈ G. (Recall that the indicator function 1 {A} of an event A maps ω ∈ A to 1, and all other elements to 0.) Since one can show that this Y is essentially unique, we denote it by E[X | G].

Moments. Our method relies on computing higher-order moments. Suppose X is a random variable with distribution d. If X takes numeric values, the kth moment of d is defined as

G(d) k E[X k ]
for k ∈ N. If X ranges over tuples, the correlations of d are defined as

G(d, {a, b}) p,q E[π a (X) p • π b (X) q ],
for p, q ∈ N, and similarly for products of three or more projections. Here, the projection π i (X) for X a tuple-valued random variable is the marginal distribution of the ith coordinate of the tuple.

Martingales. A martingale is a stochastic process with a special property: the average value of the current step is equal to the value of the previous step.

Definition 3 Let {X i } be a stochastic process adapted to filtration {F i }. We say that X is a martingale with respect to F if it satisfies the property

E[X i | F i-1 ] = X i-1 .
For a simple example, consider a symmetric random walk on the integers. Let X ∈ Z denote the current position of the walk. At each step, we flip a fair coin: if heads, we increase the position by 1, otherwise we decrease the position by 1. The sequence of positions X 0 , X 1 , . . . forms a martingale since the average position at time i is simply the position at time i -1:

E[X i | F i-1 ] = X i-1 .
Doob's decomposition. One important result in martingale theory is Doob's decomposition. Informally, it establishes that any integrable random process can be written uniquely as a sum of a martingale and a predictable process. For our purposes, it gives a constructive and purely symbolic method to extract a martingale from any arbitrary random process.

Theorem 1 (Doob's decomposition) Let X = {X i } i∈N be a stochastic process adapted to filtration {F i } i∈N where each X i has finite expected value. Then, the following process is a martingale:

M i = X 0 : i = 0 X 0 + i j=1 X j -E[X j | F j-1 ] : i > 0 If X is already a martingale, then M = X.
We will think of the stochastic process X as a seed process which generates the martingale. While the definition of the martingale involves conditional expectations, we will soon see how to automatically simplify these expectations.

Optional stopping theorem. For any martingale M , it is not hard to show that the expected value of M remains invariant at each time step. That is, for any fixed value n ∈ N, we have

E[M n ] = E[M 0 ].
The optional stopping theorem extends this equality to situations where n itself may be random, possibly even a function of the martingale.

Definition 4 Let (Ω, F, P) be a probability space with filtration {F i } i∈N . A random variable τ :

Ω → N is a stopping time if the subset {w ∈ Ω | τ (w) ≤ i} is a member of F i for each i ∈ N.
Returning to our random walk example, the first time that the position is farther than 100 from the origin is a stopping time since this time depends only on past samples. In contrast, the last time that a position is farther than 100 from the origin is not a stopping time, since this time depends on future samples. More generally, the iteration count when we exit a probabilistic loop is a stopping time since termination is a function of past samples only.

If we have a stopping time and a few mild conditions, we can apply the optional stopping theorem. 1Theorem 2 (Optional stopping) Let τ be a stopping time, and let M be a martingale. If the expected value of τ is finite, and if

|M i -M i-1 | ≤ C for all i > 0 and some constant C, then E[M τ ] = E[M 0 ].
To see this theorem in action, consider the random walk martingale S and take the stopping time τ to be the first time that |S| ≥ 100. It is possible to show that τ has finite expected value, and clearly

|S i -S i-1 | ≤ 1. So, the optional stopping theorem gives E[S τ ] = E[S 0 ] = 0.
Since we know that the position is ±100 at time τ , this immediately shows that the probability of hitting +100 is equal to the probability of hitting -100. This intuitive fact can be awkward to prove using standard probabilistic invariants, but falls out directly from a martingale analysis.

Overview of method

Now that we have seen the key technical ingredients of our approach, let us see how to combine these tools into an automated program analysis. We will take an imperative program specifying a stochastic process and a seed expression, and we will automatically synthesize a martingale and an assertion that holds at termination. We proceed in three stages: extracting a polynomial representing the stochastic process in the program, applying Doob's decomposition to the polynomial representation, and applying optional stopping. We perform symbolic manipulations to simplify the martingale and final fact.

Programs. We consider programs of the form:2 I; while e do (S; B)

where I and B are sequences of deterministic assignments (of the form X ← E), and S is a sequence of probabilistic samplings (of the form S $ ← DE). Note that we separate sample variables s ∈ S, which are the target of random samplings, from process variables x ∈ X , which are the target of deterministic assignments. This distinction will be important for our simplifications: we know the moments and correlations of sample variables, while we have less information for process variables. We require that programs assign to sample variables before assigning to process variables in each loop iteration; this restriction is essentially without loss of generality.

We take DE to be a set of standard distributions over the integers or over finite tuples of integers, to model joint distributions. For instance, we often consider the distribution Bern(1/2, {-1, 1}) ∈ DE that returns -1 and +1 with equal probability. We assume that all distributions in DE have bounded support; all moments and correlations of the primitive distributions are finite. We will also assume that distributions do not depend on the program state.

The set E of expressions is mostly standard, with a few notational conveniences for defining stochastic processes:

E ::= X | S | X [-n] process/sample/history variables | Z constants | π a (E) projections | E + E | E • E arithmetic | E < E | E ∧ E | ¬E guards
History variables X [-n] are indexed by a positive integer n and are used inside loops. The variable x[-n] refers to the value of x assigned n iterations before the current iteration. If the loop has run for fewer than n iterations, x[-n] is specified by the initialization step of our programs:

I x 1 [-n 1 ] ← e 1 ; • • • ; x k [-n k ] ← e k .
Extracting the polynomial. For programs in our fragment, each variable assigned in the loop determines a stochastic process: x is the most recent value, x[-1] is the previous value, etc. In the first stage of our analysis, we extract polynomial representations of each stochastic process from the input program.

We focus on the variables that are mutated in B-each of these variables determines a stochastic process. To keep the notation light, we will explain our process for first-order stochastic processes: we only use history variables x[-1] from the past iteration. We will also suppose that there is just one process variable and one sample variable, and only samples from the current iteration are used.

Since our expression language only has addition and multiplication as operators, we can represent the program variable x as a polynomial of other program variables:

x = P x (x[-1], s) (1)
Next, we pass to a symbolic representation in terms of (mathematical) random variables. To variable x, we associate the random variable {X i } i∈N modeling the corresponding stochastic process, and likewise for the sample variable s. By convention, i = 0 corresponds to the initialization step, and i > 0 corresponds to the stochastic process during the loop. In other words,

x[0] ← 0; while e do s $ ← d; x ← x[-1] + s desugars to x[0] ← 0; i ← 0; while e do i ← i + 1; s $ ← d; x[i] ← x[i -1] + s
in a language with arrays instead of history variables. Then, the program variable x[i] corresponds to the random variable X i . Then, Equation ( 1) and the initial conditions specified by the command I give an inductive definition for the stochastic process:

X i = P x (X i-1 , S i ).
(

) 2
Applying Doob's decomposition. The second stage of our analysis performs Doob's decomposition on the symbolic representation of the process. We know that the seed expression e must be a polynomial, so we can form the associated stochastic process {E i } i∈N by replacing program variables by their associated random variable:

E i = P e (X i , S i ). (3)
(Recall that the initial conditions X 0 and S 0 , which define E 0 , are specified by the initialization portion I of the program.) Then, Doob's decomposition produces the martingale:

M i = E 0 : i = 0 E 0 + i j=1 E j -E[E j | F j-1 ] : i > 0.
To simplify the conditional expectation, we unfold E j via Equation (3) and unroll the processes X i by one step with Equation (2). Now, we apply our simplification rules; we present a selection in Figure 1. The rules are divided into three groups (from top): linearity of expectation, conditional independence, and distribution information. The first two groups reflect basic facts about expectations and conditional expectations. The last group encodes the moments and correlations of the primitive distributions. We can pre-compute these quantities for each primitive distribution d and store the results in a table.

By the form of Equation ( 3), the simplification removes all expectations and we can give an explicit definition for the martingale:

M i = E 0 : i = 0 Q e (X i-1 , . . . , X 0 , S i , . . . , S 1 ) : i > 0, (4)
where Q e is a polynomial.

E[c • f + c • g | -] → c • E[f | -] + c • E[g | -] E[Xi-n • f | Fi-1] → Xi-n • E[f | Fi-1] (n > 0) E[Si • S i | Fi-1] → E[Si | Fi-1] • E[S i | Fi-1] (S = S ) E[S k i | Fi-1] → G(d) k (S ∼ d) E[πa(Si) p • π b (Si) q | Fi-1] → G(d a,b )p,q (S ∼ d)
Fig. 1: Selection of simplification rules Applying optional stopping. With the martingale in hand, the last stage of our analysis applies the optional stopping theorem. To meet the technical conditions of the theorem, we need two properties of the loop:

-The expected number of iterations must be finite.

-The martingale must have bounded increments.

These side conditions are non-probabilistic assertions that can already be handled using existing techniques. For instance, the first condition follows from the existence of a bounded variant [START_REF] Hart | Termination of probabilistic concurrent program[END_REF]: an integer expression v such that -0 ≤ v < K; v = 0 implies the guard is false; and the probability that v decreases is strictly bigger than throughout the loop, for and K positive constants. However in general, finding a bounded variant may be difficult; proving finite expected stopping time is an open area of research which we do not address here.

The second condition is also easy to check. For one possible approach, one can replace stochastic sampling by non-deterministic choice over the support of the distribution, and verify that the seed expression e is bounded using standard techniques [START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF][START_REF] Miné | The octagon abstract domain[END_REF][START_REF] De Moura | Z3: an efficient SMT solver[END_REF]. This suffices to show that the martingale M i has bounded increments. To see why, suppose that the seed expression is always bounded by a constant C. By Doob's decomposition, we have

|M i -M i-1 | =   i j=1 E j -E[E j | F j-1 ]   -   i-1 j=1 E j -E[E j | F j-1 ]   = |E i -E[E i | F j-1 ]| ≤ 2C,
so the martingale has bounded increments.

Thus, we can apply the optional stopping theorem to Equation (4) to conclude: Unlike the simplification step after applying Doob's decomposition, we may not be able to eliminate all expected values. For instance, there may be expected values of X at times before τ . However, if we have additional invariants about the loop, we can often simplify the fact with basic symbolic transformations.

E 0 = E[M 0 ] = E[M τ ] = E[Q e (X τ
Implementation. We have implemented our process in a Python prototype using the sympy library for handling polynomial and summation manipulations [START_REF] Joyner | Open source computer algebra systems: SymPy[END_REF].

Figure 2 shows the entire pipeline. There are three parts of the input: the program describing a stochastic process, the seed expression, and hint facts. The output is a probabilistic formula that holds at termination.

The most challenging part of our analysis is the last stage: applying OST. First, we need to meet the side conditions of the optional stopping theorem: finite expected iteration count and bounded increments. Our prototype does not verify these side conditions automatically since available termination tools are either not fully automatic [START_REF] Esparza | Proving termination of probabilistic programs using patterns[END_REF] or can only synthesize linear ranking supermartingales [START_REF] Chakarov | Probabilistic program analysis with martingales[END_REF][START_REF] Chatterjee | Algorithmic analysis of qualitative and quantitative termination problems for affine probabilistic programs[END_REF] that are insufficient for the majority of our case studies 3 . Furthermore, the final fact typically cannot be simplified without some basic information about the program state at loop termination. We include this information as hints. Hints are first-order formulae over the program variables and the loop counter (represented by the special variable t), and are used as auxiliary facts during the final simplification step. Hints can be verified using standard program verification tools since they are non-probabilistic. In our examples, we manually translate the hints and the program into the input language of EasyCrypt [START_REF] Barthe | Computer-aided security proofs for the working cryptographer[END_REF], and perform the verifications there. Automating the translation poses no technical difficulty and is left for future work.

We note that the performance of the tool is perfectly reasonable for the examples considered in the next section. For instance, it handles the "ABRA-CADABRA" example in less than 2 seconds on a modern laptop. Now, we demonstrate our approach on several classic examples of stochastic processes. In each case, we describe the code, the seed expression, and any hints needed for our tool to automatically derive the final simplified fact.

Geometric distribution. As a first application, we consider a program that generates a draw from the geometric distribution by running a sequence of coin flips.

x[0] ← 0; while (z = 0) do z $ ← Bern(p, {1, 0}); x ← x[-1] + z; end
Here, Bern(p, {1, 0}) is the distribution that returns 1 with probability p > 0, and 0 otherwise. The program simply keeps drawing until we sample 0, storing the number of times we sample 1 in x.

We wish to apply our technique to the seed expression x. First, we can the polynomial equation:

X i = X i-1 + Z i .
Applying Doob's decomposition, our tool constructs and reduces the martingale:

M i = X 0 : i = 0 X i -p • i : i > 0.
To apply optional stopping, we first need to check that the stopping time τ is integrable. This follows by taking z as a bounded variant-it remains in {0, 1} and decreases with with probability p > 0. Also, the martingale M i has bounded increments: |M i -M i-1 | should be bounded by a constant. But this is clear since we can use a loop invariant to show that |X i -X i-1 | ≤ 1, and the increment is

|M i -M i-1 | = |X i -X i-1 -p| ≤ |X i -X i-1 | + p ≤ 1 + p.
So, optional stopping shows that

0 = E[X τ -p • τ ].
With the hint x = t -1-which holds at termination-our tool replaces X τ by τ -1 and automatically derives the expected running time:

0 = E[τ -1 -p • τ ] E[τ ] = 1/(1 -p).
Gambler's ruin. Our second example is the classic Gambler's ruin process. The program models a game where a player starts with a > 0 dollars and keeps tossing a fair coin. The player wins one dollar for each head and loses one dollar for each tail. The game ends either when the player runs out of money, or reaches target of b > a dollars. We can encode this process as follows:

x[0] ← a; while (0 < x < b) do z $ ← Bern(1/2, {-1, 1}); x ← x + z; end
We will synthesize two different martingales from this program, which will yield complementary information once we apply optional stopping. For our first martingale, we use x as the seed expression. Our tool synthesizes the martingale

M i = X 0 : i = 0 X i : i > 0.
So in fact, x is already a martingale.

To apply optional stopping, we first note that x is a bounded variant: it remains in (0, b) and decreases with probability 1/2 at each iteration. Since the seed expression x is bounded, the martingale M i has bounded increments. Thus, optional stopping yields

a = E[X 0 ] = E[X τ ].
If we give the hint

(x = 0) ∨ (x = b) (5)
at termination, our prototype automatically derives

a = E[X τ • 1 {Xτ =0∨Xτ =b} ] = E[X τ • 1 {Xτ =0} ] + E[X τ • 1 {Xτ =b} ] = 0 • Pr[X τ = 0] + b • Pr[X τ = b] = b • Pr[X τ = b],
so the probability of exiting at b is exactly a/b. Now, let us take a look at a different martingale generated by the seed expression x 2 . Our prototype synthesizes the following martingale:

M i = X 2 0 : i = 0 X 2 i -i : i > 0
Again, we can apply optional stopping: x is a bounded variant, and the seed expression x 2 remains bounded in (0, b 2 ). So, we get

a 2 = E[M 0 ] = E[X 2 τ -τ ].
By using the same hint Equation ( 5), our prototype automatically derives

a 2 = E[X 2 τ • 1 {Xτ =0∨Xτ =b} ] -E[τ ] = E[X 2 τ • 1 {Xτ =0} ] + E[X 2 τ • 1 {Xτ =b} ] -E[τ ] = 0 • Pr[X τ = 0] + b 2 • Pr[X = b] -E[τ ] = b 2 • Pr[X τ = b] -E[τ ].
Since we already know that Pr[X τ = b] = a/b from the first martingale {M i } i∈N , this implies that the expected running time of the Gambler's ruin process is

E[τ ] = a(b -a).
Gambler's ruin with momentum. Our techniques extend naturally to stochastic processes that depend on variables beyond the previous iteration. To demonstrate, we'll consider a variant of Gambler's ruin process with momentum: besides just the coin flip, the gambler will also gain profit equal to the difference between the previous two dollar amounts. Concretely, we consider the following process:

x[0] ← a; x[1] ← a; while (0 < x < b) do z $ ← Bern(p, {-1, 1}); x ← x[-1] + (x[-1] -x[-2]) + z;
end Note that we must now provide the initial conditions for two steps, since the process is second-order recurrent. Given seed expression x, our tool synthesizes the following martingale:

M i = X 0 : i = 0 X 0 + X i -X i-1 : i > 0
Identical to the Gambler's ruin process, we can verify the side conditions and apply optional stopping, yielding

a = E[M 0 ] = E[M τ ] = E[X 0 + X τ -X τ -1 ].
Unfolding X 0 = a and simplifying, our tool derives the fact

E[X τ ] = E[X τ -1 ].
We are not aware of existing techniques that can prove this kind of fact-reasoning about the expected value of a variable in the iteration just prior to termination. Abracadabra. Our final example is a classic application of martingale reasoning. In this process, a monkey randomly selects a character at each time step, stopping when he has typed a special string, say "ABRACADABRA". We model this process as follows:

match0 [0] ← 1; match1 [0] ← 0; ... match11 [0] ← 0; while (match11 == 0) do s $ ← UnifMatches; match11 ← match10 [-1] * π11(s); match10 ← match9 [-1] * π10(s); ... match1 ← match0 [-1] * π1(s); end
Here, UnifMatches is a distribution over tuples that represents a uniform c draw from the letters, where the kth entry is 1 if the c matches the kth word and 0 if not. The variables match i record whether the i most recent letters match the first i letters of target word; match 0 is always 1, since we always match 0 letters. Now, we will apply Doob's decomposition. Letting L be the number of possible letters and taking the seed expression

e 1 + L • match 1 + • • • + L 11 • match 11 ,
our tool synthesizes the following martingale:

M i = 1 + L • X (1) i
+ • • • + L 11 • X (11) i
: i = 0 i j=1 L • X (1)
j + • • • + L 11 • X (11) j -L 0 • X (0) j-1 + • • • + L 10 • X (10) j-1
: i > 0,
where X (j) is the stochastic process corresponding to match j . The dependence on L is from the expectations of projections of UnifMatches, which are each 1/L-the probability of a uniformly random letter matching any fixed letter.

To apply the optional stopping theorem, note that the seed expression e is bounded in (0, L 12 ), and 1 + L + • • • + L 11 -e serves as a bounded variant: take the highest index j such that match j = 1, and there is probability 1/L that we increase the match to get match j+1 = 1, decreasing the variant. So, we have

1 = E[M 0 ] = τ j=1 E[L • X (1) j + • • • + L 11 • X (11) j ] -E[L 0 • X (0) j-1 + • • • + L 10 • X (10) j-1 ].
Our tool simplifies and uses the hints X (11) j = 0 and X (0) j = 1 for j < τ to derive

1 = L 0 • E[X (0) τ ] + • • • + L 11 • E[X (11) τ ] -E[τ ].
For the target string "ABRACADABRA", we use hints For example, if match 11 is set then the full string is matching "ABRACADABRA", so the most recently seen four characters are "ABRA". This matches the first four letters of "ABRACADABRA", so match 4 is also set. The hint can be proved from a standard loop invariant.

(
Our tool derives the expected running time:

E[τ ] = L 1 + L 4 + L 11 .
Benchmarks.

To give an idea of the efficiency of our tool, we present some benchmarks for our examples in Table 1. We measured timing on a recent laptop with a 2.4 GHz Intel Core processor with 8 GB of RAM. We did not optimize for the performance; we expect that the running time could be greatly improved with some tuning. The example miniabra is a smaller version of the ABRACADABRA example, where the alphabet is just {0, 1}, and we stop when we sample the sequence 111; fullabra is the full ABRACADABRA example.

While there is a growing body of work related to martingale techniques for program analysis (see the following section), it is not obvious how to compare benchmarks. Existing work focuses on searching for martingale expressions within some search space; this is a rather different challenge than synthesizing a single martingale from a programmer-provided seed expression. In particular, if the seed expression happens to already be a martingale by some lucky guess, our tool will simply return the seed expression after checking that it is indeed a martingale.

Related work

Martingales. Martingale theory is a central tool in fields like statistics, applied mathematics, control theory, and finance. When analyzing randomized algorithms, martingales can show tight bounds on tail events [START_REF] Motwani | Randomized Algorithms[END_REF]. In the verification community, martingales are used as invariant arguments, and as variants arguments to prove almost sure termination [START_REF] Bournez | Proving positive almost-sure termination[END_REF][START_REF] Chakarov | Probabilistic program analysis with martingales[END_REF][START_REF] Ferrer Fioriti | Probabilistic termination: Soundness, completeness, and compositionality[END_REF][START_REF] Chatterjee | Algorithmic analysis of qualitative and quantitative termination problems for affine probabilistic programs[END_REF]. Recently, martingale approaches were extended to prove more complex properties. Chakarov, Voronin, and Sankaranarayanan [START_REF] Chakarov | Deductive proof for almost sure persistence and recurrence properties[END_REF] propose proof rules for proving persistence and recurrence properties. Dimitrova, Ferrer Fioriti, Hermanns, and Majumdar [START_REF] Dimitrova | Probabilistic CTL * : The deductive way[END_REF] develop a deductive proof system for PCTL * , with proof rules based on martingales and supermartingales.

Probabilistic Hoare logic. McIver and Morgan [START_REF] Mciver | Abstraction, Refinement and Proof for Probabilistic Systems[END_REF] propose a Hoare-like logic that is quite similar to our approach of using martingales and OST. Their approach is based on weakest pre-expectations, which are an extension of Dijkstra's weakest preconditions [START_REF] Dijkstra | Guarded commands, nondeterminacy and formal derivation of programs[END_REF] based on "backward" conditional expectations. Their probabilistic invariants are similar to submartingales, as the expected value of the invariant at the beginning of the execution lower bounds the expected value of the invariant at termination. Their proof rule also requires an additional constraint to ensure soundness, but it requires a limiting argument that is more difficult to automate compared to our bounded increment condition. We could relax our condition using a weaker version of OST that generalizes their condition [START_REF] Williams | Probability with Martingales[END_REF]. Another substantial difference with our approach is that their logic supports non-deterministic choices-ours does not. It is not obvious how we can extend our synthesis approach to the non-probabilistic case as we heavily rely on the concept of filtration, not applicable in the presence of non-determinism.

Probabilistic model checking. In the last twenty years, model checking technology for probabilistic models have made impressive strides [START_REF] Kwiatkowska | PRISM 4.0: Verification of probabilistic real-time systems[END_REF][START_REF] Ciesinski | LiQuor: A tool for qualitative and quantitative linear time analysis of reactive systems[END_REF][START_REF] Katoen | The ins and outs of the probabilistic model checker MRMC[END_REF] (Baier and Katoen [START_REF] Baier | Principles of model checking[END_REF] provide overview). The main advantage of model checking is that it requires nothing from the user; our technique requires a seed expression. However, model checking techniques suffer from the state explosion problem-the time and memory consumption of the model checking algorithm depends on the number of reachable states of the program. Our approach can be used to verify infinite and parametric programs without any performance penalty, as we work purely symbolically. For example, a probabilistic model checker can find the expected running time of the gambler's ruin process for concrete values of a and b but they cannot deduce the solution for the general case, unlike our technique.

Invariant synthesis. There are several approaches for synthesizing probabilistic invariants. Katoen et al. [START_REF] Katoen | Linear-invariant generation for probabilistic programs: Automated support for proof-based methods[END_REF] propose the first complete method for the synthesis of McIver and Morgan's probabilistic linear invariants. It is an extension of the constraint solving approach by Colón, Sankaranarayanan, and Sipma [START_REF] Colón | Linear invariant generation using non-linear constraint solving[END_REF] for the synthesis of (non-probabilistic) linear invariants. Chakarov and Sankaranarayanan [START_REF] Chakarov | Probabilistic program analysis with martingales[END_REF] later extended this work to martingales and ranking supermartingales. Chakarov and Sankaranarayanan [START_REF] Chakarov | Expectation invariants for probabilistic program loops as fixed points[END_REF] propose a new notion of probabilistic invariants that generalizes the notion of supermatingales. They give a synthesis approach based on abstract interpretation, but it is not clear how their techniques can prove properties at termination time. Chen et al. [START_REF] Chen | Counterexample-guided polynomial loop invariant generation by lagrange interpolation[END_REF] propose a synthesis tool for verifying Hoare triples in the McIver and Morgan logic, using a combination of Lagrange's interpolation, sampling, and counterexample guided search. One of the novelties is that they can synthesize non-linear invariants. The main disadvantages is that one must manually check the soundness condition, and one must provide a pre-expectation. For instance, we can apply the method of Chen et al. [START_REF] Chen | Counterexample-guided polynomial loop invariant generation by lagrange interpolation[END_REF] to the gambler's ruin process only if we already know that the expected running time is a(b -a). In contrast, we can deduce E[τ ] = a(b -a) knowing only that E[τ ] is finite.

Expected running time. As the termination time of a probabilistic program is a random quantity, it is natural to measure its performance using the average running time. Rough bounds can be obtained from martingale-based termination proofs [START_REF] Ferrer Fioriti | Probabilistic termination: Soundness, completeness, and compositionality[END_REF]. Recently, Chatterjee et al. [START_REF] Chatterjee | Algorithmic analysis of qualitative and quantitative termination problems for affine probabilistic programs[END_REF] showed that arbitrary approximations can be obtained from such proofs when the program is linear. They use Azuma's inequality to obtain a tail distribution of the running time, and later they model check a finite unrolling of the loop. Monniaux [29] propose a similar approach that uses abstract interpretation to obtain the tail distribution of the running time. Kaminski, Katoen, Matheja, and Olmedo [START_REF] Kaminski | Weakest Precondition Reasoning for Expected Run-Times of Probabilistic Programs[END_REF] extend Nielson's proof system [START_REF] Nielson | A hoare-like proof system for analysing the computation time of programs[END_REF] to bound the expected running time of probabilistic programs.

Recurrence analysis. Our synthesis approach is similar to the use of recurrences relations for the synthesis of non-probabilistic invariants [START_REF] Ammarguellat | Automatic recognition of induction variables and recurrence relations by abstract interpretation[END_REF][START_REF] Rodríguez-Carbonell | Generating all polynomial invariants in simple loops[END_REF][START_REF] Kovács | Reasoning algebraically about p-solvable loops[END_REF]. The main idea is to find syntactic or semantic recurrences relations, and later simplify them using known closed forms to obtain loop invariants. In essence, we apply algebraic identities to simplify the complex martingales from Doob's decomposition. The difference is that our simplifications are more complex as we cannot always obtain a closed form but a simpler summation. However, we obtain the same closed form when we apply Doob's decomposition to inductive variables. Another difference is that we rely on the syntactic criteria to identify which values are predictable and which values are random.

Conclusion

We proposed a novel method for automatically synthesizing martingales expressions for stochastic processes. The basic idea is to transform any initial expression supplied by the user into a martingale using Doob's decomposition theorem. Our method complements the state-of-the-art synthesis approaches based on constraint solving. On one hand, we always output a martingale expression, we are able to synthesize non-inductive martingales, and since we do not rely on quantifier elimination, we can synthesize polynomial expression of very high degree. On the other hand, we do not provide any completeness result, and the shape of martingale is difficult to predict.

We considered several classical case studies from the literature, combining our synthesis method with the optional stopping theorem and non-probabilistic invariants to infer properties at termination time in a fully automatic fashion.

Future work includes extending our approach to programs with arrays and improving the tool with automated procedures for checking side-conditions. It would also be interesting to consider richer programs, say distributions with parameters that depend on program state. Another possible direction would be improving the simplification procedures; possibly, the tool could produce simpler facts. Experimenting with more advanced computer-algebra systems and designing simplification heuristics specialized to handling the conditional expectations synthesized by Doob's decomposition are both promising future directions. It would also be interesting to integrate our method as a special tool in systems for interactive reasoning about probabilistic computations.

Fig. 2 :

 2 Fig. 2: Tool pipeline

  match 11 = 1) =⇒ (match 4 = 1) (match 11 = 1) =⇒ (match 1 = 1) (match 11 = 1) =⇒ (match 0 = 1) (match 11 = 1) =⇒ (match j = 0).(for j = 0, 1, 4, 11)

Table 1 :

 1 Preliminary benchmarks.

	Example Running time (s)
	geom	0.14
	gamble	0.11
	gamble2	0.17
	miniabra	0.87
	fullabra	3.58

A basic version of the optional stopping theorem will suffice for our purposes, but there are alternative versions that don't require finite expected stopping time and bounded increments.

We focus on programs for which our method achieves full automation. For instance, we exclude conditional statements because it is difficult to fully automate the resulting simplifications. We note however that there are standard transformations for eliminating conditionals; one such transformation is if-conversion, a well-known compiler optimization[START_REF] Allen | Conversion of control dependence to data dependence[END_REF].

Although most of the ranking supermatingales needed in our case studies are nonlinear, the bounded variants are always linear.

Acknowledgments. We thank the anonymous reviewers for their helpful comments. This work was partially supported by NSF grants TWC-1513694 and CNS-1065060, and by a grant from the Simons Foundation (#360368 to Justin Hsu).


Gilles Barthe

Thomas Espitau

Benjamin Grégoire

Justin Hsu

Léo Stefanesco

Pierre-Yves Strub

Relational reasoning via probabilistic coupling

Probabilistic coupling is a powerful tool for analyzing probabilistic processes. Roughly, coupling two processes requires finding an appropriate witness process that characterizes both processes in the same probability space. Applications of coupling include reasoning about convergence of distributions, and stochastic dominance-a probabilistic version of a monotonicity property. While the mathematical definition of coupling looks rather complex and difficult to manipulate, we show that the relational program logic pRHL-the logic underlying the EasyCrypt cryptographic proof assistantinternalizes a generalization of probabilistic coupling. We demonstrate how to express and verify classic examples of couplings in pRHL, and we mechanically verifying several couplings in EasyCrypt.

Introduction

Probabilistic couplings [START_REF] Thorisson | Coupling, Stationarity, and Regeneration[END_REF][START_REF] Lindvall | Lectures on the coupling method[END_REF][START_REF] Villani | Optimal transport: old and new[END_REF] are a powerful mathematical tool for reasoning about pairs of probabilistic processes: streams of values that evolve randomly according to some rule. While the two processes may be difficult to analyze independently, a probabilistic coupling arranges processes {u i }, {v i } in the same space-for the simplest form of couplings, by viewing the pair of processes as randomly evolving pairs of values {(u i , v i )}-while ensuring certain coupling requirements. In this way, a coupling can coordinate the samples between the two processes so that the coupled process satisfies certain properties.

From the point of view of program verification, a coupling is a relational program property, which describes the relation between two programs (perhaps one program run on two different inputs, or two completely different programs). At first, it may seem that coupling is merely a single example of a relational property, but coupling is a particularly interesting property for several reasons.

Useful consequences. Couplings are known to imply many other relational properties, and are often used as a tool in mathematical proofs.

A classic use of coupling is showing that the distribution of the value of two random processes started in different locations eventually converges to the same distribution if we run the processes long enough. This property is a kind of memorylessness-or Markovian-property: The long-term behavior of the process is independent of where it started from. To prove memorylessness, the typical strategy is to couple the two processes so that their values move closer together; once the values meet, the two processes move together, yielding the same distribution.

A different use of couplings is showing that one (numeric-valued) process is, in some sense, bigger than the other. This statement has to be interpreted carefully-since both processes evolve independently, we can't guarantee that one process is always larger than the other on all traces. Stochastic domination turns out to be the right definition: for any k, we require Pr

[u ≥ k] > Pr[v ≥ k].
To verify this property, it is sufficient to find a coupling of a particular form.

Relational from non-relational. In many examples, the behavior of the second coupled process is specified by the behavior of the first; in such cases, the coupling allows us to reason just about the first process. In other words, a coupling allows us to prove certain relational properties as a property of a single program.

Compositional proofs. Not only is coupling a useful property, proofs of coupling are also intriguing. Typically, couplings are proved by considering corresponding samples of the two processes, step by step; paper proofs call this process "building a coupling", reflecting the piecewise construction of the proof. In this sense, couplings can be proved locally by considering small pieces of the programs in isolation. This feature greatly simplifies mechanical verification of couplings.

Contributions

In this paper, we apply relational program verification to probabilistic couplings. While the mathematical definition of coupling is rather involved and seemingly far from program verification technology, our primary insight is that the logic pRHL from Barthe, Grégoire, and Zanella Béguelin [START_REF] Barthe | Formal certification of code-based cryptographic proofs[END_REF] already internalizes coupling in disguise. More precisely, pRHL is built around a lifting construction, which turns a relation R on two sets A and B into a relation R † over the set of sub-distributions over A and the set of sub-distributions over B. Two programs are related by R † precisely when there exists a coupling of their output sub-distributions whose support only contains pairs of values (u, v) which satisfy R.

This observation has three immediate consequences. First, by selecting the relation R appropriately, we can model a wide variety of coupling properties, like distribution equivalence and stochastic domination. Second, by utilizing the proof system of pRHL, we have a convenient method for constructing couplings by reasoning about the programs while abstracting away the underlying details of the coupling. Finally, we can leverage EasyCrypt, a proof assistant implementing an extended version of pRHL, to mechanical verify couplings.

Preliminaries

Probabilistic coupling. We begin by giving an overview of probabilistic coupling. As we described before, a coupling places two probabilistic processes (viewed as probability distributions) in the same probabilistic space.

We will work with sub-distributions over discrete (finite or countable) sets. A sub-distribution µ over a discrete set A is a function A → [0, 1] such that a∈A µ(a) ≤ 1, and its support supp(µ) is the pre-image of (0, 1]. We let Distr(A) denote the set of sub-distributions over A. Every sub-distribution can be given a monadic structure; the unit operator maps every element a in the underlying set to its Dirac distribution δ a and the monadic composition bind(µ, F ) ∈ Distr(B) of µ ∈ Distr(A) and

F : A → Distr(B) is bind(µ, F )(b) = a∈A µ(a) × F (a)(b).
When working with sub-distributions over tuples, the probabilistic versions of the usual projections on tuples are called marginals. The first and second marginals π 1 (µ) and π 2 (µ) of a distribution µ over A×B are defined by π 1 (µ)(a) = b∈B µ(a, b) and π 2 (µ)(b) = a∈A µ(a, b). We can now formally define coupling.

Definition 1. The Frechet class F(µ 1 , µ 2 ) of two sub-distributions µ 1 and µ 2 over A and B respectively is the set of sub-distributions µ over A × B such that π 1 (µ) = µ 1 and π 2 (µ) = µ 2 . Two distributions µ 1 , µ 2 are said to be coupled with witness µ if µ ∈ F(µ 1 , µ 2 ), i.e. µ is in the Frechet class of µ 1 , µ 2 .

Lifting relations. Before introducing pRHL, we describe the lifting construction. This operation allows pRHL to make statements about pairs of (sub-)distributions, and is a generalized form of probabilistic coupling.

The idea is to define a family of couplings based on the support of the witness distribution. Given a relation R ⊆ A × B and two distributions µ 1 and µ 2 over A and B respectively, we let L R (µ 1 , µ 2 ) denote the subset of sub-distributions µ ∈ F(µ 1 , µ 2 ) such that supp(µ) ⊆ R. Given a ground relation R, we view distributions in L R as witnesses for a lifted relation on distributions.

Definition 2. The lifting of a relation

R ⊆ A × B is the relation R † ⊆ Distr(A) × Distr(B) with µ 1 R † µ 2 iff L R (µ 1 , µ 2 ) = ∅.
Before turning to the definition of pRHL, we give some intuition for why lifting is useful. Roughly, if we know two distributions are related by a lifted relation R † , we can treat two samples from the distribution as if they were related by R. In other words, the lifting machinery gives a powerful way to translate between information about distributions and information about samples. Deng and Du [START_REF] Deng | Logical, metric, and algorithmic characterisations of probabilistic bisimulation[END_REF] provide an excellent introductory exposition to lifting, and give several equivalent characterizations of lifting, including an inductive definition and a definition based on maximal flows.

A pRHL primer

We are now ready to present pRHL, a relational program logic for probabilistic computations. In its original form [START_REF] Barthe | Formal certification of code-based cryptographic proofs[END_REF], implemented in the EasyCrypt proof assistant [START_REF] Barthe | Computer-aided security proofs for the working cryptographer[END_REF], pRHL reasons about programs written in a imperative language extended with random assignments with the following syntax of commands: Assertions in the language are first-order formulae over generalized expressions. The latter are built from tagged variables x 1 and x 2 , which correspond to the interpretation of the program variable x in the first and second memories. Assertions in pRHL are deterministic and do not refer to probabilities.

c ::= x ← e | x $ ← d | if e then c else c | while e do c | skip | c; c Sample f ∈ T1 1-1 -→ T2 ∀v ∈ T1. d1(v) = d2(f v) x1 $ ← d1 ∼ x2 $ ← d2 : ∀v, Φ[v/x1, f (v)/x2] ⇒ Φ If Ψ ⇒ e1 = e2 c1 ∼ c2 : Ψ ∧ e1 ⇒ Φ c 1 ∼ c 2 : Ψ ∧ ¬e1 ⇒ Φ if e1 then c1 else c 1 ∼ if e2 then c2 else c 2 : Ψ ⇒ Φ While Φ ⇒ e1 = e2 c1 ∼ c2 : Φ ∧ e1 ⇒ Φ while e1 do c1 ∼ while e2 do c2 : Φ ⇒ Φ ∧ ¬e1
Definition 3. A pRHL judgment is a quadruple of the form c 1 ∼ c 2 : Ψ ⇒ Φ,
where Ψ and Φ are assertions, and c 1 and c 2 are separable statements, i.e. they do not have any variable in common. A judgment is valid iff for every memories m 1 and m 2 , we have

(m 1 , m 2 ) |= Ψ ⇒ ([[c 1 ]](m 1 ), [[c 2 ]](m 2 )) |= Φ † .
Judgments can be proved valid with a variety of rules.

Two-sided and one-sided rules. The pRHL logic features two-sided rules (Figure 1) and one-sided rules (Figure 2); roughly speaking, two-sided rules relate two commands with the same structure and control flow, while one-sided rules relate two commands with possibly different structure or control flow; these rules allow pRHL to express asynchronous couplings between programs that may exhibit different control flow.

We point out two rules that will be especially important for our purposes. The rule [Sample] is used for relating two sampling commands. Note that it requires an injective function f : T 1 → T 2 from the domain of the first sampling command to the domain of the second sampling command. When the two sampling commands have the same domain-as will be the case in our examples-f is simply a bijection on T = T 1 = T 2 . This bijection gives us the freedom to specify the relation between the two samples.

The rule [While] is the standard while rule adapted to pRHL. Note that we require the guard of the two commands to be equal-so in particular the two loops must make the same number of iterations-and Φ plays the role of the while loop invariant as usual.

SampleL skip ∼ c : ∀v, Ψ [v/x1] ⇒ Φ x1 $ ← d1 ∼ c : Ψ ⇒ Φ IfL c1 ∼ c : Ψ ∧ e ⇒ Φ c 1 ∼ c : Ψ ∧ ¬e1 ⇒ Φ if e1 then c1 else c 1 ∼ c : Ψ ⇒ Φ WhileL c1 ∼ c : Φ ∧ e1 ⇒ Φ while e1 do c1 lossless while e1 do c1 ∼ c : Φ ⇒ Φ Fig. 2: One-sided proof rules (selection)
Structural and program transformation rules. pRHL also features structural rules that are very similar to those of Hoare logic, including the rule of consequence and the case rule. In addition, it features a rule for program transformations, based on an equivalence relation that provides a sound approximation of semantical equivalence.

For our examples, it is sufficient that the relation models loop range splitting and biased coin splitting, as given by the following clauses:

while e do c while e ∧ e do c; while

e do c x $ ← Bern(p 1 • p 2 ) x 1 $ ← Bern(p 1 ); x 2 $ ← Bern(p 2 ); x ← x 1 ∧ x 2
Conseq c1 ∼ c2 : Ψ ⇒ Φ Ψ ⇒ Ψ Φ ⇒ Φ c1 ∼ c2 : Ψ ⇒ Φ Case c1 ∼ c2 : Ψ ∧ Ψ ⇒ Φ c1 ∼ c2 : Ψ ∧ ¬ Ψ ⇒ Φ c1 ∼ c2 : Ψ ⇒ Φ Equiv c 1 ∼ c 2 : Ψ ⇒ Φ c1 c 1 c2 c 2 c1 ∼ c2 : Ψ ⇒ Φ

From pRHL judgments to probability judgments

We will derive two kinds of program properties from the existence of an appropriate probabilistic coupling. We will first discuss the mathematical theorems, where the notation is lighter and the core idea more apparent, and then demonstrate how the mathematical version can be expressed in terms of pRHL judgments.

Total variation and coupling. The first principle bounds the distance between two distributions in terms of a probabilistic coupling. We first define the total variation distance, also known as statistical distance, on distributions.

Definition 4. Let X and X be distributions over a countable set A. The total variation (TV) distance between X and X is defined by X -X tv

1 2
a∈A |X(a) -X (a)| .

To bound the distance between two distributions, it is enough to find a coupling and bound the probability that the two coupled variables differ.

Theorem 1 (Total variation, see [START_REF] Lindvall | Lectures on the coupling method[END_REF]). Let X and X be distributions over a countable set. Then for any coupling Y = ( X, X ), we have

X -X tv ≤ Pr (x,x )∼Y [x = x ].
This theorem is useful for reasoning about convergence of distributions.

We can find a pRHL analog of this theorem. We first introduce some useful notation. For every expression e of type T and distribution µ over memories, let [[e]] µ be defined as Mlet m = µ in unit m(e); note that [[e]] µ denotes a distribution over T . Similarly, for every event E (modeled as a predicate over memories) and distribution µ over memories, let

[[E]] µ be defined as Mlet m = µ in unit E(m). Thus, [[E]] µ is the probability of event E holding in the distribution µ. Proposition 1. If c 1 ∼ c 2 : Ψ ⇒ Φ ⇒ v 1 = v 2
, where Φ exclusively refers to variables in c 1 , then for every initial memories m 1 and m 2 that satisfy the precondition, the statistical distance between

[[v 1 ]] [[c]](m1) and [[v 2 ]] [[c]](m2) is upper bounded by ¬Φ in [[¬Φ]] [[c]](m1) , i.e. [[v 1 ]] [[c1]](m1) -[[v 2 ]] [[c2]](m2) tv ≤ [[¬Φ]] [[c]](m1) .
This proposition underlies the "up-to-bad" reasoning in EasyCrypt.

Stochastic domination and coupling.

A second relational property of distributions is stochastic domination. Definition 5. Let X and X be distributions over set A with an order relation ≥. We say X stochastically dominates X , written X ≥ sd X , if for all a ∈ A,

Pr x∼X [x ≥ a] ≥ Pr x ∼X [x ≥ a].
Intuitively, stochastic domination defines a partial order on distributions over A given an order over A. Strassen's theorem shows that stochastic dominance is intimately related to coupling.

Theorem 2 (Strassen's theorem, see Lindvall [START_REF] Lindvall | Lectures on the coupling method[END_REF]). Let X and X be distributions over a countable ordered set A. Then X ≥ sd X if and only if there is a coupling Y = ( X, X ) with Y ∈ L ≥ (X, X ).

We can express the "if" direction in the following pRHL form.

Proposition 2. If c 1 ∼ c 2 : Ψ ⇒ v 1 ≥ v 2 ,
then for every initial memories m 1 and m 2 that satisfy the precondition,

[[v 1 ]] [[c]](m1) ≥ sd [[v 2 ]] [[c]](m2) .

Warming up: Random walks

We up with couplings for random walks. These numeric processes models the evolution of a token over a discrete space: at each time step the token will choose its next movement randomly. We will show that starting from any two positions, the distributions of the two positions converges as we take more steps.

The basic random walk

Our first example is a random walk on the integers. Starting at an initial position, at each step we flip a fair coin. If heads, we move one step to the right. Otherwise, we move one step to the left. The code for running process k steps is presented in the left side of Figure 4. The variable H stores the history of coin flips. While this history isn't needed for computation of the result (it is ghost code), this history will be useful for stating invariants about the process. We consider two walks that start at locations start 1 and start 2 , such that start 2start 1 = 2n ≥ 0. We want to show that the distribution on end positions in the two walks converges as k increases. From Theorem 1, it suffices to find a coupling of the two walks, i.e., a way to coordinate their random samplings.

The basic idea is to mirror the two walks. When the first process moves towards the second process, we have the second process also move closer; when the first process moves away, we have the second process move away too. When the two processes meet, we have the two processes make identical moves.

To carry out this plan, we define Σ(H) to be the number of true in H minus the number of false; in terms of the random walk, Σ(H) measures the net change in position of a process with history H. Then, we define a predicate such that P (H) holds when H contains a prefix H' with that Σ(H') = n.

Accordingly, P (H 1 ) holds when the first process has moved at least n spots to the right. Under the coupling, this means that the second process must have at least n spots to the left since the two particles are mirrored. Since the first process starts out exactly 2n to the left of the second process, P (H 1 ) is true exactly when the coupled processes have already met.

To formalize this coupling in pRHL, we aim to couple two copies of the program above, which we denote c 1 and c 2 . We relate the two while loops with rule [While] using the following invariant:

(pos 1 = pos 2 ⇒ pos 1 = i 1 + Σ(H 1 ) ∧ pos 2 = i 2 -Σ(H 1 )) ∧ (P (H 1 ) ⇒ pos 1 = pos 2 ).
The loop invariant states that before the two particles meet, their trajectories are mirrored, and that once they have met, they coincide forever.

To prove that this is an invariant, we need to relate the loop bodies. The key step is relating the two sampling operations using the rule [Sample]; note that we must provide a bijection f from booleans to booleans. We choose the bijection based on whether the two coupled walks have met or not.

More precisely, we perform a case analysis on pos 1 = pos 2 with rule [Case]. If they are equal then the walks move together, so we use the identity map for f ; this has the effect of forcing both processes to see the same sample. If the walks are at different positions, we use the negation map (¬) for f , so as to force the two processes to take opposite steps.

Putting everything together, we can prove the following judgment in pRHL:

c 1 ∼ c 2 : start 1 + 2n = start 2 ⇒ (P (H 1 ) ⇒ pos 1 = pos 2 ).
By Theorem 1, we can bound the TV distance between the final positions. If two memories m 1 , m 2 satisfy m 1 (start) + 2n = m 2 (start), we have

[[pos 1 ]] [[c1]](m1) -[[pos 2 ]] [[c2]](m2) tv ≤ [[¬P (H 1 )]] [[c1]](m1) .
Note that the right hand side depends only on the first program. In other words, proving this quantitative bound on two programs is reduced to proving a quantitative property on a single program-this is the power of coupling.

Lazy random walk on a torus

For a more interesting example of a random walk, we can consider a walk on a torus. Concretely, the position is now a d-tuple of integers in [0, k -1]. The walk first flips a fair coin; if heads it stays put, otherwise it moves. If it moves, the walk chooses uniformly in [1, d] to choose the coordinate to move, and a second fair coin to determine the direction (positive, or negative). The positions are cyclic: increasing from k -1 leads to 0, and decreasing from 0 leads to k -1.

We can simulate this walk with the program in the right side of Figure 4, where u(i) is the i-th canonical base vector in (Z/kZ) d . As before, we store the trace of the random walk in the list H. All arithmetic is done modulo k.

Like the simple random walk, we start this process at two locations start 1 and start 2 on the torus and run for k iterations. We aim to prove that the distributions of the two walks converge as k increases by coupling the two walks, iteration by iteration. Each iteration, we first choose same coordinate crd and the same direction dir in both walks. If the two positions coincide in coordinate crd, we arrange both walks to select the same direction mov, so that the walks either move together, or both stay put. If the two positions differ in crd, we arrange the walks to select opposite samples in mov so that exactly one walk moves.

As in the basic random walk, we can view our coupling as letting the first process evolve as usual, then coordinating the samples of the second process to perform the coupling. In other words, given a history H 1 of samples for the first process, the behavior of the second coupled process is completely specified.

Thus, we can define operators to extract the movements of each walk from the trace H 1 of the samplings of the first process: Σ 1 (i, H 1 ) is the drift of the ith coordinate of the first process, and Σ 2 (i, H 1 ) is the drift of the second process. Essentially, these operators encode the coupling by describing how the second process moves as a function of the first process's samples.

In pRHL, we will use the rule [While] with the following invariant:

∀i ∈ [1, d]. (Σ 1 (i, H 1 ) -Σ 2 (i, H 1 ) = ∆[i] ⇒ pos 1 [i] = pos 2 [i]) ∧ (pos 1 [i] = pos 2 [i] ⇒ pos 1 [i] = start 1 [i] + Σ 1 (i, H 1 ) ∧ pos 2 [i] = start 2 [i] + Σ 2 (i, H 1 )),
where ∆ is the vector start 2start 1 . The invariant describes the position of the two coupled processes in terms of the history H 1 . We are particularly interested in the first conjunct, which gives a condition for the two coupled processes to meet in coordinate i.

To prove that the invariant is preserved, we encode the coupling described above into pRHL, via three uses of the rule [Sample]. The first two samplesfor crd and dir-are coupled with f being identity bijections (on [1, d] and on booleans), ensuring that the processes make identical choices. When sampling mov, we inspect the history H 1 to see whether the two walks agree in position crd.

If so, we choose the identity bijection for mov; if not, we choose negation. This coupling is sufficient to verify the loop invariant.

To conclude our proof, the first conjunct in the invariant implies that we can prove the pRHL judgment c 1 ∼ c 2 : start 2start 1 = ∆ ⇒ Φ, where

Φ (∀i ∈ [1, d]. Σ 1 (i, H 1 ) -Σ 2 (i, H 1 ) = ∆[i]) ⇒ ∀i ∈ [1, d]. pos 1 [i] = pos 2 [i].
Finally, Theorem 1 implies that for any two initial memories m 1 , m 2 with m 2 (start) -m 1 (start) = ∆, we have

[[pos 1 ]] [[c1]](m1) -[[pos 2 ]] [[c2]](m2) tv ≤ [[∃i ∈ [1, d]. Σ 1 (i, H 1 )-Σ 2 (i, H 1 ) = ∆[i]]] [[c1]](m1) .
Again, proving a quantitative bound on the convergence of two distributions is reduced to proving a quantitative bound on a single program.

Combining coupling with program transformation

So far, we have seen examples where the coupling is proved directly on the two original programs c 1 and c 2 . Often, it is convenient to introduce a third program c * that is equivalent to c 1 , and then couple c * to c 2 . Applying transitivity (rule [Equiv]), this gives a coupling between c 1 and c 2 . Let's consider two examples.

Conclusion and future work

We have established the connection between relational verification of probabilistic programs using pRHL, and coupling, which plays a major role in probability theory. Furthermore, we have used the connection by verifying in pRHL several well-known examples of couplings from the literature on randomized algorithms.

More broadly, our work occupies a middle ground between the two main approaches to relational verification: (i) encoding the two programs as a single program and reasoning about the single program (techniques include include cross-products [START_REF] Zaks | Covac: Compiler validation by program analysis of the cross-product[END_REF], self-composition [START_REF] Barthe | Secure Information Flow by Self-Composition[END_REF], and product programs [START_REF] Barthe | Relational verification using product programs[END_REF]); and (ii) using a program logic to reason directly about two programs (techniques include relational Hoare logic [START_REF] Benton | Simple relational correctness proofs for static analyses and program transformations[END_REF], relational separation logic [START_REF] Yang | Relational separation logic[END_REF], and pRHL [START_REF] Barthe | Formal certification of code-based cryptographic proofs[END_REF]).

work also suggests the possibility of a practical proof technique for verifying probability results based on coupling, but a deeper understanding of the links between pRHL and coupling (and their variants thereof), and a more general framework for verification of probabilistic programs will be required.

A more general verification framework. pRHL derivations implicitly construct a sequence of couplings between the executions-modeled as sequences of distributionsof the two programs involved in the final statement. The single most important step in the construction of the coupling is performed when applying the rule [Sample] for random sampling, and more specifically when picking the bijection function f . A careful look at the rule reveals that the coupling is a deterministic coupling, as defined by Villani [START_REF] Villani | Optimal transport: old and new[END_REF]: Definition 6. A coupling (X, Y ) between two spaces X , Y is said to be deterministic if there exists a function T : X → Y such that Y = T (X).

Indeed, deterministic couplings are often sufficient in cryptographic proofs-the original motivation to develop pRHL. However, there are many examples of couplings that cannot be verified using deterministic couplings. For our examples we have worked around the difficulty by using program transformation rules. However, this approach is very specific, and does not generalize. Thus, an interesting generalization of pRHL would be to relax the requirement that the relationship between the two samplings is given by a 1-1 map, by allowing a binary relation which satisfies suitable conditions. Such a general rule could enable a more general class of couplings, and yield more principled proofs of some of the couplings studied in this paper.

Moreover, it would be interesting to extend EasyCrypt with mechanisms for handling the non-relational reasoning in couplings. To prove quantitative bounds on total variation in the random walk example, we need to bound the time it takes for a single random walk to reach a certain position. Obtaining precise bounds on the speed of convergence requires more complex reasoning, and will require an expressive program logic. Such a program logic is under development, but has not yet been integrated in EasyCrypt.

Extending to shift and path coupling. The couplings realized in the random walks are instances of exact couplings, i.e. the paths will eventually merge. A more general notion of coupling for such processes is the shift-coupling, where the merging condition is "relaxed" in the sense that the paths will merge modulo a random time shift. The general theory of path couplings provides similar-shaped inequalities as the ones in exact coupling, allowing powerful mathematical-based reasoning inside the logic with the [Conseq] rule. These coupling notions are complex, and it is not yet clear how they can be verified.

Other examples. The coupling literature is rich with further examples, from simple to challenging. We intend to verify many other examples of couplings, in particular the proof of Dynkin's card trick, and the proof of the constructive Lovasz Local Lemma, a fundamental tool used in the probabilistic method.

Fig. 1 :

 1 Fig. 1: Two-sided proof rules (selection)

Figure 3

 3 Figure 3 provides a selection of structural and program transformation rules.

Fig. 3 :

 3 Fig. 3: Structural and program transformation rules (selection)

←Fig. 4 :

 4 Fig. 4: Two random walks

Two biased coins

Consider a coin flipping process that flips a coin k times, and returns the number of heads observed. We consider this process run on two different biased coins: The first coin has probability q 1 of coming up heads, while the second coin has probability q 2 of coming up heads with q 1 ≥ q 2 . Let the distribution on the number of heads be µ 1 and µ 2 respectively.

Intuitively, it is clear that the first process is somehow bigger than the second process: it is more likely to see more heads, since the first coin is biased with a higher probability. Stochastic dominance turns out to be the proper way to formalize our intuition. To prove it, Proposition 2 implies that we just need to find an appropriate coupling of the two processes.

While it is possible to define a coupling directly by carefully coordinating the corresponding coin flips, we will give a simpler coupling that proceeds in two stages. First, we will couple a program c 1 computing µ 1 to an intermediate program c * . Then, we will show that c * is equivalent to a program c 2 computing µ 2 , thus exhibiting a coupling between µ 1 and µ 2 . Letting r = q 2 /q 1 and denoting the coin flip distribution with probability p of sampling true by Bern(p), we give the programs in Figure 5.

For the first step, we want to couple c 1 and c * . For a rough sketch, we want to use rule [While] with an appropriate loop invariant; here, n 1 ≥ n * . To show that the invariant is preserved, we need to relate the loop bodies. We use the two-sided rule [Sample] when sampling x and y (taking the bijection f to be the identity), the one-sided rule [Sample-L] to relate sampling nothing (skip) in c 1 with sampling z in c * , and the one-sided rule [IfL] to relate the two conditionals. (The one-sided rule is needed, since the two conditionals may take different branches.) Thus, we can prove the judgment c 1 ∼ c * : q 1 ≥ q 2 ∧ r = q 2 /q 1 ⇒ n 1 ≥ n * .

For the second step, we need to prove that c * is equivalent to c 2 . Here, we use a sound approximation to semantic equivalence as described in the introduction. Specifically, we have x $ ← Bern(q 1 • r) y $ ← Bern(q 1 ); r $ ← Bern(r); x ← y ∧ z for the loop bodies; showing equivalence of c * and c 2 is then straightforward. Thus, we can show c * ∼ c 2 : q 1 ≥ q 2 ∧ r = q 2 /q 1 ⇒ n * = n 2 . Applying rule [Equiv] gives the final judgment c 1 ∼ c 2 :

showing stochastic domination by Proposition 2.

Balls into bins: asynchronous coupling

The examples we have seen so far are all synchronous couplings: they relate the iterations of the while loop in lock-step. For some applications, we may want to reason asynchronously, perhaps allowing one side to progress while holding the other side fixed. One example of an asynchronous coupling is analyzing the balls into bins process. We have two bins, and a set of n balls. At each step, we throw a ball into a random bin, returning the count of both bins when we have thrown all the balls. The code is on the left side in Figure 6. Now, we would like to consider what happens when we run two processes with different numbers of balls. Intuitively, it is clear that if the first process throws

Fig. 6: Coupling balls into bins more balls than the second process, it should result in a higher load in the bins; we aim to prove that the first process stochastically dominates the second with the following coupling. Assume that the first process has more balls (n 1 ≥ n 2 ). For the first n 2 balls, we have the two process do the same thing-they choose the same bucket for their tosses. For the last n 1 -n 2 steps, the first process throws the rest of the balls. Evidently, this coupling forces the bins in the first run to have higher load than the bins in the second run.

To formalize this example, we again introduce a program c * , proving equivalence with c 1 and showing a coupling with c 2 . The code for c * is on the right side in Figure 6; we require the dummy input m to be equal to n 2 .

Proving equivalence with program c 1 is direct, using the loop range splitting transformation in EasyCrypt: while e do c while e ∧ e do c; while e do c. Once this is done, we simply need to provide a coupling between c * and c 2 . By our choice of m, we can trivially couple the first loop in c * to the (single) loop in c 2 , ensuring that Φ binA * ≥ binA 2 ∧ binB * ≥ binB 2 after the first loop.

Then, we can apply the one-sided rules to couple the second loop in c * with a skip statement in c 2 . It is straightforward to show that Φ is an invariant in [WhileL], from which we can conclude c * ∼ c 2 :

and by equivalence of c 1 and c * we have c 1 ∼ c 2 : n 1 ≥ n 2 ⇒ binA 1 ≥ binA 2 ∧ binB 1 ≥ binB 2 , enough for stochastic domination by Proposition 2.

Non-deterministic couplings: birth and death

So far, we have seen deterministic couplings, which reuse randomness from the coupled processes in the coupling; this can be seen in the [Sample] rule, when we always choose a deterministic. In this section, we will see a more sophisticated coupling that injects new randomness.

For our example, we consider a classic Markov process. Roughly speaking, a Markov process moves within a set of states each transition depending only on the current state and a fresh random sample. The random walks we saw before are classic examples of Markov processes.

A more complex Markov process is the birth and death chain. The state space is Z, and the process starts at some integer x. At every time step, if the process is at state i, the process has some probability b i of increasing by one, and some probability a i of decreasing by one. Note that a i and b i may add up to less than 1: there can be some positive probability 1 -a i -b i where the process stays fixed.

To model this process, we define a sum type Move with three elements (Left, Right and Still) which correspond to the possible moves a process can make. Then, the chains are modeled by the code in the left of Figure 7, where the distribution bd(state) is the distribution of moves from state. Just like the biased coin and balls into bins processes, we want to prove stochastic domination for two processes started at states start 1 ≥ start 2 via coupling. The difficulty is that if the processes become adjacent and they both move, the two processes may swap positions, losing stochastic domination.

The solution is to use a coupling when the two processes are on two adjacent states as in Mufa [START_REF] Mufa | Optimal markovian couplings and applications[END_REF]. Unlike the previous examples, the coupling is not deterministic: the behavior of one process is not fully determined by the randomness of the other. Our loop invariant is the usual one for stochastic domination: state 1 ≥ state 2 . To show that this invariant is preserved, we perform a case analysis on whether state 1 = state 2 , state 1 = state 2 + 1 or state 1 > state 2 + 1.

We focus on interesting case: the middle one, when the states are adjacent. Here, we perform a trick: we switch c 1 , c 7, where proj [1|2] is the first pair projection π 1 in c 1 , and the second pair projection π 2 in c 2 .

All that remains is to prove a coupling between c * 1 and c * 2 satisfying the loop invariant state 1 ≥ state 2 . With adjacent states, dcouple is given by the following function from pairs of moves to probabilities: Note that the case (Left, Right) has probability 0: this forbids the first process from skipping past the second process. Now the coupling is easy: we simply require both samples from dcouple to be the same. Since state 1 = state 2 + 1 and the distribution never returns (Left, Right), the loop invariant is trivially preserved. This shows the desired coupling, and stochastic domination by Proposition 2.

A Graphical depictions of random processes

We depict the two random walk processes in Figures 8 and9


Andrea Farina
email: andrea.farina@polimi.it

Marta Betcke

Laura Di Sieno

Andrea Bassi

Nicolas Ducros

Antonio Pifferi

Gianluca Valentini

Simon Arridge

Cosimo D'andrea

And Cosimo D'andrea

Multiple-view diffuse optical tomography system based on time-domain compressive measurements

Multiple-view diffuse optical tomography system based on time-domain compressive measurements

Compressive sensing is a powerful tool to efficiently acquire and reconstruct an image even in diffuse optical tomography (DOT) applications. In this work, a time-resolved DOT system based on structured light illumination, compressive detection, and multiple view acquisition has been proposed and experimentally validated on a biological tissue-mimicking phantom. The experimental scheme is based on two digital micromirror devices for illumination and detection modulation, in combination with a time-resolved single element detector. We fully validated the method and demonstrated both the imaging and tomographic capabilities of the system, providing state-of-the-art reconstruction quality.

multiple view measurement scheme [5,6] in order to increase the tomographic spatial resolution. Moreover, further data, such as spectral and temporal information, are crucial [7,8]. Temporal information provides three main advantages: (i) better disentanglement of absorption/scattering properties, (ii) temporal encoding of photon depth, and (iii) fluorescence lifetime quantification in the case of FMT. Spectral information (i.e., different excitation/detection wavelength) allows one to discriminate among tissue chromophores. Hence, DOT/FMT turns out to be a highly multidimensional problem with the drawback of a huge data set being generated. This represents a practical limitation of these techniques because of the extremely long acquisition and computational times, which are not typically compatible with clinical and preclinical needs. Hence, a reduction of the acquired data set by preserving the spatial resolution, or, more generally, the data set information content, is highly desirable. Following this concept, different studies have recently exploited the fact that a highly scattering medium (such as biological tissue) behaves as a low-pass filter in the spatial domain. Hence, a few illumination patterns, instead of the more typical raster scanning approach, can be adopted without losing significant spatial information [8][9][10]. This in turn leads to a reduction of the data set dimension and, consequently, of the acquisition and computational times. Recent studies have exploited such an approach in both imaging and tomographic schemes, and detection is generally carried out by a parallel detector such as a CCD, CMOS, or gated cameras [11]. Moreover, the use of a wide-field approach (such as the case of structured illumination) allows one to illuminate the sample with high power without exceeding safety limits. This improves the signal-to-noise ratio.

Recently, a patterned detection [12], following the singlepixel camera (SPC) scheme [13], has been proposed for FMT applications, as well as for PhotoAcoustics [14]. Basically, the image of the diffused light exiting the sample is modulated spatially and subsequently focused on a single element detector. This operation is equivalent to projecting the image on an element of a base set, such as Fourier, Wavelets, or Hadamard patterns. By repeated acquisitions for different base elements, In the last decade, the possibility to quantitatively reconstruct absorbing, scattering, and fluorescent inclusions within in vivo organisms attracted a great deal of interest for diagnostic purposes (e.g., tumor detection) [1], functional studies (e.g., brain oximetry) [2], and molecular imaging on small animals (e.g., pharmacological research) [3]. The general measurement scheme consists of illuminating a sample and detecting the diffused light exiting it. Then, by solving the inverse problem, based on a model of photon propagation through the biological tissue, the optical parameters in each point of the sample can be reconstructed quantitatively. These modalities are usually referred to as diffuse optical tomography (DOT) or fluorescence molecular tomography (FMT) when, respectively, the absorption/ scattering or fluorescence properties are reconstructed. The performance of DOT/FMT is mainly characterized by its capability to resolve the position and shape of inhomogeneities inside the tissue, and, consequently, improve the quantification capability of their optical parameters. Previous studies have demonstrated the importance of a dense source/detector [4]anda it is possible to recover the same image as would be measured in a conventional pixel basis. Due to the fact that a highly scattering medium acts as a low-pass filter in the spatial frequency domain, just a few frequencies are needed. This approach has the great advantage of exploiting the superior characteristics of a single detector (e.g., higher temporal resolution and larger spectral bandwidth) at a lower cost with respect to a parallel device. Moreover, compared with raster scanning, a further advantage is the acquisition speed given by a wide-field detection analogous to a structured illumination approach. Finally, it is worth mentioning that both structured illumination and detection open the possibility to acquire images and reconstructions with increasing spatial details by increasing the number of measurements.

Whereas the patterned detection approach has been successfully demonstrated for fluorescence optical tomography in a single view [12], a fully tomographic modality requires multiple views to reduce the ill-posedness in depth resolution, which leads to a more challenging experimental arrangement. In this work, we propose a multiple-view, time-domain compressed sensing DOT system exploiting Hadamard patterns in both the illumination and collection planes, and applicable to nonplanar geometries. The system has been experimentally validated on tissue phantoms with absorbing inclusions, demonstrating both imaging and tomographic capabilities.

The experimental setup is sketched schematically in Fig. 1. The sample is illuminated by a pulsed structured light while detection is carried out by either a time-resolved (TR) SPC or a continuous-wave (CW) parallel detector. The sample is placed on a rotational stage to allow different view acquisitions. By means of an acousto-optic tunable filter, light pulses at 650 nm are selected spectrally from a picosecond pulsed supercontinuum (repetition rate of 80 MHz) laser source (SuperK Extreme, NKT). Structured illumination is carried out by a digital micromirror device (DMD Discovery kit 1100, Vialux), which spatially modulates the light, and an objective lens (f 50 cm) to create an image over an area of 3 cm × 3 cm of the sample. The diffused light, exiting the sample over an area of about 2 cm × 2 cm, is imaged by a lens (f 60 cm) on a second digital micromirror device (DMD) (DMD Discovery 4100, Vialux). A flip mirror allows us to image the DMD plane on either a low-noise 16-bit cooled CCD camera (Versarray 512, Princeton Instruments) or a single element detector. The latter consists of a long-working-distance objective (10×∕0.25), which focuses the light reflected by the second DMD onto an optical fiber of 1 mm diameter. The light exiting the fiber is finally detected by a photomultiplier (PMT) (HPM-100-50, Becker and Hickl) connected to a time-correlated single-photon counting (TCSPC) board, which samples the temporal profile of the diffuse light. The system is fully computer-controlled by a homemade LabView software, enabling automated acquisition of the whole data set (illumination/detection patterns, sample rotation, and acquisition). The sample is a homogeneous cylindrical tissue mimicking a phantom (∅ 20 mm, height 50 mm) made of epoxy resin, TiO 2 (as a scatterer), and a toner (as an absorber). By means of a time-resolved spectroscopy system [15], the optical parameters were measured: absorption coefficient (μ a ) of about 0.01 mm -1 and reduced scattering coefficient (μ 0 s ) of about 1 mm -1 . Two holes (see Fig. 4), drilled into the sample (∅ 1.6 mm), allowed us to insert either solid or liquid absorbing inclusions. In particular, to better simulate a realistic perturbation, three solutions of calibrated ink and Intralipid have been prepared [16], giving μ 0 s ∼ 1 mm -1 (the same as that of the background) and μ a of 0.05∕0.05 mm -1 (Exp 1), 0.05∕0.1 mm -1 (Exp 2), and 0.05∕0.03 mm -1 (Exp 3) for inclusions A and B.

Initially, images have been acquired by means of the CCD camera on the detection side and a low-cost camera on the illumination side, to register the illumination/detection area over the sample. Then 360 shadows of the object (every 1°) have been acquired to create the mesh [6]. It is worth emphasizing that precise calibration is critical to achieving an accurate simulation of the forward problem, which in turn is a prerequisite to obtain a high-quality tomographic reconstruction.

Measurements have been performed on the phantom with and without the absorbing inclusions. The acquisition procedure is carried out by a complete 360°rotation of the sample with steps of 45°(eight views). For each view, 8 × 8 ordered Walsh-Hadamard (WH) patterns, covering an area of 1.3 cm × 1.3 cm on the sample, have been used for both illumination and detection. Each WH pattern consists of two states (-1 to 1). Hence, two positive patterns (ranging from 0 to 1), complementary to one another, have been acquired and properly subtracted to obtain the desired WH pattern. The acquisition time for each pattern is 1 s with 800 kHz as maximum count rate to fulfill the single-photon statistics. This last parameter is, indeed, the limiting factor on the overall acquisition time, which is about 25 min. A full-pixel image can be recovered by applying a fast WH inverse transform to the detected data [17].

For the reconstruction of the absorption map in the volume, the following objective function has been constructed:

Ψx 1 2 X n y n -f n x f 0n x 2 τRx; ( 1
)
where x is the absorption coefficient in every mesh element, y n is the measurement performed with the SPC, f n and f 0n are the forward models, heterogeneous and homogeneous, respectively, τ is the hyper-parameter, R is a regularization functional, and n is the measurement index. The software TOAST, a finiteelement-based solver [18], has been used to calculate f n x. In order to minimize the objective function in Eq. ( 1), a damped Gauss-Newton method based on a one-dimensional line-search algorithm [19] has been implemented. A total variation regularization functional has been used. In the calculation of both the forward model and the Jacobian, the instrumental response function (IRF) has been taken into account by convolution in time. First, measurements have been carried out using black solid rods as inclusions to demonstrate the imaging capability of the system and to estimate the number of patterns to be used in the tomographic reconstruction. In particular, time-resolved data acquired by the SPC have been integrated over time to obtain CW data and compared with the CCD images. An example of the images acquired by the CCD and the ones based on the SPC (by spatially modulating the detection) is shown in Fig. 2. We observe a good agreement between the two images, which is improved by increasing the number of adopted patterns as reported in Fig. 2, where the root-mean-square error (RMSE) is reported as a function of the WH pattern order. In particular, we do not observe a significant improvement for a WH pattern order higher than 8. Moreover, it is possible to observe that the RMSE plot for the inhomogeneous phantom presents a higher variability among the different views with respect to the homogeneous case. It is worth stressing that the number of required patterns strongly depends on the optical parameters/shape of the sample and the position/dimension of the inclusions.

In order to explore the imaging capability of the proposed method, the relative contrast, calculated as the difference between heterogeneous and homogeneous images divided by the homogeneous one, provided by one solid inclusion, in both the CCD and SPC images, is reported, for eight different views, in Fig. 3. In particular, three cases are reported: (i) the sample is illuminated with 16 × 16 ordered WH patterns while the detection side has a uniform square pattern (Fig. 3, first row); (ii) the sample is illuminated with a uniform square pattern while the detection side is spatially modulated by 16 × 16 ordered WH patterns (Fig. 3, second row); (iii) CCD images by using a uniform square illumination pattern are also reported (Fig. 3, third row). Cases (ii) and (iii) show good agreement; in particular, the presence of the absorbing inclusion can be clearly observed when it is located, during sample rotation, closer to the detector (see the vertical blue bar in the image at 0°or 45°f or OUT). On the contrary, for the other views, the inclusion cannot be observed because of the scattering. In case (i), there is no correspondence between the images acquired by the CCD and SPC. In particular, we observe that, by modulating the illumination, we can better observe the inclusion for views where it is closer to the illumination source (see the vertical blue bar in the image at 180°or 225°for IN). In fact, the SPC approach measures the integral of the signal; then the imaging capability is not influenced by the scattering events followed by photons after impinging on the inclusion as occurs for the CCD [20]. These examples demonstrate the imaging capability of the proposed method and, in particular, the importance of the choice of illumination/detection patterns according to the view, the sample (shape and optical parameters) and inclusions, for both imaging and reconstruction.

As a first demonstration of the tomographic capability of the proposed system, three tomographic reconstructions by using an early gate of the time-resolved profile have been carried out. For each view, a single constant illumination pattern was used and detection was performed with 8 × 8 WH patterns while the temporal gate has been chosen corresponding to the rising edge of the TR signal, here resulting in a time window of 500 ps length. The homogeneous measurements have been used to scale the inhomogeneous phantom data to match the magnitude of the forward temporal point-spread functions (TPSFs). The mesh used for the forward problem has 120,000 elements and 1016 TPSFs have been generated (127 WH patterns for eight views) and sampled in 156 temporal steps of 8 ps length. The computational time for the forward problem is about 10 s on a machine with 10 Dual Intel Xeon processors of 2.3 GHz each, an Nvidia Tesla K-40 GPU, and 64 GB of RAM memory. First, the eight TR measurements of the homogeneous phantom with planar illumination have been used to retrieve the background optical properties. For this purpose, a Levenberg-Marquardt fitting procedure with TOAST as a forward solver has been used, obtaining μ a 0.01 mm -1 and μ 0 s 0.88 mm -1 showing a good agreement with the properties measured using the spectroscopy system mentioned above. The reconstruction has been carried out on a regular grid of 85,731 points containing the whole cylindrical mesh with a voxel size of 0.5 mm × 0.5 mm × 1 mm. The Jacobian has been calculated using the adjoint method [21] deploying the fast Fourier transform for fast implementation of the temporal convolution. The Gauss-Newton algorithm has been terminated after three iterations, after which the reconstruction ceased to improve. The overall reconstruction time was about 3 h.

Figures 4(a), 4(c), and 4(e) show the tomographic reconstructions of μ a at different vertical slices. Due to the limited field of view of both illumination and detection, only a part of the cylinder can be reconstructed (about 16 mm centered at about 14 mm from the top). We observe good reconstruction quality pertaining to both the localization and relative contrast of the two inclusions. By fitting the reconstructed inclusions with a three-dimensional Gaussian function, we obtained a total contrast of about 0.3 times the truth for all cases. Moreover, in Figs. 4(b), 4(d), and 4(f ), normalized line profiles across the inclusions at the z 13 mm plane are shown for all the experiments. We observe a worsening of the localization for inclusion B of Exp 3 probably due to the reduced contrast produced by the low-absorbing solution poured in it. Finally, in order to quantify the localization capability of the reconstruction, the center of mass (COM) for each inclusion has been computed on a region twice larger than the inclusions (see Table 1).

In conclusion, a fully tomographic time-resolved DOT system based on sampling in the spatial frequency domain (both illumination and detection spaces) and multiple view acquisition has been proposed and validated on a tissue-mimicking phantom, demonstrating state-of-the-art reconstruction quality. Moreover, the imaging capability of the system has been validated in CW by comparing SPC with standard CCD acquisition, showing the importance of the choice of illumination/detection patterns for imaging purposes. Future work will be devoted to the optimization of the data set (choice of illumination/detection patterns, number of views, and temporal gates) and system improvements (detection efficiency, calibration procedure) in order to strongly reduce the acquisition time, while preserving or even increasing the information content. In particular, adaptive basis scan approaches will be investigated [22].

Funding. Royal Society (IE140259); Fondazione Cariplo (2013-0615); Laserlab-Europe EU-H2020 654148.

Fig. 1 .

 1 Fig. 1. Experimental setup.

Fig. 3 .

 3 Fig. 3. Reconstruction of μ a using an early gate. Slices are shown from z 19 mm (top of the cylinder) to z 4 mm with a step of 1 mm. The z 0 mm plane is at the middle of the cylinder. Panels (a), (c), and (e) show reconstructions (values are in mm -1 ) for Exp 1, Exp 2, Exp 3, respectively. Panels (b), (d), and (f) show the normalized profiles of the absorption perturbations along line Γ connecting the centers of the two inclusions on the plane at z 13 mm.

Fig. 2 .

 2 Fig. 2. Relative contrast for different rotation views (columns) in the case of modulated illumination (first row), modulated detection (second row), and CCD images by using uniform square illumination (third row). Images in rows 1 and 2 were obtained by a WH inverse transform.

Fig. 4 .

 4 Fig. 4. Comparison between CCD and SPC images outcoming from the homogeneous (first row) and heterogeneous (second row) phantoms. First column: image measured with the CCD. Second column: image measured by the SPC using 16 × 16 WH patterns for a single view. Third column: RMSE as a function of WH order for eight views.

Table 1 .

 1 COM Coordinates in Millimeters

		x A	y A	z A	x B	y B	z B
	TRUE	4.1	3.0	12.5	-4.3	-0.2	12.5
	Exp 1	3.5	2.4	13.1	-4.5	0.1	13.2
	Exp 2	4.6	3.4	12.8	-3.9	-0.5	12.3
	Exp 3	3.15	2.4	13.0	-6.2	0.2	15.1

Acknowledgment. We acknowledge the support of NVIDIA with the donation of the Tesla K40 GPU used for this research.