Loïc Le Folgoc

Hervé Delingette

Antonio Criminisi

Nicholas Ayache

Quantifying Registration Uncertainty with Sparse Bayesian Modelling

Keywords: Registration, Sparse Bayesian Learning, Uncertainty Quantification, MCMC, Reversible Jump, Automatic Relevance Determination

We investigate uncertainty quantification under a sparse Bayesian model of medical image registration. Bayesian modelling has proven powerful to automate the tuning of registration hyperparameters, such as the trade-off between the data and regularization functionals. Sparsity-inducing priors have recently been used to render the parametrization itself adaptive and data-driven. The sparse prior on transformation parameters effectively favors the use of coarse basis functions to capture the global trends in the visible motion while finer, highly localized bases are introduced only in the presence of coherent image information and motion. In earlier work, approximate inference under the sparse Bayesian model was tackled in an efficient Variational Bayes (VB) framework. In this paper we are interested in the theoretical and empirical quality of uncertainty estimates derived under this approximate scheme vs. under the exact model. We implement an (asymptotically) exact inference scheme based on reversible jump Markov Chain Monte Carlo (MCMC) sampling to characterize the posterior distribution of the transformation and compare the predictions of the VB and MCMC based methods. The true posterior distribution under the sparse Bayesian model is found to be meaningful: orders of magnitude for the estimated uncertainty are quantitatively reasonable, the uncertainty is higher in textureless regions and lower in the direction of strong intensity gradients.

I. INTRODUCTION

Non-rigid image registration is an ill-posed task that supplements limited, noisy data with 'inexact but useful' prior knowledge to infer an optimal deformation between images of interest [START_REF] Sotiras | Deformable medical image registration: A survey[END_REF]. As a standard processing step in many pipelines for medical imaging, for computational anatomy & physiology, registration would benefit from the development of principled strategies to analyze its output and subsequently re-evaluate model assumptions. Bayesian modelling provides a framework to explicitly incorporate prior assumptions and reassess their relevance in retrospect. We focus here on another expected benefit of Bayesian approaches that is, the possibility to quantify uncertainty in the optimal solution.

Probabilistic approaches to registration and uncertainty quantification are not yet widespread in the literature. Gee and Bajcsy [START_REF] Gee | Elastic matching: Continuum mechanical and probabilistic analysis[END_REF] laid the groundwork for a Bayesian interpretation of registration, extending the mechanical formulation of Broit [START_REF] Broit | Optimal registration of deformed images[END_REF]. Exploiting the Gaussian Markov random L. Le Folgoc is with the Asclepios Research Project, Inria Sophia Antipolis and with the Microsoft-Inria Joint Centre, France. H. Delingette and N. Ayache are with the Asclepios Research Project, Inria Sophia Antipolis, 2004 route des Lucioles BP 93, 06902, Sophia Antipolis, France.

A. Criminisi is with the Machine Learning and Perception Group, Microsoft Research Cambridge, United Kingdom. field structure inherited from a finite-element discretization of the domain, they characterize the posterior distribution of displacements by Gibbs sampling. Risholm et al [START_REF] Risholm | Bayesian characterization of uncertainty in intra-subject non-rigid registration[END_REF] extend the approach to the case of unknown confidence on the observed data and on model priors respectively, aiming to address the critical issue of finding an objective trade-off between data fit and regularity-inducing priors. The so-called temperature hyperparameters are treated as latent variables and approximately marginalized over, while a Markov chain with full dimensional Metropolis-Hastings transitions traverses the space of transformation parameters. The aforementioned authors proceed in the framework of Markov Chain Monte Carlo (MCMC) sampling to explore the posterior distribution of model parameters. MCMC sampling yields an arbitrarily good characterization of the posterior provided that enough samples can be drawn within the available computational budget -inference becomes exact in non finite time. In practice, the computational burden and the technicality of the Markov chain implementation quickly become limiting factors. As an alternative, Variational Bayes (VB) inference provides tools to efficiently approximate the (true) posterior distribution on a chosen family of variational (approximate) posteriors. The choice of variational posterior realizes a trade-off between the computational burden and the quality of the estimates. Using a parametric (FFD) representation of the displacement field, Simpson et al [START_REF] Simpson | Probabilistic inference of regularisation in non-rigid registration[END_REF], [START_REF] Simpson | Probabilistic nonlinear registration with spatially adaptive regularisation[END_REF] approximate the posterior distribution within a 'convenient' family for which transformation parameters and model hyperparameters factorize. The variational factorization renders the approach applicable to real scale registration tasks. As a drawback estimates of uncertainty quantify variability in the displacement field conditionally to the inferred hyperparameters, but disregard uncertainty induced by hyperparameter variability. Although uncertainty quantification is peripheral to their work, Richard et al [START_REF] Richard | A SAEM algorithm for the estimation of template and deformation parameters in medical image sequences[END_REF] develop for the related task of atlas building a mixed SAEM and MCMC approach where nodes of the finite-element mesh are updated via Metropolis-Hastings-Within-Gibbs transitions; Zhang et al [START_REF] Zhang | Bayesian estimation of regularization and atlas building in diffeomorphic image registration[END_REF] implement a mixed SAEM and Hybrid Monte Carlo approach for a Bayesian MAP estimation of the template and of temperature hyperparameters in a diffeomorphic setting.

In this paper, we compare the approximate posterior returned by a Variational Bayes method with an MCMC method based on the same underlying model. We focus on the Bayesian model of registration developed in earlier work [START_REF] Le Folgoc | Sparse bayesian registration of medical images for self-tuning of parameters and spatially adaptive parametrization of displacements[END_REF]. The main goal of this model is to allow not only for the automatic determination of registration parameters (such as the trade-off between image similarity and regularization functionals), but also for a data-driven, multiscale, spatially adaptive parametrization of deformations via the recourse to a sparsity-inducing prior on transformation parameters.

Our contribution is twofold. Firstly, the complexity of the model renders inference non trivial. While in our previous work approximate inference was conducted on the grounds of Variational Bayes, we adopt here an exact MCMC-based approach. At a high level, the space of transformation parameters is explored by a reversible jump Markov chain [START_REF] Green | Reversible jump markov chain monte carlo computation and bayesian model determination[END_REF]. It provides a principled mechanism to elegantly jump between competing parametrizations of the displacement field, regardless of their dimensionalities, without the prohibitely expensive computation of so-called Bayes factors. This allows to seamlessly refine the parametrization of the transformation, adapting the granularity of the parametrization to the granularity of the underlying motion and the local informativeness of the image, all the while exploring the most likely deformations. At a lower level, we capitalize on closed form marginalization of most nuisance variables, and integrate second-order knowledge of the posterior distribution in proposal kernels. This yields an algorithm that reliably and consistently traverses the parameter space towards the most likely deformations in spite of the model intricacies.

Secondly, we compare the expectation and uncertainty predicted by both the fast (approximate) Variational Bayes inference and the (asymptotically) exact MCMC inference scheme both on empirical and theoretical grounds. We found that the expectation is typically well approximated by the VB inference, but that the uncertainty is underestimated. We exhibit two mechanisms that explain this behaviour. Furthermore we show that uncertainties predicted by the exact model are consistent with intuition: the orders of magnitude are sound, the uncertainty is higher in textureless regions and lower in the direction of strong intensity gradients.

The article unfolds as follows. In part II we describe the sparse Bayesian model of registration and devise a principled strategy for exact inference. The proposed design of the Markov chain exploits insight gained about the model to bypass standard impediments of MCMC schemes. Hyperparameter uncertainty is fully accounted for by marginalization of the nuisance variables. In part III we review breakdown scenarii in which the approximate posterior significantly departs from the true posterior, leading to poor approximate predictive uncertainty. In part IV we conduct preliminary experiments to assess the validity of MCMC uncertainty estimates.

II. STATISTICAL MODEL AND INFERENCE

Registration infers, from prior knowledge and limited data D, a transformation of space Ψ that pairs homologous features in objects of interests (e.g. organs or vessels, in a medical setting). The section starts with a succinct description of the registration model, and offers insight into its mechanisms. Fig. 1 provides a graphical representation thereof. An MCMC approach for systematic characterization of the posterior distribution is then devised. For pairwise registration of a fixed image J and a moving image I, a mixture-of-Gaussians model (GMM) of intensity residuals is adopted here as a flexible and robust variant of the widespread sum of squared differences (SSD). Fig. 2 summarizes this model of data in graphical form, making explicit the nodes P, H P , D of Fig. 1. The multiple components of the GMM naturally cope with the fact that intensity residuals may rightfully take high outlier values for an undeterminate fraction of voxels, because of acquisition artefacts, heteroscedastic noise and model inaccuracies. At voxel center v i , the intensity residual r i = J(v i ) -I[Ψ -1 (v i )] is assigned to the lth component of the mixture, 1 ≤ l ≤ L, if the L-way categorical variable c i ∈ {1 • • • L} takes value l. If so the residual r i follows a normal distribution N (0, β -1 l ). The component assignment c i follows a categorical distribution and takes value l with probability ρ l , normalized such that L l=1 ρ l = 1. For distinct voxels v i and v j , residuals r i and r j (resp. component assignments c i and c j ) are assumed to be independent. The corresponding GMM energy E D (D, Ψ; β, ρ) is given by Eq.

A. Bayesian Model of Registration

E D (D, Ψ; β) = β 2 N i=1 T i -Ψ(t i ) 2 . (1)
(2), with Z l = 2π/β l a normalizing constant:

- N i=1 log L l=1 ρ l Z l exp - β l 2 J[v i ] -I[Ψ -1 (v i )] 2 (2)
Fig. 3 shows the typical profile of the GMM energy comparatively with the SSD. The assumption of independence of voxelwise residuals is known not to hold (see e.g. [START_REF] Simpson | Probabilistic inference of regularisation in non-rigid registration[END_REF], [START_REF] Le Folgoc | Sparse bayesian registration of medical images for self-tuning of parameters and spatially adaptive parametrization of displacements[END_REF]) and to affect the outcome of the probabilistic registration. Since a proper probabilistic account of correlations in intensity residuals is both beyond the scope of this work and irrelevant to the ensuing developments, the Virtual Decimation scheme of [START_REF] Simpson | Probabilistic inference of regularisation in non-rigid registration[END_REF] is reproduced instead for simplicity.

2) Transformation parametrization: A small deformation standpoint is adopted for convenience. The displacement field

u : x ∈ Ω ⊂ R d → u(x) = Ψ -1 (x) -x ∈ R d is parametrized by a linear combination of M basis functions φ k (.) with associated weight w k ∈ R d : u(x) = 1≤k≤M φ k (x)w k = φ(x) w .
(3)

φ(x) = φ 1 (x) • • • φ M (x) and w = w 1 • • • w M are respectively the concatenation, for k = 1 • • • M , of φ k (x)
and w k . Arbitrary choices of basis functions φ k are possible. B-splines (e.g. [START_REF] Rueckert | Nonrigid registration using free-form deformations: application to breast MR images[END_REF]) present desirable properties in terms of smoothness and interpolation. Here the φ k 's instead consist of multiscale Gaussian radial basis functions (RBFs) whose centers lie on a regular grid of points (typically, decimated voxel centers). Multiscale Gaussian RBFs possess attractive analytical and computational properties.

3) Transformation priors: The weights w are endowed with a generalized Spike-&-Slab prior that favours both smoothness of the resulting displacement field and sparsity in its parametrization. The properties of this prior are central to the proposed 'sparse Bayesian' modelling and to our analysis thereof. Each basis φ k is assigned a distinct activation variable z k that controls its inclusion in the active parametrization (or exclusion therefrom). If z k = 0 the basis φ k is pruned out of the active parametrization. We do so by designing p(w k |z k = 0) as a Dirac distribution centered at 0. If z k = 1 the basis φ k is included in the parametrization. The prior on such bases is designed as a joint, structured Gaussian distribution that penalizes lack of smoothness in the induced displacement field [START_REF] Ashburner | A fast diffeomorphic image registration algorithm[END_REF]. Let us denote by S the set of such indices k for which z k = 1 and by w S the concatenation of the corresponding subset of weights {w k , k ∈ S}. For an arbitrary linear differential operator D, we wish to penalize high values of the quadratic energy Du 2 = w S R S w S , where R S is the |S| × |S| matrix whose k, l-th coefficient is Dφ k |Dφ l . The Gaussian distribution N (w S 0, {λ d|S|} -1 R -1 S ) is a natural choice of prior for p(w S |S), that we adopt henceforth. Note that the covariance normalization by d|S|, where d is the image dimension, departs from that of [START_REF] Le Folgoc | Sparse bayesian registration of medical images for self-tuning of parameters and spatially adaptive parametrization of displacements[END_REF]. Under this prior λ d|S| • w S R S w S is χ 2 (d|S|) distributed so that λ immediately relates to the expectation of the energy:

E p(w S |S) ( Du 2 ) = λ -1 and E p(w) ( Du 2 ) = λ -1 .
The prior over all weights w conditioned on the state of the gate variables z = z 1 • • • z M is best summarized in the form of Eq. ( 4), where -S is the complement of S:

p(w|z, λ) = N (w S 0, 1 λ d|S| R -1 S ) • N (w -S 0, 0) . (4)
4) Hyperpriors: Parameters introduced in the specification of priors are in turn treated as latent variables. λ is endowed with a Gamma prior Γ(λ|a 0 , b 0 ) that is conjugate to p(w|z, λ). The parameters β l (resp. β) involved in the likelihood model for image (resp. landmark) registration are endowed with independent Gamma priors Γ(β l |γ 0 , δ 0 ). The noise mixture

proportions ρ = {ρ 1 • • • ρ L } are assigned a Dirichlet prior Dir(ρ|κ), with κ = (κ 1 • • • κ L ).
Independent Bernoulli priors B(z k |π k ) on each z k constitute a natural, conjugate hyperprior specification for the activation variables z. The positive mass 1 -π k concentrated at w k = 0 as a result explicitly encodes sparsity. Assuming all π k = π 0 to be equal, all parametrizations using the same number of active bases |S| are a priori equally probable. In addition the cost of including a new basis in the active parametrization is independent of the current number of active bases. However, we opt instead for a stronger prior, p(z) ∝ Γ( d|S| 2 ) -1 . The Gamma function Γ(•) is a natural extension of the (integer) factorial to real values, yielding a prior that increasingly penalizes each new inclusion. This prior was found to perform better w.r.t. sparsity, as can be theoretically argued from the analysis of the marginal prior p(w|z).

B. Model analysis

1) Marginal prior and marginal likelihood: Critical insight into the statistical model can be gained by considering the prior p(w|z, H) and likelihood p(D|w, c, H) with so called temperature parameters λ and β marginalized over, e.g.:

p(w|z, H) = R+ p(w|z, λ, H)p(λ|H)dλ . (5)
The multivariate Student distribution t ν (•|µ, Λ) with location parameter µ, inverse scale matrix Λ and ν degrees of freedom naturally appears in analytic derivations, yielding the following expressions for the prior and likelihood:

p(w|z, H) = N w -S 0, 0 t ν λ w S 0, a 0 b 0 d|S| R S (6)
p(D|w, c, H) = L l=1 t ν l I l • Ψ -1 w J l , γ 0 δ 0 I (7)
where ν λ = 2a 0 , ν l = 2γ 0 , S is the set of active bases and

|S| = k z k its cardinal. J l = • • • J[v i ] • • • i|ci=l is the
vector of voxel values in image J, for those voxels assigned to component l, and

I l • Ψ -1 w = • • • I[Ψ -1 w (v i )] • • • i|ci=l is similarly defined for the warped image I • Ψ -1
w . For a fixed choice of active bases z, the posterior distribution of the weights p(w|z, c, D, H) is proportional to the product of the prior Eq. ( 6) and likelihood Eq. (7). In the limit of uninformative hyperpriors a 0 , γ 0 → 0, β 0 , δ 0 → 0 and assuming L = 1 for the sake of illustration,

p(w|z, c, D, H) ∝ N w -S 0, 0 1 χ lik [w] N 1 χ pr [w] d|S| . (8)
where χ lik [w] 2 is the data error and χ pr [w] 2 = Du w 2 the regularizing energy. In particular the posterior distribution is invariant to rescaling of the data error, and hence to rescaling of the intensity profile, after marginalizing over temperature parameters. Note also that, for a fixed parametrization z, the ratio of posterior probabilities of two distinct parameter sets w 1 and w 2 may become arbitrarily overwhelmed by the prior as the number of bases in the parametrization grows (|S| N ). If not for sparsity, this might render MCMC characterization of the posterior unreliable (using e.g. Metropolis Hastings transitions), potentially making its outcome dependent on the size of the parametrization. Fortunately the proposed sparse model has a clear mechanism to prevent overparametrization and render overlapping bases largely mutually exclusive, as discussed next.

2) Prior probability of basis inclusion: Interactions between overlapping bases can be better understood by looking at the probability p(z k |w -k , z -k , H) of inclusion of a new basis z k given a known configuration z -k for the other bases and their associated weights w -k . The state w -k of other bases informs us about the expected regularity of the signal u w , introducing dependencies between z k and z -k conditionally to w -k . Denoting by z (resp. z) the state with z k = 1 (resp. z k = 0), we see from Bayes' rule that:

p(z k = 1|w -k , z -k ) p(z k = 0|w -k , z -k ) = p(w -k | z) p(w -k |z) p( z) p(z) (9)
where the dependence on hyperparameters is made implicit for convenience of notations. Leaving details of derivations aside, we note that in the limit of uninformative values, the ratio of Eq. ( 9) takes the form of

p( z) p(z) |κ k | |R k,k | 1/2 1 - µ k pr R k,k µ k pr w -k R S w -k - d|S| 2 ( 10
)
where S is the set of active bases (excluding k),

µ k pr = -R -1 k,k R k w -k and κ k = R k,k -R k R -1 S R k .
The middle factor penalizes the inclusion of basis k if it overlaps with bases in the active set S, in the sense of the metric induced by R. κ k is a measure of overlap of basis k with all bases in the active set S and is null if basis k is perfectly collinear to S. The rightmost factor favors the inclusion of basis k if it is a priori expected to yield a significant increase in regularity.

C. Posterior Exploration by MCMC Sampling

For any set of points X = {x 1 • • • x n } in the admissible domain Ω, consider the vector of displacements u X = u(x 1 ) • • • u(x n ) . We wish to characterize the joint posterior distribution p(u X |D, H) of any such vector of displacements for any discrete set X. To that aim we merely need to characterize the posterior distribution p(w|D, H) of the weights w involved in the parametrization of the transformation Ψ -1 sufficiently well.

1) Related work: MCMC methods are tools of predilection to explore arbitrarily complex distributions in a principled manner. Gibbs sampling [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the bayesian restoration of images[END_REF] cycles between latent variables, sampling from their conditional distributions in turn while other model variables remain fixed. It is attractive when conditional distributions are known in closed form whereas the joint distribution is untractable or computationally costly to sample. When the conditional cannot be sampled directly, a component-wise proposal may be used instead within a Metropolis-Hastings (MH) step (Metropolis-Within-Gibbs). Unfortunately, Gibbs sampling of temperature parameters is prone to failure, with the chain drifting away from regions of high probability for the duration of any finite MCMC run. Collapsing temperature parameters λ, β when sampling regressor variables w is highly opportune. In the context of registration, Risholm et al. [START_REF] Risholm | Bayesian characterization of uncertainty in intra-subject non-rigid registration[END_REF] propose a MH scheme where marginalizing over temperature parameters induces the expensive computation of partition functions, for which an intricate procedure based on Laplace approximations is designed. In the proposed model, the computation of partition functions (specifically, marginal likelihoods, a.k.a. evidences) may arise as well when sampling gate variables z k . Selecting a specific configuration z can be interpreted as a choice between competing models of varying complexity and dimensionality. The problem of estimating the evidence for a model is well studied in the statistical literature. A variety of methods exist, ranging from the straightforward Laplace approximation to more principled approaches typically exploiting samples from the (possibly augmented) posterior, including Chib's method [START_REF] Chib | Marginal likelihood from the gibbs output[END_REF], importance sampling, bridge sampling, path sampling (see e.g. [START_REF] Gelman | Simulating normalizing constants: From importance sampling to bridge sampling to path sampling[END_REF]) and reversible jump MCMC [START_REF] Green | Reversible jump markov chain monte carlo computation and bayesian model determination[END_REF]. The latter approach is in fact primarily concerned with sampling from a posterior distribution involving competing models (freely jumping between models in the process) and merely obtains evidence ratios as a byproduct. Reversible jump MCMC is appealing in our setting where competing models z are organized in series of nested models of increasing complexity, rendering its machinery mostly invisible. Reversible jump MCMC proceeds in the general framework of Metropolis-Hastings, hence a sound proposal must be crafted. We derive a sensible family of proposals from a modal analysis of the posterior distribution.

2) Modal analysis of the posterior & proposal: For the model described in II-A, the Laplace approximation of the (conditional) posterior p(w|D, z, c, H) around its mode w * = arg max w p(w|D, z, c, H) takes the following form:

-log p(w|D, z, c, H) ≈ 1 2 γ 0 + N/2 δ 0 + χ 2 lik /2 (T * -Φw) H * (T * -Φw)+ 1 2 a 0 + |S|/2 b 0 + χ 2 pr |S|/2 d|S| w R S w + const . ( 11
)
where for the sake of illustration we take a single component mixture (L = 1, c i = 1 for all i, β = β). χ 2 pr = w * R S w * is the energy in the displacement field, χ 2 lik is the data error

χ 2 lik = N i=1 (J[v i ]-I[Ψ -1 * (v i )]
) 2 and we discard higher order terms in b 0 , δ 0 .

T * = T 1 * • • • T N * is
β * = γ 0 + N/2 δ 0 + χ 2 lik /2 , λ * = a 0 + |S|/2 b 0 + χ 2 pr |S|/2 (12)
are commensurable to temperature parameters. The approximation of the conditional posterior is Gaussian (Eq. ( 11) is quadratic) and admits the more obvious canonical form

N (µ, Σ), with µ = ΣΦ (β * H * )T * and Σ = (Φ β * H * Φ + λ * |S|R S ) -1 .
The Laplace approximation provides a reasonable approximation of the posterior and a judicious starting point to design proposals. Component-wise proposals that leave most of the activation variables z l and the corresponding weights w l unchanged will be of particular interest to us (cf. section II-C3). A natural idea is to use the conditionals wk ∼ N (µ k pos , Σ k ) of the Laplace approximation N (µ, Σ) as proposal distributions. Because they neither require the actual computation of µ and Σ nor involve inner products φ k (β * H * )φ l , these 'Gibbs-like' proposals are computationally appealing. As a final tweak to alleviate modal assumptions, we reintroduce dependency on the current value of w k , yielding the following component-wise proposal instead, with 0 ≤ r HMALA ≤ 1 and s ≥ 1:

q k (w k → wk ) = N ( wk | m k (w k ), sΣ k ) (13) m k (w k ) = (1 -r HMALA ) w k + r HMALA µ k pos ( 14
)
If not set to 1, the factor s accounts for potentially fatter tails of the true conditional posterior in the proposal. µ k pos and Σ k depend on H * and T * , which in the formal reasoning based on the Laplace approximation are computed around

Ψ -1 * (•) = Id+φ(•) w * .
In fact T * and H * can be replaced by T w and H w computed from a (local) quadratic approximation of p(w|D, z, c, λ * , β * ) around the current Ψ -1 (•) = Id + φ(•) w. In that case Eq. ( 13), ( 14) exactly coincide with a component-wise Hessian preconditioned Metropolis Adjusted Langevin Algorithm (HMALA) [START_REF] Roberts | Exponential convergence of langevin distributions and their discrete approximations[END_REF]- [START_REF] Zhang | Quasi-newton methods for markov chain monte carlo[END_REF], which exploits first and second order local information about the target distribution for increased efficiency. However the local approximation generates additional computations at each step and offers little gain if we expect the posterior to be unimodal. Given our experimental settings, we use the global approximation with adaptation during the burn-in phase (at that stage λ * , β * , T * and H * are recomputed every few iterations from statistics |S|, χ 2 pr , χ 2 lik averaged with decaying weights over past samples). 3) Reversible jump MCMC scheme: The groundwork for this scheme was laid in sections II-B1, II-B2, II-C2. The reversible jump procedure itself lets us generate samples of the joint posterior p(w, z, c|D, H) with temperature parameters marginalized over. Dropping irrelevant variables in the generated samples, we obtain samples of the marginals of interest, e.g. p(w|D, H). The reversible jump scheme simply proposes to move from a current state w, z, c to a new state w, z, c and computes a Metropolis-Hastings acceptance ratio for the proposal, leading to acceptance or rejection of the new state. For the sake of simplicity, proposals for a new state of w, z may be made separately from those of c. For the latter, the most natural proposal exactly results in collapsed Gibbs sampling of each c i , see e.g. [START_REF] Murphy | Machine learning: a probabilistic perspective[END_REF] 1 . For w, z we design basic moves that -when combined -allow to add, remove or switch active bases as well as update several components of w. These basic moves are combined to craft proposal distributions Q(w, z → w, z) for which the probability of a move w, z → w, z has direct symmetries with that of the reverse move w, z → w, z, so that the acceptance ratio min 1, p( w, z, c|D, H) p(w, z, c|D, H)

Q( w, z, c → w, z, c) Q(w, z, c → w, z, c) (15)
becomes particularly straightforward to compute. The basic moves are: a) Basis removal. For a basis k such that z k = 1, set zk = 0 and wk = 0. The symmetric move is the basis addition.

b) Component-wise update. For a basis k such that z k = 1, propose a new wk ∼ q k (w k → wk ) according to Eq. ( 13), if on-off or update then

Algorithm 1: Proposal Q k (w, z, reverse traversal → w, z, reverse traversal * ).
• For update: set I = {k}.

• For on-off : Set I ⊂ S\{k} to a list of nneighb active bases, favouring bases well-aligned with φ k . • If reverse traversal = 1, reverse the ordering of I. for l ∈ I do wnew l ∼ q l ( w → wnew ) and wl = wnew l if φ k inactive and on-off then Set zk = 1 and wk ∼ q k (w k → •), using rHMALA = 1. else if φ k inactive and exchange then Set zk * = 1 and wk * ∼ q k * (w k * → wk * ) (rHMALA = 1). if on-off or exchange then Switch the state of the binary variable reverse traversal. [START_REF] Chib | Marginal likelihood from the gibbs output[END_REF] with a fixed 0 ≤ r HMALA ≤ 1. This move is its own symmetric (using the reverse update). c) Basis addition. For a basis k such that z k = 0, set zk = 1 and propose a new wk according to Eq. ( 13), ( 14) with r HMALA = 1. The symmetric move is the basis removal.

The family of proposals Q k (•) that we design combines these basic moves in such a way that when reversed, the sequence of moves induced by the proposal Q k (w, z, c → w, z, c) coincides exactly with the sequence of moves induced by Q k ( w, z, c → w, z, c). The proposal and reverse proposal travel along the same path in opposite directions, drastically reducing the computational load when evaluating Eq. [START_REF] Gelman | Simulating normalizing constants: From importance sampling to bridge sampling to path sampling[END_REF]. Each proposal Q k revolves primarily around the corresponding basis φ k and is defined as per Algorithm 1 (where we introduced a binary variable reverse traversal to address technicalities). Using Q k , we define a transition kernel P k conventionally: given the current state w t , z t , c t , we propose a new state c = c t , w, z ∼ Q k (w t , z t → •). The state is accepted with probability given by Eq. ( 15), in which case we set (w t+1 , z t+1 , c t+1 ) = ( w, z, c); otherwise we stay at the current state and (w t+1 , z t+1 , c t+1 ) = (w t , z t , c t ). Computation of the acceptance ratio is relatively straightforward by construction, since the ratio of posterior probabilities involved in Eq. ( 15) can be rewritten as:

p(D| w, z, c, H) p(D|w, z, c, H) • p( w| z, c, H)p( z|H) p(w|z, c, H)p(z|H) (16)
The leftmost factor is a ratio of likelihoods and need only be evaluated once for a proposed transition. As the denominator is known from the previous iteration, only the numerator need be evaluated. In the context of registration, this part corresponds to the image term and would involve costly computations if evaluated repeatedly. Note also that for basis functions with compact support (or approximately so), only part of the image term need be updated to evaluate the ratio. The ratio on the right-hand side and the ratio of proposals are simply decomposed over the sequence of previously defined basic moves, then efficiently evaluated using Eq. ( 6), ( 13), ( 14) and expressions similar to Eq. ( 9), [START_REF] Green | Reversible jump markov chain monte carlo computation and bayesian model determination[END_REF]. For the latter, statistics κ k are kept up to date (for all bases) using efficient rank one updates derived in [START_REF] Le Folgoc | Sparse bayesian registration of medical images for self-tuning of parameters and spatially adaptive parametrization of displacements[END_REF]. Alternatively, the necessary statistic κ k can be recomputed from scratch only for the bases under consideration. This is usually much more efficient (cf. algorithmic complexity in II-C5). Each transition kernel P k satisfies a detailed balance condition. In terms of these transition kernels, the MCMC chain proceeds as follows. Random variables k 1 , k 2 , . . . taking values in {1, 2, . . . , M } are chosen according to some scheme and the corresponding transition kernel P kt is used at time t. Conventional schemes include the random-scan, where the {k t } are i.i.d uniform, and the deterministic scan that cycles through {1, 2, . . . , M } in natural order (see e.g. [START_REF] Roberts | Harris recurrence of metropoliswithin-gibbs and trans-dimensional markov chains[END_REF]). For the random scan, the global transition kernel also satisfies detailed balance conditions. For both schemes, the MCMC chain has stationary distribution p(w, z, c|D, H) after incorporating collapsed Gibbs updates of c. Highlights of the MCMC scheme main constituents are summarized in Fig. 10.

4) Markov chain mixing improvement: Similarly to Gibbs sampling of temperature parameters, Gibbs sampling of voxel GMM assignments c within updates separated from those of w, z potentially hampers the mixing of the Markov chain for any finite, practical duration of the MCMC run. If at any point in time, a data point that should be regarded as an outlier (e.g. an image artifact), or a group of such points, is assigned to a 'non-outlier' mixture component, the disjoint sampling generally causes the chain to remain stuck in the vicinity of the corresponding local mode of the posterior p(w, z, c|D, H): the desired reverse assignment move virtually occurs with probability zero after readjustment of w, z. This defect is critical as such failure scenarii happen with overwhelming probability. Fortunately, joint proposals for w, z, c can be designed at little cost, even more so after noting that the component-wise proposals for w k (Eq. ( 13), ( 14)) and z k only indirectly depend on c. The transition Q k (w, z, c → w, z, c) proceeds in two steps. First, w, z is proposed as per Algorithm 1. Then, c is sampled by component-wise collapsed Gibbs sampling of each ci ∼ p(c i | cj<i , c j>i , w, D, H) in turn. For efficiency, only the subset of voxels in the support of updated basis functions is sampled, and voxel assignments are updated only once in case of overlapping supports. The twostep move is accepted or rejected based on the acceptance ratio [START_REF] Gelman | Simulating normalizing constants: From importance sampling to bridge sampling to path sampling[END_REF], replacing c by c where necessary. The order of voxel traversal is reversed according to the state of reverse traversal. Sampling c and computing its contribution to the acceptance ratio exclusively involves the residuals r i and ri of updated voxels prior and after the update w, z → w, z, which were already required to compute the likelihood change in Eq. ( 16).

5) Algorithmic complexity: The algorithmic complexity associated to a transition kernel P k (proposal and acceptancereject) is O(|S| • |I + | + l∈I+ V l + L C), noting I + the set of updated bases, V l the number of voxels in the support of basis φ l and C the number of voxels whose assignments ci are resampled. The first term includes part of the cost of the proposal w, z → w, z and its impact on the ratio of prior probabilities. The second term is replicated three times and can be heavily parallelized in each case: once to compute µ l pos , Σ l in Eq. ( 13), ( 14) for l ∈ I + , twice to evaluate and store differences in the displacement fields (resp. residual images) over the support of basis functions in I + following their update. The last term accounts for all computations related to resampled voxel GMM assignments c. When a move that involves the inclusion or removal of a basis function from the active set is accepted, an additional O(|S| 2 + M • |S|) cost is involved to maintain statistics κ k over all bases in the dictionary, with the right-hand term being parallelizable into M disjoint O(|S|) operations. The O(M • |S|) cost upon inclusion or deletion of a basis can be replaced by a O(|S| 2 ) cost per proposed move, which is usually more efficient.

6) Initialization: The chain is initialized from the output of the deterministic algorithm presented in [START_REF] Le Folgoc | Sparse bayesian registration of medical images for self-tuning of parameters and spatially adaptive parametrization of displacements[END_REF] which progresses greedily in the space of parameters {z, λ, P } towards a local maximum of their joint posterior. We comment, however, that any registration algorithm could reasonably be used to initialize the chain.

III. PREDICTIVE UNCERTAINTIES: MARGINAL LIKELIHOOD MAXIMIZATION VS. EXACT INFERENCE

The 'sparse Bayesian' model presented in Fig. 1 is inspired by the Spike-&-Slab model of Mitchell and Beauchamp [START_REF] Mitchell | Bayesian variable selection in linear regression[END_REF] and the Relevance Vector Machine (RVM) proposed by Tipping [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF] for tasks of regression and classification. In the latter work, the author approaches the problem of inferring an optimal sparse regression function from the standpoint of Automatic Relevance Determination (ARD). Point estimates of the hyperparameters that govern basis selection (and in fact of all hyperparameters) are sought in a first step by maximizing the marginal likelihood or evidence as per Eq. ( 17):

θ * = arg max θ p(D|θ, H) = arg max θ w p(D|w, θ, H)p(w|θ, H)dw ( 17
)
where θ = {z, P, λ} using our notations. If non-uniform, proper hyperpriors on θ are assumed, θ * maximizes the posterior p(θ|D, H) ∝ p(D|θ, H)p(θ|H) instead. In a second step, the distribution of weights w k is characterized conditionally to the selected model,

p(w|D, H) ≈ p(w|θ * , D, H) . ( 18
)
This strategy is typically successful in reaching strongly sparse solutions with good predictive power but, above all else, is motivated by its computational efficiency. Dedicated schemes relying on linear algebra and rank one updates make it possible to efficiently, iteratively build the set |S| of relevant basis functions φ k from scratch. See for instance [START_REF] Tipping | Fast marginal likelihood maximisation for sparse bayesian models[END_REF], and [START_REF] Le Folgoc | Sparse bayesian registration of medical images for self-tuning of parameters and spatially adaptive parametrization of displacements[END_REF] for an extension to the wider family of priors required for registration tasks. The approximation of Eq. ( 18) is justified by observing that the full posterior p(w|D, H) is obtained by summing over all conditional posteriors p(w|θ, D, H), conditioned on the value θ, weighted by the posterior probability p(θ|D, H) for this value:

p(w|D, H) = θ p(w|θ, D, H)p(θ|D, H)dθ . (19)
Now if the available data D is informative enough, p(θ|D, H) will be sharply peaked around its mode(s). In the limit case where p(θ|D, H) is a Dirac centered at its single mode θ * , Eq. 18 is retrieved exactly, and the two-step scheme outlined in Eq. ( 17), ( 18) is justified. Moreover in the case of sparsity governing parameters z = z 1 • • • z M , Tipping [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF] argues that, even if several combinations of parameters are highly probable due to the presence of redundant functions φ k in the dictionary of bases, they should roughly lead to the same optimal solution u * and an approximate mode (or the expectation) of p(u|D, H) should still be correctly evaluated. Regardless, we now demonstrate why this evidence-based approximation will typically fail to properly approximate higher order moments of the full posterior, resulting for instance in poor approximation of the real predictive uncertainty. There are two main breakdown situations for the evidence-based approximation of the full posterior assumed in Eq. 18.

Firstly in absence of data, the assumption that the posterior distribution p(θ|D, H) of hyperparameters is well approximated by a Dirac collapses. Indeed the posterior then resembles the prior distribution p(θ|H), which is typically flat. This scenario is relevant in the case of basis selection parameters z k , since associated basis functions φ k have a local support over which reliable data may be missing. Away from data and without strong incentive to include the basis to increase the deformation regularity, the probability of basis inclusion (resp. exclusion) is π k (resp. 1 -π k ), and for neutral values of π k , the choice of excluding the basis is arbitrary. Secondly and even in the presence of data, many combinations of active bases could have quasi-identical probability. When using radial basis functions for instance, the location of basis centers can be slightly perturbed without significantly affecting the posterior probability of the new configuration. The optimal value of basis weights w under two such perturbations will slightly differ however, as well as the resulting transformation Ψ. The evidence-based approximation of Eq. ( 18) relies on a single -perhaps only marginally superiorconfiguration, whereas the true posterior sums over all such configurations, as seen from Eq. ( 19). As it turns out, 'basis wiggling' accounts for a significant part of the uncertainty.

IV. EXPERIMENTS AND RESULTS

The following experiments aim to qualitatively evaluate the consistency of posterior distributions inferred by the Variational Bayes approximate inference scheme and the MCMC asymptotically exact inference scheme.

A. Material & Experimental Setting

We focus on the 2D registration example of Fig. 5. For the approximate-based inference, the methodology of [START_REF] Le Folgoc | Sparse bayesian registration of medical images for self-tuning of parameters and spatially adaptive parametrization of displacements[END_REF] is used without change. The multiscale dictionary hence uses Gaussian RBFs at three different scales (isotropic, σ = 6mm, 12mm and 24mm), for a total of approximately 7 • 10 4 basis functions, of which no more than 50 -100 are typically active at a time (both with VB and MCMC approaches). We set the differential operator D to the Laplacian of the displacement field. The Gaussian Mixture model of intensity residuals has L = 5 components: hyperparameter values β * for the proposal distribution learned during the burn-in phase (sec. II-C2) indicate that 2 or 3 components would suffice. No strong dependence of the results on the number of components was observed. All hyperpriors use small uninformative values a 0 = b 0 = γ 0 = δ 0 = 10 -10 , and κ 0 = 0.5 for a noninformative (Jeffreys) Dirichlet prior. The MCMC chain was run for roughly 7 • 10 5 transitions and 500 samples were regularly extracted. Approximately 7 • 10 4 additional samples were discarded as part of the burn-in phase, during which the parameters of the proposal distribution were fine-tuned (cf. section II-C2). The tuning relies on a set of sufficient statistics, such as the average energy and the average voxelwise square intensity residuals per sample. The averages are computed using a scheme that downweights the early samples: a fixed learning rate is initially applied before reverting to a classical (inverse linear) weighting, drawing inspiration from the SAEM scheme of e.g. [START_REF] Richard | A SAEM algorithm for the estimation of template and deformation parameters in medical image sequences[END_REF]. The free parameter s controlling the spread of proposals compared to the second-order approximation of the posterior (section II-C2) was set to 1 (spread unchanged). The observed acceptance rate varied between 20-45% under sensible variations of the experimental setting, and between 27-34% during the run of interest with the settings described above. Examples of samples are reported in Fig. 6. As an order of magnitude, the run takes 10 minutes on a standard laptop with a naive implementation.

Finally to gain more insight into the behaviour of VB and MCMC approaches, we also experiment with an MCMC chain that proceeds as described above except for the choice of active basis functions which, instead of exploring various configurations, is fixed to that of the Variational Bayes approach. We refer to this experiment as Fixed Basis MCMC (FBMCMC). The setting is entirely identical to that of the full MCMC, but the only transitions proposed are component-wise updates as opposed to exchange, addition or removal of basis functions.

B. Results

1) Naive alternated sampling vs. joint sampling: Fig. 7 demonstrates the benefit of a careful design of the Markov chain. The left-most figure displays the estimated mean displacement, under the aforementioned experimental setting, if moves in the space of transformation parameters are done separately from the resampling of voxelwise assignments to components of the noise mixture instead of jointly (right-most figure). In this example, a local discrepancy in the intensity profiles of the fixed and moving images induces a spurious maximum in the joint posterior distribution of transformation parameters and voxel labels (cf. section II-C4). A systematic drift towards this mode was observed in all runs where the sampling was performed in an alternated manner, for the whole duration of the run, whereas systematic recovery was observed under the improved scheme. Similar observations were made in experiments where temperature parameters were treated by Gibbs sampling instead of analytically marginalized over.

2) VB vs. MCMC -estimated displacement: Fig. 8 reports the mean displacement reported respectively by the evidencebased inference scheme and by the MCMC inference scheme. As anticipated from the discussion of section III, very good agreement between the evidence-based and MCMC-based estimates of the displacement is observed. Upon close inspection, Fig. 6. Three example samples returned by the MCMC run. minor differences are noted in some areas with flat intensity profiles or otherwise low confidence (such as that resulting from artefacts, or disagreeing intensities in the fixed and moving image). Their magnitude is lower than the level of uncertainty in the output of registration, as estimated from the MCMC scheme.

3) VB vs. MCMC -uncertainty estimates: Fig. 9 compares the estimates of uncertainty obtained from the MCMC characterization of the posterior and those obtained from the Variational Bayes inference. For the MCMC inference, relevant statistics are estimated from the set of samples returned by the run. To study the spatial localization of uncertainty, we visualize at each voxel center x i the 2×2 covariance matrix of the posterior distribution p(u(x i )|I, J, H) of the corresponding displacement vector u(x i ). This is reasonable under the assumption that the posterior on displacements is approximately mono-modal and Gaussian. The voxelwise covariance matrix, or its square root (homogeneous to a standard deviation), can be visualized as a 2D tensor that encodes uncertainty at this point along any direction. Fig. 9 displays the resulting tensor map (bottom row) and a scalar summary (upper row).

On the one hand, the order of magnitude of uncertainties under the true posterior (typically ∼ 1mm for a 95% confidence interval), as estimated by MCMC sampling, is consistent with both the magnitude of the underlying motion (no more than 5mm, see fig. 8) and the resolution (voxel dimensions:

1.25mm × 1.25mm). As expected, uncertainty is higher in regions with little structured content (no intensity gradients) and in the direction of contours. On the other hand, the VB scheme does not appear to reliably approximate the true uncertainty. Its order of magnitude is generally underestimated. Moreover, VB-based uncertainty may lack spatial coherence in regions that are textureless, with a flat intensity profile (e.g. in the right ventricle on Fig. 5). This hints at the fact that, when relying on the evidence-based (VB) scheme, regions of high uncertainty are localized nearby the inferred (unique) set of active basis functions. 4) Fixed Basis MCMC: Fig. 8 (second row) and Fig. 9 (middle column) report the estimates of the mean and uncertainty for the Fixed Basis MCMC scheme. The estimated mean displacement is in good agreement with both approaches. Moreover the magnitude of the difference between FBMCMC and VB (resp. FBMCMC and MCMC) is generally below that of the residual displacement between VB and MCMC. The FBMCMC approach, similarly to the VB approach, underestimates uncertainty in regions of flat intensity (e.g. bottom right of the image) and displays small localized uncertainty peaks. The magnitude of the predicted uncertainty is globally consistent with that of the VB scheme (Fig. 9, similar tensor sizes in the first and second rows), albeit sometimes slightly superior, typically nearby active basis functions (e.g. in the anterior part of the right ventricle).

V. DISCUSSION

A. Markov Chain design for efficient and reliable inference

The proposed model of registration copes with various unknowns in the image and transformation model: the noise level and its spatial variability, the regularity of the hidden motion, and the optimal parametrization of the displacement field itself. This renders inference challenging and special care has to be taken in the design of the Markov chain. Typically several joint configurations of the regularization hyperparameters, the noise levels (and voxelwise mixture assignments) and the displacement field constitute local maxima in the joint parameter space. To prevent the chain from remaining stuck around poor local maxima, it proved useful to analytically marginalize over nuisance variables (noise and regularization levels, mixture proportions) as well as to jointly sample transformation parameters and voxel assignments to mixture components (as opposed to alternate between sampling one or the other).

The reversible jump MCMC framework that we implement has strong connection with the jump diffusion process of Grenander and Miller [START_REF] Grenander | Representations of knowledge in complex systems[END_REF] and the birth-and-death kernel framework. It allows to move freely in the space of transformation parameters but also and conccurently, in the space of admissible parametrizations. By circumventing the costly computation of Bayes factors (ratios of evidence for competing parametrizations), it effectively renders MCMC inference tractable for the sparse Bayesian model of registration, even with large dictionaries of basis functions (∼ 10 5 in our experiments).

Full-dimensional moves over the space of transformation parameters were not implemented, as calibrating such transitions calls for the particularly expensive computation of large (non-diagonal) Hessian matrices. This renders them inefficient unless e.g., exploiting dedicated procedures inspired from limited memory quasi-Newton methods [START_REF] Zhang | Quasi-newton methods for markov chain monte carlo[END_REF]. Component-wise transitions are also particularly suitable provided that the set of active bases must be jointly explored.

B. Variational Bayes and MCMC inference

Experimental results point towards a good empirical correspondence between the mean estimates of displacement returned by the approximate VB inference and (asymptotically) exact MCMC inference, particularly in the presence of informative data. Unfortunately they also evidence limitations of the approximate VB scheme for purposes of uncertainty quantification. This defect is offset by a significantly faster running time for the VB scheme (one order of magnitude).

As shown in section III, VB inferences selects a single parametrization by means of marginal likelihood maximization, although this optimal parametrization will often be only marginally so. Discarding all marginally sub-optimal parametrizations results in generally unreliable estimates of uncertainty. This is also evidenced by the lower magnitude of uncertainty predicted by the fixed basis MCMC scheme. Approximate schemes that circumvent this issue can likely be devised, for instance by keeping track of several sets of relevant explanatory variables [START_REF] Schniter | Fast bayesian matching pursuit[END_REF]. The uncertainty on optimal basis locations could be accounted for in the VB scheme, either in an ad-hoc manner by local perturbations of the basis centers when sampling a transformation from the variational posterior, or in a more formal way by regarding basis centers as random variables whose associated variational posteriors are to be estimated.

The RVM * basis augmentation strategy of Rasmussen and Candela [START_REF] Rasmussen | Healing the relevance vector machine through augmentation[END_REF] partially addresses the second issue of uncertainty underestimation in absence of data. It is applicable only when voxelwise estimates of uncertainty are expected, as opposed to full transformation samples. Another strategy would be to relax the form of the variational posterior family so that it can better approximate the true posterior away from data, with the constraint that the computational burden remain suitably low under such a factorization. Alternatively we note instead the high potential for parallelization of the proposed MCMC approach, which could make it more amenable to routine use on real data.

Finally VB makes parametric assumptions about the form of the posterior distribution, and infers hyperparameter distributions whereas the proposed MCMC scheme generally marginalizes over such hyperparameters. This is likely to account for some of the minor differences observed between the VB and Fixed Basis MCMC approaches (sec. IV-B4).

C. Underlying assumptions of the Sparse Bayesian registration model

The validity of model assumptions may affect the quality of uncertainty estimates. Possible biases stem on the one hand from the inexactness of the generative model of images (modelling the intensity residual as a mixture of Gaussians, discarding spatial correlations between residuals), on the other hand from inexactness of the transformation model (the parametrization of the transformation as well as the choice of prior). Their impact was not thoroughly explored so far, but this work provides the methodological framework to do so.

The assumption that source and target image intensities coincide up to spatially varying noise mostly holds in the context of mono-modal registration. For multi-modal registration, a mapping function between source and target image intensities should be used (as in e.g. [START_REF] Leventon | Multi-modal volume registration using joint intensity distributions[END_REF], [START_REF] Guimond | Three-dimensional multimodal brain warping using the demons algorithm and adaptive intensity corrections[END_REF]) and can be regressed within a probabilistic framework [START_REF] Janoos | Bayesian characterization of uncertainty in multi-modal image registration[END_REF].

VI. CONCLUSION

In this article we explored the properties of the proposed sparse Bayesian model of registration for the purpose of uncertainty quantification. We emphasize the distinction between the Bayesian model itself and inference schemes used to estimate posterior distributions under this model. In previous work [START_REF] Le Folgoc | Sparse bayesian registration of medical images for self-tuning of parameters and spatially adaptive parametrization of displacements[END_REF] an efficient but approximate inference scheme was developed, based on Variational Bayesian arguments and the principle of marginal likelihood maximization. In the present work we design a reversible jump Markov chain that characterizes the exact posterior arbitrarily well (provided that enough samples can be drawn) and answer the two following questions. Firstly, does the fast approximate scheme provide faithful estimates of expectation and uncertainty? Secondly, is the sparse Bayesian model of registration useful for the purpose of uncertainty quantification? We evidence limitations of the approximate inference scheme for uncertainty quantification, but show that the true posterior distribution itself is meaningful: orders of magnitude for the true uncertainty (as characterized by MCMC sampling) are quantitatively reasonable, the uncertainty is higher in textureless regions and lower in the direction of strong intensity gradients.

Fig. 1 .

 1 Fig. 1. Graphical model of registration. The generative model of data D involves a transformation Ψ of space, and noise governed by a set of underlying parameters P . Hyperpriors (with hyperparameters H P ) are in turn imposed over the noise parameters. The transformation is parametrized as a linear combination of predefined basis functions {φ k , k = 1 • • • M } with associated weights w k . Priors on the transformation smoothness and on the relevance of individual bases introduce additional parameters (λ and z k respectively). Random variables are circled, hyperparameters are squared. Arrows capture conditional dependencies. Shaded nodes are observed variables or fixed hyperparameters. The transformation Ψ is fully determined by its parent nodes (the φ k and w k ), hence the doubly circled node. The content of plates is replicated (M times).

1 )

 1 Likelihood model: The generative model of data makes explicit the relationship between the data D and the spatial mapping Ψ. It is specified by a likelihood model p(D|Ψ; P ) (often conditioned on a set of hyperparameters P ) that typically assumes the form of a Boltzmann distribution p(D|Ψ; P ) ∝ exp -E D (D, Ψ; P ). For landmark registration, a transformation that approximately maps corresponding key points {t i } and {T i }, i = 1 • • • N , between a template object and a target object is sought. A standard choice of energy is the sum of squared distances between pairings, up to multiplicative factor:

Fig. 2 .

 2 Fig. 2. Graphical representation of the generative data model (using the same graphical symbols as in Fig. 1). Residuals between the fixed image J and the warped image I • Ψ -1 are assumed to be distributed according to a mixture of L Gaussian components whose parameters ρ l (probability of falling in the lth component) and β l (inverse variance a.k.a. precision parameter for the lth Gaussian component) are regarded as latent variables. c n ∈ {1 • • • L} assigns the corresponding voxel to one of the L mixture components.

Fig. 3 .

 3 Fig. 3. Energies corresponding to example GMM (black) and SSD (dotted) models. The voxelwise penalty is shown as a function of the intensity residual. The effective soft threshold on the penalty incurred for large intensity residuals results in increased robustness of the GMM.

  nneighb is an integer fixed in advance. Set w = w, z = z.Draw one of3 competing events: on-off , exchange, update. if φ k inactive and update then Exit. No action to implement (as wk = 0). if exchange and φ k active then Draw an inactive basis φ k * to replace φ k . A proposal that favours well-aligned bases is designed. else if exchange and φ k inactive then Draw an active basis φ k * to replace φ k . A proposal that favours well-aligned bases is designed. if φ k active and on-off or exchange then Set zk = 0 and wk = 0.

Fig. 4 .

 4 Fig. 4. Comparison of approximate evidence-based inference and faithful MCMC inference for the sparse Bayesian model, on a 1D regression task. Data points (black dots) are sampled with additive i.i.d. Gaussian noise from the true signal (dashed green line). The consistence of the fast and faithful estimates of the regressor function (black lines) is satisfactory (w.r.t.uncertainty levels), even more so in the presence of data. Esimates of uncertainty (grey ribbon), however, can be inconsistent.

Fig. 5 .

 5 Fig. 5. Registration setting: (left) fixed image, and (right) moving image, at resolution 1.25mm × 1.25mm.

Fig. 7 .

 7 Fig. 7. Comparison of estimates of the posterior mean returned by MCMC characterization. (Left) Sampling activation variables z and corresponding weights w, alternatively with voxel mixture assignments c. (Right) Joint sampling, as per the approach proposed in section II-C4.

Fig. 8 .

 8 Fig. 8. (Top Row) Comparison of the posterior mean displacement returned by VB (left) vs. MCMC (right), and difference between the two (middle). (Bottom Row) Mean displacement returned by Fixed Basis MCMC (middle) and the difference with the VB (resp. MCMC) estimate (left, resp. right).

Fig. 9 .

 9 Fig. 9. Estimates of uncertainty obtained by characterizing the posterior distribution of the sparse Bayesian model by (Left column) Variational Bayes (Right column) MCMC sampling (Middle column) Fixed Basis MCMC sampling. (Second row) Tensor visualization of the displacement uncertainty: each tensor encodes the square root of the 2 × 2 covariance displacement matrix at this location. Tensor elongation along a direction indicates higher uncertainty along that direction. The color scheme encodes the direction of the first eigenvector. (First row) Trace of the square root covariance.

Fig. 10 .

 10 Fig. 10. Main constituents of the MCMC scheme.

  a set of virtual pairings whose value does not depend on β, λ. H * is a block diagonal matrix whose ith diagonal block H * i is the d × d precision matrix associated to the ith virtual pairing T i * . The factors stemming from the marginalization:

A complete and concise summary of the relevant derivations and schemes is given in http://www.kamperh.com/notes/kamper bayesgmm13.pdf

ACKNOWLEDGMENT Part of this work was funded by the European Research Council through the ERC Advanced Grant MedYMA (2011-291080) on Biophysical Modeling and Analysis of Dynamic Medical Images. The first author was funded by the Microsoft Research -Inria Joint Centre.


Bruno Guillaume
email: bruno.guillaume@inria.fr

Nicolas Lefebvre
email: nicolas.lefebvre@inria.fr

Crowdsourcing Complex Language Resources: Playing to Annotate Dependency Syntax

This article presents the results we obtained on a complex annotation task (that of dependency syntax) using a specifically designed Game with a Purpose, ZombiLingo. 1 We show that with suitable mechanisms (decomposition of the task, training of the players and regular control of the annotation quality during the game), it is possible to obtain annotations whose quality is significantly higher than that obtainable with a parser, provided that enough players participate. The source code of the game and the resulting annotated corpora (for French) are freely available.

Introduction

For better or worse, 2 the overwhelming domination of machine-learning systems in natural language processing (NLP) over the past 25 years and the increasing use of evaluation campaigns and shared tasks have put human-annotated corpora at the heart of the field. Manual annotation is now the place where linguistics hides.

The availability of manually annotated corpora of high quality (or, at least, reliability) is therefore key to the development of the field in any given language. However, the creation of such resources is notoriously costly, especially when complex annotations, e.g. for dependency syntax, are at issue. For example, the cost of the Prague Dependency Treebank was estimated at $600,000 in [START_REF] Böhmová | [END_REF].

Over the years, many solutions have been investigated in the attempt to lower manual annotation costs. One obvious avenue is to use an appropriate annotation tool, as shown for example in [START_REF] Dandapat | [END_REF]. As a complementary aid, NLP systems can be used to reduce the annotation burdeneither beforehand, for example with tag dictionaries [START_REF][END_REF] and more generally with preannotation [START_REF] Skjaerholt | Influence of preprocessing on dependency syntax annotation: speed and agreement[END_REF], or more iteratively during the annotation process, for instance through active learning [START_REF] Baldridge | Active learning and the total cost of annotation[END_REF]. These solutions have proved efficient and have indeed helped reduce the annotation cost, but the creation of a large annotated corpus in the traditional manner remains very expensive.

Another way to address the issue is simply to limit the amount paid to the human annotators. This is the case with microworking crowdsourcing, especially through the use of platforms like Amazon Mechanical Turk, via which the workers are (micro)paid to perform simplified tasks [START_REF] Snow | [END_REF]. Apart from the ethical issues raised by these platforms (detailed in (Fort et al., 2011)), microworking platforms do not support true training of the workers: tests can be set up to select them, but they have to perform hidden work to train themselves (as shown in [START_REF] Gupta | [END_REF]). Consequently, tasks must be simplified to be manageable by workers: for example, a task related to recognition of textual entailment across several sentences (which might actually reflect entailment, neutrality, or contradiction) could be simplified by presenting only one pair of sentences and a binary response to a single question, such as "Would most people say that if the first sentence is true, then the second sentence must be true?" [START_REF] Bowman | [END_REF]. To paraphrase [START_REF] Dandapat | [END_REF], it seems that there is no easy escape from the high cost of complex linguistic annotation.

We present here an on-line game that enables the production of quality annotations of complex phenomena (here, dependency syntax), tested on French. The produced corpus is constantly growing and is designed to be (i) completely free and available, and (ii) of sufficient quality. We first examine the existing treebanks for French and their limitations and detail some previous experiments with Games with a Purpose (GWAPs). Then we present the game we developed and the evaluation we performed on the produced annotations. We finally conclude by discussing the potentials and limits of GWAPs for the creation of complex language resources.

Previous Work

Treebanking for French

As surprising as it may seem, French has long been relatively low-resourced, primarily due to legal issues: lexicons and annotated corpora existed but could not be used or redistributed freely. This is in particular the case for the French Treebank (FTB or corpus arboré de Paris 7) [START_REF] Abeillé | [END_REF], which is available for research purposes only and cannot be freely redistributed. 3 Several versions are reported in the literature, with the size varying from 12,351 sentences and 350,947 words (for the FTB-UC) to 18,535 sentences and 557,149 words4 (for the FTB-SPRML).

To try and circumvent this restriction, Candito and Seddah created the Sequoia corpus [START_REF] Candito | Le corpus Sequoia : annotation syntaxique et exploitation pour l'adaptation d'analyseur par pont lexical[END_REF], which is freely available5 under a LGPL-LR license, but is limited to 67,038 tokens. The same authors developed an additional question bank with 23,236 tokens [START_REF] Seddah | Hard time parsing questions: Building a questionbank for french[END_REF]. Both corpora use the same annotation guide and set of relations as the FTB.

A Universal Dependency corpus (McDonald et al., 2013) was created for French and is freely available under a CC BY-NC-SA license. In version 1.3, released in May 2016, it contains 401,960 tokens, but it "has not been manually corrected systematically". 6 Moreover, the annotation format for Universal Dependencies suffers from certain drawbacks, for example, it does not distinguish between arguments and modifiers for nominal complements of verbs.

Other treebanks exist for French, but they either concern spoken language, as with Rhapsodie (Lacheret et al., 2014) or the oral Treebank (Abeillé and Crabbé, 2013), or specific language types, like the Social Media Treebank [START_REF] Seddah | [END_REF].

This situation (in which references exist, but a large, fully available, manually annotated corpus is lacking) makes dependency syntax for French an appropriate candidate for testing a new paradigm for complex linguistic resource development: the use of on-line Games with a Purpose.

Playing to Create Language Resources

Games with a Purpose are games in which participants, knowingly or not, create data by playing. They are not serious games as such, as their main purpose is not to train people, but to produce data (such as annotations, lexicon entries, image labels, etc). GWAPs for NLP are broadly surveyed in [START_REF] Lafourcade | [END_REF], and a detailed analysis of their performance in comparison with other means of language resource production is provided in [START_REF] Chamberlain | [END_REF]. Thus, we will focus here on the complexity of the resource to be produced, rather than on giving an exhaustive list of games.

Most GWAPs rely on the players' knowledge of the world and innate capabilities to enable creation of data, in our case language resources. This reliance was seen for example in the very well-known ESP Game [START_REF] Von Ahn | Labeling images with a computer game[END_REF], in which participants played by labeling images. Another early GWAP is JeuxDeMots7 [START_REF] Mathieu Lafourcade | Making people play for lexical acquisition[END_REF], in which players created a lexical network of more than 47 million relations between more than 900,000 entries (terms and named entities), simply by entering terms and named entities associated in a more or less specified way with another entry. The video games presented in [START_REF] Jurgens | It's All Fun and Games until Someone Annotates: Video Games with a Purpose for Linguistic Annotation[END_REF]) also seem 8 to target the users' intuition and knowledge of the world to perform word sense disambiguation and a mapping from WordNet senses to images. A more limited gamified interface is used in the Language Quiz 9 , which asks the participants to perform sentence or tweet sentiment analysis, without any training. 10 Other GWAPs benefit from the players' school knowledge. One such game is Phrase Detectives 11 [START_REF] Poesio | [END_REF], in which participants are asked to identify the antecedent of a noun phrase in a text. The task is more complex than labeling images or associating ideas, so the game includes a mandatory short training phase, allowing the players to become familiar with it. Almost 8,000 players participated in the game, which enabled the annotation of anaphora relations in a 162,000 word corpus. The agreement between players and experts was evaluated overall at 84%.

The wordrobe website 12 presents a family of eight gamified tasks 13 whose results help to improve the Groningen Meaning Bank [START_REF] Venhuizen | [END_REF]. The proposed tasks all relate to semantic disambiguation (noun vs verb, co-reference identification, named entity annotation, etc) and while some are relatively easy -like Play Twins (noun vs verb) or Play Names (named entity annotation) -most require some more advanced (at least school-level) knowledge. However, the interface does not provide a training phase and the only help available is a short guide to the task.

In a very different domain (biochemistry), the creators of FoldIt 14 demonstrated that people with no specific knowledge of the subject could be trained to perform complex protein folding tasks with impressive results: players helped find the solution to a long-standing problem concerning protein crystal structure (Khatib et al., 2011). To achieve these results, players are trained within "introductory levels", solving puzzles of varying complexity "[...] introducing the player to new problem related concepts as though they are the rules of a game." (Cooper et al., 2010).

This approach seemed particularly suited to our task (dependency annotation with nearly twenty relations to annotate), so we decided to build upon it, asking the player to annotate one relation at a time, with a specific training process for each. We report here the results we obtained, which show that a GWAP can be used for a complex linguistic annotation task that is not directly linked to world knowledge or to human intuition. To our knowledge, no other such experiment has been carried out to date. The model of dependency syntax that we use is well-known in linguistics and NLP but not in the French educational system, so we can assume that people without a linguistic or NLP background have no knowledge of it.

It should be noted that an un-gamified crowdsourcing-based method has been proposed in [START_REF] Hana | Getting more data -schoolkids as annotators[END_REF][START_REF] Hladká | [END_REF] for the dependency syntax annotation of Czech. (In the Czech Republic, the presentation of syntax in school is very close to dependency syntax.) However, the authors report in [START_REF] Hana | Getting more data -schoolkids as annotators[END_REF] that the accuracy of the annotations they obtained is significantly lower than that of their parser. Moreover, they evaluate the results on a small set of 100 sentences selected or created by the authors. In their more recent publication, they report only the tree editing distance between the produced annotation and the reference annotation; hence, we are unable to compare their results to those presented here.

Designing a Game for a Complex Task

We describe below the mechanisms used in the game to take into account the complexity of the task and the game's avoidance of direct reliance on intuition. We have chosen to rely on the Sequoia annotation guidelines (which largely follow the FTB guidelines). This is a natural choice if we want to use available resources for training players and for evaluation. Moreover, existing parsers for French have been 8 The games do not seem to be running as of mid-July 2016, so we could not test them. 9 See: http://quiz.ucomp.eu 10 Two of the authors tested it in September 2016 and were at that time the only participants. 11 See: https://anawiki.essex.ac.uk/phrasedetectives/. 12 See: http://wordrobe.housing.rug.nl. 13 The game features are reduced to a leader board and a betting option, we are therefore reluctant to call it a game, although there is obviously a continuum here.

14 See: http://fold.it.

developed with the same schema.

Decomposing the Complexity of the Task

Since the analysis of a whole sentence is complex, 15 we decided to decompose it: the full annotation of a sentence is split into atomic tasks, where each task is linked to one type of dependency relation. Thus, the player is trained on one type of relation at a time, instead of twenty, and must focus on only one type of dependency relation across several sentences. S/he therefore has fewer elements of information to remember at a time. Moreover, each relation was awarded a level which reflects its difficulty (from the point of view of the game developers) and corresponds to a level in the game. The player can therefore choose amongst more and more relations as s/he progresses in the game. We decided not to include certain relations in the current version of the game. This concerns punct, which is not consistently annotated in reference, reldep, which is an underspecified relation used in different types of contexts, and relations which are not processed by Talismane. Note that the results given in Section 4 is based on all the relations except punct.

For each relation, explanations and examples are provided with the different contexts where it may appear. 16In addition to explanations, it is also important to give players examples of real utterances taken from a corpus. Reference examples are taken from the Sequoia corpus and are used in the two game mechanisms, TRAINING and CONTROL, described below.

We used only a part of the corpus in order to keep sentences aside for the evaluation of the produced annotations (see Section 4). In

Playing the Game

The organization of ZombiLingo is illustrated in Figure 1, which presents how the TRAINING and play phases are articulated with both the sub-corpora described in Table 1 and the raw corpora to be annotated (usually extracted from Wikipedia). The reference corpus is used in three different phases, TRAINING, CONTROL and EVAL, as explained below.

TRAINING phase

Following the decomposition of the task, for each dependency relation, a specific TRAINING phase is required before entering the game. During this phase, sentences from the REF T rain&Control corpus are presented to the player and feedback is given in case of error.

Figure 2 illustrates the feedback given to the players in the TRAINING phase: the player was asked to find the second conjunct of a coordination et (and) and s/he wrongly answered Europe. A skull and crossbones flashes and a message revealing the right answer parvient (reaches) are displayed.

An advantage of offering a separate TRAINING for each relation is that the player does not have to wait long before starting the game. Once s/he is connected to the game, only a few minutes are required before starting actual play and production of annotations.

Play phase

In the general mode, the player chooses a relation from those available and a sequence of ten questions is proposed. The questions vary as follows: • For relations where the dependent element tends to be unique for a given governor (e.g. a verb has one subject, a noun has one determiner, and so on), the governor is given and the players must find the dependent;

• For relations involving several dependent elements (e.g. a verb has several prepositional phrases modifying it, a noun may be modified by several adjectives, and so on), the dependent element is given and the players must find the governor.

The player's answer is registered in the database and points are awarded if the same answer was previously given by a parser or by another player. 17 The number of points (in the game, they are represented as brains) that can be won on a certain relation depends on a global measure of the sentence complexity (measured as the maximum number of nested dependency relations) and on the level of the relation in the game.

CONTROL mechanism

TRAINING is not sufficient to ensure that players annotate data consistently: if they did the TRAINING on a certain relation a long time ago or if they play a lot, their ability to annotate this relation may decrease. To deal with this problem, we added a CONTROL mechanism during the game. The player is sometimes asked to annotate a relation in a sentence taken from the REF T rain&Control corpus. If the player fails to find the right answer, feedback including the solution is displayed (as in Figure 2). After a given number of failures on the same relation, the player has to redo the corresponding TRAINING. Using this CONTROL mechanism, we can also estimate our degree of confidence in a specific player on a specific relation (see Section 4) and take this into account to weigh his or her answer.

Behind the Curtain

Preprocessing Data

In the first version of the game, when new sentences were added, a parser was used to pre-annotate the sentences. Then, the players were asked to confirm or correct the parser's predictions. The main drawback of this arrangement is that the player would have a large number of very easy annotations to decide on and would usually agree with the parser. Modifying this arrangement, we now use two parsers to pre-annotate the corpora and ask the participants to play items for which the parsers give different annotations. The two parsers used are Talismane [START_REF] Urieli | Robust French syntax analysis: reconciling statistical methods and linguistic knowledge in the Talismane toolkit[END_REF] and FRDEP-PARSE (Guillaume and [START_REF] Guillaume | Dependency Parsing with Graph Rewriting[END_REF].

Talismane is a statistical parser trained on the training part of the FTB. It was evaluated on the dev and test part of the FTB with LAS scores ranging from 86.8% to 88.5%. As for FRDEP-PARSE, it is a parser which combines statistical methods for POS-tagging (using MElt [START_REF] Denis | Coupling an annotated corpus and a lexicon for state-of-the-art pos tagging[END_REF]) and symbolic methods for dependency parsing (with graph rewriting). FRDEP-PARSE was evaluated on Sequoia with a LAS score of 76.04% and a precision of 85.96%. (The system returns partial dependency syntax structures.)

It is important that the two parsers are based on different paradigms (statistical and symbolic): we can hope that the two tools are complementary and will produce different types of errors.

The next section will provide the two parsers' detailed results for the REF Eval corpus.

Exporting Data

The players' annotations are stored in a database and each annotation receives a score. When a player is asked a question, we consider the set of possible answers in the database and adjust the score as follows:

• If the player's answer belongs to the set, the score of the answer is increased and the scores of its competing annotations (the rest of the set) are decreased.

• If the player gives an answer not in the set, a new annotation is created in the database with a default score and the scores of the answers in the set are decreased.

The positive or negative score adjustments are weighted by the level of the player, we thus award higher confidence to heavy players (who have usually reached higher levels) than to beginners. When a corpus is exported, for each token (lexical unit), we consider all the annotations in the database for which it is a dependent element and select the one with the highest score. Thus, each token receives exactly one governor with one relation and we can ensure that the exported corpus contains well-formed dependency trees.

4 Quantitative and Qualitative Evaluations

Participation and Production

If a crowdsourcing approach to annotation is to succeed, enough participants must be induced to participate and create data. As of 2016, July 10, there were 647 players registered for our game. They have collectively produced 107,719 annotations.

As a reminder, Table 2 shows a quantitative comparison of the existing corpora for French. What this table does not show is that the corpus produced through the game is still growing, and will grow as long as we support and advertise the application. On the ZombiLingo corpus the average number of players' annotations by tokens (called the density) is 0.84 (107,719 / 128,046). Figure 3 shows the density of annotations per relation. The coverage is clearly not homogeneous: the density on the relation aff (affix) is greater than 6, whereas for some relations it remains below 1. This means that we need to attract at least some players to play a wider variety of relations types or to direct them towards some under-played relations, such as mod (modifier) or coord (coordination).

Qualitative Evaluation Methodology

In the first version of the game, we used the CONTROL mechanism to evaluate the quality of annotations produced by participants. We observed that 85.4% of CONTROL items were correctly annotated.

While this evaluation showed that it is possible to train players for a difficult task, it remains a partial evaluation, with severe drawbacks. It measures the performance of the participants on reference sentences but it gives no information about the annotations produced on pre-annotated sentences: if such a sentence contains a dependency relation not predicted by the parser, the player will never play this relation in this sentence. We call this omission silence in the produced data. On the other hand, if a relation is predicted by the parser but cannot be found in the sentence, the player will be asked to give an opinion on a question for which there is no sensible answer: we call such nonsensical situation noise in the data.

During the game, the user is supposed to click on a specific image (crossed bones) if a question has no answer but, although we train the participants on this alternative, we know that it is underused. For instance, when a player is tasked with finding the object of a verb, s/he tends to search for a word which looks like an object without realizing that this verb may actually have no object at all.

In order to obtain a precise evaluation of the produced annotations, we entered a part of our reference corpus (named REF Eval ) into the game (see Figure 1), as if it were a raw corpus. We used it in the following way:

• the raw text corresponding to the REF Eval corpus is put into the general pipeline (i.e., it is parsed with the two parsers, the differences in the annotation being proposed to the players; this is the EVAL mechanism in the Figure 1);

• we report the scores obtained by each parser (i.e., recall/precision/F-measure) and the score of the corpus EXP Eval , as exported by the game.

Results

The REF Eval corpus contains 12,660 relations which can be played. On this subset of relations, the two parsers give the same output in 10,134 of the cases. Their analysis is correct in 96.84% of the cases (i.e. in 9,814 cases out of 10,134) and it is wrong in only 320 cases, which represents 3.16% of the output (i.e. 320 cases out of 10,134). However, even in cases where the two parsers provide two different wrong answers, the players will be asked to give their opinion on these cases and, hopefully, will produce a correct final annotation. In the end, the number of wrong annotations that will never be proposed to the players represents only 2.53% (i.e. 320 cases out of 12,660) of the relations which can be played.

In Table 3, we report the scores in three settings: the Talismane parser, the FRDEP-PARSE parser, and data exported by the game (EXP Eval ). In all of these setting, we may have partial dependency structures, so we report recall, precision and F-measure. We consider all relations (except punctuation which is not consistently annotated in the Sequoia corpus). Note that the score of the annotations produced by the game is lower than that of the second parser. This difference comes about because some relations were not played by enough players, as already observed in Subsection 4.1. To explore in more detail how the players influence the corpus score, we must look at these values relation by relation.

First, 3,575 game items correspond to the REF Eval corpus and only 2,644 game actions were performed by the player on these items (so the average density is 0.74). Moreover, as already observed in Figure 3, the distribution of the game actions is not at all homogeneous. If we discard relations for which the number of occurrences in the REF Eval corpus is too low18 to obtain significant results, the average density of game actions on game items ranges from 0.04 on mod (modifier) to 4.65 on aux.tps (tense auxiliary).

In Figure 4, we report the density of annotations on the corpus REF Eval and Figure 5 gives, for the same relations in the same order, the values of the F-measure in the three settings (with Talismane in light gray, FRDEP-PARSE in a darker gray and finally the game export in dark gray).

For the relations where the density is higher than 1 (seen in the left hand part of the figures), we observe that the F-measure computed on the corpus from the game is always higher than that of the two parsers. From these experiments, we can conclude that we manage to obtain annotations with a quality significantly higher than that obtainable with a parser, provided that enough players participate. It is worth noting that amongst the relations which are densely played, some are considered as complex, such as dep.coord (the player has to find the head of the second conjunct of a coordination). Figure 2 shows an example of such a difficult case: the player has to read and to understand the whole sentence to give the right answer (it is a long distance dependency).

The next challenge is to induce players to annotate a wider variety of relations so as to increase the quality of the whole corpus. Relations with a very low density are either relations with a high number of occurrences (the REF Eval corpus contains 1,874 mod relations) or relations which are less intuitive for the players. We will take this last factor into account and improve the documentation on these relations.

Conclusion

We have presented ZombiLingo, a game designed for a complex linguistic annotation task, namely dependency syntax. A first prototype of the game was released in July 2014 and an engineer started working on a production version in October 2015, completing the first version by the end of that year. As of July 2016, the game had enabled the production of more than 100,000 annotations for French, with a precision of 0.93. These results are very promising, especially for low-resourced languages. We still need to fill some annotation gaps (e.g. for relations insufficiently played) and have developed for this purpose a new duel mode, which will allow senior players to compete with each other in the annotation of whole sentences.

The most difficult factor when using GWAPs is to find ways of attracting and keeping participants [START_REF] Poesio | [END_REF], which requires a continuing communication effort. We have experienced "waves" of players following specific events we participated in or challenges we organized. We also attracted new players when advertising the game on social networks, but this effort must be regularly maintained and renewed.

The game source code is freely available on GitHub 19 under an CeCILL open-source license 20 . The code is designed to be easily adaptable to any human language: all messages are isolated from the code, which is Unicode compliant. We plan to adapt it to English and to a less-resourced language in the coming months.

The resource created for French is directly and freely available under a CC BY-NC-SA license from the game website 21 and is updated every night.

Figure 1 :

 1 Figure 1: Organization of the different mechanisms and corpora

Figure 2 :

 2 Figure 2: The main interface of the game during the TRAINING phase.

Figure 3 :

 3 Figure 3: Density of answers for each relation.

Figure 4 :

 4 Figure 4: Density of answers for each relation on the REF Eval corpus (a whitespace separates relations where density is greater than 1 from the others).

  Figure 5: F-measures for the two parsers and the game, on each relation. (As in Figure4, relations are split into two groups, with density greater than 1 on the left and density lower than 1 on the right.)

Table 1 ,

 1 we show how the reference corpus is split for different uses.

	REF T rain&Control	REF Eval	Unused
	50%	25%	25%
	1,549 sentences	776 sentences	774 sentences

Table 1 :

 1 Uses of the 3,099 sentences of the Sequoia reference corpus.

Table 2 :

 2 ZombiLingo corpus size, as compared to other existing French corpora annotated with dependency syntax.

		Sequoia 7.0 UD-French 1.3 FTB-UC FTB-SPMRL ZombiLingo
	Sentences	3,099	16,448	12,351	18,535	5,221
	Tokens	67,038	401,960	350,947	557,149	128,046
	Tokens/sentence	21.6	24.4	28.4	30.1	24.5

Table 3 :

 3 Evaluation of the data produced by the game on the REF Eval corpus (except punctuation).

		Talismane FRDEP-PARSE Game
	LAS / Recall	0.745	0.743	0.674
	Precision	0.759	0.852	0.932
	F-measure	0.752	0.794	0.782

See http://www.llf.cnrs.fr/fr/Gens/Abeille/French-Treebank-fr.php.

These numbers correspond to the version found on: http://gforge.inria.fr/projects/fdtb-v1.

See here: http://deep-sequoia.inria.fr/.

See http://universaldependencies.org/fr/overview/introduction.html.

The game is still running online: http://www.jeuxdemots.org.

Insight concerning the complexity of syntactic annotation for the Penn Treebank is provided in[START_REF] Marcus | [END_REF], where the learning curve for syntax was estimated to be twice that for part-of-speech (two months vs one).

This is much like in an annotation guide but with simpler vocabulary and fewer details.

Unfortunately, if a player is the first to give a right answer, which is not given by the pre-annotation, s/he receives no points. A mechanism to give points off-line in this situation is planned but not yet implemented.

We have discarded relations with less than 25 occurrences in REF Eval , namely aux.caus (causative auxiliary, 4 occurrences), arg (specific relation linking two parts of a range, 2 occurrences) and dis (dislocation, 1 occurrence).

Acknowledgements

The development of ZombiLingo is funded by Inria, through an ADT grant. The French Ministry of culture has also helped us to communicate and publish results concerning ZombiLingo through two grants.

Above all, we want to thank all of the ZombiLingo players, without whom the resource would simply not exist! Special thanks to players JYA, Chouchou and Xohwohxo, who gave us valuable feedback on the game. We also wish to thank Djamé Seddah, for his help in building a solid state of the art and Guy Perrier, who, as Professor Frankenperrier, helped us to write the guidelines and to check errors in the reference. Finally, we want to thank the reviewers, who gave us interesting suggestions for improving this article.