David Lopez-Paz

Maxime Oquab
email: maxime.oquab@inria.fr

REVISITING CLASSIFIER TWO-SAMPLE TESTS

The goal of two-sample tests is to assess whether two samples, S P ∼ P n and S Q ∼ Q m , are drawn from the same distribution. Perhaps intriguingly, one relatively unexplored method to build two-sample tests is the use of binary classifiers. In particular, construct a dataset by pairing the n examples in S P with a positive label, and by pairing the m examples in S Q with a negative label. If the null hypothesis "P = Q" is true, then the classification accuracy of a binary classifier on a held-out subset of this dataset should remain near chance-level. As we will show, such Classifier Two-Sample Tests (C2ST) learn a suitable representation of the data on the fly, return test statistics in interpretable units, have a simple null distribution, and their predictive uncertainty allow to interpret where P and Q differ. The goal of this paper is to establish the properties, performance, and uses of C2ST. First, we analyze their main theoretical properties. Second, we compare their performance against a variety of state-of-the-art alternatives. Third, we propose their use to evaluate the sample quality of generative models with intractable likelihoods, such as Generative Adversarial Networks (GANs). Fourth, we showcase the novel application of GANs together with C2ST for causal discovery.

INTRODUCTION

One of the most fundamental problems in statistics is to assess whether two samples, S P ∼ P n and S Q ∼ Q m , are drawn from the same probability distribution. To this end, two-sample tests (?) summarize the differences between the two samples into a real-valued test statistic, and then use the value of such statistic to accept1 or reject the null hypothesis "P = Q". The development of powerful two-sample tests is instrumental in a myriad of applications, including the evaluation and comparison of generative models. Over the last century, statisticians have nurtured a wide variety of two-sample tests. However, most of these tests are only applicable to one-dimensional examples, require the prescription of a fixed representation of the data, return test statistics in units that are difficult to interpret, or do not explain how the two samples under comparison differ.

Intriguingly, there exists a relatively unexplored strategy to build two-sample tests that overcome the aforementioned issues: training a binary classifier to distinguish between the examples in S P and the examples in S Q . Intuitively, if P = Q, the test accuracy of such binary classifier should remain near chance-level. Otherwise, if P = Q and the binary classifier is able to unveil some of the distributional differences between S P and S Q , its test accuracy should depart from chance-level. As we will show, such Classifier Two-Sample Tests (C2ST) learn a suitable representation of the data on the fly, return test statistics in interpretable units, have simple asymptotic distributions, and their learned features and predictive uncertainty provide interpretation on how P and Q differ. In such a way, this work brings together the communities of statistical testing and representation learning.

The goal of this paper is to establish the theoretical properties and evaluate the practical uses of C2ST. To this end, our contributions are:

• We review the basics of two-sample tests in Section ??, as well as their common applications to measure statistical dependence and evaluate generative models.

• We analyze the attractive properties of C2ST (Section ??) including an analysis of their exact asymptotic distributions, testing power, and interpretability.

• We evaluate C2ST on a wide variety of synthetic and real data (Section ??), and compare their performance against multiple state-of-the-art alternatives. Furthermore, we provide examples to illustrate how C2ST can interpret the differences between pairs of samples.

• In Section ??, we propose the use of classifier two-sample tests to evaluate the sample quality of generative models with intractable likelihoods, such as Generative Adversarial Networks (?), also known as GANs.

• As a novel application of the synergy between C2ST and GANs, Section ?? proposes the use of these methods for causal discovery.

TWO-SAMPLE TESTING

The goal of two-sample tests is to assess whether two samples, denoted by S P ∼ P n and S Q ∼ Q m , are drawn from the same distribution (?). More specifically, two-sample tests either accept or reject the null hypothesis, often denoted by H 0 , which stands for "P = Q". When rejecting H 0 , we say that the two-sample test favors the alternative hypothesis, often denoted by H 1 , which stands for "P = Q". To accept or reject H 0 , two-sample tests summarize the differences between the two samples (sets of identically and independently distributed examples):

S P := {x 1 , . . . , x n } ∼ P n (X) and S Q := {y 1 , . . . , y m } ∼ Q m (Y ) (1)
into a statistic t ∈ R. Without loss of generality, we assume that the two-sample test returns a small statistic when the null hypothesis "P = Q" is true, and a large statistic otherwise. Then, for a sufficiently small statistic, the two-sample test will accept H 0 . Conversely, for a sufficiently large statistic, the two-sample test will reject H 0 in favour of H 1 .

More formally, the statistician performs a two-sample test in four steps. First, decide a significance level α ∈ [0, 1], which is an input to the two-sample test. Second, compute the two-sample test statistic t. Third, compute the p-value p = P (T ≥ t|H 0 ), the probability of the two-sample test returning a statistic as large as t when H 0 is true. Fourth, reject H 0 if p < α, and accept it otherwise.

Inevitably, two-sample tests can fail in two different ways. First, to make a type-I error is to reject the null hypothesis when it is true (a "false positive"). By the definition of p-value, the probability of making a type-I error is upper-bounded by the significance level α. Second, to make a type-II error is to accept the null hypothesis when it is false (a "false negative"). We denote the probability of making a type-II error by β, and refer to the quantity π = 1β as the power of a test. Usually, the statistician uses domain-specific knowledge to evaluate the consequences of a type-I error, and thus prescribe an appropriate significance level α. Within the prescribed significance level α, the statistician prefers the two-sample test with maximum power π.

Among others, two-sample tests serve two other uses. First, two-sample tests can measure statistical dependence (?). In particular, testing the independence null hypothesis "the random variables X and Y are independent" is testing the two-sample null hypothesis "P (X, Y ) = P (X)P (Y )". In practice, the two-sample test would compare the sample

S = {(x i , y i )} n i=1 ∼ P (X, Y ) n to a sample S σ = {(x i , y σ(i) )} n i=1 ∼ (P (X)P (Y )) n
, where σ is a random permutation of the set of indices {1, . . . , n}. This approach is consistent when considering all possible random permutations. However, since independence testing is a subset of two-sample testing, specialized independence tests may exhibit higher power for this task (?).

Second, two-sample tests can evaluate the sample quality of generative models with intractable likelihoods, but tractable sampling procedures. Intuitively, a generative model produces good samples Ŝ = {x i } n i=1 if these are indistinguishable from the real data S = {x i } n i=1 that they model. Thus, the two-sample test statistic between Ŝ and S measures the fidelity of the samples Ŝ produced by the generative model. The use of two-sample tests to evaluate the sample quality of generative models include the pioneering work of ?, the use of Maximum Mean Discrepancy (MMD) criterion (?????), and the connections to density-ratio estimation (????).

Over the last century, statisticians have nurtured a wide variety of two-sample tests. Classical two-sample tests include the t-test (?), which tests for the difference in means of two samples; the Wilcoxon-Mann-Whitney test (??), which tests for the difference in rank means of two samples; and the Kolmogorov-Smirnov tests (??) and their variants (?), which test for the difference in the empirical cumulative distributions of two samples. However, these classical tests are only efficient when applied to one-dimensional data. Recently, the use of kernel methods (?) enabled the development of twosample tests applicable to multidimensional data. Examples of these tests include the MMD test (?), which looks for differences in the empirical kernel mean embeddings of two samples, and the Mean Embedding test or ME (??), which looks for differences in the empirical kernel mean embeddings of two samples at optimized locations. However, kernel two-sample tests require the prescription of a manually-engineered representation of the data under study, and return values in units that are difficult to interpret. Finally, only the ME test provides a mechanism to interpret how P and Q differ.

Next, we discuss a simple but relatively unexplored strategy to build two-sample tests that overcome these issues: the use of binary classifiers.

CLASSIFIER TWO-SAMPLE TESTS (C2ST)

Without loss of generality, we assume access to the two samples S P and S Q defined in (??), where x i , y j ∈ X , for all i = 1, . . . , n and j = 1, . . . , m, and m = n. To test whether the null hypothesis H 0 : P = Q is true, we proceed in five steps. First, construct the dataset

D = {(x i , 0)} n i=1 ∪ {(y i , 1)} n i=1 =: {(z i , l i )} 2n i=1 .
Second, shuffle D at random, and split it into the disjoint training and testing subsets D tr and D te , where D = D tr ∪ D te and n te := |D te |. Third, train a binary classifier f : X → [0, 1] on D tr ; in the following, we assume that f (z i ) is an estimate of the conditional probability distribution p(l i = 1|z i ). Fourth, return the classification accuracy on D te :

t = 1 n te (zi,li)∈Dte I I f (z i ) > 1 2 = l i (2)
as our C2ST statistic, where I is the indicator function. The intuition here is that if P = Q, the test accuracy (??) should remain near chance-level. In opposition, if P = Q and the binary classifier unveils distributional differences between the two samples, the test classification accuracy (??) should be greater than chance-level. Fifth, to accept or reject the null hypothesis, compute a p-value using the null distribution of the C2ST, as discussed next.

NULL AND ALTERNATIVE DISTRIBUTIONS

Each term

I [I(f (z i ) > 1/2) = l i ]
appearing in (??) is an independent Bernoulli(p i ) random variable, where p i is the probability of classifying correctly the example z i in D te .

First, under the null hypothesis H 0 : P = Q, the samples S P ∼ P n and S Q ∼ Q m follow the same distribution, leading to an impossible binary classification problem. In that case, n te t follows a Binomial(n te , p = 1 2 ) distribution. Therefore, for large n te , we can use the central limit theorem to approximate the null distribution of (??) by N ( 1 2 , 1 4nte ). Second, under the alternative hypothesis H 1 : P = Q, the statistic n te t follows a Poisson Binomial distribution, since the constituent Bernoulli random variables may not be identically distributed. In the following, we will approximate such Poisson Binomial distribution by the Binomial(n, p) distribution, where p = 1 n n i=1 p i (?). Therefore, we can use the central limit theorem to approximate the alternative distribution of (??) by N (p, p(1-p) nte ).

TESTING POWER

To analyze the power (probability of correctly rejecting false null hypothesis) of C2ST, we assume that the our classifier has an expected (unknown) accuracy of H 0 : t = 1 2 under the null hypothesis "P = Q", and an expected accuracy of H 1 : t = 1 2 + under the alternative hypothesis "P = Q", where ∈ (0, 1 2 ) is the effect size distinguishing P from Q. Let Φ be the Normal cdf, n te the number of samples available for testing, and α the significance level. Then, Theorem 1. Given the conditions described in the previous paragraph, the approximate power of the statistic

(??) is Φ √ nte-Φ -1 (1-α)/2 √ 1 4 -2 .
See Appendix ?? for a proof. The power bound in Theorem ?? has an optimal order of magnitude for multi-dimensional problems (???). These are problems with fixed d and n → ∞, where the power bounds do not depend on d.

Remark 1. We leave for future work the study of quadratic-time C2ST with optimal power in highdimensional problems (?). These are problems where the ratio n/d → c ∈ [0, 1], and the power bounds depend on d. One possible line of research in this direction is to investigate the power and asymptotic distributions of quadratic-time C2ST statistics

1 nte(nte-1) i =j I[I(f (z i , z j ) > 1 2 ) = l i ],
where the classifier f (z, z ) predicts if the examples (z, z ) come from the same sample.

Theorem ?? also illustrates that maximizing the power of a C2ST is a trade-off between two competing objectives: choosing a classifier that maximizes the test accuracy and maximizing the size of the test set n te . This relates to the well known bias-variance trade-off in machine learning. Indeed, simple classifiers will miss more nonlinear patterns in the data (leading to smaller test accuracy), but call for less training data (leading to larger test set sizes). On the other hand, flexible classifiers will miss less nonlinear patterns in the data (leading to higher test accuracy), but call for more training data (leading to smaller test sizes). Formally, the relationship between the test accuracy, sample size, and the flexibility of a classifier depends on capacity measures such as the VC-Dimension (?). Note that there is no restriction to perform model selection (such as cross-validation) on D tr .

Remark 2. We have focused on test statistics (??) built on top of the zero-one loss 0-1 (y, y ) = I[y = y ] ∈ {0, 1}. These statistics give rise to Bernoulli random variables, which can exhibit high variance. However, our arguments are readily extended to real-valued binary classification losses. Then, the variance of such real-valued losses would describe the norm of the decision function of the classifier two-sample test, appear in the power expression from Theorem ??, and serve as a hyper-parameter to maximize power as in (?, Section 3).2

INTERPRETABILITY

There are three ways to interpret the result of a C2ST. First, recall that the classifier predictions f (z i ) are estimates of the conditional probabilities p(l i = 1|z i ) for each of the samples z i in the test set. Inspecting these probabilities together with the true labels l i determines which examples were correctly or wrongly labeled by the classifier, with the least or the most confidence. Therefore, the values f (z i ) explain where the two distributions differ. Second, C2ST inherit the interpretability of their classifiers to explain which features are most important to distinguish distributions, in the same way as the ME test (?). Examples of interpretable features include the filters of the first layer of a neural network, the feature importance of random forests, the weights of a generalized linear model, and so on. Third, C2ST return statistics t in interpretable units: these relate to the percentage of samples correctly distinguishable between the two distributions. These interpretable numbers can complement the use of p-values.

PRIOR USES

The reduction of two-sample testing to binary classification was introduced in (?), studied within the context of information theory in (??), discussed in (??), and analyzed (for the case of linear discriminant analysis) in (?). The use of binary classifiers for two-sample testing is increasingly common in neuroscience: see (??) and the references therein. Implicitly, binary classifiers also perform two-sample tests in algorithms that discriminate data from noise, such as unsupervised-assupervised learning (?), noise contrastive estimation (?), negative sampling (?), and GANs (?). a significant improvement in performance when increasing the flexibility of these classifiers (e.g., increasing the number of hidden neurons or decreasing the number of nearest neighbors). When analyzing one-dimensional data, we compare the performance of C2ST-NN and C2ST-KNN against the Wilcoxon-Mann-Whitney test (??), the Kolmogorov-Smirnov test (??), and the Kuiper test (?). In all cases, we also compare the performance of C2ST-NN and C2ST-KNN against the linear-time estimate of the Maximum Mean Discrepancy (MMD) criterion (?), the ME test (?), and the SCF test (?). We use a significance level α = 0.05 across all experiments and tests, unless stated otherwise. We use Gaussian approximations to compute the null distributions of C2ST-NN and C2ST-KNN. We use the implementations of the MMD, ME, and SCF tests gracefully provided by ?, the scikitlearn implementation of the Kolmogorov-Smirnov and Wilcoxon tests, and the implementation from https://github.com/aarchiba/kuiper of the Kuiper test. The implementation of our experiments is available at https://github.com/lopezpaz/classifier_tests.

EXPERIMENTS ON TWO-SAMPLE TESTING

EXPERIMENTS ON TWO-SAMPLE TESTING

Control of type-I errors We start by evaluating the correctness of all the considered two-sample tests by examining if the prescribed significance level α = 0.05 upper-bounds their type-I error.

To do so, we draw x 1 , . . . , x n , y 1 , . . . , y n ∼ N (0, 1), and run each two-sample test on the two samples {x i } n i=1 and {y i } n i=1 . In this setup, a type-I error would be to reject the true null hypothesis. Figure ??(a) shows that the type-I error of all tests is upper-bounded by the prescribed significance level, for all n ∈ {25, 50, 100, 500, 1000, 5000, 10000} and 100 random repetitions. Thus, all tests control their type-I error as expected, up to random variations due to finite experiments.

Gaussian versus Student

We consider distinguishing between samples drawn from a Normal distribution and samples drawn from a Student's t-distribution with ν degrees of freedom. We shift and scale both samples to exhibit zero-mean and unit-variance. Since the Student's t distribution approaches the Normal distribution as ν increases, a two-sample test must focus on the peaks of the distributions to distinguish one from another. Figure ??(b,c) shows the percentage of type-II errors made by all tests as we vary separately n and ν, over 100 trials (random samples). We set n = 2000 when ν varies, and let ν = 3 when n varies. The Wilcoxon-Mann-Whitney exhibits the worst performance, as expected (since the ranks mean of the Gaussian and Student's t distributions coincide) in this experiment. The best performing method is the the one-dimensional Kuiper test, followed closely by the multi-dimensional tests C2ST-NN and ME.

Independence testing on sinusoids For completeness, we showcase the use two-sample tests to measure statistical dependence. This can be done, as described in Section ??, by performing a two-sample test between the observed data {(x i , y i )} n i=1 and {(x i , y σ(i) )} n i=1 , where σ is a random

Distinguishing between NIPS articles

We consider the problem of distinguishing between some of the categories of the 5903 articles published in the Neural Information Processing Systems (NIPS) conference from 1988 to 2015, as discussed in ?. We consider articles on Bayesian inference (Bayes), neuroscience (Neuro), deep learning (Deep), and statistical learning theory (Learn). Table ?? shows the type-I errors (Bayes-Bayes row) and powers (rest of rows) for the tests reported in (?), together with C2ST-NN, at a significance level α = 0.01, when averaged over 500 trials. In these experiments, C2ST-NN achieves maximum power, while upper-bounding its type-I error by α.

Distinguishing between facial expressions Finally, we apply C2ST-NN to the problem of distinguishing between positive (happy, neutral, surprised) and negative (afraid, angry, disgusted) facial expressions from the Karolinska Directed Emotional Faces dataset, as discussed in (?). See the fourth plot of Figure ??, first two-rows, for one example of each of these six emotions. Table ?? shows the type-I errors (± vs ± row) and the powers (+ vsrow) for the tests reported in (?), together with C2ST-NN, at α = 0.01, averaged over 500 trials. C2ST-NN achieves a near-optimal power, only marginally behind the perfect results of SCF-full and MMD-quad.

EXPERIMENTS ON GENERATIVE ADVERSARIAL NETWORK EVALUATION

Since effective generative models will produce examples barely distinguishable from real data, twosample tests arise as a natural alternative to evaluate generative models. Particularly, our interest is to evaluate the sample quality of generative models with intractable likelihoods, such as GANs (?). GANs implement the adversarial game

min g max d E x∼P (X) [log(d(x))] + E z∼P (Z) [log(1 -d(g(z)))] , (3)
where d(x) depicts the probability of the example x following the data distribution P (X) versus being synthesized by the generator. This is according to a trainable discriminator function d. In the adversarial game, the generator g plays to fool the discriminator d by transforming noise vectors z ∼ P (Z) into real-looking examples g(z). On the opposite side, the discriminator plays to distinguish between real examples x and synthesized examples g(z). To approximate the solution to (??), alternate the optimization of the two losses (?) given by Under the formalization (??), the adversarial game reduces to the sequential minimization of L d (d) and L g (g), and reveals the true goal of the discriminator: to be the C2ST that best distinguishes data examples x ∼ P and synthesized examples x ∼ P , where P is the probability distribution induced by sampling z ∼ P (Z) and computing x = g(z). The formalization (??) unveils the existence of an arbitrary binary classification loss function (See Remark ??), which in turn decides the divergence minimized between the real and fake data distributions (?).

L d (d) = E x [ (d(x), 1)] + E z [ (d(g(z)), 0)] , L g (g) = E x [ (d(x), 0)] + E z [ (d(g(z)), 1)] . (
Unfortunately, the evaluation of the log-likelihood of a GANs is intractable. Therefore, we will employ a two-sample test to evaluate the quality of the fake examples x = g(z). In simple terms, evaluating a GAN in this manner amounts to withhold some real data from the training process, and use it later in a two-sample test against the same amount of synthesized data. When the two-sample test is a binary classifier (as discussed in Section ??), this procedure is simply training a fresh discriminator on a fresh set of data. Since we train and test this fresh discriminator on held-out examples, it may differ from the discriminator trained along the GAN. In particular, the discriminator trained along with the GAN may have over-fitted to particular artifacts produced by the generator, thus becoming a poor C2ST.

We evaluate the use of two-sample tests for model selection in GANs. To this end, we train a number of DCGANs (?) on the bedroom class of LSUN (?) and the Labeled Faces in the Wild (LFW) dataset (?). We reused the Torch7 code of ? to train a set of DCGANs for {1, 10, 50, 100, 200} epochs, where the generator and discriminator networks are convolutional neural networks (?) with {1, 2, 4, 8} × gf and {1, 2, 4, 8} × df filters per layer, respectively. We evaluate each DCGAN on 10, 000 held-out examples using the fastest multi-dimensional two-sample tests: MMD, C2ST-NN, and C2ST-KNN.

Our first experiments revealed an interesting result. When performing two-sample tests directly on pixels, all tests obtain near-perfect test accuracy when distinguishing between real and synthesized (fake) examples. Such near-perfect accuracy happens consistently across DCGANs, regardless of the visual quality of their examples. This is because, albeit visually appealing, the fake examples contain checkerboard-like artifacts that are sufficient for the tests to consistently differentiate between real and fake examples. ? discovered this phenomenon concurrently with us.

On a second series of experiments, we featurize all images (both real and fake) using a deep convolutional ResNet (?) pre-trained on ImageNet, a large dataset of natural images (?). In particular, we use the resnet-34 model from ?. Reusing a model pre-trained on natural images ensures that the test will distinguish between real and fake examples based only on natural image statistics, such as Gabor filters, edge detectors, and so on. Such a strategy is similar to perceptual losses (?) and inception scores (?). In short, in order to evaluate how natural the images synthesized by a DCGAN look, one must employ a "natural discriminator". Table ?? shows three GANs producing poor samples and three GANs producing good samples for the LSUN and LFW datasets, according to the MMD, C2ST-KNN, C2ST-NN tests on top of ResNet features. See Appendix ?? for the full list of results. Although it is challenging to provide with an objective evaluation of our results, we believe that the rankings provided by two-sample tests could serve for efficient early stopping and model selection.

Remark 3 (How good is my GAN? Is it overfitting?). Evaluating generative models is a delicate issue (?), but two-sample tests may offer some guidance. In particular, good (non-overfitting) generative models should produce similar two-sample test statistics when comparing their generated samples to both the train-set and the test-set samples.3 As a general recipe, prefer the smallest (in number of parameters) generative model that achieves the same and small two-sample test statistic when comparing their generated samples to both the train-set and test-set samples.

We have seen that GANs of different quality may lead to the same (perfect) C2ST statistic. To allow a finer comparison between generative models, we recommend implementing C2ST using a margin classifier with finite norm, or using as statistic the whole area under the C2ST training curve (on train-set or test-set samples).

EXPERIMENTS ON INTERPRETABILITY

We illustrate the interpretability power of C2ST. First, the predictive uncertainty of C2ST sheds light on where the two samples under consideration agree or differ. In the context of GANs, this interpretability is useful to locate captured or dropped modes. In the first plot of Figure ??, a C2ST-NN separates two bivariate Gaussian distributions with different means. When performing this separation, the C2ST-NN provides an explicit decision boundary that illustrates where the two distributions separate from each other. In the second plot of Figure ??, a C2ST-NN separates a Gaussian distribution from a Student's t distribution with ν = 3, after scaling both to zero-mean and unit-variance. The plot reveals that the peaks of the distributions are their most differentiating feature. Finally, the third plot of Second, the features learnt by the classifier of a C2ST are also a mechanism to understand the differences between the two samples under study. The third plot of Figure ?? shows six examples from the Karolinska Directed Emotional Faces dataset, analyzed in Section ??. In that same figure, we arrange the weights of the first linear layer of C2ST-NN into the feature most activated at positive examples (bottom left, positive facial expressions), the feature most activated at negative examples (bottom middle, negative facial expressions), and the "discriminative feature", obtained by substracting these two features (bottom right). The discriminative feature of C2ST-NN agrees with the one found by (?): positive and negative facial expressions are best distinguished at the eyebrows, smile lines, and lips. A similar analysis ? on the C2ST-NN features in the NIPS article classification problem (Section ??) reveals that the features most activated for the "statistical learning theory" category are those associated to the words inequ, tight, power, sign, hypothesi, norm, hilbert. The features most activated for the "Bayesian inference" category are those associated to the words infer, markov, graphic, conjug, carlo, automat, laplac.

Method

EXPERIMENTS ON CONDITIONAL GANS FOR CAUSAL DISCOVERY

In causal discovery, we study the causal structure underlying a set of d random variables X 1 , . . . , X d .

In particular, we assume that the random variables X 1 , . . . , X d share a causal structure described by a collection of Structural Equations, or SEs (?). More specifically, we assume that the random variable X i takes values as described by the SE X i = g i (Pa(X i , G), N i ), for all i = 1, . . . , d. In the previous, G is a Directed Acyclic Graph (DAG) with vertices associated to each of the random variables X 1 , . . . , X d . Also in the same equation, Pa(X i , G) denotes the set of random variables which are parents of X i in the graph G, and N i is an independent noise random variable that follows the probability distribution P (N i ). Then, we say that X i → X j if X i ∈ Pa(X j ), since a change in X i will cause a change in X j , as described by the i-th SE. The goal of causal discovery is to infer the causal graph G given a sample from P (X 1 , . . . , X d ). For the sake of simplicity, we focus on the discovery of causal relations between two random variables, denoted by X and Y . That is, given the sample D = {(x i , y i )} n i=1 ∼ P n (X, Y ), our goal is to conclude whether "X causes Y ", or "Y causes X". We call this problem cause-effect discovery (?). In the case where X → Y , we can write the cause-effect relationship as:

x ∼ P (X), n ∼ P (N ), y ← g(x, n).

(5)

The current state-of-the-art in the cause-effect discovery is the family of Additive Noise Models, or ANM (?). These methods assume that the SE (??) allow the expression y ← g(x) + n, and exploit the independence assumption between the cause random variable X and the noise random variable N to analyze the distribution of nonlinear regression residuals, in both causal directions. Unfortunately, assuming independent additive noise is often too simplistic (for instance, the noise could be heteroskedastic or multiplicative). Because of this reason, we propose to use Conditional Generative Adversarial Networks, or CGANs (?) to address the problem of cause-effect discovery.

Our motivation is the shocking resemblance between the generator of a CGAN and the SE (??): the random variable X is the conditioning variable input to the generator, the random variable N is the noise variable input to the generator, and the random variable Y is the variable synthesized by the generator. Furthermore, CGANs respect the independence between the cause X and the noise N by construction, since n ∼ P (N ) is independent from all other variables. This way, CGANs bypass the additive noise assumption naturally, and allow arbitrary interactions g(X, N ) between the cause variable X and the noise variable N .

To implement our cause-effect inference algorithm in practice, recall that training a CGAN from X to Y minimizes the two following objectives in alternation: [ (d(x, y), 1)] + E x,z [ (d(x, g(x, z)), 0)] , L g (g) = E x,y [ (d(x, y), 0)] + E x,z [ (d(x, g(x, z)), 1)] . Our recipe for cause-effect is to learn two CGANs: one with a generator g y from X to Y to synthesize the dataset D X→Y = {(x i , g y (x i , z i ))} n i=1 , and one with a generator g x from Y to X to synthesize the dataset D X←Y = {(g x (y i , z i ), y i )} n i=1 . Then, we prefer the causal direction X → Y if the two-sample test statistic between the real sample D and D X→Y is smaller than the one between D and D Y →X . Thus, our method is Occam's razor at play: declare the simplest direction (in terms of conditional generative modeling) as the true causal direction.

L d (d) = E x,y
Table ?? summarizes the performance of this procedure when applied to the 99 Tübingen cause-effect pairs dataset, version August 2016 (?). RCC is the Randomized Causation Coefficient of (?). The Ensemble-CGAN-C2ST trains 100 CGANs, and decides the causal direction by comparing the top generator obtained in each causal direction, as told by C2ST-KNN. The need to ensemble is a remainder of the unstable behaviour of generative adversarial training, but also highlights the promise of such models for causal discovery.

CONCLUSION

Our take-home message is that modern binary classifiers can be easily turned into powerful two-sample tests. We have shown that these classifier two-sample tests set a new state-of-the-art in performance, and enjoy unique attractive properties: they are easy to implement, learn a representation of the data on the fly, have simple asymptotic distributions, and allow different ways to interpret how the two samples under study differ. Looking into the future, the use of binary classifiers as two-sample tests provides a flexible and scalable approach for the evaluation and comparison of generative models (such as GANs), and opens the door to novel applications of these methods, such as causal discovery.

Figure 1 :

 1 Figure 1: Results (type-I and type-II errors) of our synthetic two-sample test experiments.

Figure 2 :

 2 Figure 2: Interpretability of C2ST. The color map corresponds to the value of p(l = 1|z).

  Figure ?? displays, for the LFW and LSUN datasets, five examples classified as real with high uncertainty (first row, better looking examples), and five examples classified as fake with high certainty (second row, worse looking examples).

Table 2 :

 2 Type-I errors (first row) and powers (second row) in distinguishing facial expressions. permutation. Since the distributions P (X)P (Y ) and P (X, Y ) are bivariate, only the C2ST-NN, C2ST-KNN, MMD, and ME tests compete in this task. We draw (x i , y i ) according to the generative model x i ∼ N (0, 1), i ∼ N (0, γ 2 ), and y i ∼ cos(δx i ) + i . Here, x i are iid examples from the random variable X, and y i are iid examples from the random variable Y . Thus, the statistical dependence between X and Y weakens as we increase the frequency δ of the sinusoid, or increase the variance γ 2 of the additive noise.

	Problem Bayes-Bayes 215 n te ME-full ME-grid SCF-full SCF-grid MMD-quad MMD-lin C2ST-NN .012 .018 .012 .004 .022 .008 .002
	Bayes-Deep Bayes-Learn 138 216 Bayes-Neuro 394 Learn-Deep 149 Learn-Neuro 146	.954 .990 1.00 .956 .960	.034 .774 .300 .052 .572	.688 .836 .828 .656 .590	.180 .534 .500 .138 .360	.906 1.00 .952 .876 1.00	.262 .238 .972 .500 .538	1.00 1.00 1.00 1.00 1.00
	Table 1: Type-I errors (first row) and powers (rest of rows) in distinguishing NIPS papers categories.
	Problem n te ME-full ME-grid SCF-full SCF-grid MMD-quad MMD-lin C2ST-NN .010 .012 .014 .002 .018 .008 .002 ± vs. ± 201 + vs. -201 .998 .656 1.00 .750 1.00 .578 .997

Figure ??(d,e,f)

shows the percentage of type-II errors made by C2ST-NN, C2ST-KNN, MMD, and ME as we vary separately n, δ, and γ over 100 trials. We let n = 2000, δ = 1, γ = 0.25 when fixed. Figure

??(d,e,f

) reveals that among all tests, C2ST-NN is the most efficient in terms of sample size, C2ST-KNN is the most robust with respect to high-frequency variations, and that C2ST-NN and ME are the most robust with respect to additive noise.

Table 3 :

 3

Results on GAN evaluation. Lower test statistics are best. Full results in Appendix ??.

Table 4 :

 4 Results on cause-effect discovery on the Tübingen pairs experiment.

	Accuracy	ANM-HSIC IGCI RCC CGAN-C2ST Ensemble C2ST type 67% 71% 76% 73% 82% KNN 70% 73% NN 58% 65% MMD

For clarity, we abuse statistical language and write "accept" to mean "fail to reject".

For a related discussion on this issue, we recommend the insightful comment by Arthur Gretton and Wittawat Jitkrittum, available at https://openreview.net/forum?id=SJkXfE5xx.

As discussed with Arthur Gretton, if the generative model memorizes the train-set samples, a sufficiently large set of generated samples would reveal such memorization to the two-sample test. This is because some unique samples would appear multiple times in the set of generated samples, but not in the test-set of samples.

B PROOF OF THEOREM ??

Our statistic is a random variable T ∼ N 1 2 , 1 4nte under the null hypothesis, and T ∼ N 1 2 + , n -1 te 1 4 -2 under the alternative hypothesis. Furthermore, at a significance level α, the threshold of our statistic is

; under this threshold we would accept the null hypothesis. Then, the probability of making a type-II error is

Therefore, the power of the test is

which concludes the proof.

C ACKNOWLEDGEMENTS

We are thankful to L. Bottou, B. Graham, D. Kiela, M. Rojas-Carulla, I. Tolstikhin, and M. Tygert for their help in improving the quality of this manuscript. This work was partly supported by ERC grant LEAP (no. 336845) and CIFAR Learning in Machines & Brains program.


Aurélie Pirayre

Yuling Zheng

Laurent Duval

Jean-Christophe Pesquet

HOGMep: VARIATIONAL BAYES AND HIGHER-ORDER GRAPHICAL MODELS APPLIED TO JOINT IMAGE RECOVERY AND SEGMENTATION

Keywords: variational Bayes, higher-order graphical models, multicomponent images, segmentation, restoration

Variational Bayesian approaches have been successfully applied to image segmentation. They usually rely on a Potts model for the hidden label variables and a Gaussian assumption on pixel intensities within a given class. Such models may however be limited, especially in the case of multicomponent images. We overcome this limitation with HOGMep, a Bayesian formulation based on a higher-order graphical model (HOGM) on labels and a Multivariate Exponential Power (MEP) prior for intensities in a class. Then, we develop an efficient statistical estimation method to solve the associated problem. Its flexibility accommodates to a broad range of applications, demonstrated on multicomponent image segmentation and restoration.

INTRODUCTION

Variational Bayesian Approximation (VBA) has been widely used in various applications such as in graphical model learning [START_REF] Jordan | An introduction to variational methods for graphical models[END_REF], image processing [START_REF] Likas | A variational approach for Bayesian blind image deconvolution[END_REF][START_REF] Chantas | Variational Bayesian image restoration based on a product of tdistributions image prior[END_REF][START_REF] Chen | Variational Bayesian methods for multimedia problems[END_REF][START_REF] Zheng | Efficient variational Bayesian approximation method based on subspace optimization[END_REF], source separation [START_REF] Choudrey | Variational Methods for Bayesian Independent Component Analysis[END_REF] or super-resolution [START_REF] Babacan | Variational Bayesian super resolution[END_REF]. Its popularity is mainly due to its ability to generate efficient approximations of the posterior distribution of the variables to be estimated. In image segmentation, VBA-based approaches traditionally uses a Potts (resp. Gaussian) model for hidden label variables (resp. pixel intensities corresponding to a given label) [START_REF] Ayasso | Joint NDT image restoration and segmentation using Gauss-Markov-Potts prior models and variational Bayesian computation[END_REF][START_REF] Chaari | Variational solution to the joint detection estimation of brain activity in fMRI[END_REF][START_REF] Mcgrory | Variational Bayes for estimating the parameters of a hidden Potts model[END_REF].

We propose a more generic Bayesian model, combining two main priors on the sought data. The first one concerns the latent label variables for which an HOGM is employed for clustering or classification [START_REF] Marinari | Cluster algorithms for the generalized 3d, 3q Potts model[END_REF][START_REF] Zheleva | Higher-order graphical models for classification in social and affiliation networks[END_REF]. The second one adopts a MEP distribution [START_REF] Gómez-Sánchez-Manzano | Sequences of elliptical distributions and mixtures of normal distributions[END_REF][START_REF] Marnissi | Generalized multivariate exponential power prior for wavelet-based multichannel image restoration[END_REF] for pixel intensities in a given class. To the best of our knowledge, these priors have not been jointly used for image recovery and segmentation tasks. The novel model is suitable for applications dealing with multicomponent or multivariate images (e.g. multi/hyperspectral data) [START_REF] Chaux | A nonlinear Stein based estimator for multichannel image denoising[END_REF][START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized Gaussian distributions with an application to multicomponent texture discrimination[END_REF][START_REF] Chierchia | A nonlocal structure tensor-based approach for multicomponent image recovery problems[END_REF]. One of the advantages of the proposed approach is that it allows us to estimate the associated hyperparameters by defining suitable priors. While alternatives could solve segmentation [START_REF] Pereyra | Estimating the granularity coefficient of a Potts-Markov random field within a Markov chain Monte Carlo algorithm[END_REF][START_REF] Sodjo | Joint segmentation of multiple images with shared classes: a Bayesian nonparametrics approach[END_REF][START_REF] Bioucas-Dias | Alternating direction optimization for image segmentation using hidden Markov measure field models[END_REF] or reconstruction [START_REF] Chouzenoux | A majorize-minimize subspace approach for 2-0 image regularization[END_REF][START_REF] Briceño-Arias | Proximal algorithms for multicomponent image processing[END_REF] problems, VBA benefits from low computational cost while providing high quality results.

Section 2 introduces the HOGMep Bayesian formulation for joint restoration and segmentation. The VBA methodology solving the related problem is described in Section 3. The proposed joint recovery and segmentation approach is favorably evaluated and compared within a deconvolution context in Section 4.

HOGMep BAYESIAN FORMULATION

Inverse problem model

We concentrate on the standard inverse problem consisting of recovering an unknown signal x from a degraded one y. We consider a linear model with additive noise formulated as y = Hx + n, where, for (M, N, B) ∈ (N * ) 3 , y ∈ R M is the observed data, x ∈ R N B corresponds to the unknown signal to be recovered, H ∈ R M ×N B represents a linear degradation operator and n is a noise (statistically independent of x). We are interested in B-component images where x = [x 1 , . . . , x N ] and, for every pixel i ∈ {1, . . . , N }, xi = (xi,1, . . . , xi,B) . Assuming a zero-mean white Gaussian noise with inverse variance (or precision) γ, the likelihood p(y|x, γ) is given by the normal distribution N (Hx, γ -1 I). In conjunction with the recovery task, our objective is to perform a classification of the components of x. For this purpose, L being the number of expected classes, the label field is encoded by a vector of hidden variables z ∈ {1, . . . , L} N .

Prior model on the unknown data and hyperpriors

HOGM for label variables [START_REF] Marinari | Cluster algorithms for the generalized 3d, 3q Potts model[END_REF][START_REF] Zheleva | Higher-order graphical models for classification in social and affiliation networks[END_REF] deals with cliques of arbitrary size instead of being limited to pairwise interactions between variables. The chosen prior is:

p(z) ∝ exp   S s=1 (i 1 ,...,is)∈Ns Vs(zi 1 , . . . , zi s )   , (1)
where S is the size of the maximal clique and, for every s ∈ {1, . . . , S}, the function Vs is a potential function of order s, and Ns is the set of cliques of size s. The model contains a prior weighting parameter λ, not explicitly written in [START_REF] Jordan | An introduction to variational methods for graphical models[END_REF].

In addition to a prior on hidden label variables, we assume that the conditional distribution of x given label variables z is a Multivariate Exponential Power (MEP) distribution denoted by M. Such a prior is well suited to multicomponent images [START_REF] Marnissi | An auxiliary variable method for Langevin based MCMC algorithms[END_REF]. We denote by zi ∈ {1, . . . , L} the label of the class to which the pixel indexed by i belongs, and by z = (z1, . . . , zN ) the label vector for all pixels of the image. Furthermore, for every class labeled by l ∈ {1, . . . , L}, m l , Ω l and β l denote the parameters of the MEP distribution associated with label value l. By assuming that β l ∈ (0, 1], the MEP distribution can be expressed as a Gaussian Scale Mixture (GSM) [START_REF] Gómez-Sánchez-Manzano | Multivariate exponential power distributions as mixtures of normal distributions with Bayesian applications[END_REF]:

M(xi; m l , Ω l , β l ) = R + N (xi; m l , u -1 i Ω -1 l )p(ui|β l )dui,
where p(ui|β l ) denotes the probability density function (pdf) of ui with shape parameter β l and u = (u1, . . . , uN ) gathers all the latent variables. When β l < 1, this pdf can be expressed as a function of a positive alpha-stable distribution. When β l = 1, it degenerates into a Dirac distribution [START_REF] Gómez-Sánchez-Manzano | Multivariate exponential power distributions as mixtures of normal distributions with Bayesian applications[END_REF]. As a result, we have

p(x|z, m, Ω, β) = N i=1 p(xi|zi, m, Ω, β),
where p(xi|zi = l, m, Ω, β) = M(xi; m l , Ω l , β l ). Our Bayesian model involves four types of hyper-parameters: inverse noise variance γ, mean variables (m l ) 1≤l≤L , inverse covariance matrices (Ω l ) 1≤l≤L and shape parameters (β l ) 1≤l≤L . In this work, we assume that, for all the classes l ∈ {1, . . . , L}, the shape parameters of the MEP distribution are identical, i.e. β1 = • • • = βL = β. This single parameter is fixed in advance, according to prior knowledge. If G and W denote Gamma and Wishart distributions, to jointly estimate all the remaining parameters, we define a hyperprior on each of them:

p(γ) = G(ā, b), p(m l ) = N ( μ, Λ), p(Ω l ) = W( Γ, ν),
Figure 1 summarizes the dependency relationships between variables involved in this hierarchical graphical model.

Joint posterior distribution

Using Bayes' rule, the joint posterior distribution reads

p(x, u, z, γ, m, Ω|y, β) (2)
∝ p(y|x, γ)

N i=1 p(xi|zi, ui, m, Ω)p(ui|β) p(z)p(γ) × L l=1 p(m l )p(Ω l )
Its intricate form is due to the dependence between the unknown variables. VBA can circumvent this issue [START_REF] Ayasso | Joint NDT image restoration and segmentation using Gauss-Markov-Potts prior models and variational Bayesian computation[END_REF][START_REF] Chaari | Variational solution to the joint detection estimation of brain activity in fMRI[END_REF][START_REF] Mcgrory | Variational Bayes for estimating the parameters of a hidden Potts model[END_REF], and is indeed shown in the following section to provide an elegant solution to the Bayesian inference.

y x u m Ω z γ β μ, Λ Γ, ν λ ā, b
Fig. 1: Dependency relationships between variables.

VARIATIONAL BAYESIAN APPROACH

Let Θ = (Θj) 1≤j≤J be a vector containing all the variables (x, u, z, γ, m, Ω) to be estimated. As previously mentioned, we want to estimate the true posterior distribution p(Θ|y). The core idea of VBA is to approximate p(Θ|y) by a separable distribution, denoted by q(Θ), such that q(Θ) = J j=1 q(Θj). The minimization of the Kullback-Leibler (KL) divergence between q(Θ) and p(Θ|y) allows us to determine the optimal separable approximation:

KL(q(Θ)||p(Θ|y)) = q(Θ) ln q(Θ) p(Θ|y) dΘ. (3)
The optimal approximate distribution can be computed [START_REF] Šmídl | The Variational Bayes Method in Signal Processing[END_REF]:

(∀j ∈ {1, . . . , J}) q(Θj) ∝ exp ln p(y, Θ) q -Θ j , (

where q-Θ j = i =j q(Θi), and . q denotes the expectation value with respect to a probability distribution q. Implicit relations between pdfs q(Θj) 1≤j≤J generally prevent analytical expressions for q(Θ). These distributions are determined in an iterative way, by updating one of the separable components q(Θj) 1≤j≤J while fixing the others. In this work, the approximation takes the form:

q(Θ) = q(γ) N i=1 q(xi, zi)q(ui) L l=1 q(m l )q(Ω l ) , (5)
with q(xi, zi) = q(xi|zi)q(zi). Hence, using (4), for every i ∈ {1, . . . , N } and l ∈ {1, . . . , L}, the optimal solutions for q(xi|zi), q(zi), q(m l ), q(Ω l ) and q(γ) are such that

q(xi|zi = l) = N (η i,l , Ξ i,l ), q(zi = l) = π i,l , q(m l ) = N (µ l , Λ l ), q(Ω l ) = W(Γ l , ν l ), q(γ) = G(a, b). (6)
Except for q(ui), the optimization of these distributions can be performed iteratively, since they belong to known parametrized distribution families. The MEP prior distribution can be expressed as a GSM. Hence, the mean value of ui can be determined without an analytical expression for q(ui) [START_REF] Palmer | Variational EM algorithms for non-Gaussian latent variable models[END_REF][START_REF] Zheng | Wavelet based unsupervised variational Bayesian image reconstruction approach[END_REF]. This allows us to derive an approximate distribution of xi, which depends on the mean value of ui.

In the following, assuming that k ∈ N designates the iteration number, we describe how to update these distributions by deriving closed form expressions for their parameters.

Determination of model pdf q(xi, zi)

As mentioned in [START_REF] Choudrey | Variational Methods for Bayesian Independent Component Analysis[END_REF], for l ∈ {1, . . . , L}, q k+1 (xi|zi = l) is a Gaussian distribution. Using (4), its covariance matrix Ξ i,l and mean η i,l are given, at the iteration k + 1, by

Ξ k+1 i,l = γ k H i Hi + u k i Ω k l -1 (7)
η k+1 i,l = Ξ k+1 i,l γ k H i (y - j<i Hj x k+1 j - j>i Hj x k j ) + u k i Ω k l µ k l . (8)
Hi ∈ R M ×B is the columnwise decomposition of H. Furthermore, for an arbitrary variable a, a k is its expectation at iteration k with respect to the pdf q k (a). We can then derive the expressions of the probabilities q k+1 (zi = l):

q k+1 (zi = l) = π k+1 i,l ∝ |Ξ k+1 i,l ||Γ k l | 1/2 exp 1 2 (η k+1 i,l ) (Ξ k+1 i,l ) -1 η k+1 i,l + 1 2 B b=1 ψ ν k l + 1 -b 2 - 1 2 u k i tr[(Λ k l + µ k l (µ k l ) ) Ω k l ] + Ṽ k l , ( 9
)
where ψ is the digamma function and

V k l = V1(l) + S-1 s=1
Vs+1(l, zi 1 , . . . , zi s )

j =i q k (z j )
.

For every j ∈ {1, . . . , N }, we have x k+1 j = L l=1 π k+1 j,l η k+1 j,l . Instead of computing the q(ui) distribution, we focus on the expectation of ui which can be expressed as follows, thanks to the integral form of GSM [START_REF] Palmer | Variational EM algorithms for non-Gaussian latent variable models[END_REF]:

u k+1 i = β L l=1 π k+1 i,l tr[A k Ω k l ] β-1 , ( 10
)
where

A k = (η k+1 i,l -µ k l )(η k+1 i,l -µ k l ) + Ξ k+1 i,l + Λ k l .
3.2. Finding of hyperprior pdfs q(m l ), q(Ω l ) and q(γ)

As shown by ( 6), q k+1 (m l ) is a Gaussian distribution with mean µ k+1 l and covariance matrix Λ k+1 l expressed as

Λ k+1 l = Λ-1 + Ω k l N i=1 π k+1 i,l u k+1 i -1 , (11)
µ k+1 l = Λ k+1 l Λ-1 μ + Ω k l N i=1 π k+1 i,l u k+1 i η k+1 i,l . (12)
Then, the mean value of q k+1 (m l ) is m k+1 l = µ k+1 l . The Wishart distribution q k+1 (Ω l ) is parametrized by

ν k+1 l = ν + N i=1 π k+1 i,l , (13)
Γ k+1 l = N i=1 π k+1 i,l u k+1 i Ãk + Γ-1 -1 , (14)
where

Ãk = (η k+1 i,l -µ k+1 l )(η k+1 i,l -µ k+1 l ) + Ξ k+1 i,l + Λ k+1 l . The mean value of q k+1 (Ω l ) is Ω k+1 l = ν k+1 l Γ k+1 l . (15)
The q k+1 (γ) distribution is a Gamma distribution with parameters a k+1 = ā + M 2 and

b k+1 = b + 1 2 y-H x k+1 2 + N i=1 L l=1 π k+1 i,l tr[H i HiΞ k+1 i,l ]. (16)
From standard properties of the Gamma distribution, the expectation of γ at the iteration k + 1 is equal to

γ k+1 = a k+1 b k+1 = 2ā + M 2b k+1 . (17)

Resulting algorithm

Altogether, our algorithm can be summed up as follows:

Set initial values:

η 0 i,l , Ξ 0 i,l , u 0 i , π 0 i,l , µ 0 l , Λ 0 l , Γ 0 l , ν 0 l , b 0 , and set a k ≡ ā + M 2 . Compute x 0 i = L l=1 π 0 i,l η 0 i,l , Ω 0 l = ν 0 l Γ 0 l , and γ 0 = 1/b 0 (ā + M/2).
For k = 0, 1, . . .

Update parameters Ξ k+1

i,l and η k+1 i,l of q k+1 (xi|zi = l) using ( 7) and [START_REF] Ayasso | Joint NDT image restoration and segmentation using Gauss-Markov-Potts prior models and variational Bayesian computation[END_REF]. Compute π k+1 i,l from (9).

Update mean values u k+1

i of q k+1 (ui) using (10). 3. Update parameters Λ k+1 l and µ k+1 l of q k+1 (m l ) using ( 11) and ( 12). 4. Update parameters ν k+1 l and Γ k+1 l of q k+1 (Ω l ) using ( 13) and [START_REF] Marnissi | Generalized multivariate exponential power prior for wavelet-based multichannel image restoration[END_REF]. Compute Ω k+1 l from (15). 5. Update parameter b k+1 of q k+1 (γ) using [START_REF] Verdoolaege | Geodesics on the manifold of multivariate generalized Gaussian distributions with an application to multicomponent texture discrimination[END_REF]. Compute γ k+1 from (17).

SIMULATION RESULTS

The performance of HOGMep is assessed through segmentation and recovery tasks in a deconvolution context. Experiments are performed on both synthetic and benchmark color images. Channels are degraded by distinct blur operators and additive Gaussian noises.

A natural comparison is made with VB-MIG [START_REF] Ayasso | Joint NDT image restoration and segmentation using Gauss-Markov-Potts prior models and variational Bayesian computation[END_REF], in which our MEP prior is restricted to a Gaussian one. Restoration results are also compared to those obtained with the recent state-of-the-art variational approach 3MG [START_REF] Chouzenoux | A majorize-minimize subspace approach for 2-0 image regularization[END_REF][START_REF] Chouzenoux | A Majorize-Minimize memory gradient algorithm applied to Xray tomography[END_REF] and segmentation with the iterated conditional mode (ICM) algorithm [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF]. For the 'synth' image with ground truth, results in SNR are given in Segmentation performance, summarized in Table 1b, is numerically assessed through variation of information (VI, [START_REF] Meilȃ | Comparing clusterings by the variation of information[END_REF]) and Rand index (RI, [START_REF] Rand | Objective criteria for the evaluation of clustering methods[END_REF]) measures. Segmented images and their binary difference to the original segmentation are displayed in Figure 4. HOG-Mep offers a very low rate of wrongly assigned pixels, and globally correctly labels all regions. For the 'peppers' image, Table 2 gathers objective SNR and SSIM performance, while restoration is depicted in Figure 5, and segmentation results in Figure 6. Through the obtained results, we conclude that HOGMep outperforms state-of-theart methods in both recovery and segmentation problems.

ICM VB-MIG HOGMep

Fig. 6: Segmentation results for 'peppers' and the corresponding silhouette. Silhouette image: blackish pixels reflect silhouette indices close to -1 while whitish pixels correspond to scores near 1. Using this color code, brighter zones reflect satisfying clustering.

CONCLUSIONS

In this work, we have performed a joint image recovery and segmentation via a Bayesian approach. This formulation combines two main priors on the sought data: an HOGM for the latent label variables and a MEP distribution for the unknown signal within a given class. We have proposed an original VBA method for solving the related problem where most of the hyper-parameters are automatically estimated. The provided simulation results show that HOGMep delivers good results for both recovery and segmentation tasks.

Fig. 2 :Fig. 3 :Fig. 4 :

 234 Fig. 2: RGB channels ('synth'): (a) original (b) blurred and noisy. Restored with (c) 3MG, (d) VB-MIG, (e) HOGMep. Degraded 3MG HOGMep

Fig. 5 :

 5 Fig. 5: Restoration results for the color images ('peppers').

Table 1a

 1a

					, and displayed in
	Figure 2 (restored channels) and 3 (color image). HOGMep consis-
	tently yields the best objective results, with crisper images.
	Channel	Initial 3MG VB-MIG HOGMep
	Red	10.73 11.21	22.62	24.25
	Green	12.42 15.17	16.41	17.14
	Blue	3.92	2.60	18.39	19.75
	Color	6.98	6.08	19.41	20.63
	(a) Channel and color restoration results in terms of SNR. Best per-
	formers are in bold.
			ICM	VB-MIG HOGMep
	VI		0.344	0.13	0.006
	RI		92.72	98.24	99.98
	% misclass.	41.10	4.70	0.05
	(b) Segmentation: variation of information [30], Rand index [31].

Table 1 :

 1 Objective restoration and segmentation performance ('synth').

Table 2 :

 2 Objective restoration performance ('peppers'). Best performers are in bold.

	Channel		Initial 3MG VB-MIG HOGMep
	Red	SNR 24.30 28.99 SSIM 0.938 0.972	19.94 0.939	33.43 0.972
	Green	SNR 18.36 22.82 SSIM 0.824 0.903	21.13 0.875	23.96 0.904
	Blue	SNR 13.77 14.19 SSIM 0.688 0.748	14.12 0.779	14.30 0.792
	Color	SNR 19.70 22.25	21.77	23.00
	Degraded	3MG	HOGMep
			VB-MIG		Original