81 research outputs found
Accelerating proximal Markov chain Monte Carlo by using an explicit stabilised method
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 that combines several gradient evaluations to
significantly accelerate its convergence speed, similarly to accelerated
gradient optimisation 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 -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.Comment: 28 pages, 13 figures. Accepted for publication in SIAM Journal on
Imaging Sciences (SIIMS
Revisiting maximum-a-posteriori estimation in log-concave models
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 optimisation. Despite its success and wide
adoption, MAP estimation is not theoretically well understood yet. The
prevalent view in the community is that MAP estimation is not proper Bayesian
estimation in a decision-theoretic sense because it does not minimise a
meaningful expected loss function (unlike the minimum mean squared error (MMSE)
estimator that minimises the mean squared loss). This paper addresses this
theoretical gap by presenting a 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 is the Bregman divergence associated with the negative log posterior
density. We then show that the MAP estimator is the only Bayesian estimator
that minimises the expected canonical loss, and that the posterior mean or MMSE
estimator minimises the dual canonical loss. We also study the question of MAP
and MSSE estimation performance in large scales and establish a universal bound
on the expected canonical error as a function of dimension, offering new
insights into the good performance observed in convex problems. These results
provide a new understanding of MAP and MMSE estimation in log-concave settings,
and of the multiple roles that convex geometry plays in imaging problems.Comment: Accepted for publication in SIAM Imaging Science
Seismic Tomography Using Variational Inference Methods
Seismic tomography is a methodology to image the interior of solid or fluid
media, and is often used to map properties in the subsurface of the Earth. In
order to better interpret the resulting images it is important to assess
imaging uncertainties. Since tomography is significantly nonlinear, Monte Carlo
sampling methods are often used for this purpose, but they are generally
computationally intractable for large datasets and high-dimensional parameter
spaces. To extend uncertainty analysis to larger systems we use variational
inference methods to conduct seismic tomography. In contrast to Monte Carlo
sampling, variational methods solve the Bayesian inference problem as an
optimization problem, yet still provide probabilistic results. In this study,
we applied two variational methods, automatic differential variational
inference (ADVI) and Stein variational gradient descent (SVGD), to 2D seismic
tomography problems using both synthetic and real data and we compare the
results to those from two different Monte Carlo sampling methods. The results
show that variational inference methods can produce accurate approximations to
the results of Monte Carlo sampling methods at significantly lower
computational cost, provided that gradients of parameters with respect to data
can be calculated efficiently. We expect that the methods can be applied
fruitfully to many other types of geophysical inverse problems.Comment: 26 pages, 14 figure
Accelerating Bayesian computation in imaging
The dimensionality and ill-posedness often encountered in imaging inverse problems are a challenge for Bayesian computational methods, particularly for state-of-the-art sampling alternatives based on the Euler-Maruyama discretisation of the Langevin diffusion process. In this thesis, we address this difficulty and propose alternatives to accelerate Bayesian computation in imaging inverse problems, focusing on its computational aspects.
We introduce, as our first contribution, a highly efficient proximal Markov chain Monte Carlo (MCMC) methodology, based on a state-of-the-art approximation known as the proximal stochastic orthogonal Runge-Kutta-Chebyshev (SK-ROCK) method. It has the advantage of cleverly combining multiple gradient evaluations to significantly speed up convergence, similar to accelerated gradient optimisation techniques. We rigorously demonstrate the acceleration of the Markov chains in the 2-Wasserstein distance for Gaussian models as a function of the condition number κ.
In our second contribution, we propose a more sophisticated MCMC sampler, based on the careful integration of two advanced proximal Langevin MCMC methods, SK-ROCK and split Gibbs sampling (SGS), each of which uses a unique approach to accelerate convergence. More precisely, we show how to integrate the proximal SK-ROCK sampler with the model augmentation and relaxation method used by SGS at the level of the Langevin diffusion process, to speed up Bayesian computation at the expense of asymptotic bias. This leads to a new, faster proximal SK-ROCK sampler that combines the accelerated quality of the original sampler with the computational advantages of augmentation and relaxation.
Additionally, we propose the augmented and relaxed model to be considered a generalisation of the target model rather than an approximation that situates relaxation in a bias-variance trade-off. As a result, we can carefully calibrate the amount of relaxation to boost both model accuracy (as determined by model evidence) and sampler convergence speed. To achieve this, we derive an empirical Bayesian method that automatically estimates the appropriate level of relaxation via maximum marginal likelihood estimation.
The proposed methodologies are demonstrated in several numerical experiments related to image deblurring, hyperspectral unmixing, tomographic reconstruction and inpainting. Comparisons with Euler-type proximal Monte Carlo approaches confirm that the Markov chains generated with our methods exhibit significantly faster convergence speeds, achieve larger effective sample sizes, and produce lower mean square estimation errors with the same computational budget
Consistent Inversion of Noisy Non-Abelian X-Ray Transforms
For a simple surface, the non-linear statistical inverse problem of
recovering a matrix field from discrete, noisy
measurements of the -valued scattering data of a solution of a
matrix ODE is considered (). Injectivity of the map was established by [Paternain, Salo, Uhlmann; Geom.Funct.Anal. 2012].
A statistical algorithm for the solution of this inverse problem based on
Gaussian process priors is proposed, and it is shown how it can be implemented
by infinite-dimensional MCMC methods. It is further shown that as the number
of measurements of point-evaluations of increases, the statistical
error in the recovery of converges to zero in -distance at a
rate that is algebraic in , and approaches for smooth matrix
fields . The proof relies, among other things, on a new stability
estimate for the inverse map .
Key applications of our results are discussed in the case to
polarimetric neutron tomography, see [Desai et al., Nature Sc.Rep. 2018] and
[Hilger et al., Nature Comm. 2018
Uncertainty quantification in medical image synthesis
Machine learning approaches to medical image synthesis have shown
outstanding performance, but often do not convey uncertainty information. In this chapter, we survey uncertainty quantification methods in
medical image synthesis and advocate the use of uncertainty for improving clinicians’ trust in machine learning solutions. First, we describe basic
concepts in uncertainty quantification and discuss its potential benefits in
downstream applications. We then review computational strategies that
facilitate inference, and identify the main technical and clinical challenges.
We provide a first comprehensive review to inform how to quantify, communicate and use uncertainty in medical synthesis applications
- …