{'fraction_non_alphanumeric': 0.042139691918126186, 'fraction_numerical': 0.009432369698248576, 'mean_word_length': 4.095258574346844, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stochastic Algorithm For Parameter Estimation For Dense Deformable Template Mixture Model

Estimating probabilistic deformable template models is a new approach in the fields of computer vision and probabilistic atlases in computational anatomy. A first coherent statistical framework modelling the variability as a hidden random variable has been given by Allassonni\`ere, Amit and Trouv\'e in [1] in simple and mixture of deformable template models. A consistent stochastic algorithm has been introduced in [2] to face the problem encountered in [1] for the convergence of the estimation algorithm for the one component model in the presence of noise. We propose here to go on in this direction of using some"SAEM-like"algorithm to approximate the MAP estimator in the general Bayesian setting of mixture of deformable template model. We also prove the convergence of this algorithm toward a critical point of the penalised likelihood of the observations and illustrate this with handwritten digit images.

1. Introduction. The issue of representing and analysing some geometrical structures upon which some deformations can act is a challenging question in applied mathematics as well as in Computational Anatomy. One central point is the modelisation of varying objects, and the quantification of this variability with respect to one or several reference models which will be called templates. This is known as "Deformable Templates" [12]. To our best knowledge, the problem of constructing probabilistic models of variable shapes in order to statistically quantify this variability has not been successfully addressed yet in spite of its importance. For example, modelling the anatomical variability of organs around an ideal shape is of a crucial interest in the medical domain in order to find some characteristic differences between populations (pathological and control), or to exhibit some pathological kind of deformations or shapes of an organ.
Many solutions have been proposed to face the problem of the template definition. They go from some generalised Procruste's means with a variational [11] or statistical [10] point of view to some statistical models like Active Appearance Model [6] or Minimum Description Length methods [16]. Unfortunately, all these methods only focus on the template whereas the geometrical variability is computed afterwards (using PCA). This contradicts with the fact that a metric is required to compute the template through the computation of deformations. Moreover, they do not really differ from the variational point of view since they consider the deformations as some nuisance parameters which have to be estimated and not as some unobserved random variables.
The main goal of this paper is to propose a coherent estimation of both photometric model and geometrical distribution in a given population. Another issue addressed here is the clustering problem. Given a set of images, the statistical estimation of the component weights and of the image labels is usually supervised, at least the number of components is fixed. The templates of each component and the label are estimated iteratively (for example in methods like K-means) but the geometry, and related to this the metric used to compute the distances between elements, is still fixed. Moreover, the label, which is not observed, is, as the deformations, considered as a parameter and not as a hidden random variable. These methods do not lead to a statistical coherent framework for the understanding of deformable template estimation and all these iterative algorithms derived from those approaches do not have a statistical interpretation as the parameter optimisation of a generative model describing the data.
In this paper we consider the statistical framework for dense deformable templates developed by Allassonniere, Amit and Trouve in [1] in the generalised case of mixture model for multicomponent estimation. Each image taken from a database is supposed to be generated from a noisy random deformation of a template image picked randomly among a given set of possible templates. 1 All the templates are assumed to be drawn from a common prior distribution on the template image space. To propose a generative model, each deformation and each image label have to be considered as hidden variables. The templates, the parameters of the deformation laws and the components weights are the parameters of interest. This generative model allows to automatically decompose the database into components and, at the same time, estimates the parameters corresponding to each component while increasing the likelihood of the observations. Given this parametric statistical Bayesian model, the parameter estimation is performed in [1] by a Maximum A Posteriori (MAP). The authors carry out this estimation problem using a deterministic and iterative scheme based on the EM (Expectation Maximisation) algorithm where the posterior distribution is approximated by a Dirac measure on its mode. Unfortunately, this gives an algorithm whose convergence toward the MAP estimator cannot be proved. Moreover, as shown by the experiments in that paper, the convergence is lost within a noisy setting.
Our goal in this paper is to propose some stochastic iterative method to reach the MAP estimator for which we will be able to get a convergence result as already done for the one component case in [2]. We propose to use a stochastic version of the EM algorithm to reach the maximum of the posterior distribution.We use the Stochastic Approximation EM (SAEM) algorithm introduced by Delyon et al in [7] coupled with a Monte Carlo Markov Chain (MCMC) method. This coupling algorithm has been introduced by Kuhn and Lavielle in [14] in the case where the missing variables had a compact support. Contrary to the one component model where we can couple the iteration of the SAEM algorithm with the Markov chain evolution (cf. [2]), we show here that it cannot be driven numerically. We need to consider an alternative method. We propose to simulate the hidden variables using some auxiliary Markov chains, one per component, to approach the posterior distribution. We prove the convergence of our algorithm for a non compact setting by adapting Delyon's theorem about general stochastic approximations and introducing truncation on random boundaries as in [5].
The paper is organised as follows: in Section 2 we first recall the observation mixture model proposed by Allassonniere, Amit and Trouve in [1]. In Section 3, we describe the stochastic algorithm used in our particular setting. Section 4 is devoted to the convergence theorem. Experiments on 2D real data sets are presented in Section 5. The proofs of the convergence of the algorithm are postponed in Section 6 whereas Conclusion and Discussion are given in Section 7.
2. The observation model. We consider the framework of multicomponent model introduced in [1]. Given a sample of gray level images (y i ) 1≤i≤n observed on a grid of pixels {v u ∈ D ⊂ R 2 , u ∈ Λ} where D is a continuous domain and Λ the pixel network, we are looking for some template images which explain the population. Each of these images is a real function I 0 : R 2 → R defined on the whole plane. An observation y is supposed to be a discretisation on Λ of a deformation of one of the templates plus an independent additive noise. This leads to assume the existence of an unobserved deformation field z : R 2 → R 2 such that for u ∈ Λ : where ǫ denotes an additive noise.
2.1. Models for templates and deformations. We use the same framework as chosen in [1] to describe both the templates I 0 and the deformation fields z. Our model takes into account two complementary sides: photometric -indexed by pcorresponding to the templates and the noise variances, and geometric -indexed by gcorresponding to the deformations. The templates I 0 and the deformations z are assumed to belong to some finite dimensional subspaces of two reproducing kernels Hilbert spaces V p and V g (determined by their respective kernels K p and K g ). We choose a representation of both of them by finite linear combinations of the kernels centred at some fixed landmark points in the domain D: (v p,j ) 1≤j≤kp respectively (v g,j ) 1≤j≤kg . They are therefore parametrised by the coefficients α ∈ R kp and β ∈ (R kg ) 2 which yield: ∀v ∈ D, 2.2. Parametrical model. In this paper, we consider a mixture of the deformable template models which allows a fixed number τ m of components in each training set. This means that the data will be separated in τ m (at most) different components by the algorithm.
Therefore, for each observation y i , we consider the pair (β i , τ i ) of unobserved variables which correspond respectively to the deformation field and to the label of image i. We denote below by y t (y t 1 , . . . , y t n ), by β t (β t 1 , . . . , β t n ) and by τ t (τ 1 , . . . , τ n ). The generative model is: is the action of the deformation on the template I α , for u in Λ and δ t is the Dirac function on t. The parameters of interest are the vectors (α t ) 1≤t≤τm coding the templates, the variances (σ 2 t ) 1≤t≤τm of the additive noises, the covariance matrices (Γ g,t ) 1≤t≤τm of the deformation fields and the component weights (ρ t ) 1≤t≤τm . We denote by (θ t , ρ t ) 1≤t≤τm the parameters so that θ t corresponds to the parameters composed of the photometric part (α t , σ 2 t ) and the geometric part Γ g,t for component t. We assume that for all 1 ≤ t ≤ τ m , the parameter θ t = (α t , σ 2 t , Γ g,t ) belongs to the open space Θ defined as Θ = { (α, σ 2 , Γ g ) | α ∈ R kp , |α| < R, σ > 0, Γ g ∈ Sym + 2kg , * (R) } , where R is an arbitrary positive constant and Sym + 2kg , * (R) is the set of strictly positive symmetric matrices. Concerning the weights (ρ t ) 1≤t≤τm , we assume that they belong to the set ̺ = (ρ t ) 1≤t≤τm ∈]0, 1[ τm | τm t=1 ρ t = 1 .

Remark 1.
This yields a generative model: given the parameters of the model, to get a realisation of an image, we first draw a label τ with respect to the probability law τm t=1 ρ t δ t . Then, we simulate a deformation field β using the covariance matrix corresponding to component τ according to N (0, Γ g,τ ). We apply it to the template of the τ th component. Last, we add an independent Gaussian noise of variance σ 2 τ .
We choose a normal distribution for the unobserved deformation variable because of the background we have in image analysis. Indeed, the registration problem is an issue that has been studied deeply for the past two decades. The goal is, given two images, to find the best deformation that will match one image close to the other. Such methods require to choose the kind of deformations that are allowed (smooth, diffeomorphic, etc). These conditions are equivalent, for some of these methods, to choose a covariance matrix that enables to define an inner product between two deformations coded by a vector β (cf. [18,4]). The regularisation term of the matching energy in the small deformation framework treated in this paper can be written as : β t Γ −1 g β. This looks like the logarithm of the density of a Gaussian distribution on β with 0 mean and a covariance matrix Γ g . The link between these two points of view has been given in [1]; the mode of the posterior distribution equals the solution of a general matching problem. This is why we therefore set on the deformation vector β such a distribution. Moreover, many experiments have been run using a large variety of such a matrix which gives us now a good initial guess for our parameter. This leads us to consider a Bayesian approach with a weakly informative prior.
2.3. The Bayesian approach. The information given by the image analysis background is here introduced mathematically in terms of prior laws on the parameters of model (2.1). As already mentioned in the previous paragraph, this background knowledge enables to determine a good initial guess for the laws and the values of the hyper-parameters. As well as for the covariance matrix Γ g , the same arguments are true for the noise variance σ 2 . In the registration viewpoint, this variance is the tradeoff between the deformation cost and the data attachment term that compose the energy to minimise. An empirical good initial guess is therefore known as well.
On another hand, the high dimensionality of the parameters can lead to degenerated maximum likelihood estimator when the training sample is small. While introducing prior distributions, the estimation with small samples is still possible. The importance of these prior distributions in the estimation problem has been shown in [1]. The solution of the estimation equation can be interpreted as barycenters between the hyper-parameters of the priors and the empirical values. This ensures easy computations and other theoretical properties as for example, the invertibility of the covariance matrix Γ g . The role of the other hyper-parameters are discussed in the experiments.
We use a generative model which includes natural standard conjugate prior distributions with fixed hyper-parameters. These distributions are an inverse-Wishart priors on each Γ g,t and σ 2 t and a normal prior on each α t , for all 1 ≤ t ≤ τ m . All priors are assumed independent. Then, is the scalar product of two matrices A and B and tr stands for the trace.
For the prior law ν ρ , we choose the Dirichlet distribution, D(a ρ ), with density , with fixed parameter a ρ .
3. Parameter estimation using a stochastic version of the EM algorithm. For the sake of simplicity, let us denote by N 2nk g and by T {1, . . . , τ m } n so that the missing deformation variables take their values in R N and the missing labels in T . We also introduce the following notations: η = (θ, ρ) with θ = (θ t ) 1≤t≤τm and ρ = (ρ t ) 1≤t≤τm .
In our Bayesian framework, we choose the MAP estimator to estimate the parameters: where q B (η|y) denotes the distribution of η conditionally to y.
Remark 2. Even if we are working in a Bayesian framework, we do not want to estimate the distributions of our parameters. Knowing the distribution of the template image and its possible deformations is not of great interest from an image analysis point of view. Indeed, people are more interested, in particular in the medical imaging community, in an atlas which characterises the populations of shapes that they consider rather than its distribution. Moreover, the distribution of the deformation law makes even less sense. This is the reason why we focus on the MAP.
In practice, to reach this estimator, we maximise this posterior distribution using a Stochastic Approximation EM (SAEM) algorithm coupled with a Monte Carlo Markov Chain (MCMC) method. Indeed, due to the intractable computation of the E step of the EM algorithm introduced by [8] encountered in this complex non linear setting, we follow a stochastic variation called SAEM proposed in [7]. However, again due to the expression of our model, the simulation required in this algorithm cannot be performed directly. Therefore, we propose to use some MCMC methods to reach this simulation as proposed by Kuhn and Lavielle in [14] and done for the one component model in [2]. Unfortunately, the direct generalisation of the algorithm presented in [2] paper turns out to be of no use in practice because of some trapping state problems (cf. Subsection 3.2.). This suggests to go back to some other extension of the SAEM procedure.
3.1. The SAEM algorithm using MCMC methods. Let us first recall the SAEM algorithm. It generates a sequence of estimated parameters (η k ) k which converges towards a critical point of η → log q(y, η) under some mild assumptions (cf. [7]). These critical points coincide with the critical points of η → log q B (η|y). The k th iteration consists in three steps: Simulation step : the missing data, here the deformation parameters and the labels, (β, τ ), are drawn with respect to the distribution of (β, τ ) conditionally to y denoted by π η , using the current parameter η k−1

2)
Stochastic approximation step : given (∆ k ) k a decreasing sequence of positive step-sizes, a stochastic approximation is done on the quantity log q(y, β, τ , η), using the simulated value of the missing data :

3)
Maximisation step : the parameters are updated in the M-step, Initial values Q 0 and η 0 are arbitrarily chosen. We notice that the density function of the model proposed in paragraphs 2.2 and 2.3 belongs to the curved exponential family. That is to say that the complete likelihood can be written as: , where the sufficient statistic S is a Borel function on R N PT taking its values in an open subset S of R m and ψ, φ two Borel functions on ΘP̺. (Note that S, φ and ψ may depend also on y, but since y will stay fixed in the sequel, we omit this dependency). Thanks to this property of our model, it is equivalent to do the stochastic approximation on the complete log-likelihood as well as on the sufficient statistics. This yields equation (3.3) to be replaced by the following stochastic approximation s of the sufficient statistics S : We now introduce the following function: L : SPΘP̺ → R as L(s; η) = −ψ(η) + s, φ(η) . It has been proved in [1] that there exists a critical functionη : S → ΘP̺ which is a zero of ∇L. It is straightforward to prove that this function satisfies: ∀η ∈ ΘP̺, ∀s ∈ S, L(s;η(s)) ≥ L(s; η) so that the maximisation step (3.4) becomes: Concerning the simulation step, in our model, the simulation of the missing variables with respect to the conditional distribution π η cannot be carried out. Indeed, its probability density function (pdf) has a close form but rather complicated ; it does not correspond to some usual pdf. One solution proposed in [14] for such cases is to couple the SAEM algorithm with Monte Carlo Markov Chain (MCMC) method. However, we do not fit exactly into their requirements since the missing variable β does not have a compact support. We introduce an ergodic Markov chain whose stationary distribution is the conditional distribution π η . We denote its transition kernel by Π η . The simulation step (3.2) is thus replaced by the following step: The most common choice of kernel is an accept-reject step which is carried out through a Metropolis-Hastings algorithm. Unfortunately, in our particular setting, we deal with large dimensions for the missing variables. This made us move to some other kind of MCMC methods, like a Gibbs sampler, to simulate our missing variables.
3.2. The transition of the MCMC method: a hybrid Gibbs sampler. If we consider the full vector (β, τ ) as a single vector of missing data, we can use the hybrid Gibbs sampler on R N PT as follows. For any b ∈ R and 1 ≤ j ≤ N , let us denote by β b→j the unique configuration which is equal to β everywhere except the coordinate j where β j b→j = b and by β −j the vector β without the coordinate j. Each coordinate of the deformation field β j is updated using a Metropolis-Hastings step where the proposal is given by the conditional distribution of β j |β −j , τ coming from the current Gaussian distribution with the corresponding parameters (pointed by τ ). Then, the last coordinates corresponding to the missing variable τ are drawn with respect to q(τ |β, y, η).
Even if this procedure provides an estimated parameter sequence which would theoretically converge toward the MAP estimator, in practice, as mentioned in [19], it would take a quite long time to reach its limit because of the trapping state problem: when a small number of observations are assigned to a component, the estimation of the component parameters is hardly concentrated and the probability of changing the label of an image to this component or from this component to another is really small (most of the time under the computer precision).
We can interpret this from an image analysis viewpoint: the first iteration of the algorithm gives a random label to the training set and computes the corresponding maximiser η = (θ, ρ). Then, for each image, according to its current label, it simulates a deformation field which only takes into account the parameters of this given component. Indeed, the simulation of β through the Gibbs sampler involves a proposal whose corresponding Markov chain has q(β|τ , y, η) as stationary distribution. Therefore, the deformation tries to match y to the deformed template of the given component τ . The deformation field tries to get a better connection between the component parameters and the observation, and there is only small probability that the observation given this deformation field will be closer to another component. The update of the label τ is therefore conditional to this deformation which would not leave much chance to switch component.
To overcome the trapping state problem, we will simulate the optimal label, using as many Markov chains in β as the number of components so that each component has a corresponding deformation which "computes" its distance to the observation. Then we can simulate the optimal deformation corresponding to that optimal label.
Since we aim to simulate (β, τ ) through a transition kernel that has q(β, τ |y, η) as stationary distribution, we simulate τ with a kernel whose stationary distribution is q(τ |y, η) and then β through a transition kernel that has q(β|τ , y, η) as stationary distribution.
For the first step, we need to compute the weights q(t|y i , η) ∝ q(t, y i |η) for all 1 ≤ t ≤ τ m and all 1 ≤ i ≤ n which cannot be easily reached. However, for any density function f , for any image y i and for any 1 ≤ t ≤ τ m , we have (3.7) Obviously the computation of this expectation w.r.t. the posterior distribution is not tractable either but we can approximate it by a Monte Carlo sum. However, we cannot easily simulate variables through the posterior distribution q(·|y i , t, η) as well, so we use some realisations of an ergodic Markov chain having q(·|y i , t, η) as stationary distribution instead of some independent realisations of this distribution. The solution we propose is the following: suppose we are at the k th iteration of the algorithm and let η be the current parameters. Given any initial deformation field ξ 0 ∈ R 2kg , we run, for where the normalisation is done such that their sum over t equals one, involving the dependence on all the auxiliary variables ξ i . The ergodic theorem ensures the convergence of our approximation toward the expected value. We then simulate Concerning the second step, we update β by re-running J times the hybrid Gibbs sampler Π η,τ on R N starting from a random initial point β 0 in a compact subset of R N . The size of J will depend on the iteration k of the SAEM algorithm in a sense that will be precised later, thus we now index it by k.
The density function f involved in the Monte Carlo sum above needs to be specified to get the convergence result proved in the last section of this paper. We show that using the prior on the deformation field enables to get the sufficient conditions for convergence. This density is the Gaussian density function and depends on the component we are working with: Algorithm 2 shows the detailed iteration. Remark 3. The use of one simulation of β per component is a point that was already used in [1] while computing the best matching, β * , for all components by minimising the corresponding energies. This gives as many β * as components for each image. Then, according to these best matchings, it computed the best component which therefore pointed the matching to consider.
3.3. Truncation on random boundaries. Since our missing data have a non-compact support, some of the convergence assumptions of such algorithms [14] are not satisfied. This leads to consider a truncation algorithm as suggested in [7] and extended in [2].
Let (K q ) q≥0 be an increasing sequence of compact subsets of S such as ∪ q≥0 K q = S and K q ⊂ int(K q+1 ), ∀q ≥ 0. Let K be a compact subset of R N . Let Π η be a transition kernel of an ergodic Markov chain on R N having π η as stationary distribution. We construct the homogeneous Markov chain ((β k , τ k , s k , κ k )) k≥0 as explained in Algorithm 1. As long as the stochastic approximation does not wander out the current compact set, we run our "SAEM-MCMC" algorithm. As soon as this previous condition is not satisfied, we reinitialise the sequences of s and (β, τ ) using a projection (for more details see [7]). The current compact is then enlarge. To point toward the current compact, we use a counter sequence (κ k ) k which remains unchanged when the previous condition is satisfied and increases to point toward a bigger compact when re-projecting.
4. Convergence theorem of the multicomponent procedure. In this particular section the variances of the components (σ 2 t ) 1≤t≤τm are fixed. Alleviating this condition is not straightforward and is an issue of our future work.
To prove the convergence of our parameter estimate toward the MAP, we have to go back to a convergence theorem which deals with general stochastic approximations. Indeed, the SAEM-MCMC algorithm introduced and detailed above is a Robbins-Monro type stochastic approximation procedure. One common tool to prove the w.p.1 convergence of such a stochastic approximation has been introduced by Kushner and Clark in [15]. However, some of the assumptions they require are intractable with our procedure (in particular concerning the mean field defined below). This leads us to slightly adapt the convergence theorem for stochastic approximations given in [7].
end for Algorithm 2 Transition step k → k + 1 using a hybrid Gibbs sampler on (β, τ ) where (e k ) k≥1 and (r k ) k≥1 are random processes defined on the same probability space taking their values in an open subset S of R ns ; h is referred to as the mean field of the algorithm; (r k ) k≥1 is a remainder term and (e k ) k≥1 is a stochastic excitation.
To be able to get a convergence result, we consider the truncated sequence (s k ) k defined as follow: let S a ⊂ S ands k = s k−1 + ∆ k h(s k−1 ) + ∆ k e k + ∆ k r k , where The projections k can be made through different ways (cf. [7]). We will use Delyon's theorem which gives sufficient conditions for the sequence (s k ) k≥0 truncated on random boundaries to converge with probability one. The theorem we state here is a generalisation of the theorem presented in [7]. Indeed, we have add the existence of an absorbing set for the stochastic approximation sequence. The proof of this theorem can be carried out the same way Delyon et al. do theirs adding the absorbing set. This is why it is not detailed here.
Theorem 4.1. We consider the sequence (s k ) k≥0 given by the truncated procedure (4.2). Assume that : (SA0') There exists a closed convex set S a ⊂ S such that for all k ≥ 0, s k ∈ S a w.p.1. Assumption (SA2), which involves a Lyapunov function w, replaces the usual condition that, w.p.1, the sequence comes back infinitely often in a compact set which is not easy to check in practice. In addition, assumptions (STAB1') and (STAB2) give a recurrent condition introducing a Lyapunov function W which controls the excursions outside the compact sets. The two Lyapunov functions w and W do not need to be the same. Another interesting point is that the truncation does not change the mean field and therefore the stationary points of the sequence.
This theorem does not ensure the convergence of the sequence to a maximum of the likelihood but to one of its critical points. To ensure that the critical point reached is a maximum, we would have to satisfy two other conditions (called (LOC1-2) in [7]) which are typical conditions. That is to say, it requires that the critical points are isolated and for every stationary points s * ∈ L ′ the Hessian matrix of the observed log-likelihood is negative definite.
We now want to apply this theorem to prove the convergence of our "SAEM like" procedure where the missing variables are not simulated through the posterior density function but by a kernel which can be as close as wanted -increasing J k -to this posterior law (generalising Theorem 3 in [7]).
Let us consider the following stochastic approximation: (β k , τ k ) are simulated by the transition kernel described in the previous section and which can be connected to the Robbins-Monro procedure using the notations introduced in [7]: let F = (F k ) k≥1 be the filtration where F k is the σ−algebra generated by the random variables (S 0 , β 1 , . . . , β k , τ 1 , . . . , τ k ), E π η is the expectation with respect to the posterior distribution π η and  Let (s k ) k≥0 be the truncated sequence defined in equation (4.2), K a compact set of R N and K 0 ⊂ S(R N ) a compact subset of S. Then, for all β 0 ∈ K, τ 0 ∈ T and s 0 ∈ K 0 , we have whereP β 0 ,τ 0 ,s0,0 is the probability measure associated with the chain (Z k = (β k , τ k , s k , κ k )) k≥0 starting at (β 0 , τ 0 , s 0 , 0).
The first assumption which concerns the two sequences involved in the algorithm, is not restrictive at all since these sequences can be chosen arbitrarily.
The second assumption, however, is more complex. This is required to satisfy the assumptions of Theorem 4.1. This is a condition we have not proved yet and is part of our future work.
Proof. The proof of this theorem is given in Section 6. We give here a quick sketch to emphasise the main difficulties and differences between our proof and the convergence proof of the SAEM algorithm given in [7].
Even if the only algorithmic difference between our algorithm and the SAEM algorithm is the simulation of the missing data which is not done with respect to the posterior law q(β, τ |y, η) but through an approximation which can be arbitrarily close, this yields a very different proof. Indeed, whereas for the SAEM algorithm, the stochastic approximation leads to a Robbins-Monro type equation (4.1) with no residual term r k , our method induces one. The first difficulty is therefore to prove that this residual term tends to 0 while the number of iterations k tends to infinity. Our proof is decomposed into two part, the first one concerning the deformation variable β and the second one the label τ . The first term requires to prove the geometric ergodicity of the Markov chain in β generated through our kernel. For this purpose, we prove some typical sufficient conditions which include the existence of a small set for the transition kernel and a drift condition. Then, we use for the second term some concentration inequalities for non stationary Markov chains to prove that the kernel associated with the label distribution converges toward the posterior distribution q(τ |y, η).
Proof. This is a direct consequence of the smoothness of the function s →η(s) on S and of Lemma 2 of [7] which links the sets L and L ′ .

USPS database.
To illustrate the previous algorithm for the deformable template model, we are considering handwritten digit images. For each digit, referred as class later, we learn two templates, the corresponding noise variances and the geometric covariance matrices. We use the USPS database which contains a training set of around 7000 images. Each picture is a (16P16) gray level image with intensity in [0, 2] where 0 corresponds to the black background. In Figure 5.1 we show some of the training images used for the statistical estimation. 5.1.1. General setting of the algorithm. A natural choice for the prior laws on α and Γ g is to set 0 for the mean on α and to induce the two covariance matrices by the metric of the spaces V p and V g involving the correlation between the landmarks through the kernels: define the square matrices M p (k, k ′ ) = K p (v p,j , v p,j ′ ) ∀1 ≤ k, k ′ ≤ k p , and M g (k, k ′ ) = K g (v g,j , v g,j ′ ) ∀1 ≤ k, k ′ ≤ k g . Then Σ p = M −1 p and Σ g = M −1 g . In our experiments, we have chosen Gaussian kernels for both K p and K g , where the standard deviations are fixed: σ p = 0.2 and σ g = 0.12 (for an estimation on [−1.5, 1.5] 2 and [−1, 1] respectively).
For the stochastic approximation step-size, we allow a heating period k h which corresponds to the absence of memory for the first iterations. This allows the Markov chain to reach an area of interest in the posterior probability density function q(β, τ |y, η) before exploring this particular region. In the experiments presented, the heating time k h lasts up to 150 iterations and the whole algorithm is stopped at, at most, 200 iterations depending on the data set (noisy or not). This number of iterations corresponds to a point when the convergence is reached. We choose, as suggested in [7] the step-size sequence as ∆ k = 1 for all 1 ≤ k ≤ k h and ∆ k = (k − k h ) −0.6 otherwise.
The multicomponent case has to face the problem of its computational time. Indeed, as we have to approximate the posterior density by running J k elements of τ m independent Markov chains, the computation time increases linearly with J k . In our experiments, we have chosen a fixed J for every EM iteration, J = 50. This is enough to get a good approximation thanks to the coupling between the iterations of the SAEM algorithm and the iterations of the Markov chains. Indeed, even if the parameter η is modified along the SAEM iterations, its successive jumps are small enough to ensure the convergence of the MCMC method. Intuitively speaking, it is equivalent to consider not only 50 iterations of the MCMC method but 50 times the number of SAEM iterations.

The estimated templates.
We are showing here the results of the statistical learning algorithm for our generative model. To avoid the problems shown in [2], we choose the same initialisation of the template parameter α as they did, that is to say, we set the initial value of α such that the corresponding I α is the mean of the gray-level training images.  In Figure 5.2, we show the two estimated templates obtained by the multicomponent procedure with 40 training examples per class. It appears that, as for the mode approximation algorithm which results are presented on this database in [1], the two components reached are meaningful, such as the 2 with and without loop or American and European 7. They even look alike.
In Figure 5.3, we show a second run of the algorithm with a different database, the training images are randomly selected in the whole USPS training set. We can see that there are some variability, in particular for digit 7 where there were no European 7 in this training set. This generates two different clusters still relevant for this digit. The other digits are quite stable, in particular the strongly constrained ones (like 3, 5, 8 or 9).

The photometric noise variance.
Even if we prove the convergence result for fixed component noise variances, we still try to learn them in the experiments. The same behaviour for our stochastic EM as for the mode approximation EM algorithm done in [1] is observed for the noise variances. Indeed, allowing the decomposition of the class into components enables the model to better fit the data yielding a lower residual noise. In addition, the stochastic algorithm enables to look around the whole posterior distribution and not only focusing on its mode which increases the accuracy of the geometric covariance and the template estimation. This yields lower noise required to explain the gap between the model and the truth. The evolution of the estimated variances for the two components of each digits are presented in Figure 5.4.
The convergence of this variance for some very constrained digits like digit 1 is faster. This is due to the well defined templates and geometric variability in the class which can be easily captured. Therefore, a very low level of noise is required very quickly. On the other hand, some very variable digits like digit 2 are slower to converge. The huge geometric variability adding to very complex shapes for the templates lead to a more difficult estimation and therefore more iterations before convergence.
Last point that can be noticed is the convergence of the European 7 which looks slower than the other component (American 7). The reason of this behaviour is that there are only two images of such a 7 in the training set and it takes a longer time for the algorithm to put together and only together these two shapes so that the clustering is better with respect to the likelihood. The other 7 does not suffer from this problem and converges faster.

The estimated geometric distribution.
To be able to compare the learnt geometry, we draw some synthetic examples using the mixture model with the learnt parameters. Even when the templates look similar, the separation between two components can be justified by the different geometry distributions. To show the effects of the geometry on the components, we have drawn some "2" with their respective parameters in the four top rows of Figure 5.5. For each component, we have drawn the deformation given by the variable β and its opposite −β since, as soon as one is learnt, because of the symmetry of the centred Gaussian distribution, the opposite deformation is learnt at the same time. This is why sometimes, one of the two looks strange whereas the other looks like some element of the training set.
The simulation is done using a common standard Gaussian distribution which is then multiplied by a square root of the covariance matrix we want to apply. We can see the effects of the covariance matrix on both templates and the large variability learnt. This has to be compared with the bottom rows of Figure 5.5, where the two samples are drawn on the one template but with the covariance matrix of the other one. Even if these six lines represent some "2"s, the bottom ones suffer from the geometrical tendency of the other cluster and do not look as natural. This shows the variability of the models into classes.

Medical images.
We also test the algorithm on a database which consists in 47 2D medical images. Each of them represents the splenium (back of the corpus calosum) and a part of the cerebellum. Some of the training images are shown in Figure 5.6 first row. The results of the estimation are presented in Figure 5.6, second row. The first picture presented, (a), is the gray level mean of the 47 images. The second one, (b), shows the estimated template computed with the Fast Approximation with Mode Algorithm presented in [1] for a single component model. This algorithm is an EM-like algorithm where the E step is simplified. The posterior distribution of the hidden variable is approximated by a Dirac distribution on its mode. This yields a deterministic algorithm, quite simple to implement but with no theoretical convergence properties. It shows a well contrasted splenium whereas the cerebellum remains a little bit blurry (note that it is still much better that the simple mean (a)).
This picture has to be compared with picture (c) which gives the estimated template computed with our algorithm with τ m = 1. The great improvement from the gray level mean of the images (a) or the FAM estimation (b) to our estimations is obvious. In particular, the splenium is still very contrasted, better localised and the cerebellum is reconstructed with several branches. The background presents several structures whereas the other estimates are blurry. The two anatomical shapes are relevant representants of the ones observed in the training set.
The estimation has been done while enabling the decomposition of the database into two components with our SAEM-MCMC algorithm presented here. The two estimated templates are shown in Figure 5.6 (d) and (e). The differences can be seen in particular on the shape of the splenium where the boundaries are more or less curved. The thickness of the splenium varies as well between the two estimates. The position of the fornix is also different, being closer to the boundary of the image. The number of branches in the two cerebella also tends to be different from one template to the other (4 in the first component and 5 in the second one).
The estimation suffers from the small number of images we have. This can be seen in the estimation of the background which is blurry in both images. To be able to explain the huge variability of the two anatomical shapes, more components would be interesting but at the same time more images required so that the components will not end up empty.
6. Proof of Theorem 4.2. We recall that in this section the variances of the components are fixed. This reduces the parameters θ t to (α t , Γ g,t ) for all 1 ≤ t ≤ τ m .
First let exhibit sufficient statistics for the model. The complete log-likelihood equals: where K β p α = z β I α and . denotes the Euclidean norm. This emphasises five sufficient statistics given in their matricial form for all 1 ≤ t ≤ τ m , Thus we apply the stochastic approximation at iteration k of the algorithm leading to: for 0 ≤ m ≤ 4 and rewrite the maximisation step. The weights and the covariance matrix are updated as follows: The photometric parameters are solution of the following system: s k,0,τ |Λ|+ap s k,0,τ (s k,4,τ + (α τ,k ) t s k,2,τ α τ,k − 2(α τ,k ) t s k,1,τ ) + a p σ 2 0 , which can be solved iteratively for each component τ starting with the previous values. In this part, all σ 2 τ,k are fixed and this leads to an explicit form for the parameters α τ,k .
We will now apply Theorem 4.1 to prove Theorem 4.2.
(SA0') is satisfied with the set S a defined by Thanks to the convexity of this set, the new value s k defined as a barycenter remains in S a . Assumption (SA1) is trivially satisfied since we can choose our step-size sequence (∆ k ) k .
(SA2) holds as already proved in [2] for the one component case with w(s) = −l(η(s)) such as (STAB1'(ii)) with the same function W (s) = −l(η(s)). These conditions imply the contraction property of the Lyapunov function w and the convergence of the stochastic approximation under some conditions on the perturbations.
We need to suppose, like in the one component case [2], that the critical points of our model are in a compact subset of S which stands for (STAB1'(i)). This is an assumption which has to be considered in a future work.
We will now focus on (STAB2) which is the assumption which gives the control of the perturbations required for the convergence.
We first show the convergence to zero of the remainder term |r k |1 W (s k−1 )≤M for any positive integer M . We denote by π k = πη (s k ) for any k ≥ 0. We have t,i and by We can now rewrite Denoting M η k−1 = max τ R N |S(β, τ )|q(β|τ , y, η k−1 )dβ, we obtain finally We will first show that the Gibbs sampler kernel Π η,τ satisfies a lower bound condition and a Drift condition (MDRI) to get its geometric ergodicity (as it has been done in [2]).
(MDRI) For any s ∈ S and any τ ∈ T , Πη (s),τ is irreducible and aperiodic. In addition there exists a function V : R N → [1, ∞[ such that for any p ≥ 1 and any compact subset K ⊂ S, there exist a set C, an integer m, constants 0 < κ < 1, B > 0 , δ > 0 and a probability measure ν such that Notation 1. Let (e j ) 1≤j≤N be the canonical basis of the β-space and for any 1 ≤ j ≤ N , let E η,τ ,j { β ∈ R N | β, e j η,τ = 0} be the orthogonal of Span{e j } and p η,τ ,j be the orthogonal projection on E η,τ ,j i.e.
Inequality (6.6) is equivalent to the existence of a small set C for the kernel Πη (s),τ independent of s ∈ K. We recall here the definition of a small set: Definition 1. (cf. [17]) A set E ∈ B(R N ) is called a small set for the kernel Π if there exist an integer m > 0 and a non trivial measure ν m on B(R N ), such that for all β ∈ E, B ∈ B(R N ), Π m (β, B) ≥ ν m (B).
When this holds, we say that E is ν m -small.
We now prove the following lemma: Lemma 6.1. Let E be a compact subset of R N and K be a compact subset of S, then E is a small set of R N for (Πη (s),τ ) s∈K,τ ∈T .
Proof. The transition probability kernel of our Markov chain on β is defined as follows : for coordinate j, the kernel is Then note that there exists a c > 0 such that for any η ∈ ΘP̺, any β ∈ R N and any b ∈ R, the acceptance rate r j (β j , b; β −j , η, τ ) is uniformly lower bounded by a c so that for any 1 ≤ j ≤ N and any non-negative function f , where g 0,1 is the probability density function of the standard N (0, 1). By induction, we have where p η,τ ,q,r = p η,τ ,r • p η,τ ,r−1 • · · · • p η,τ ,q for any integer q ≤ r and p η,τ ,N +1,N = Id R N . Let A η,τ ∈ L(R N ) be the linear mapping on z N 1 = (z 1 , · · · , z N ) defined by z j p η,τ ,j+1,N (e j )/ e j η,τ .
One easily checks that for any 1 ≤ k ≤ N , Span{ p η,τ ,j+1,N (e j ), k ≤ j ≤ N } = Span{e j , k ≤ j ≤ N } so that A η,τ is an invertible mapping. By a change of variable, we get where g µ,Σ stands for the probability density function of the normal law N (µ, Σ).
This property also implies the φ-irreducibility of the Markov Chain generated by Π η,τ . Moreover, the existence of a ν 1 -small set implies the aperiodicity of the chain (cf. [17]). Now consider the Drift condition (6.7).
This finishes the proof of (6.7) and in the same time the (MDRI).
Thus, applying the Proposition 1 and Lemma 6.4 to the Drift conditions of Lemmas 6.2 and 6.3, we get that each Gibbs sampler kernel Π η,τ is geometrically ergodic.
Let us now go back to the convergence of the first part of the residual term (6.4) towards 0.
Since for each τ the function β → S(β, τ ) belongs to L V , since we have proved that each transition kernel Π η,τ is geometrically ergodic and since the set V ℓ is compact, we can deduce that the first term (6.4) converges to zero as J k tends to infinity.
We now consider the second term (6.5). We first need to prove that M η k 1 W (s k−1 )≤M is uniformly bounded that is to say the integral of the sufficient statistics are uniformly bounded on {W (s k−1 ) ≤ M } ; we only need to focus on the sufficient statistic which is not bounded itself: let (j, m) ∈ {1, ..., 2k g } 2 : |β j β m |q(β|τ , y, η k−1 )dβ1 η k−1 ∈V ℓ ≤ |β j β m | q(β, τ , y, η k−1 ) q(τ , y, η k−1 ) where C(V ℓ ) is a constant depending only on the set V ℓ ,Γ g,τ is the diagonal block matrix with all the Γ g,τi given by the label vector τ and we have changed the variable in the last inequality and Q is a quadratic form in β whose coefficients are continuous functions of elements of the matrix Γ g . So we obtain that for all M > 0 there exists ℓ > 0 such that for all integer k we have: We now prove the convergence to 0 of the second term of the product involved in (6.5). Let us denote by R τ ,y,k for the term |R J k (τ |y, η k−1 ) − q(τ |y, η k−1 )|. Thus we have: S J k (τ i , y i |ξ τi,i , η k−1 ) s S J k (s, y i |ξ s,i , η k−1 ) − q(τ i , y i |η k−1 ) q(y i |η k−1 ) Q(ξ i )dξ i , where we denote by S J (t, y i |ξ t,i , η) the quantity 1 We write each term of this sum as follows: S J k (s, y i |ξ s,i , η k−1 )) q(y i |η k−1 ) τm s=1 S J k (s, y i |ξ s,i , η k−1 ) + (S J k (τ i , y i |ξ τi,i , η k−1 ) − q(τ i , y i |η k−1 )) τm s=1 S J k (s, y i |ξ s,i , η k−1 ) q(y i |η k−1 ) τm s=1 S J k (s, y i |ξ s,i , η k−1 ) .
Using the first inequality of Theorem 2 of [9], we get: P (A k,i,t ) ≤ c 1 exp −c 2 J k ζ 2 k q(t,yi|η k−1 ) 4 , where c 1 and c 2 are independent of k since (η k ) only moves in a compact set V ℓ thanks to the condition 1 W (s k−1 ≤M) . This yields: sup ξ,s S J k (s, y i |ξ s,i , η k−1 ) q(y i |η k−1 )