Abstract. This paper studies plug-and-play (PnP) Langevin sampling strategies for Bayesian inference in low-photon Poisson imaging problems, a challenging class of problems with significant applications in astronomy, medicine, and biology. PnP Langevin sampling offers a powerful framework for Bayesian image restoration, enabling accurate point estimation as well as advanced inference tasks, including uncertainty quantification and visualization analyses, and empirical Bayesian inference for automatic model parameter tuning. Herein, we leverage and adapt recent developments in this framework to tackle challenging imaging problems involving weakly informative Poisson data. Existing PnP Langevin algorithms are not well suited for low-photon Poisson imaging due to high solution uncertainty and poor regularity properties, such as exploding gradients and nonnegativity constraints. To address these challenges, we explore two strategies for extending Langevin PnP sampling to Poisson imaging models: (i) an accelerated PnP Langevin method that incorporates boundary reflections and a Poisson likelihood approximation, and (ii) a mirror sampling algorithm that leverages a Riemannian geometry to handle the constraints and the poor regularity of the likelihood without approximations. The effectiveness of these approaches is evaluated and contrasted through extensive numerical experiments and comparisons with state-of-the-art methods. The source code accompanying this paper is available at https://github.com/freyyia/pnp-langevin-poisson .

Text-to-image latent diffusion models (LDMs) have recently emerged as powerful generative models with great potential for solving inverse problems in imaging. However, leveraging such models in a Plug & Play (PnP), zero-shot manner remains challenging because it requires identifying a suitable text prompt for the unknown image of interest. Also, existing text-to-image PnP approaches are highly computationally expensive. We herein address these challenges by proposing a novel PnP inference paradigm specifically designed for embedding generative models within stochastic inverse solvers, with special attention to Latent Consistency Models (LCMs), which distill LDMs into fast generators. We leverage our framework to propose LAtent consisTency INverse sOlver (LATINO), the first zero-shot PnP framework to solve inverse problems with priors encoded by LCMs. Our conditioning mechanism avoids automatic differentiation and reaches SOTA quality in as little as 8 neural function evaluations. As a result, LATINO delivers remarkably accurate solutions and is significantly more memory and computationally efficient than previous approaches. We then embed LATINO within an empirical Bayesian framework that automatically calibrates the text prompt from the observed measurements by marginal maximum likelihood estimation. Extensive experiments show that prompt self-calibration greatly improves estimation, allowing LATINO with PRompt Optimization to define new SOTAs in image reconstruction quality and computational efficiency. The code is available at https://latino-pro.github.io

Abstract. Image restoration tasks often admit a wide range of solutions that are equally consistent with the observed data. Quantifying this uncertainty is crucial for the robust interpretation of restored images and their reliable use in science and decision-making. We consider the classification of images reconstructed from noisy and corrupted measurements, with special attention to the quantification of uncertainty in the delivered classification results. We address this problem by constructing a Bayesian statistical approach that combines learning-based image priors and image classifiers with explicit image observation models specified during inference. Following the manifold hypothesis, we assume that the image prior is supported on a submanifold of the ambient space—which we learn from uncorrupted training data using a variational autoencoder—and use as a classifier a support vector machine operating in this low-dimensional representation. The observation model is incorporated during inference time through its likelihood function. Bayesian computation is then efficiently performed by leveraging variants of the unadjusted Langevin algorithm that operate directly on the submanifold and are robust to multimodality. This results in a robust image classification method that provides uncertainty estimates that are provably well-posed, derived from Bayesian decision theory rigorously and transparently, and which incorporate physical and instrumental aspects of the data acquisition process through Bayes’ theorem. We demonstrate the effectiveness of the proposed approach through experiments on the MNIST and CelebA datasets, where we achieve accurate uncertainty estimates, as measured by the expected calibration error, even in challenging image restoration problems with significant inherent uncertainty.

Scientific imaging problems are often severely ill-posed and hence have significant intrinsic uncertainty. Accurately quantifying the uncertainty in the solutions to such problems is therefore critical for the rigorous interpretation of experimental results as well as for reliably using the reconstructed images as scientific evidence. Unfortunately, existing imaging methods are unable to quantify the uncertainty in the reconstructed images in a way that is robust to experiment replications. This paper presents a new uncertainty quantification methodology based on an equivariant formulation of the parametric bootstrap algorithm that leverages symmetries and invariance properties commonly encountered in imaging problems. Additionally, the proposed methodology is general and can be easily applied with any image reconstruction technique, including unsupervised training strategies that can be trained from observed data alone, thus enabling uncertainty quantification in situations where there is no ground truth data available. We demonstrate the proposed approach with a series of experiments and comparisons with alternative state-of-the-art uncertainty quantification strategies. In all our experiments, the proposed equivariant bootstrap delivers remarkably accurate high-dimensional confidence regions and outperforms the competing approaches in terms of estimation accuracy, uncertainty quantification accuracy, and computing time. These empirical findings are supported by a detailed theoretical analysis of equivariant bootstrap for linear estimators.

Abstract. This paper presents a new accelerated proximal Markov chain Monte Carlo methodology to perform Bayesian inference in imaging inverse problems with an underlying convex geometry. The proposed strategy takes the form of a stochastic relaxed proximal-point iteration that admits two complementary interpretations. For models that are smooth or regularized by Moreau–Yosida smoothing, the algorithm is equivalent to an implicit midpoint discretization of an overdamped Langevin diffusion targeting the posterior distribution of interest. This discretization is asymptotically unbiased for Gaussian targets and shown to converge in an accelerated manner for any target that is [Formula: see text]-strongly log-concave (i.e., requiring in the order of [Formula: see text] iterations to converge, similar to accelerated optimization schemes), comparing favorably to Pereyra, Vargas Mieles, and Zygalakis [ SIAM J. Imaging Sci., 13 (2020), pp. 905–935], which is only provably accelerated for Gaussian targets and has bias. For models that are not smooth, the algorithm is equivalent to a Leimkuhler–Matthews discretization of a Langevin diffusion targeting a Moreau–Yosida approximation of the posterior distribution of interest and hence achieves a significantly lower bias than conventional unadjusted Langevin strategies based on the Euler–Maruyama discretization. For targets that are [Formula: see text]-strongly log-concave, the provided nonasymptotic convergence analysis also identifies the optimal time step, which maximizes the convergence speed. The proposed methodology is demonstrated through a range of experiments related to image deconvolution with Gaussian and Poisson noise with assumption-driven and data-driven convex priors. Source codes for the numerical experiments of this paper are available from https://github.com/MI2G/accelerated-langevin-imla .

Abstract. This paper presents a novel stochastic optimization methodology to perform empirical Bayesian inference in semi-blind image deconvolution problems. Given a blurred image and a parametric class of possible operators, the proposed optimization approach automatically calibrates the parameters of the blur model by maximum marginal likelihood estimation, followed by (non-blind) image deconvolution by maximum a posteriori estimation conditionally to the estimated model parameters. In addition to the blur model, the proposed approach also automatically calibrates the noise level as well as any regularization parameters. The marginal likelihood of the blur, noise, and regularization parameters is generally computationally intractable, as it requires calculating several integrals over the entire solution space. Our approach addresses this difficulty by using a stochastic approximation proximal gradient optimization scheme, which iteratively solves such integrals by using a Moreau–Yosida regularized unadjusted Langevin Markov chain Monte Carlo algorithm. This optimization strategy can be easily and efficiently applied to any model that is log-concave and by using the same gradient and proximal operators that are required to compute the maximum a posteriori solution by convex optimization. We provide convergence guarantees for the proposed optimization scheme under realistic and easily verifiable conditions and subsequently demonstrate the effectiveness of the approach with a series of deconvolution experiments and comparisons with alternative strategies from the state of the art

We propose a new methodology for leveraging deep generative priors for Bayesian inference in imaging inverse problems. Modern Bayesian imaging often relies on score-based diffusion generative priors, which deliver remarkable point estimates but significantly underestimate uncertainty. Push-forward models such as variational auto-encoders and generative adversarial networks provide a robust alternative, leading to Bayesian models that are provably well-posed and which produce accurate uncertainty quantification results for small problems. However, push-forward models scale poorly to large problems because of issues related to bias, mode collapse and multimodality. We propose to address this difficulty by embedding a conditional deep generative prior within an empirical Bayesian framework. We consider generative priors with a super-resolution architecture, and perform inference by using a Bayesian computation strategy that simultaneously computes the maximum marginal likelihood estimate (MMLE) of the low-resolution image of interest, and draws Monte Carlo samples from the posterior distribution of the high-resolution image, conditionally to the observed data and the MMLE. The methodology is demonstrated with an image deblurring experiment and comparisons with the state-of-the-art.

Abstract. This paper studies a new and highly efficient Markov chain Monte Carlo (MCMC) methodology to perform Bayesian inference in low-photon imaging problems, with particular attention given to situations involving observation noise processes that deviate significantly from Gaussian noise, such as binomial, geometric, and low-intensity Poisson noise. These problems are challenging for many reasons. From an inferential viewpoint, low-photon numbers lead to severe identifiability issues, poor stability, and high uncertainty about the solution. Moreover, low-photon models often exhibit poor regularity properties that make efficient Bayesian computation difficult, e.g., hard nonnegativity constraints, nonsmooth priors, and log-likelihood terms with exploding gradients. More precisely, the lack of suitable regularity properties hinders the use of state-of-the-art Monte Carlo methods based on numerical approximations of the Langevin stochastic differential equation (SDE), as both the SDE and its numerical approximations behave poorly. We address this difficulty by proposing an MCMC methodology based on a reflected and regularized Langevin SDE, which is shown to be well-posed and exponentially ergodic under mild and easily verifiable conditions. This then allows us to derive four reflected proximal Langevin MCMC algorithms to perform Bayesian computation in low-photon imaging problems. The proposed approach is demonstrated with a range of experiments related to image deblurring, denoising, and inpainting under binomial, geometric, and Poisson noise.

Abstract. Developing efficient Bayesian computation algorithms for imaging inverse problems is challenging due to the dimensionality involved and because Bayesian imaging models are often not smooth. Current state-of-the-art methods often address these difficulties by replacing the posterior density with a smooth approximation that is amenable to efficient exploration by using Langevin Markov chain Monte Carlo (MCMC) methods. Such methods rely on gradient or proximal operators to exploit geometric information about the target posterior density and scale efficiently to large problems. An alternative approach is based on data augmentation and relaxation, where auxiliary variables are introduced in order to construct an approximate augmented posterior distribution that is amenable to efficient exploration by Gibbs sampling. This paper proposes a new accelerated proximal MCMC method called latent space SK-ROCK (ls-SK-ROCK), which tightly combines the benefits of the two aforementioned strategies. Additionally, instead of viewing the augmented posterior distribution as an approximation of the original model, we propose to consider it as a generalization of this model. Following on from this, we empirically show that there is a range of values for the relaxation parameter for which the accuracy of the model improves and propose a stochastic optimization algorithm to automatically identify the optimal amount of relaxation for a given problem. In this regime, ls-SK-ROCK converges faster than competing approaches from the state of the art, and it also achieves better accuracy since the underlying augmented Bayesian model has a higher Bayesian evidence. The proposed methodology is demonstrated with a range of numerical experiments related to image deblurring and inpainting, as well as with comparisons with alternative approaches from the state of the art. An open-source implementation of the proposed MCMC methods is available from https://github.com/luisvargasmieles/ls-MCMC .

This paper proposes a new methodology for performing Bayesian inference in imaging inverse problems where the prior knowledge is available in the form of training data. Following the manifold hypothesis, we adopt a data-driven prior that is supported on a submanifold of the ambient space, which we can learn from the training data using a generative model, such as a variational autoencoder or generative adversarial network. We establish the existence and well-posedness of the associated posterior distribution and posterior moments under easily verifiable conditions, providing a rigorous underpinning for Bayesian estimators and uncertainty quantification analyses. Bayesian computation is performed using a parallel tempered version of the pCN algorithm on the manifold, which is shown to be ergodic and robust to the nonconvex nature of these data-driven models. In addition to point estimators and uncertainty quantification analyses, we derive a model misspecification test to automatically detect situations where the data-driven prior is unreliable, and we explain how to identify the dimension of the latent space directly from the training data. The proposed approach is illustrated with a range of experiments with the MNIST dataset and is compared with some variational and message passing image reconstruction approaches from the state of the art that also use data-driven regularization. A model accuracy analysis suggests that the Bayesian probabilities reported by the proposed data-driven models are also accurate under a frequentist definition of probability, suggesting that the learnt prior is close to the true marginal distribution of the unknown image.

Since the seminal work of Venkatakrishnan et al. in 2013, Plug & Play (PnP) methods have become ubiquitous in Bayesian imaging. These methods derive Minimum Mean Square Error (MMSE) or Maximum A Posteriori (MAP) estimators for inverse problems in imaging by combining an explicit likelihood function with a prior that is implicitly defined by an image denoising algorithm. The PnP algorithms proposed in the literature mainly differ in the iterative schemes they use for optimisation or for sampling. In the case of optimisation schemes, some recent works guarantee the convergence to a fixed point, albeit not necessarily a MAP estimate. In the case of sampling schemes, to the best of our knowledge, there is no known proof of convergence. There also remain important open questions regarding whether the underlying Bayesian models and estimators are well defined, well-posed, and have the basic regularity properties required to support these numerical schemes. To address these limitations, this paper develops theory, methods, and provably convergent algorithms for performing Bayesian inference with PnP priors. We introduce two algorithms: 1) PnP-ULA (Unadjusted Langevin Algorithm) for Monte Carlo sampling and MMSE inference; and 2) PnP-SGD (Stochastic Gradient Descent) for MAP inference. Using recent results on the quantitative convergence of Markov chains, we establish detailed convergence guarantees for these two algorithms under realistic assumptions on the denoising operators used, with special attention to denoisers based on deep neural networks. We also show that these algorithms approximately target a decision-theoretically optimal Bayesian model that is well-posed. The proposed algorithms are demonstrated on several canonical problems such as image deblurring, inpainting, and denoising, where they are used for point estimation as well as for uncertainty visualisation and quantification.

Discretized Langevin diffusions are efficient Monte Carlo methods for sampling from high dimensional target densities that are log-Lipschitz-smooth and (strongly) log-concave. In particular, the Euclidean Langevin Monte Carlo sampling algorithm has received much attention lately, leading to a detailed understanding of its non-asymptotic convergence properties and of the role that smoothness and log-concavity play in the convergence rate. Distributions that do not possess these regularity properties can be addressed by considering a Riemannian Langevin diffusion with a metric capturing the local geometry of the log-density. However, the Monte Carlo algorithms derived from discretizations of such Riemannian Langevin diffusions are notoriously difficult to analyze. In this paper, we consider Langevin diffusions on a Hessian-type manifold and study a discretization that is closely related to the mirror-descent scheme. We establish for the first time a non-asymptotic upper-bound on the sampling error of the resulting Hessian Riemannian Langevin Monte Carlo algorithm. This bound is measured according to a Wasserstein distance induced by a Riemannian metric ground cost capturing the squared Hessian structure and closely related to a self-concordance-like condition. The upper-bound implies, for instance, that the iterates contract toward a Wasserstein ball around the target density whose radius is made explicit. Our theory recovers existing Euclidean results and can cope with a wide variety of Hessian metrics related to highly non-flat geometries.

This paper presents a detailed theoretical analysis of the three stochastic approximation proximal gradient algorithms proposed in our companion paper [A. F. Vidal et al., SIAM J. Imaging Sci., 13 (2020), pp. 1945--1989] to set regularization parameters by marginal maximum likelihood estimation. We prove the convergence of a more general stochastic approximation scheme that includes the three algorithms of [A. F. Vidal et al., SIAM J. Imaging Sci., 13 (2020), pp. 1945--1989] as special cases. This includes asymptotic and nonasymptotic convergence results with natural and easily verifiable conditions, as well as explicit bounds on the convergence rates. Importantly, the theory is also general in that it can be applied to other intractable optimization problems. A main novelty of the work is that the stochastic gradient estimates of our scheme are constructed from inexact proximal Markov chain Monte Carlo samplers. This allows the use of samplers that scale efficiently to large problems and for which we have precise theoretical guarantees.

We present a highly efficient proximal Markov chain Monte Carlo methodology to perform Bayesian computation in imaging problems. Similarly to previous proximal Monte Carlo approaches, the proposed method is derived from an approximation of the Langevin diffusion. However, instead of the conventional Euler--Maruyama approximation that underpins existing proximal Monte Carlo methods, here we use a state-of-the-art orthogonal Runge--Kutta--Chebyshev stochastic approximation [A. Abdulle, I. Aimuslimani, and G. Vilmart, SIAM/ASA J. Uncertain. Quantif., 6 (2018), pp. 937--964] that combines several gradient evaluations to significantly accelerate its convergence speed, similarly to accelerated gradient optimization methods. The proposed methodology is demonstrated via a range of numerical experiments, including non-blind image deconvolution, hyperspectral unmixing, and tomographic reconstruction, with total-variation and $\ell_1$-type priors. Comparisons with Euler-type proximal Monte Carlo methods confirm that the Markov chains generated with our method exhibit significantly faster convergence speeds, achieve larger effective sample sizes, and produce lower mean-square estimation errors at equal computational budget.

Many imaging problems require solving an inverse problem that is ill-conditioned or ill-posed. Imaging methods typically address this difficulty by regularizing the estimation problem to make it well- posed. This often requires setting the value of the so-called regularization parameters that control the amount of regularization enforced. These parameters are notoriously difficult to set a priori and can have a dramatic impact on the recovered estimates. In this work, we propose a general empirical Bayesian method for setting regularization parameters in imaging problems that are convex w.r.t. the unknown image. Our method calibrates regularization parameters directly from the observed data by maximum marginal likelihood estimation and can simultaneously estimate multiple regularization parameters. Furthermore, the proposed algorithm uses the same basic operators as proximal optimization algorithms, namely gradient and proximal operators, and it is therefore straightforward to apply to problems that are currently solved by using proximal optimization techniques. Our methodology is demonstrated with a range of experiments and comparisons with alternative approaches from the literature. The considered experiments include image denoising, nonblind image deconvolution, and hyperspectral unmixing, using synthesis and analysis priors involving the $\ell_1$, total-variation, total-variation and $\ell_1$, and total-generalized-variation pseudonorms. A detailed theoretical analysis of the proposed method is presented in our companion paper [V. De Bortoli et al., SIAM J. Imaging Sci., 13 (2020), pp. 1990--2028].

Maximum-a-posteriori (MAP) estimation is the main Bayesian estimation methodology in imaging sciences, where high dimensionality is often addressed by using Bayesian models that are log-concave and whose posterior mode can be computed efficiently by convex optimization. However, despite its success and wide adoption, MAP estimation is not theoretically well understood yet. In particular, the prevalent view in the community is that MAP estimation is not proper Bayesian estimation in the sense of Bayesian decision theory because it does not minimize a meaningful expected loss function (unlike the minimum mean squared error (MMSE) estimator that minimizes the mean squared loss). This paper addresses this theoretical gap by presenting a general decision-theoretic derivation of MAP estimation in Bayesian models that are log-concave. A main novelty is that our analysis is based on differential geometry and proceeds as follows. First, we use the underlying convex geometry of the Bayesian model to induce a Riemannian geometry on the parameter space. We then use differential geometry to identify the so-called natural or canonical loss function to perform Bayesian point estimation in that Riemannian manifold. For log-concave models, this canonical loss coincides with the Bregman divergence associated with the negative log posterior density. Following on from this, we show that the MAP estimator is the only Bayesian estimator that minimizes the expected canonical loss, and that the posterior mean or MMSE estimator minimizes the dual canonical loss. We then study the question of MAP and MMSE estimation performance in high dimensions. Precisely, we establish a universal bound on the expected canonical error as a function of image dimension, providing new insights on the good empirical performance observed in convex problems. Together, these results provide a new understanding of MAP and MMSE estimation in log-concave settings, and of the multiple beneficial roles that convex geometry plays in imaging problems. Finally, we illustrate this new theory by analyzing the regularization-by-denoising Bayesian models, a class of state-of-the-art imaging models where priors are defined implicitly through image denoising algorithms, and an image denoising model with a wavelet shrinkage prior.

We propose a Bayesian uncertainty quantification method for large-scale imaging inverse problems. Our method applies to all Bayesian models that are log-concave, where maximum a posteriori (MAP) estimation is a convex optimization problem. The method is a framework to analyze the confidence in specific structures observed in MAP estimates (e.g., lesions in medical imaging, celestial sources in astronomical imaging), to enable using them as evidence to inform decisions and conclusions. Precisely, following Bayesian decision theory, we seek to assert the structures under scrutiny by performing a Bayesian hypothesis test that proceeds as follows: first, it postulates that the structures are not present in the true image, and then seeks to use the data and prior knowledge to reject this null hypothesis with high probability. Computing such tests for imaging problems is generally very difficult because of the high dimensionality involved. A main feature of this work is to leverage probability concentration phenomena and the underlying convex geometry to formulate the Bayesian hypothesis test as a convex problem, which we then efficiently solve by using scalable optimization algorithms. This allows scaling to high-resolution and high-sensitivity imaging problems that are computationally unaffordable for other Bayesian computation approaches. We illustrate our methodology, dubbed BUQO (Bayesian Uncertainty Quantification by Optimization), on a range of challenging Fourier imaging problems arising in astronomy and medicine. MATLAB code for the proposed uncertainty quantification method is available on GitHub.

This paper presents an empirical Bayesian method to estimate regularisation parameters in imaging inverse problems. The method calibrates regularisation parameters directly from the observed data by maximum marginal likelihood estimation, and is useful for inverse problems that are convex. A main novelty is that maximum likelihood estimation is performed efficiently by using a stochastic proximal gradient algorithm that is driven by two proximal Markov chain Monte Carlo samplers, intimately combining modern optimisation and sampling techniques. The proposed methodology is illustrated with an application to total-variation image deconvolution, where it compares favourably to alternative Bayesian and non-Bayesian approaches from the state of the art.

Modern imaging methods rely strongly on Bayesian inference techniques to solve challenging imaging problems. Currently, the predominant Bayesian computation approach is convex optimization, which scales very efficiently to high-dimensional image models and delivers accurate point estimation results. However, in order to perform more complex analyses, for example, image uncertainty quantification or model selection, it is necessary to use more computationally intensive Bayesian computation techniques such as Markov chain Monte Carlo methods. This paper presents a new and highly efficient Markov chain Monte Carlo methodology to perform Bayesian computation for high-dimensional models that are log-concave and nonsmooth, a class of models that is central in imaging sciences. The methodology is based on a regularized unadjusted Langevin algorithm that exploits tools from convex analysis, namely, Moreau--Yoshida envelopes and proximal operators, to construct Markov chains with favorable convergence properties. In addition to scaling efficiently to high-dimensions, the method is straightforward to apply to models that are currently solved by using proximal optimization algorithms. We provide a detailed theoretical analysis of the proposed methodology, including asymptotic and nonasymptotic convergence results with easily verifiable conditions, and explicit bounds on the convergence rates. The proposed methodology is demonstrated with four experiments related to image deconvolution and tomographic reconstruction with total-variation and $\ell_1$ priors, where we conduct a range of challenging Bayesian analyses related to uncertainty quantification, hypothesis testing, and model selection in the absence of ground truth.

This paper presents a detailed theoretical analysis of the Langevin Monte Carlo sampling algorithm recently introduced in Durmus et al. (Efficient Bayesian computation by proximal Markov chain Monte Carlo: when Langevin meets Moreau, 2016) when applied to log-concave probability distributions that are restricted to a convex body $K$. This method relies on a regularisation procedure involving the Moreau-Yosida envelope of the indicator function associated with $K$. Explicit convergence bounds in total variation norm and in Wasserstein distance of order $1$ are established. In particular, we show that the complexity of this algorithm given a first order oracle is polynomial in the dimension of the state space. Finally, some numerical experiments are presented to compare our method with competing MCMC approaches from the literature.

Solutions to inverse problems that are ill-conditioned or ill-posed may have significant intrinsic uncertainty. Unfortunately, analyzing and quantifying this uncertainty is very challenging, particularly in high-dimensional problems. As a result, while most modern mathematical imaging methods produce impressive point estimation results, they are generally unable to quantify the uncertainty in the solutions delivered. This paper presents a new general methodology for approximating Bayesian high-posterior-density credibility regions in inverse problems that are convex and potentially very high-dimensional. The approximations are derived by using recent concentration of measure results related to information theory for log-concave random vectors. A remarkable property of the approximations is that they can be computed very efficiently, even in large-scale problems, by using standard convex optimization techniques. In particular, they are available as a by-product in problems solved by maximum-a-posteriori estimation. The approximations also have favorable theoretical properties, namely, they outer-bound the true high-posterior-density credibility regions, and they are stable with respect to model dimension. The proposed methodology is illustrated on two high-dimensional imaging inverse problems related to tomographic reconstruction and sparse deconvolution, where the approximations are used to perform Bayesian hypothesis tests and explore the uncertainty about the solutions, and where proximal Markov chain Monte Carlo algorithms are used as a benchmark to compute exact credible regions and measure the approximation error.

This paper presents a new Bayesian estimation technique for hidden Potts-Markov random fields with unknown regularisation parameters, with application to fast unsupervised K -class image segmentation. The technique is derived by first removing the regularisation parameter from the Bayesian model by marginalisation, followed by a small-variance-asymptotic (SVA) analysis in which the spatial regularisation and the integer-constrained terms of the Potts model are decoupled. The evaluation of this SVA Bayesian estimator is then relaxed into a problem that can be computed efficiently by iteratively solving a convex total-variation denoising problem and a least-squares clustering ( K -means) problem, both of which can be solved straightforwardly, even in high-dimensions, and with parallel computing techniques. This leads to a fast fully unsupervised Bayesian image segmentation methodology in which the strength of the spatial regularisation is adapted automatically to the observed image during the inference procedure, and that can be easily applied in large 2D and 3D scenarios or in applications requiring low computing times. Experimental results on synthetic and real images, as well as extensive comparisons with state-of-the-art algorithms, confirm that the proposed methodology offer extremely fast convergence and produces accurate segmentation results, with the important additional advantage of self-adjusting regularisation parameters.

This paper presents a new Bayesian collaborative sparse regression method for linear unmixing of hyperspectral images. Our contribution is twofold; first, we propose a new Bayesian model for structured sparse regression in which the supports of the sparse abundance vectors are a priori spatially correlated across pixels (i.e., materials are spatially organized rather than randomly distributed at a pixel level). This prior information is encoded in the model through a truncated multivariate Ising Markov random field, which also takes into consideration the facts that pixels cannot be empty (i.e., there is at least one material present in each pixel), and that different materials may exhibit different degrees of spatial regularity. Second, we propose an advanced Markov chain Monte Carlo algorithm to estimate the posterior probabilities that materials are present or absent in each pixel, and, conditionally to the maximum marginal a posteriori configuration of the support, compute the minimum mean squared error estimates of the abundance vectors. A remarkable property of this algorithm is that it self-adjusts the values of the parameters of the Markov random field, thus relieving practitioners from setting regularization parameters by cross-validation. The performance of the proposed methodology is finally demonstrated through a series of experiments with synthetic and real data and comparisons with other algorithms from the literature.

This paper presents a fast converging Riemannian steepest descent method for nonparametric statistical active contour models, with application to image segmentation. Unlike other fast algorithms, the proposed method is general and can be applied to any statistical active contour model from the exponential family, which comprises most of the models considered in the literature. This is achieved by first identifying the intrinsic statistical manifold associated with this class of active contours, and then constructing a steepest descent on that manifold. A key contribution of this paper is to derive a general and tractable closed-form analytic expression for the manifold's Riemannian metric tensor, which allows computing discrete gradient flows efficiently. The proposed methodology is demonstrated empirically and compared with other state of the art approaches on several standard test images, a phantom positron-emission-tomography scan and a B-mode echography of in-vivo human dermis.

This letter considers the problem of computing the Cramer–Rao bound for the parameters of a Markov random field. Computation of the exact bound is not feasible for most fields of interest because their likelihoods are intractable and have intractable derivatives. We show here how it is possible to formulate the computation of the bound as a statistical inference problem that can be solve approximately, but with arbitrarily high accuracy, by using a Monte Carlo method. The proposed methodology is successfully applied on the Ising and the Potts models.

This paper addresses the problem of estimating the Potts parameter β jointly with the unknown parameters of a Bayesian model within a Markov chain Monte Carlo (MCMC) algorithm. Standard MCMC methods cannot be applied to this problem because performing inference on β requires computing the intractable normalizing constant of the Potts model. In the proposed MCMC method, the estimation of β is conducted using a likelihood-free Metropolis-Hastings algorithm. Experimental results obtained for synthetic data show that estimating β jointly with the other unknown parameters leads to estimation results that are as good as those obtained with the actual value of β. On the other hand, choosing an incorrect value of β can degrade estimation performance significantly. To illustrate the interest of this method, the proposed algorithm is successfully applied to real bidimensional SAR and tridimensional ultrasound images.

This paper addresses the problem of jointly estimating the statistical distribution and segmenting lesions in multiple-tissue high-frequency skin ultrasound images. The distribution of multiple-tissue images is modeled as a spatially coherent finite mixture of heavy-tailed Rayleigh distributions. Spatial coherence inherent to biological tissues is modeled by enforcing local dependence between the mixture components. An original Bayesian algorithm combined with a Markov chain Monte Carlo method is then proposed to jointly estimate the mixture parameters and a label-vector associating each voxel to a tissue. More precisely, a hybrid Metropolis-within-Gibbs sampler is used to draw samples that are asymptotically distributed according to the posterior distribution of the Bayesian model. The Bayesian estimators of the model parameters are then computed from the generated samples. Simulation results are conducted on synthetic data to illustrate the performance of the proposed estimation strategy. The method is then successfully applied to the segmentation of in vivo skin tumors in high-frequency 2-D and 3-D ultrasound images.

This paper addresses the problem of estimating the statistical distribution of multiple-tissue non-stationary ultrasound images of skin. The distribution of multiple-tissue images is modeled as a finite mixture of Heavy-Tailed Rayleigh distributions. An original Bayesian algorithm combined with a Markov chain Monte Carlo method is then derived to jointly estimate the mixture parameters and a label vector associating each voxel to a tissue. Precisely, a hybrid Metropolis-within-Gibbs sampler is proposed to draw samples that are asymptotically distributed according to the posterior distribution of the Bayesian model. These samples are then used to compute the Bayesian estimators of the model parameters. Simulation results are conducted on synthetic data to illustrate the performance of the proposed estimation strategy. The method is then successfully applied to the detection of an in-vivo skin lesion in a high frequency 3D ultrasound image.