Wojciech Zaremba
email: woj.zaremba@gmail.com

Arthur Gretton
email: arthur.gretton@gmail.com

Matthew Blaschko
email: matthew.blaschko@inria.fr

Équipe Galen

B-tests: Low Variance Kernel Two-Sample Tests

A family of maximum mean discrepancy (MMD) kernel two-sample tests is introduced. Members of the test family are called Block-tests or B-tests, since the test statistic is an average over MMDs computed on subsets of the samples. The choice of block size allows control over the tradeoff between test power and computation time. In this respect, the B-test family combines favorable properties of previously proposed MMD two-sample tests: B-tests are more powerful than a linear time test where blocks are just pairs of samples, yet they are more computationally efficient than a quadratic time test where a single large block incorporating all the samples is used to compute a U-statistic. A further important advantage of the B-tests is their asymptotically Normal null distribution: this is by contrast with the U-statistic, which is degenerate under the null hypothesis, and for which estimates of the null distribution are computationally demanding. Recent results on kernel selection for hypothesis testing transfer seamlessly to the B-tests, yielding a means to optimize test power via kernel choice.

Introduction

Given two samples {x i } n i=1 where x i ∼ P i.i.d., and {y i } n i=1 , where y i ∼ Q i.i.d, the two sample problem consists in testing whether to accept or reject the null hypothesis H 0 that P = Q, vs the alternative hypothesis H A that P and Q are different. This problem has recently been addressed using measures of similarity computed in a reproducing kernel Hilbert space (RKHS), which apply in very general settings where P and Q might be distributions over high dimensional data or structured objects. Kernel test statistics include the maximum mean discrepancy [START_REF] Gretton | A kernel two-sample test[END_REF][START_REF] Fromont | Kernels based tests with nonasymptotic bootstrap approaches for two-sample problems[END_REF] (of which the energy distance is an example [START_REF] Sejdinovic | Hypothesis testing using pairwise distances and associated kernels[END_REF][START_REF] Baringhaus | On a new multivariate two-sample test[END_REF][START_REF] Székely | Testing for equal distributions in high dimension[END_REF]), which is the distance between expected features of P and Q in the RKHS; the kernel Fisher discriminant [START_REF] Harchaoui | Testing for homogeneity with kernel Fisher discriminant analysis[END_REF], which is the distance between expected feature maps normalized by the feature space covariance; and density ratio estimates [START_REF] Yamada | Relative density-ratio estimation for robust distribution comparison[END_REF]. When used in testing, it is necessary to determine whether the empirical estimate of the relevant similarity measure is sufficiently large as to give the hypothesis P = Q low probability; i.e., below a user-defined threshold α, denoted the test level. The test power denotes the probability of correctly rejecting the null hypothesis, given that P = Q.

The minimum variance unbiased estimator MMD u of the maximum mean discrepancy, on the basis of n samples observed from each of P and Q, is a U-statistic, costing O(n 2 ) to compute. Unfortunately, this statistic is degenerate under the null hypothesis H 0 that P = Q, and its asymptotic distribution takes the form of an infinite weighted sum of independent χ 2 variables (it is asymptotically Gaussian under the alternative hypothesis H A that P = Q). Two methods for empirically estimating the null distribution in a consistent way have been proposed: the bootstrap [START_REF] Gretton | A kernel two-sample test[END_REF], and a method requiring an eigendecomposition of the kernel matrices computed on the merged samples from P and Q [START_REF] Gretton | A fast, consistent kernel twosample test[END_REF]. Unfortunately, both procedures are computationally demanding: the former costs O(n 2 ), with a large constant (the MMD must be computed repeatedly over random assignments of the pooled data); the latter costs O(n 3 ), but with a smaller constant, hence can in practice be faster than the bootstrap. Another approach is to approximate the null distribution by a member of a simpler parametric family (for instance, a Pearson curve approximation), however this has no consistency guarantees.

More recently, an O(n) unbiased estimate MMD l of the maximum mean discrepancy has been proposed [START_REF] Gretton | A kernel two-sample test[END_REF]Section 6], which is simply a running average over independent pairs of samples from P and Q. While this has much greater variance than the U-statistic, it also has a simpler null distribution: being an average over i.i.d. terms, the central limit theorem gives an asymptotically Normal distribution, under both H 0 and H A . It is shown in [START_REF] Gretton | Optimal kernel choice for large-scale two-sample tests[END_REF] that this simple asymptotic distribution makes it easy to optimize the Hodges and Lehmann asymptotic relative efficiency [START_REF] Serfling | Approximation Theorems of Mathematical Statistics[END_REF] over the family of kernels that define the statistic: in other words, to choose the kernel which gives the lowest Type II error (probability of wrongly accepting H 0 ) for a given Type I error (probability of wrongly rejecting H 0 ). Kernel selection for the U-statistic is a much harder question due to the complex form of the null distribution, and remains an open problem.

It appears that MMD u and MMD l fall at two extremes of a spectrum: the former has the lowest variance of any n-sample estimator, and should be used in limited data regimes; the latter is the estimator requiring the least computation while still looking at each of the samples, and usually achieves better Type II error than MMD u at a given computational cost, albeit by looking at much more data (the "limited time, unlimited data" scenario). A major reason MMD l is faster is that its null distribution is straightforward to compute, since it is Gaussian and its variance can be calculated at the same cost as the test statistic. A reasonable next step would be to find a compromise between these two extremes: to construct a statistic with a lower variance than MMD l , while retaining an asymptotically Gaussian null distribution (hence remaining faster than tests based on MMD u ). We study a family of such test statistics, where we split the data into blocks of size B, compute the quadratic-time MMD u on each block, and then average the resulting statistics. We call the resulting tests B-tests. As long as we choose the size B of blocks such that n/B → ∞, we are still guaranteed asymptotic Normality by the central limit theorem, and the null distribution can be computed at the same cost as the test statistic. For a given sample size n, however, the power of the test can increase dramatically over the MMD l test, even for moderate block sizes B, making much better use of the available data with only a small increase in computation.

The block averaging scheme was originally proposed in [START_REF] Ho | Two-stage U-statistics for hypothesis testing[END_REF], as an instance of a two-stage Ustatistic, to be applied when the degree of degeneracy of the U-statistic is indeterminate. Differences with respect to our method are that Ho and Shieh compute the block statistics by sampling with replacement [13, (b) p. 863], and propose to obtain the variance of the test statistic via Monte Carlo, jackknife, or bootstrap techniques, whereas we use closed form expressions. Ho and Shieh further suggest an alternative two-stage U-statistic in the event that the degree of degeneracy is known; we return to this point in the discussion. While we confine ourselves to the MMD in this paper, we emphasize that the block approach applies to a much broader variety of test situations where the null distribution cannot easily be computed, including the energy distance and distance covariance [START_REF] Sejdinovic | Hypothesis testing using pairwise distances and associated kernels[END_REF][START_REF] Baringhaus | On a new multivariate two-sample test[END_REF][START_REF] Székely | Testing for equal distributions in high dimension[END_REF] and Fisher statistic [START_REF] Harchaoui | Testing for homogeneity with kernel Fisher discriminant analysis[END_REF] in the case of two-sample testing, and the Hilbert-Schmidt Independence Criterion [START_REF] Gretton | A kernel statistical test of independence[END_REF] and distance covariance [START_REF] Székely | Measuring and testing dependence by correlation of distances[END_REF] for independence testing. Finally, the kernel learning approach of [START_REF] Gretton | Optimal kernel choice for large-scale two-sample tests[END_REF] applies straightforwardly, allowing us to maximize test power over a given kernel family. Code is available at http://github.com/wojzaremba/btest.

Theory

In this section we describe the mathematical foundations of the B-test. We begin with a brief review of kernel methods, and of the maximum mean discrepancy. We then present our block-based average MMD statistic, and derive its distribution under the H 0 (P = Q) and H A (P = Q) hypotheses. The central idea employed in the construction of the B-test is to generate a low variance MMD estimate by averaging multiple low variance kernel statistics computed over blocks of samples. We show simple sufficient conditions on the block size for consistency of the estimator. Furthermore, we analyze the properties of the finite sample estimate, and propose a consistent strategy for setting the block size as a function of the number of samples.

Definition and asymptotics of the block-MMD

Let F k be an RKHS defined on a topological space X with reproducing kernel k, and P a Borel probability measure on X . The mean embedding of P in F k , written µ k (p) ∈ F k is defined such In both plots, the number of samples is fixed at 500. As we vary B, we trade off the quality of the finite sample Gaussian approximation to the null distribution, as in Theorem 2.3, with the variances of the H0 and HA distributions, as outlined in Section 2.1. In (b) the distribution under H0 does not resemble a Gaussian (it does not pass a level 0.05 Kolmogorov-Smirnov (KS) normality test [START_REF] Kolmogorov | Sulla determinazione empirica di una legge di distribuzione[END_REF][START_REF] Smirnov | Table for estimating the goodness of fit of empirical distributions[END_REF]), and a Gaussian approximation results in a conservative test threshold (vertical green line). The remaining empirical distributions all pass a KS normality test.

that

E x∼p f (x) = f, µ k (p) F k for all f ∈ F k ,
and exists for all Borel probability measures when k is bounded and continuous [START_REF] Berlinet | Reproducing Kernel Hilbert Spaces in Probability and Statistics[END_REF][START_REF] Gretton | A kernel two-sample test[END_REF]. The maximum mean discrepancy (MMD) between a Borel probability measure P and a second Borel probability measure Q is the squared RKHS distance between their respective mean embeddings,

η k (P, Q) = µ k (P ) -µ k (Q) 2 F k = E xx k(x, x ) + E yy k(y, y ) -2E xy k(x, y), (1)
where x denotes an independent copy of x [START_REF] Gretton | A kernel method for the two-sample-problem[END_REF]. Introducing the notation z = (x, y), we write

η k (P, Q) = E zz h k (z, z ), h(z, z ) = k(x, x ) + k(y, y ) -k(x, y ) -k(x , y). (2)
When the kernel k is characteristic, then η k (P, Q) = 0 iff P = Q [START_REF] Sriperumbudur | Hilbert space embeddings and metrics on probability measures[END_REF]. Clearly, the minimum variance unbiased estimate MMD u of η k (P, Q) is a U-statistic.

By analogy with MMD u , we make use of averages of h(x, y, x , y ) to construct our two-sample test. We denote by ηk (i) the ith empirical estimate MMD u based on a subsample of size B, where 1 ≤ i ≤ n B (for notational purposes, we will index samples as though they are presented in a random fixed order). More precisely,

ηk (i) = 1 B(B -1) iB a=(i-1)B+1 iB b=(i-1)B+1,b =a h(z a , z b ). (3)
The B-test statistic is an MMD estimate obtained by averaging the ηk (i). Each ηk (i) under H 0 converges to an infinite sum of weighted χ 2 variables [START_REF] Gretton | A fast, consistent kernel twosample test[END_REF]. Although setting B = n would lead to the lowest variance estimate of the MMD, computing sound thresholds for a given p-value is expensive, involving repeated bootstrap sampling [START_REF] Efron | An Introduction to the Bootstrap[END_REF][START_REF] Lloyd | Continuous univariate distributions[END_REF], or computing the eigenvalues of a Gram matrix [START_REF] Gretton | A fast, consistent kernel twosample test[END_REF].

In contrast, we note that ηk (i) i=1,..., n B are i.i.d. variables, and averaging them allows us to apply the central limit theorem in order to estimate p-values from a normal distribution. We denote the average of the ηk (i) by ηk ,

ηk = B n n B i=1 ηk (i). (4)
We would like to apply the central limit theorem to variables ηk (i) i=1,..., n

B

. It remains for us to derive the distribution of ηk under H 0 and under H A . We rely on the result from [START_REF] Gretton | A kernel method for the two-sample-problem[END_REF]Theorem 8] for H A . According to our notation, for every i, Theorem 2.1 Assume 0 < E(h 2 ) < ∞, then under H A , ηk converges in distribution to a Gaussian according to

B 1 2 (η k (i) -MMD 2 ) D → N (0, σ 2 u ), (5)
where

σ 2 u = 4 E z [(E z h(z, z )) 2 -E z,z (h(z, z ))] 2 .
This in turn implies that

ηk (i) D → N (MMD 2 , σ 2 u B -1 ). (6)
For an average of {η k (i)} i=1,..., n B , the central limit theorem implies that under

H A , ηk D → N MMD 2 , σ 2 u (Bn/B) -1 = N MMD 2 , σ 2 u n -1 . (7)
This result shows that the distribution of H A is asymptotically independent of the block size, B.

Turning to the null hypothesis, [START_REF] Gretton | A kernel method for the two-sample-problem[END_REF]Theorem 8] additionally implies that under H 0 for every i,

Theorem 2.2 B ηk (i) D → ∞ l=1 λ l [z 2 l -2], (8)
where z l ∼ N (0, 2) 2 i.i.d, λ l are the solutions to the eigenvalue equation

X k(x, x )ψ l (x)dp(x) = λ l ψ l (x ), (9)
and k(

x i , x j ) := k(x i , x j ) -E x k(x i , x) -E x k(x, x j ) + E x,x k(x, x ) is the centered RKHS kernel.
As a consequence, under H 0 , ηk (i) has expected variance 2B -2 ∞ l=1 λ 2 . We will denote this variance by CB -2 . The central limit theorem implies that under H 0 ,

ηk D → N 0, C B 2 n/B -1 = N 0, C(nB) -1 (10)
The asymptotic distributions for ηk under H 0 and H A are Gaussian, and consequently it is easy to calculate the distribution quantiles and test thresholds. Asymptotically, it is always beneficial to increase B, as the distributions for η under H 0 and H A will be better separated. For consistency, it is sufficient to ensure that n/B → ∞.

A related strategy of averaging over data blocks to deal with large sample sizes has recently been developed in [START_REF] Kleiner | A scalable bootstrap for massive data[END_REF], with the goal of efficiently computing bootstrapped estimates of statistics of interest (e.g. quantiles or biases). Briefly, the approach splits the data (of size n) into s subsamples each of size B, computes an estimate of the n-fold bootstrap on each block, and averages these estimates. The difference with respect to our approach is that we use the asymptotic distribution of the average over block statistics to determine a threshold for a hypothesis test, whereas [START_REF] Kleiner | A scalable bootstrap for massive data[END_REF] is concerned with proving the consistency of a statistic obtained by averaging over bootstrap estimates on blocks.

Convergence of Moments

In this section, we analyze the convergence of the moments of the B-test statistic, and comment on potential sources of bias.

The central limit theorem implies that the empirical mean of {η k (i)} i=1,..., n B converges to E(η k (i)). Moreover it states that the variance {η k (i)} i=1,..., n B converges to

E(η k (i)) 2 -E(η k (i) 2 )
. Finally, all remaining moments tend to zero, where the rate of convergence for the jth moment is of the order n B j+1 2 [START_REF] Von | On the convergence of moments in the central limit theorem[END_REF]. This indicates that the skewness dominates the difference of the distribution from a Gaussian.

Under both H 0 and H A , thresholds computed from normal distribution tables are asymptotically unbiased. For finite samples sizes, however, the bias under H 0 can be more severe. From Equation [START_REF] Gretton | A kernel statistical test of independence[END_REF] we have that under H 0 , the summands, ηk (i), converge in distribution to infinite weighted sums of χ 2 distributions. Every unweighted term of this infinite sum has distribution N (0, 2) 2 , which has finite skewness equal to 8. The skewness for the entire sum is finite and positive,

C = ∞ l=1 8λ 3 l , (11)
as λ l ≥ 0 for all l due to the positive definiteness of the kernel k. The skew for the mean of the ηk (i) converges to 0 and is positively biased. At smaller sample sizes, test thresholds obtained from the standard Normal table may therefore be inaccurate, as they do not account for this skew. In our experiments, this bias caused the tests to be overly conservative, with lower Type I error than the design level required (Figures 2 and5).

Finite Sample Case

In the finite sample case, we apply the Berry-Esséen theorem, which gives conservative bounds on the ∞ convergence of a series of finite sample random variables to a Gaussian distribution [START_REF] Andrew | The accuracy of the gaussian approximation to the sum of independent variates[END_REF].

Theorem 2.3 Let X 1 , X 2 , . . . , X n be i.i.d. variables. E(X 1 ) = 0, E(X 2 1 ) = σ 2 > 0, and

E(|X 1 | 3 ) = ρ < ∞. Let F n be a cumulative distribution of n i=1
Xi √ nσ , and let Φ denote the standard normal distribution. Then for every x,

|F n (x) -Φ(x)| ≤ Cρσ -3 n -1/2 , ( 12
)
where C < 1.

This result allows us to ensure fast point-wise convergence of the B-test. We have that ρ(η k ) = O(1), i.e., it is dependent only on the underlying distributions of the samples and not on the sample size. The number of i.i.d. samples is nB -1 . Based on Theorem 2.3, the point-wise error can be upper bounded by

O(1) O(B -1 ) 3 2 √ n B = O( B 2 √ n ) under H A .
Under H 0 , the error can be bounded by

O(1) O(B -2 ) 3 2 √ n B = O( B 3.5 √ n ).
While the asymptotic results indicate that convergence to an optimal predictor is fastest for larger B, the finite sample results support decreasing the size of B in order to have a sufficient number of samples for application of the central limit theorem. As long as B → ∞ and n B → ∞, the assumptions of the B-test are fulfilled.

By varying B, we make a fundamental tradeoff in the construction of our two sample test. When B is small, we have many samples, hence the null distribution is close to the asymptotic limit provided by the central limit theorem, and the Type I error is estimated accurately. The disadvantage of a small B is a lower test power for a given sample size. Conversely, if we increase B, we will have a lower variance empirical distribution for H 0 , hence higher test power, but we may have a poor estimate of the number of Type I errors (Figure 1). A sensible family of heuristics therefore is to set

B = [n γ ] ( 13
)
for some 0 < γ < 1, where we round to the nearest integer. In this setting the number of samples available for application of the central limit theorem will be [n (1-γ) ]. For given γ computational complexity of the B-test is O n 1+γ . We note that any value of γ ∈ (0, 1) yields a consistent estimator. We have chosen γ = 1 2 in the experimental results section, with resulting complexity O n 1.5 : we emphasize that this is a heuristic, and just one choice that fulfils our assumptions.

Experiments

We have conducted experiments on challenging synthetic and real datasets in order to empirically measure (i) sample complexity, (ii) computation time, and (iii) Type I / Type II errors. We evaluate B-test performance in comparison to the MMD l and MMD u estimators, where for the latter we compare across different strategies for null distribution quantile estimation. 1: Sample complexity for tests on the distributions described in Figure 3. The fourth column indicates the minimum number of samples necessary to achieve Type I and Type II errors of 5%. The fifth column is the computation time required for 2000 samples, and is not presented for settings that have unsatisfactory sample complexity.

Synthetic data

Following previous work on kernel hypothesis testing [START_REF] Gretton | Optimal kernel choice for large-scale two-sample tests[END_REF], our synthetic distributions are 5 × 5 grids of 2D Gaussians. We specify two distributions, P and Q. For distribution P each Gaussian has identity covariance matrix, while for distribution Q the covariance is non-spherical. Samples drawn from P and Q are presented in Figure 3. These distributions have proved to be very challenging for existing non-parametric two-sample tests [START_REF] Gretton | Optimal kernel choice for large-scale two-sample tests[END_REF]. We employed three different kernel selection strategies in the hypothesis test. First, we used a Gaussian kernel with σ = 1, which approximately matches the scale of the variance of each Gaussian in mixture P . While this is a somewhat arbitrary default choice, we selected it as it performs well in practice (given the lengthscale of the data), and we treat it as a baseline. Next, we set σ equal to the median pairwise distance over the training data, which is a standard way to choose the Gaussian kernel bandwidth [START_REF] Schölkopf | Support vector learning[END_REF], although it is likewise arbitrary in this context. Finally, we applied a kernel learning strategy, in which the kernel was optimized to maximize the test power for the alternative P = Q [START_REF] Gretton | Optimal kernel choice for large-scale two-sample tests[END_REF]. This approach returned a non-negative linear combination combination of base kernels, where half the data were used in learning the kernel weights (these data were excluded from the testing phase).

The base kernels in our experiments were chosen to be Gaussian, with bandwidths in the set σ ∈ {2 -15 , 2 -14 , . . . , 2 10 }. Testing was conducted using the remaining half of the data. As B grows, the Type II error drops quickly when the kernel is appropriately chosen. The kernel selection method is described in [START_REF] Gretton | Optimal kernel choice for large-scale two-sample tests[END_REF], and closely approximates the baseline performance of the well-informed user choice of σ = 1.

For comparison with the quadratic time Ustatistic MMD u [START_REF] Gretton | A fast, consistent kernel twosample test[END_REF][START_REF] Gretton | A kernel two-sample test[END_REF], we evaluated four null distribution estimates: (i) Pearson curves, (ii) gamma approximation, (iii) Gram matrix spectrum, and (iv) bootstrap. For methods using Pearson curves and the Gram matrix spectrum, we drew 500 samples from the null distribution estimates to obtain the 1 -α quantiles, for a test of level α. For the bootstrap, we fixed the number of shuffles to 1000. We note that Pearson curves and the gamma approximation are not statistically consistent. We considered only the setting with σ = 1 and σ set to the median pairwise distance, as kernel selection is not yet solved for tests using MMD u [START_REF] Gretton | Optimal kernel choice for large-scale two-sample tests[END_REF].

In the first experiment we set the Type I error to be 5%, and we recorded the Type II error. We conducted these experiments on 2000 samples over 1000 repetitions, with varying block size, B. Figure 4 presents results for different kernel choice strategies, as a function of B. The median heuristic performs extremely poorly in this experiment. As discussed in [9, Section 5], the reason for this failure is that the lengthscale of the difference between the distributions P and Q differs from the lengthscale of the main data variation as captured by the median, which gives too broad a kernel for the data.

In the second experiment, our aim was to compare the empirical sample complexity of the various methods. We again fixed the same Type I error for all methods, but this time we also fixed a Type II error of 5%, increasing the number of samples until the latter error rate was achieved. Column four of Table 1 shows the number of samples required in each setting to achieve these error rates. We additionally compared the computational efficiency of the various methods. The computation time for each method with a fixed sample size of 2000 is presented in column five of Table 1. All experiments were run on a single 2.4 GHz core.

Finally, we evaluated the empirical Type I error for α = 5% and increasing B. Figure 2 displays the empirical Type I error, where we note the location of the γ = 0.5 heuristic in Equation [START_REF] Ho | Two-stage U-statistics for hypothesis testing[END_REF]. For the user-chosen kernel (σ = 1, Figure 2(a)), the number of Type I errors closely matches the targeted test level. When median heuristic is used, however, the test is overly conservative, and makes fewer Type I errors than required (Figure 2(b)). This indicates that for this choice of σ, we are not in the asymptotic regime, and our Gaussian null distribution approximation is inaccurate. Kernel selection via the strategy of [START_REF] Gretton | Optimal kernel choice for large-scale two-sample tests[END_REF] alleviates this problem (Figure 2(c)). This setting coincides with a block size substantially larger than 2 (MMD l ), and therefore achieves lower Type II errors while retaining the targeted Type I error.

Musical experiments

In this set of experiments, two amplitude modulated Rammstein songs were compared (Sehnsucht vs. Engel, from the album Sehnsucht). Following the experimental setting in [9, Section 5], samples from P and Q were extracts from AM signals of time duration 8.3 × 10 -3 seconds in the original audio. Feature extraction was identical to [START_REF] Gretton | Optimal kernel choice for large-scale two-sample tests[END_REF], except that the amplitude scaling parameter was set to 0.3 instead of 0.5. As the feature vector had size 1000 we set the block size B = √ 1000 = 32. Table 2 summarizes the empirical Type I and Type II errors over 1000 repetitions, and the average computation times. Figure 5 shows the average number of Type I errors as a function of B: in this case, all kernel selection strategies result in conservative tests (lower Type I error than required), indicating that more samples are needed to reach the asymptotic regime. Figure 1 shows the empirical H 0 and H A distributions for different B.

Discussion

We have presented experimental results both on a difficult synthetic problem, and on real-world data from amplitude modulated audio recordings. The results show that the B-test has a much better

B = 2 .Figure 1 :

 21 Figure1: Empirical distributions under H0 and HA for different regimes of B for the music experiment (Section 3.2). In both plots, the number of samples is fixed at 500. As we vary B, we trade off the quality of the finite sample Gaussian approximation to the null distribution, as in Theorem 2.3, with the variances of the H0 and HA distributions, as outlined in Section 2.1. In (b) the distribution under H0 does not resemble a Gaussian (it does not pass a level 0.05 Kolmogorov-Smirnov (KS) normality test[START_REF] Kolmogorov | Sulla determinazione empirica di una legge di distribuzione[END_REF][START_REF] Smirnov | Table for estimating the goodness of fit of empirical distributions[END_REF]), and a Gaussian approximation results in a conservative test threshold (vertical green line). The remaining empirical distributions all pass a KS normality test.

Figure 2 :

 2 Figure 2: Type I errors on the distributions shown in Figure 3 for α = 5%: (a) MMD, single kernel, σ = 1, (b) MMD, single kernel, σ set to the median pairwise distance, and (c) MMD, non-negative linear combination of multiple kernels. The experiment was repeated 30000 times. Error bars are not visible at this scale.

Figure 3 :

 3 Figure 3: Synthetic data distributions P and Q. Samples belonging to these classes are difficult to distinguish.

Figure 4 :

 4 Figure4: Synthetic experiment: number of Type II errors vs B, given a fixed probability α of Type I errors. As B grows, the Type II error drops quickly when the kernel is appropriately chosen. The kernel selection method is described in[START_REF] Gretton | Optimal kernel choice for large-scale two-sample tests[END_REF], and closely approximates the baseline performance of the well-informed user choice of σ = 1.

Acknowledgments We thank Mladen Kolar for helpful discussions. This work is partially funded by ERC Grant 259112, and by the Royal Academy of Engineering through the Newton Alumni Scheme.

sample complexity than MMD l over all tested kernel selection strategies. Moreover, it is an order of magnitude faster than any test that consistently estimates the null distribution for MMD u (i.e., the Gram matrix eigenspectrum and bootstrap estimates): these estimates are impractical at large sample sizes, due to their computational complexity. Additionally, the B-test remains statistically consistent, with the best convergence rates achieved for large B. The B-test combines the best features of MMD l and MMD u based two-sample tests: consistency, high statistical efficiency, and high computational efficiency.

A number of further interesting experimental trends may be seen in these results. First, we have observed that the empirical Type I error rate is often conservative, and is less than the 5% targeted by the threshold based on a Gaussian null distribution assumption (Figures 2 and5). In spite of this conservatism, the Type II performance remains strong (Tables 1 and2), as the gains in statistical power of the B-tests improve the testing performance (cf. Figure 1). Equation [START_REF] Gretton | A fast, consistent kernel twosample test[END_REF] implies that the size of B does not influence the asymptotic variance under H A , however we observe in Figure 1 that the empirical variance of H A drops with larger B. This is because, for these P and Q and small B, the null and alternative distributions have considerable overlap. Hence, given the distributions are effectively indistinguishable at these sample sizes n, the variance of the alternative distribution as a function of B behaves more like that of H 0 (cf. Equation [START_REF] Gretton | A kernel two-sample test[END_REF]). This effect will vanish as n grows.

Finally, [START_REF] Ho | Two-stage U-statistics for hypothesis testing[END_REF] propose an alternative approach for U-statistic based testing when the degree of degeneracy is known: a new U-statistic (the TU-statistic) is written in terms of products of centred U-statistics computed on the individual blocks, and a test is formulated using this TU-statistic. Ho and Shieh show that a TU-statistic based test can be asymptotically more powerful than a test using a single U-statistic on the whole sample, when the latter is degenerate under H 0 , and nondegenerate under H A . It is of interest to apply this technique to MMD-based two-sample testing.


come

Mechanics of the discrete and continuous systems in the light of the non-standard analysis


Boris Lesner

Bruno Zanuttini

Efficient Policy Construction for MDPs Represented in Probabilistic PDDL

We present a novel dynamic programming approach to computing optimal policies for Markov Decision Processes compactly represented in grounded Probabilistic PDDL. Unlike other approaches, which use an intermediate representation as Dynamic Bayesian Networks, we directly exploit the PPDDL description by introducing dedicated backup rules. This provides an alternative approach to DBNs, especially when actions have highly correlated effects on variables. Indeed, we show significant improvements on several planning domains from the International Planning Competition. Finally, we exploit the incremental flavor of our backup rules for designing promising approaches to policy revision.

Introduction

Markov Decision Processes have become a standard framework for modelling decision-theoretic planning tasks. While classical approches rely on dynamic programming algorithms such as Value Iteration, the AI community quickly focused on factored approches to solve MDPs, since practical problems involves millons or even billions of states. Central to such methods are avoiding to enumerate all states and exploiting the structure present in most practical problems. Factored MDP algorithms represent states and actions in a symbolic way. States are described using valuations over a set of variables and actions describe changes on variables.

The Probabilistic Planning Domain Description Language [START_REF] Younes | PPDDL1.0: An Extension to PDDL for Expressiong Planning Domains with Probabilistic Effects[END_REF]) is a widely adopted language for specifying probabilistic and decision-theoretic planning problems. However, state-of-the-art approaches for computing optimal policies for MDPs rely on a representation by Dynamic Bayesian Networks (DBNs), and therefore, though very efficient on most problems, they require a translation from PPDDL.

We propose a new approach for computing optimal policies for MDPs, which directly uses the PPDDL specification. Our approach relies on classical dynamic programming, but performs action backup using the new notion of a "frameless" action-value. Using such backups in classical value iteration yields a new algorithm, which turns out to be complementary to other factored approaches. We Preliminaries Markov Decision Processes A Markov Decision Process M = (S, A, T, R) involves finite sets of states S and actions A. For s, s ′ ∈ S and a ∈ A, the transition probability T (s ′ |s, a) denotes the probability of reaching state s ′ after taking action a in s. For states s, s ′ and action a, the reward function R(s, a, s ′ ) determines the payoff of taking action a in state s and ending up in state s ′ .

Given an MDP M , one can compute an optimal policy π h : S → A which maximizes expected reward at some given horizon h. π h can be derived from the optimal value function at h, written V h (s) = max a∈A Q a V h-1 (s) with

Q a V (s) = s ′ ∈S
T (s ′ |s, a) (R(s, a, s ′ ) + γV (s ′ )) (1) and V 0 (s) = 0 for all s ∈ S. The optimal policy π h is then given by π h (s) = argmax a Q a V h-1 (s). The discount factor γ ∈ [0, 1], reduces the importance of future rewards. Infinite-horizon value functions can be approximated arbitrarily closely when γ < 1 : let V * = lim h→∞ V h denote the infinite horizon value function, then as soon as

||V h -V h-1 || ∞ ≤ ǫ(1-γ) 2γ holds, ||V h -V * || ∞ ≤ ǫ holds.
PPDDL The Probabilistic Planning Domain Description Language [START_REF] Younes | PPDDL1.0: An Extension to PDDL for Expressiong Planning Domains with Probabilistic Effects[END_REF] allows to model probabilistic and decision-theoretic planning problems. We present here its components, domains and problems.

Planning domains consist of predicates and action schemata. Predicates encode Boolean state variables1 . Action schemata represent both the transition and reward functions. Each action schemata a consists in a precondition, that is, a formula characterizing the states where a is applicable, and an effect. The effect describes the changes that may be applied to a state when taking action a. Importantly, PPDDL actions obey the frame assumption: variables not explicitely modified by an action remains the same after taking the action. PPDDL1.0 defines the following types of effects: A planning problem, on a given domain, defines a set of objects used to instanciate predicates and action schemata. In turn, instanciated predicates encode a set X of propositional state variables and instanciated action schemata define a set of actions A, each action a having precondition p a and effect e a defined over X. A problem may also define a formula Γ characterizing the set of goal states associated with a reward R Γ . If the problem has no goal (like classical MDPs), we assume Γ = ⊥. Goal-oriented problems can be easily dealt with in the MDP framework, by ensuring the following for any goal state g |= Γ : no action can be taken in g and for any horizon h, V h (g) = R Γ . A problem may also define an initial state, but since we focus on the more general policy construction problem, the initial state is ignored. As an example, Figure 1 presents the PPDDL description of the move action of the well-known COFFEE problem.

For t, ℓ two sets of literals, write t ⊲ ℓ for the set of literals matching t, and completed as in ℓ over other variables, that is, t ⊲ ℓ = t ∪ (ℓ \ {p | p ∈ t}); e.g., xy ⊲ xȳzt = xyzt. Given a state s (a full instanciation of states variables X), each effect e described by a PPDDL problem induces a probability distribution D(e, s) over a set of pairs (t, r) (Younes and Littman 2004, Sec. 4.1). Each such pair consists in a "basic effet" t, that is, a consistent set of literals modified in state s to obtain next state s ′ = t ⊲ s, and an immediate reward r. Writing (t, r) : p for "effect (t, r) occurs with probability p", D(e, s) is recursively defined as follows:

• D(x, s) = {(x, 0) : 1} • D(not x, s) = {(x, 0) : 1} • D(increase (reward) r, s) = {(∅, r) : 1} • D(decrease (reward) r, s) = {(∅, -r) : 1} • D(when φ e, s) = D(e, s) if s |= φ otherwise{(∅, 0) : 1} • D(prob p 1 e 1 . . . , s) = i {(t, r) : p i p | (t, r) : p ∈ D(e i , s)} • D(and , s) = {(∅, 0) : 1} • D(and e 1 e 2 . . . , s) = {(t ′ ∪ t, r ′ + r) : p ′ p | (t ′ , r ′ ) : p ′ ∈ D(e 1 , s
. Let V 1 , V 2 : {0, 1} n → R, r ∈ R, s ∈ {0, 1} n ,
φ a formula on {x 1 , . . . , x n } and t a conjunction of literals:

(V 1 op V 2 )(s) = V 1 (s) op V 2 (s) for op ∈ {+, ×, -} max(V 1 , V 2 )(s) = max(V 1 (s), V 2 (s)) (V 1 + r)(s) = V 1 (s) + r (φ × V 1 )(s) = V 1 (s) if s |= φ, 0 otherwise ITE(φ, V 1 , V 2 )(s) = V 1 (s) if s |= φ, V 2 (s) otherwise (V 1 ) [t] (s) = V 1 (t ⊲ s) (∃x i V 1 )(s) = (V 1 ) [xi] (s) + (V 1 ) [xi] (s)
In the factored MDP litterature (see below), ADDs revealed to be candidate of choice for representing structured transitions matrices and value functions, since in practical problems many states have the same value, allowing memory savings and fast value function manipulation.

Related Work

Most work on computing policies for factored MDPs has focused on DBN representations. Such representations use conditional probability tables (CPTs), which give the probability of each x ′ being made true or false by each action, depending on the values of its parents in the DBN. PPDDL descriptions can be converted to DBNs [START_REF] Younes | PPDDL1.0: An Extension to PDDL for Expressiong Planning Domains with Probabilistic Effects[END_REF], and CPTs can be compactly represented using ADDs [START_REF] Hoey | SPUDD: Stochastic planning using decision diagrams[END_REF][START_REF] St-Aubin | APRICODD: Approximate policy construction using decision diagrams[END_REF][START_REF] Feng | Symbolic heuristic search for factored Markov decision processes[END_REF][START_REF] Feng | Symbolic generalization for on-line planning[END_REF][START_REF] Guestrin | Effcient Solution Algorithms for Factored MDPs[END_REF]. Then action regression is performed by simple ADD operations.

Though this representation is natural when variables are independent, correlated effects such as (prob 0.5 xy 0.5 xȳ) introduce dependencies between postaction values. The usual trick is to introduce auxiliary variables whose values encode which effect occurred, and deterministic dependencies from them to post-action variables (Younes and Littman 2004, Sec. 5), but naturally, this increases the size of diagrams manipulated. Another problem is independent conditions on the same variable, e.g., (and (when x 1 (prob p 1 x)) (when x 2 (prob p 2 x))), which make CPTs blow up if auxiliary variables are not (:action move :effect (and (when (in-office) (probabilistic 0.9 (not (in-office)))) (when (not (in-office)) (probabilistic 0.9 (in-office))) (when (and (raining) (not (has-umbrella))) (probabilistic 0.9 (is-wet) 0.1 (when (not (is-wet)) (increase (reward) 0.2)))) (when (or (not (raining)) (has-umbrella)) (when (not (is-wet)) (increase (reward) 0.2))) (when (user-has-coffee) (increase (reward) 0.8))))

Figure 1: PPDDL descritpion of the move action of the COFFEE problem. The three first "when" clauses read: with probability 0.9 the agent will change location; concurrently, starting from a state where it rains and the agent has no umbrella, with probability 0.9 the agent will be wet in next state; otherwise, if it is not already wet, it will remain dry and receive a 0.2 reward. used. Nevertheless, algorithms using such representations can be very efficient, both in the exact [START_REF] Hoey | SPUDD: Stochastic planning using decision diagrams[END_REF] and approximative (St-Aubin, [START_REF] St-Aubin | APRICODD: Approximate policy construction using decision diagrams[END_REF][START_REF] Guestrin | Effcient Solution Algorithms for Factored MDPs[END_REF] cases. Some work has also been done at the relational level [START_REF] Lang | Planning with Noisy Probabilistic Relational Rules[END_REF].

A somehow intermediate representation between PPDDL and DBNs is Probabilistic STRIPS Operators [START_REF] Kushmerick | An algorithm for probabilistic planning[END_REF]. Under this representation, actions have mutually exclusive conditions, each associated with a probabilistic effect of the form

(prob p 1 e 1 . . . p k e k )
where each e i is a conjunction of literals. As far as we know, no work has focused on computing optimal MDP policies for such representations in the grounded case. Nevertheless, they are directly used by algorithms which work at the relational level [START_REF] Boutilier | Practical Solution Techniques for First-Order MDPs[END_REF].

When using ADDs, DBN representation can be exponentially more succinct than PPDDL. Conversely, because of nested effects, PPDDL may be exponentially more succinct than DBNs or PSOs. For more details about these issues, we refer the reader to [START_REF] Boutilier | Decision theoretic planning: Structural assumptions and computational leverage[END_REF][START_REF] Rintanen | Expressive Equivalence of Formalism for Planning with Sensing[END_REF][START_REF] Younes | PPDDL1.0: An Extension to PDDL for Expressiong Planning Domains with Probabilistic Effects[END_REF].

Updating Value Functions

We now come to our new approach for exact backup of value functions in MDPs, which we integrate in value iteration.

F-values

The process of action backup refers to the task of computing the value Q a V from a value function V . We propose to perform this task incrementally with respect to a PPDDL description of actions, using an intermediate representation of values as frameless action-values. As other approaches, we use ADDs for efficient storage and manipulation of real-valued functions of Boolean variables. Definition 1 Let V : S → R be a value function and e be an effect. The frameless action-value (F -value for short) of e with respect to V is the function F e V : S ×2 X → R which for s ∈ S and ℓ ′ ∈ 2 X , gives the expected reward of applying e in s and forcing the values of ℓ ′ on unaffected variables:

F e V (s, ℓ ′ ) = E (t,r)∼D(e,s) [r + V (t ⊲ ℓ ′ )]
Intuitively, F -values are just like usual action-values Q a V , but they make no frame assumption on the unmodified variables. As we will see, this allows to handle conjunctive effects in an incremental fashion.

Example 1 Let x be a variable and

V (x) = 0, V (x) = 1. Then the F -value F e V of e = (when x (not x)) wrt V is given by F e V (x, •) = 1, F e V (x, x ′ ) = 0, F e V (x, x′ ) = 1.
In words, for s = x, the F -value leaves the possibility that x persists, but also that it changes to x ′ . This is to be compared with the action-value Q e V (•) = 1. Another example is given by the right branch of the F -value depicted on Figure 2 (d) for e = (when x (prob 0.6 x 0.4 y)): for s satisfying x and ℓ ′ = x′ ȳ′ , the expected value is 5.2. Indeed, in s, either:

• effect x occurs, and y is not affected; hence the value of y ′ is taken from ℓ ′ , resulting in state x ′ ȳ′ , which has value 2 (Figure 2 upper left); this occurs with probability 0.6; • or y occurs, resulting in state x′ y ′ with value 10; this occurs with probability 0.4.

Hence the F -value of e in s, with unaffected variables modified as in ℓ ′ , is 0.6 × 2 + 0.4 × 10 = 5.2.

An important feature of F -values is that they embed actionvalues. Indeed, when discounting and enforcing the frame assumption we get (see the proof in Theorem 2):

Q a V (s) = F ea γV (s, s) (2)
Example 2 Considering Figure 2 (d) and s = xȳ, Q e V (s) is retrieved by evaluating the ADD in xȳx ′ ȳ′ , yielding 3.2, which is indeed the expected value of e in s (probability 0.6 of setting x to true, resulting in state x ′ ȳ′ with value 2, and 0.4 of setting y to true, resulting in x ′ y ′ with value 5).

Backup rules Our rules take as input an F -value V ′ and an effect e described in PPDDL, and produce a new F -value B V ′ (e) which accounts for the application of e on V ′ . Since F -values are functions of S × 2 X , we introduce a second set of primed variables X ′ = {x ′ | x ∈ X} to represent them with ADDs. B V ′ (e) is recursively defined as follows:

R1.1 B V ′ (x) = V ′ [x ′ ] R1.2 B V ′ (not x) = V ′ [x ′ ] R2.1 B V ′ (increase (reward) r) = V ′ + r R2.2 B V ′ (decrease (reward) r) = V ′ -r R3 B V ′ (when φ e) = ITE(φ, B V ′ (e), V ′ ) R4 B V ′ (prob p 1 e 1 . . . p k e k ) = k i=1 p i × B V ′ (e i ) R5.1 B V ′ (and) = V ′ R5.2 B V ′ (and e 1 e 2 . . . e k ) = B B V ′ (e1) (and e 2 . . . e k ) Theorem 1 Given a, V , let V ′ (•, s) = V (s) be the primed version of V . Then the backup rule B V ′ (e a ) computes F ea V . Proof.
We show by induction on the PPDDL description of e that for any function W

: S × 2 X → R, B W (e)(s, ℓ ′ ) is E(t,r)∈D(e,s) [r + W (s, t ⊲ ℓ ′ )].
The result follows because V ′ depends on its second argument only, hence V ′ (s, t⊲ℓ ′ ) = V (t⊲ℓ ′ ). For brevity, we write E(t,r) [•] for E(t,r)∼D(e,s) [•].

Let e = x.

Then B W (e)(s, ℓ ′ ) = W [x ′ ] (s, ℓ ′ ) (Rule 1.1). Because D(e, s) = {(x, 0) : 1}, this is W (s, {x ′ } ⊲ ℓ ′ ) = E(t,r) [r + W (s, t ⊲ ℓ ′ )], as desired. Case (not x) is dual. For e = (increase (reward) r), B W (e)(s, ℓ ′ ) is r + W (s, ℓ ′ ). Since D(e, s) = {(∅, r) : 1}, this is indeed E(t,r) [r + W (s, t ⊲ ℓ ′ )]. Case decrease is dual. Now let e = (when φ e ′ ). If s |= φ, then B W (e)(s, ℓ ′ ) is B W (e ′ )(s, ℓ ′ ), which is correct by the induction hypothesis (IH). If s |= φ, we get B W (e)(s, ℓ ′ ) = W (s, ℓ ′ ), which is E(t,r)∼D(e,s) [r+W (s, t ⊲ ℓ ′ )] since D(e, s) = {(∅, 0) : 1}.
For e = (prob p 1 e 1 . . . p k e k ), we have

B W (e)(s, ℓ ′ ) = i p i B W (e i )(s, ℓ ′ ). This is p i E(t i ,ri) [r i + W (s, t i ⊲ ℓ ′ )]
by IH, and correctness follows by linearity of expectations.

The (and) case is obvious. Finally, for e = (and e 1 e 2 ) we have

B W (e)(s, ℓ ′ ) = B B W (e1) (e 2 )(s, ℓ ′ ). This is E(t 2 ,r2) [r 2 + B W (e 1 )(s, t 2 ⊲ ℓ ′ )] by IH on e 2 , and E(t 2 ,r2) r 2 + E(t 1 ,r1) [r 1 + W (s, t 1 ⊲ (t 2 ⊲ ℓ ′ ))] by IH on e 1 . Now because effects are consistent we have t 1 ⊲ (t 2 ⊲ ℓ ′ ) = (t 1 ∪ t 2 ) ⊲ ℓ ′ , hence as desired we get E(t 1 ,r1),(t2,r2) [r 1 + r 2 + W (s, (t 1 ∪ t 2 ) ⊲ ℓ ′ )].
At this point, we described how to compute F -values F e V using ADDs. For a complete action backup, it remains to compute Q a V from F ea V . Equation 2 suggests to enforce the frame assumption on the ADD B V (e a ) by computing the ADD Persist(B V (e a )), where Persist is defined by:

Persist(F ) = ∃X ′ [x 1 ↔ x ′ 1 × • • • × x n ↔ x ′ n × F ] (3) Finally, given V ′ (•, s) = γV (s) (discounted,
primed version of value function V ), and accounting for the precondition p a and the goal Γ, we get:

Q a V = ITE(p a × ¬Γ, Persist(B V ′ (e a )
), -∞) Assigning -∞ to states not satisfying p a and to goal states, rules this action out when maximizing over actions. A complete value function update example is given on figure 2.

Putting it all together, we propose the algorithm RBAB (for "Rule-Based Action Backup"), whose pseudocode is given in Algorithm 1. RBAB performs symbolic value iteration to compute value functions and eventually a policy.

Theorem 2 RBAB computes an ǫ-optimal infinite-horizon policy.

Proof.

We show by induction on h that after the h-th iteration of the main loop, W is the optimal hsteps-to-go value function V h . Let a be an action. By

Algorithm 1: Rule-Based Action Backup (RBAB) V ← Γ × R Γ repeat W ← -∞; V ′ ← γ × prime variables(V ) foreach action a do Q ← ITE(p a × ¬Γ, Persist(B V ′ (e a )), -∞) W ← max(W, Q) W ← ITE(Γ, R Γ , W ) // Goal states → R Γ converged ← ||V -W || ∞ ≤ ǫ(1-γ) 2γ V ← W until converged // Extract policy π W ← -∞; V ′ ← γ × prime variables(V ) π ← 0 foreach action a i , i = 1, . . . , |A| do Q ← ITE(p a × ¬Γ, Persist(B V ′ (e a )), -∞) G ← Q > W // ∀s, G(s) ⇔ Q(s) > W (s) π ← max(π, G × i) W ← ITE(G, Q, W ) return π × ¬Γ ; // Goal states → action 0
Algorithm 2: Persist Input: N , an ADD with primed and unprimed variables Input: A, a set of literals to be forced (initially

∅) if N is a constant ADD then return N v ← var (N ) // Force literals. if v ∈ A then return Persist(Then(N ), A \ {v}) if v ∈ A then return Persist(Else(N ), A \ {v}) if (N, A) → R is in cache then return R x ←swap(v)// swap(v ′ )=v;swap(v)=v ′ T ← Persist(Then(N ), A ∪ {x}) E ← Persist(Else(N ), A ∪ {x}) // Merge over unprimed variable. if v is primed then u ← x else u ← v R ← ITE(u, T, E) insert (N, A) → R into cache return R
Theorem 1 and the induction hypothesis,

B V ′ (e a )(s, ℓ ′ ) is E(t,r)∼D(e a ,s) r + γV ′ h-1 (s, t ⊲ ℓ ′ )) . Since Persist forces x ′ ↔ x for all primed variables in B V ′ (e a ), we get V ′ h-1 (s, t ⊲ ℓ ′ ) = V h-1 (t ⊲ s).
Moreover, Persist abstracts primed variables away, hence the function Q computed by the algorithm depends on unprimed variables only, so that when s satisfies p a × ¬Γ we have

Q(s) = E(t,r)∼D(e a ,s) [r + γV h-1 (t ⊲ s)] = Q a V h-1 (
s), as desired. Observing that the remainder of the algorithm is classical Value Iteration concludes the proof.

Optimizations A critical step of Algorithm 1 is the computation of Persist. While Equation 3 directly maps to ADD operations, this requires several traversals of the ADD. To circumvent this, we use the dedicated implementation given as Algorithm 2, which requires a single pass.

One major drawback to ADDs is their sensitivity to variable ordering. Indeed, the size of the diagrams may increase dramatically if a wrong ordering is choosen, resulting in poor performance. Reordering algorithms provided by ADD packages may take some time to find a good, let alone optimal, variable ordering which reduces diagrams sizes, hence they must be used with caution. Experiments have shown that forcing x, x ′ to stay adjacent and using approximate reordering at the end of each iteration provides the best results.

Finally, note that the backup of probabilistic effects may benefit from a dedicated (possibly n-ary) multiply-and-add operator, but we haven't evaluated this empirically.

Empirical Results

RBAB was implemented as a client for the MDPsim server, allowing us to use domains from the International Planning Competition (IPC) and compare our approach with SPUDD.

Along with MDPsim is shipped MTBDDclient, an implementation of SPUDD which first translates PPDDL to DBNs. Since MTBDDclient is intended to be a demonstration program, we optimized it in two ways: first, a policy is built only in a last iteration; second, we introduced variable reordering heuristics to enhance its performance further. We tested two versions of SPUDD: the one called "1 by 1" considers variables one at a time (backup and abstraction) at each iteration, and "matrix" precomputes the "complete action diagram" (Hoey et al. 1999, Sec. 4.2).

We empirically compared the running times of all three solvers on domains taken from the 5th and 6th IPC (probabilistic track). Obviously enough, since we focused on the problem of finding optimal control policies, the hardest instances of the domains were out of reach for the solvers. We kept the problems where at least one solver could produce an ǫ-optimal discounted infinite-horizon policy with ǫ = 0.1 and γ = 0.9 in less than one hour.

Figure 3 shows the results on domains pitchcatch, rectange-tireworld, search-and-rescue, drive, and driveunrolled. Clearly, these results show the complementarity of the different solvers. Indeed, no algorithm is always better (or worse, for that matter) than any other.

On the pitchcatch domain, RBAB behaves better than "matrix", but "1 by 1" is still better. On the other hand, on the rectangle-tireworld domain RBAB outperforms both. The most conclusive results are on search-and-rescue, where RBAB appears to be exponentially faster than both. Moreover this problem is one of the two we tested (along with pitchcatch) for which the number of non exclusive conditions under an and statement, was greater than 1 and increasing linearily with the instance size. This suggests that this is a good indicator to expect good RBAB performance.

The results for the drive domains (drive, drive-unrolled and drive-unrolled2) are very interesting. Indeed, all three domains express the same traffic intersection problem, but in the "unrolled" version, actions are decomposed into a set of smaller, equivalent ones. On the drive instances RBAB is the worst algorithm, but its behavior clearly improves on the "unrolled" version, exploiting a smaller PPDDL descrip-tion of actions. Domain drive-unrolled2 is another variant which we designed, where we merged almost identical actions (look at light east with look at light west and look at light north with look at light south) into more general, but equivalent ones. The effects of these actions were also simplified by moving up redundant effects in the effect hierarchy. The resulting domain exploits as much as possible the succinctness of PPDDL. As depicted in Figure 3, RBAB outperforms other solvers on this variant. This confirms that our approach is able to exploit the succinctness of the problem specification in PPDDL.

For exhaustivity, the table on Figure 3 summarizes the results for domains where only a few instances were solved. Interestingly, on these domains the "1 by 1" solver reaches at least as many instances as the others, whereas it was almost always outperformed in the other domains. Apart from this, each solver clearly outperforms both others on one of these hard domains.

Policy revision using F-values

So far, F -values have been used as a tool for computing usual action-values. We now describe how they can be used for some forms of policy revision.

By "revision", we mean here the general problem of computing the new value of an action a, given its "old" value and a modification of its description. More precisely, we write a for an action, F a V for its F -value with respect to some value function V , and a ′ for another action intended to be a modified version of a. Our general goal is to compute F a ′ V given F a V and the PPDDL descriptions of a, a ′ . For instance, assuming V is the optimal infinite-horizon value function for the problem at hand, with π an associated optimal policy, we have that F a V is the optimal (frameless) action-value function for a. Then the revised function F a ′ V gives the expected value of performing a ′ , then following the old policy π from the resulting state on. In equations, for any state s, F a ′ V (s, s) = E(t,r)∼D(e a ′ ,s) [r + V (t ⊲ s))]. Notably, this revision problem arises in model-based reinforcement learning (RL), where action descriptions and rewards are incrementally learnt by the agent, and hence the corresponding action values need to be constantly revised. In particular, Algorithms RTDP-RMAX and RTDP-IE [START_REF] Strehl | Incremental model-based learners with formal learning-time guarantees[END_REF]) perform a single backup step, as above, when rewards and transitions are updated. For instance, in an RL scenario, after it has gathered some further experience an agent may want to revise the relative probabilities of effects x and y in Figure 2, from 0.6/0.4 to 0.8/0.2, resulting in the diagram depicted on Figure 2 (j).

A natural way to solve revision problems is to store the value function V and, when revision is required, to compute Q a ′ V from scratch. We propose instead to store the F -value of a (recall that at any time, Q a V can be retrieved efficiently by

Q a V (s) = F a V (s, s)).
As we show next, using the richness of information of F -values, we can save a substantial amount of computation for revision. The counterpart is increased memory needs, since F -values are bigger than usual actionvalues. How this tradeoff speed/memory is ideally resolved depends on the application at hand.

0.6B V (x) y ′ 3 1.2 (a) 7 5.2 y ′ x ′ 5 y ′ 3.2 B V (probabilistic 0.6x 0.4y) 10 y ′ x ′ 5 2 B V (when x (probabilistic 0.6x 0.4y)) x ↔ x ′ × y ↔ y ′ = 0.4B V (y) (d) (e) (f) (c) (b)
x

′ x x ′ y y ′ y ′ 10 y 5 y ′ 0 y ′ 3.2 + (j) y ′ (g)
x ′

(h) + -4.17 Adding an effect The first, simplest case is when the revision of a consists in adding an effect e, namely, when e a ′ = (and e a e) (we assume that e is consistent with a).

′ p ′ p × F [x ′ ] α × F [x ′ ∧y ′ ] F a ′ V 1-p ′ 1-p × F [y ′ ]
Then it follows directly from Theorem 1 that F a ′ V can be computed from F a V simply as B F a V (e). Observe that no knowledge of V is needed, and that only the "difference" between the revised and the old action has to be treated. Moreover, we wish to emphasize that F -values are necessary for this to be possible, in the sense that usual action-values do not carry enough information.

Example 3 As shown in Example 1, for a = (when x (not x)) and V (x) = 0, V (x) = 1, we have Q a V (•) = 1.
Hence when revising a into a ′ = (and (when x (not x)) (when (not x) x)), usual action values alone cannot distinguish between cases V (x) = 0 as above, and, say, V (x) = 10, whereas revision should yield Q a ′ V (x) = 0 in the former case, but Q a ′ V (x) = 10 in the latter.

Revising rewards/costs It is straightforward to see that revising the reward or cost of an action in a group of states S can be handled by simply adding an effect of the form (when φ (increase (reward) r)), where φ characterizes S and r is the revision amount, and revising as in the previous paragraph. We however wish to note that this particular case can be handled as well using usual action-values only (in both cases, simply add φ × r to the old value function).

Revising probabilities We now consider the case when the probabilities of some effects are revised, namely, when (prob p e . . .) is revised into (prob p ′ e . . .). We restrict ourselves to the case where this clause occurs at the highest level in the description of the action, namely, when a is of the form (and e 1 . . . e k (when φ (prob p e . . .))) (wlog, since and is commutative). We also assume that e is a conjunction of simple effects (literals), as in PSO representations. However, we place no restriction on other e i 's.

For clarity, we start with the case when there is only one probabilistic effect, namely, when we want to revise a = (and e 1 . . . e k (when φ (prob p e))) to a ′ = (and e 1 . . . e k (when φ (prob p ′ e))). The next proposition shows that again, revision can be performed with the sole knowledge of F a V . We omit the proof because it is a special case of Proposition 2 (with f = ∅), but the intuition in this particular case is that revising p to p ′ amounts to add the effect (when φ (prob δ e)) to a, because this means for e not to occur with probability (1p)(1δ) = (1p ′ ) (Proposition 2 shows correctness for negative δ). Proposition 1 Let a be an action of the form (and e 1 . . . e k (when φ (prob p e))), where e is a conjunction of simple effects and p = 1. Let V be a value function, let p ′ = p, and let a ′ be the action (and e 1 . . . e k (when φ (prob p ′ e))). Then, with

δ = 1 -1-p ′ 1-p , F a ′ V can be computed from F a V as ITE(φ, δ × B F a V (e) + (1 -δ)F a V , F a V )
We now turn to the more complex case where a clause of the form (prob p e (1p) f ) is revised. Again, we assume that it occurs at the highest level of a and that e, f are conjunctions of simple effects. Moreover, we assume that e, f are mutually consistent (contain no opposite literals).

So assume p is revised to p ′ . We cannot simply "add" new clauses, say (prob δ e e) and (prob δ f f ), as in the case of a single outcome, since this would incorrectly introduce possible outcomes with e and f both occurring. Once again however, it turns out that the F -value of a can be revised without recomputing it from scratch. Proposition 2 Let a be an action of the form (and e 1 . . . e k (when φ (prob p e (1p) f ))), where e, f are mutually consistent conjunctions of simple effects and p, 1p = 0.

Let V be a value function, let p ′ = p, and let a ′ be the action (and e 1 . . . e k (when φ (prob

p ′ e (1 -p ′ ) f ))). Then F a ′ V can be computed from F a V as ITE(φ, p ′ p × F [e ′ ] + 1 -p ′ 1 -p × F [f ′ ] -α × F [e ′ ∪f ′ ] , F a V ) with α = p ′ (1-p) p + p(1-p ′ ) 1-p , F [e ′ ] = (F a V ) [e ′
] with e ′ the primed version of e, and similarly for

F [f ′ ] , F [e ′ ∪f ′ ] .
Proof. Write W for B V ′ (and e 1 . . . e k ). By Theorem 1 and because e, f are conjunctions of simple effects, the correct value is given by

F a ′ V = ITE(φ, p ′ × W [e ′ ] + (1 -p ′ ) × W [f ′ ] , W ) (4
) Moreover, for the same reasons we have

F a V = ITE(φ, p × W [e ′ ] + (1 -p) × W [f ′ ] , W ) (5)
Clearly the result holds for states not satisfying φ. Now for states satisfying φ we claim that

p ′ p × F [e ′ ] + 1 -p ′ 1 -p × F [f ′ ] -α × F [e ′ ∪f ′ ] (6) is the correct function. Indeed, observing (W [e ′ ] ) [e ′ ] = W [e ′ ] , (W [f ′ ] ) [f ′ ] = W [f ′ ] , and (W [e ′ ] ) [f ′ ] = (W [f ′ ] ) [e ′ ] = W [e ′ ∪f ′ ]
(since e ′ , f ′ are consistent together), we get from Equation 5:

p ′ p × F [e ′ ] = p ′ × W [e ′ ] + p ′ (1 -p) p × W [e ′ ∪f ′ ] 1 -p ′ 1 -p × F [f ′ ] = p(1 -p ′ ) (1 -p) × W [e ′ ∪f ′ ] + (1 -p ′ ) × W [f ′ ] α × F [e ′ ∪f ′ ] = p ′ (1 -p) p + p(1 -p ′ ) (1 -p) W [e ′ ∪f ′ ]
Substituting these values in Equation 6shows equivalence with the "then" argument in Equation 4, as desired.

When there are more than 2 effects or effects are not consistent together, it can be shown that probabilities still can be revised from old F -values alone, albeit with a potential combinatorial explosion (depending on how effects interact). For space reasons however, we do not elaborate on this here. It can also be easily seen that this contruction extends to any deterministic effect e (possibly with nesting), albeit with corrective terms proportional to the reward associated to e (which is 0 for conjunctions of simple effects).

The construction of Proposition 2 is illustrated on the right part of Figure 2 for revising (prob 0.6 x 0.4 y) to (prob 0.8 x 0.2 y).

Using F-values for control As we already insisted, Fvalues represent the possible post-action combinations of variables, taking into account explicit effects, but not enforcing any frame assumption on other variables. We now briefly discuss an application to control.

Assume you want to embark a policy for some MDP, and the evolution of some variables is stochastic, independent of your actions, but might be known one step in advance. To make things concrete, assume you know the probability for rain to fall (rain) on each day, but you also have the opportunity, at execution time, to know whether or not it will rain in the next state (by watching weather forecast, say).

Precisely, assume that you have two deterministic actions, namely taking an umbrella (take-umb), with cost 1 in all states, and doing nothing (noop), with cost 0 in all states. Assume moreover that in any state, rain falls with probability 0.2, and does not with probability 0.8, independently of whether it was falling before. Finally, in any state, if it rains and you do not carry your umbrella (umb) you incur a penalty of -10, while in other cases you incur no penalty. Now let s = rain umb. Under the usual semantics, action take-umb has a (1 step-to-go) value of -1 + 0.2 × 0 + 0.8 × 0 = -1 in s, and noop has a value of 0 + 0.2 × -10 + 0.8 × 0 = -5. Hence the best action in s is to take the umbrella. Now assume that at execution time, in state s you learn that it will not rain in the next state (say, tomorrow), that is, you learn that rain ′ is false. Then action-values will not tell you what the best action is with this new information, since they do not make explicit the role of rain in expected values. Constrastingly, assume you have F -values available, namely, values F a V (a ∈ {take-umb,noop}), where V is an optimal value function and F a V is computed wrt V without taking into account the probability of rain. Then you can simply evaluate them in s with ℓ ′ = rain ⊲ s. This will tell you exactly the expected value of actions knowing that it will not rain tomorrow (and using probabilities 0.2/0.8 from tomorrow on). Precisely, in this case, the F -value for take-umb will still be -1, but that for noop will become 0, making noop now the best action.

We insist that F -values embed action-values. In this case, if you miss the weather forecast at execution time, you can simply add the effect (prob 0.2 rain 0.8 rain) to both actions, update their F -values, and retrieve the action-values.

Another case where such control is particularly relevant is in multi-agent systems. Here you may know at execution time whether some other agents will help you by setting the values of some variables which you do not control. We leave a detailed investigation of such applications to multi-agent systems and cooperation for future work, but we think that they are very promising.

As a final note, we remark that this problem can alternatively be solved by introducing specific fluents with an epistemic flavor, say, two fluents coding whether you know the forecast, and what it is in case you know. Then action descriptions can be conditional on their values, and here usual action-values are enough. Nevertheless, it is clear that this makes the state space blow up, in particular if there are many such "short-term predictable" fluents.

Conclusion

We presented RBAB, an algorithm for computing optimal policies for MDPs described in PPDDL. To that end, we introduced the notion of a frameless action-value, and showed that using them for action backup in a value iteration scheme yields an algorithm which is complementary to previous approaches, and interesting in particular for correlated action effects and compact PPDDL descriptions. We also discussed some promising applications of this notion to some policy revision and control problems.

Our future work will focus on goal-oriented problems, especially using reachability analysis for computing partial policies over the set of reachable states, as in the LAO* family of algorithms. Also promising is the adaptation of our techniques to PPDDL at the relational level, for which there has been a lot of excitement recently [START_REF] Boutilier | Practical Solution Techniques for First-Order MDPs[END_REF][START_REF] Lang | Planning with Noisy Probabilistic Relational Rules[END_REF]. To that end we plan to use First-Order ADDs [START_REF] Wang | First order decision diagrams for relational MDPs[END_REF]. On the short term, we also plan to enrich our approach with the use of Affine ADDs [START_REF] Sanner | Affine algebraic decision diagrams (AADDs) and their application to structured probabilistic inference[END_REF], which we believe will improve the behaviour of our algorithm just as it does for SPUDD and other solvers.

  ′ , y ′ ) over primed value function V

Figure 2 :Figure 3 :

 23 Figure 2: (Left) Application of the backup rules on a value function V for the effect (when x (prob 0.6 x 0.4 y)). Diagrams (a) and (b) present the backup of simple effects (weighted by their probabilities), (c) is the diagram for the probabilistic effect and (d) for the conditional one. Diagram (e) is the intermediate diagram of Equation 3 and (f) is the final action-value function. (Right) Revision of the F -value for the prob clause with 0.8/0.2 instead of 0.6/0.4.

•

  Simple Effects specify the truth value update for a predicate p; (p) (resp. (not p)) indicates that p will be true (resp. false) after execution of the effect.• Update Effects specify how the reward function is affected. They are either (increase (reward) r) or (decrease (reward) r) with obvious meaning. • Conditional Effects of the form (when φ e), where φ is a formula over the domain predicates and e an effect, mean that e occurs only in states satifying φ. • Probabilistic Effects of the form (probabilistic p 1 e 1 . . . p k e k ) mean that effect e i occurs with probability p i . We write (prob p 1 e 1 . . . p k e k ) for short. In case i p i < 1, the empty effect implicitely occurs with probability 1i p i . • Conjunctive Effects of the form (and e 1 . . . e k ), make all e i 's occur synchronously. Effects e i must be consistent together, i.e., not contain (x) and (not x) (whatever the combination of outcomes if e i 's are probabilistic).

We omit functions, which represent numeric state-variables.