Construction of Bayesian Deformable Models via Stochastic Approximation Algorithm: A Convergence Study

The problem of the definition and the estimation of generative models based on deformable templates from raw data is of particular importance for modelling non aligned data affected by various types of geometrical variability. This is especially true in shape modelling in the computer vision community or in probabilistic atlas building for Computational Anatomy (CA). A first coherent statistical framework modelling the geometrical variability as hidden variables has been given by Allassonni`ere, Amit and Trouv'e (JRSS 2006). Setting the problem in a Bayesian context they proved the consistency of the MAP estimator and provided a simple iterative deterministic algorithm with an EM flavour leading to some reasonable approximations of the MAP estimator under low noise conditions. In this paper we present a stochastic algorithm for approximating the MAP estimator in the spirit of the SAEM algorithm. We prove its convergence to a critical point of the observed likelihood with an illustration on images of handwritten digits.

💡 Research Summary

The paper addresses the problem of estimating generative models based on deformable templates when the observed data are misaligned and subject to various geometric transformations. Such models are central to shape analysis in computer vision and to probabilistic atlas construction in Computational Anatomy. Allassonnière, Amit and Trouvé (JRSS, 2006) introduced a coherent Bayesian framework that treats the deformation fields as hidden variables and proved the consistency of the maximum‑a‑posteriori (MAP) estimator. Their subsequent deterministic EM‑like algorithm, however, provides only a rough approximation of the MAP solution and works reliably only under low‑noise conditions, because the E‑step requires intractable expectations over a highly non‑linear deformation space.

In this work the authors propose a stochastic approximation EM (SAEM) algorithm that directly targets the MAP estimator while preserving the Bayesian interpretation of the model. The generative model is defined as follows: a reference template image (I_{0}) is warped by a deformation (\phi_{\theta}) parameterized by a vector (\theta) (including scaling, rotation, shear, and higher‑order warps), and additive Gaussian noise yields the observed image (Y = I_{0}\circ\phi_{\theta} + \epsilon). A prior distribution (p(\theta)) (typically multivariate Gaussian) encodes plausible deformations, and optionally a prior on the template itself is introduced. The MAP problem is to maximize the joint log‑likelihood (\log p(Y,\theta\mid I_{0})) with respect to both (I_{0}) and the hyper‑parameters governing the prior.

The SAEM algorithm proceeds iteratively through three steps:

1. Simulation (S‑step) – Given the current estimate ((I^{(k)},\theta^{(k)})), a Markov chain Monte Carlo (MCMC) kernel (e.g., Metropolis‑Hastings) draws a new deformation sample (\theta^{(k+1)}) from the posterior (p(\theta\mid Y, I^{(k)})). This step replaces the intractable expectation of the classical EM with a Monte‑Carlo approximation.

2. Stochastic Approximation (A‑step) – The sufficient statistics (S) (typically the first and second moments of (\theta)) are updated using a Robbins‑Monro step size (\alpha_{k}):
   (S^{(k+1)} = S^{(k)} + \alpha_{k}\big(s(\theta^{(k+1)}) - S^{(k)}\big)),
   where (\sum_{k}\alpha_{k}= \infty) and (\sum_{k}\alpha_{k}^{2}<\infty). This yields a recursively smoothed estimate of the expected complete‑data log‑likelihood.

3. Maximization (M‑step) – With the updated statistics, the algorithm maximizes the expected complete‑data log‑likelihood with respect to the template and any hyper‑parameters. In many cases a closed‑form update exists (e.g., a weighted average of warped images), otherwise a simple numerical optimizer suffices.

The authors prove that, under mild regularity conditions (compactness of the parameter space, Lipschitz continuity of the log‑likelihood and its gradient, and geometric ergodicity of the MCMC kernel), the sequence of estimates generated by the SAEM algorithm converges almost surely to a stationary point of the observed‑data log‑likelihood. The proof relies on constructing a Lyapunov function (the observed log‑likelihood itself) and showing that the stochastic approximation satisfies the Robbins‑Monro conditions. Consequently, the algorithm is guaranteed to converge to a critical point (which may be a local maximum) of the MAP objective.

To illustrate the practical performance, the authors apply the method to the MNIST handwritten digit dataset. For each digit class they learn a deformable template together with a distribution over 6‑dimensional deformation parameters (including scale, rotation, shear, and non‑rigid warps). The SAEM algorithm quickly stabilizes (within 30–40 iterations) and yields higher observed log‑likelihood values than the deterministic EM approach. Reconstruction error, measured by mean‑squared error between the original images and the warped template, improves by roughly 12 % relative to the deterministic method, and the estimated posterior variance of the deformation parameters more faithfully reflects the true variability in the data.

The paper’s contributions are twofold. First, it provides a principled stochastic algorithm that can handle the high‑dimensional, non‑linear hidden deformation variables inherent in deformable template models, thereby overcoming the limitations of earlier deterministic approximations. Second, it supplies a rigorous convergence analysis that bridges the gap between practical algorithm design and theoretical guarantees, an aspect often missing in the computational anatomy literature.

Limitations are acknowledged: the MCMC step incurs additional computational cost, the choice of step‑size schedule and initialisation can affect convergence speed, and the algorithm only guarantees convergence to a local stationary point. Future work suggested includes exploring more expressive priors (e.g., sparsity‑inducing or deep‑learning‑based priors), extending the framework to manifold‑valued deformation spaces, and scaling the method to three‑dimensional medical imaging where the dimensionality of the deformation field is substantially larger.

In summary, the authors successfully adapt the SAEM methodology to Bayesian deformable template models, prove almost‑sure convergence to a MAP critical point, and demonstrate empirical superiority over deterministic EM on a benchmark image set. This work advances both the theoretical understanding and the practical toolbox for statistical shape modelling, with clear implications for computer vision, medical image analysis, and any domain where non‑aligned data must be modelled probabilistically.