PaperlayerView exact revision on arXiv

Source-linked AI summary

Efficient Bayesian computation by proximal Markov chain Monte Carlo: when Langevin meets Moreau

Alain Durmus, Eric Moulines, Marcelo Pereyra

arXiv:1612.07471v1stat.COstat.ME

TL;DR

High-dimensional Bayesian imaging needs methods beyond efficient convex optimisation when uncertainty quantification, hypothesis testing, or model selection is required. The paper develops a regularised Langevin MCMC method using Moreau-Yoshida envelopes and proximal operators, with convergence analysis and imaging experiments. The resulting method is designed for log-concave, non-smooth models and can be applied to models handled by proximal convex optimisation.

  • Problem

    Convex optimisation efficiently provides point estimates for high-dimensional imaging but does not support key analyses such as uncertainty quantification, hypothesis testing, and model selection.

  • Method

    The paper regularises non-smooth log-concave targets with a Moreau-Yoshida envelope and applies unadjusted Langevin MCMC using proximal optimisation components.

  • Results

    ULA offers non-asymptotic total-variation guarantees with explicit dependence on step size and dimension, while the proposed regularised target is log-concave, Lipschitz continuously differentiable, and arbitrarily close to the original target.

  • Takeaways & Limitations

    The methodology extends efficient proximal convex-optimisation workflows to Bayesian computation for high-dimensional, non-smooth imaging models.

  • Takeaways & Limitations

    The method targets convex, log-concave posterior models and its samples are initially distributed differently from π, requiring λ and γ to be chosen appropriately for closeness.

Abstract

from arXiv · show

Modern imaging methods rely strongly on Bayesian inference techniques to solve challenging imaging problems. Currently, the predominant Bayesian computation approach is convex optimisation, 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 non-smooth, a class of models that is central in imaging sciences. The methodology is based on a regularised unadjusted Langevin algorithm that exploits tools from convex analysis, namely Moreau-Yoshida envelopes and proximal operators, to construct Markov chains with favourable 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 optimisation algorithms. We provide a detailed theoretical analysis of the proposed methodology, including asymptotic and non-asymptotic 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.

1. Introduction

Bayesian imaging balances efficient convex optimisation for point estimation against MCMC methods needed for uncertainty quantification, hypothesis testing, and model selection. The paper proposes a proximal MCMC methodology designed for high-dimensional, non-smooth models, with convergence guarantees and imaging applications.

  • Image estimation problems span denoising, deconvolution, compressive sensing, super-resolution, tomography, inpainting, source separation, fusion, and phase retrieval.
  • Bayesian imaging commonly uses log-concave posteriors whose MAP estimates can be computed efficiently with high-dimensional convex optimisation.Convex optimisation is theoretically well understood and broadly applicable, but does not provide all Bayesian analyses.
  • MCMC methods address hierarchical and empirical Bayesian models and enable hypothesis testing and model selection beyond optimisation-based analyses.
  • Convex optimisation and MCMC have complementary strengths, motivating workflows that combine efficient full-dataset analysis with targeted in-depth simulation.
  • The proposed proximal MCMC methodology applies to many models solved by convex optimisation, including models addressed by forward-backward splitting.The paper also provides simple convergence conditions and convergence-rate bounds.
  • Four experiments study image deconvolution and tomographic reconstruction with total-variation and ℓ1 sparse priors for model comparison and uncertainty quantification.

2. Bayesian analysis and computation

This section formulates convex Bayesian imaging inverse problems and explains why standard Langevin methods struggle with high-dimensional, non-smooth posteriors. It presents ULA as an efficient smooth-target method while motivating proximal regularisation for the imaging setting.

  • 2.2. Imaging inverse problems.: The posterior π represents knowledge about an unknown image after combining observations with prior information in a convex inverse problem.
  • 2.2. Imaging inverse problems.: The model class uses U = f + g, with smooth convex f and proper, convex, lower-semicontinuous g, covering non-smooth imaging regularisers and constraints.
  • 2.2. Imaging inverse problems.: MAP estimation is often computationally efficient in large problems through proximal convex optimisation, whereas high-dimensional integration with respect to π is more expensive.
  • 2.2. Imaging inverse problems.: Optimisation-based imaging generally cannot assess solution uncertainty or intrinsically compare alternative models without ground truth.
  • 2.3. Bayesian computation: unadjusted and Metropolis-adjusted Langevin algorithms.: ULA discretises overdamped Langevin diffusions using a step size γ and independent standard Gaussian variables to generate a Markov chain.
  • 2.3. Bayesian computation: unadjusted and Metropolis-adjusted Langevin algorithms.: ULA has non-asymptotic total-variation bounds with explicit dependence on step size γ and dimension d, and its efficiency deteriorates at most linearly with d under strong convexity.
  • 2.3. Bayesian computation: unadjusted and Metropolis-adjusted Langevin algorithms.: ULA and MALA are not well defined for non-smooth targets, while insufficient regularity can make ULA explosive and MALA non-geometrically ergodic.

3. Proximal MCMC: Moreau-Yosida regularised Unadjusted Langevin Algorithm

The method smooths non-smooth log-concave targets with a Moreau-Yosida approximation, then applies unadjusted Langevin sampling using proximal operators. The approximation preserves useful regularity while remaining arbitrarily close to the original target, and the resulting MYULA chains admit convergence guarantees and practical bias-control mechanisms.

  • Regularisation: Moreau-Yosida smoothing replaces the non-smooth potential with Uλ, yielding a stable Langevin approximation while allowing πλ to approach π as λ decreases.The regularisation parameter controls the trade-off between smoothness and approximation error.
  • Regularisation: The envelope gλ preserves convexity, is continuously differentiable, and is gradient-Lipschitz even when g is non-smooth.Its gradient-Lipschitz constant is controlled by λ, enabling standard Langevin analysis.
  • Surrogate target: The surrogate density πλ is proper, log-concave, continuously differentiable, and converges to π in total variation as λ ↓ 0.These properties follow under the stated H1 and H2 assumptions.
  • Illustrations: For Laplace and uniform examples, the smooth approximations converge toward the target as λ decreases, with the Laplace total variation error vanishing quadratically in λ.The stated general linear bound does not apply to the uniform density.
  • MYULA: MYULA applies Euler-Maruyama discretisation to the Langevin diffusion targeting πλ, using the proximal operator of g in its update.The proximal operator is already widely used and efficiently computed in imaging optimisation methods.
  • Bias correction: The MYULA stationary distribution differs from πλ because of discretisation, while importance sampling or Metropolis-Hastings correction can address associated bias.The paper focuses on MYULA without these corrections; MYMALA is described as a separate method under development.
  • Convergence analysis: Theoretical analysis establishes geometric convergence to an approximation controlled by λ and γ, non-asymptotic finite-iteration error bounds, and explicit dimension dependence.The analysis also provides practical guidance for selecting λ and γ.

H 3. There exist a minimizer x⋆of Uλ, ηc > 0 and Rc ≥0 such that for all x ∈Rd, ∥x −x⋆∥≥Rc,

The analysis establishes geometric convergence and finite-iteration error bounds for MYULA under tractable assumptions. It also gives iteration-complexity results and practical guidance for balancing regularization and discretization bias against sampling efficiency.

  • Theoretical guarantees: Under H1 and H4, the Euler-Maruyama Markov kernel satisfies geometric ergodicity under suitable step-size conditions.The kernel is irreducible, strongly aperiodic, and satisfies a Foster-Lyapunov drift condition under the stated assumptions.
  • Theoretical guarantees: The unadjusted Langevin algorithm generates samples close to πλ when γ is chosen sufficiently small, although its invariant distribution differs from πλ.The approximation error decreases with the step size, while πλ itself approximates the target π through regularization.
  • Theoretical guarantees: MYULA admits a non-asymptotic total-variation bound between π and the marginal laws of its finite-iteration samples.The bound explicitly accounts for the regularization and discretization parameters and supports finite-sample analysis.
  • Iteration complexity: The worst-case iteration count to reach precision ε is of order d^5 log^2(ε^-1) ε^-2 for the general model class.More precise bounds are available when Uλ is strongly convex outside a bounded region.
  • Parameter selection: The step size γ must lie in (0, λ/(Lfλ+1)] for Euler-Maruyama stability, creating a bias-variance trade-off.Larger γ accelerates the chain and lowers estimation variance but can increase asymptotic bias; smaller γ has the opposite effect.
  • Parameter selection: The experiments use λ = 1/Lf and γ ∈ [L^-1f/4] and obtain estimation errors of approximately 1%.The authors recommend relatively large γ values for efficient computation in high-dimensional imaging settings.

4. Experimental results

The experiments apply MYULA to image deconvolution, tomography, and microscopy for model selection and posterior uncertainty analyses beyond MAP estimation. Across these tasks, MYULA agrees closely with Px-MALA while offering substantially higher computational efficiency.

  • Image deconvolution with total-variation prior: MYULA performs Bayesian model selection without ground truth, assigning posterior probabilities 0.964, 0.036, and <0.001 to M1, M2, and M3, respectively.These probabilities agree with PSNR-based rankings and provide evidence favouring M1.
  • Image deconvolution with total-variation prior: MYULA reaches the posterior typical set in around 10^2 iterations, compared with 10^4 for Px-MALA.The comparison uses chains initialized from a common starting point under model M1.
  • Image deconvolution with wavelet frame: In the wavelet-frame experiment, M2 receives the highest posterior probability, 0.42, ahead of M1 at 0.32 and M3 at 0.26.The ground-truth-free result agrees with PSNR-based evaluation and identifies M2 as the most appropriate model for the data.
  • Computational comparisons: MYULA is approximately one to two orders of magnitude more efficient than Px-MALA in the reported deconvolution and tomography efficiency analyses.The wavelet-frame comparison reports an order-of-magnitude advantage per iteration, while the tomography comparison reports two orders of magnitude in integrated autocorrelation time.
  • Tomographic image reconstruction: The tomography credible region contains an alternative image lacking the structure of interest, so the data cannot confidently establish that structure’s presence.The counterexample has U(x†) = 1.27 × 10^4 and lies in the 90% credible region because η0.10 = 2.34 × 10^4.
  • Microscopy experiment: Microscopy uncertainty quantification estimates 99% positional uncertainty of ±78nm vertically and ±125nm horizontally, close to an independently reported average precision of about 80nm.MYULA’s threshold approximation error is reported as approximately 0.1% relative to Px-MALA.

5. Conclusion

The paper introduces a proximal MCMC methodology for Bayesian computation in high-dimensional, log-concave, non-smooth imaging models. It combines Moreau-Yoshida regularisation with unadjusted Langevin sampling and provides convergence guarantees with explicit dimension dependence.

  • The methodology targets convex, non-smooth imaging inverse problems that are mainly addressed using convex optimisation.
  • Moreau-Yoshida regularisation produces a log-concave, Lipschitz continuously differentiable approximation suitable for unadjusted Langevin MCMC.
  • The analysis establishes asymptotic and non-asymptotic convergence results with explicit convergence-rate dependence on model dimension.

Appendix A. Proof of Proposition 1

The appendix proves integrability and coercivity properties for convex functions used in the methodology, then derives corresponding properties for the Moreau-Yoshida regularisation. The proof proceeds through convex geometry, growth bounds, and monotone convergence.

  • The proof begins by establishing that a lower semicontinuous convex function is finite on a non-empty open subset of R^d.It selects d + 1 points with linearly independent differences and uses their convex hull, which has non-empty interior.
  • Convexity bounds the function on the selected convex hull by M_co, while lower boundedness makes M_co finite.
  • The sublevel set {g ≤ M_co + 1} is shown to be bounded by contradiction using the volume growth of convex hulls containing points arbitrarily far from v_0.
  • The bounded sublevel-set argument yields a linear growth condition outside a ball centered at x_g.
  • Under the proposition's assumptions, the Moreau-Yoshida envelope preserves the required integrability and convergence properties, with the limit obtained through monotone convergence.

Appendix B. Model selection using improper priors

The appendix explains how to perform Bayesian model selection when alternative models share the same improper prior. Shared prior structure allows marginal posterior probabilities to be defined despite the individual joint densities being undefined.

  • Using improper priors in model selection can make each model's joint density undefined.
  • When alternative models share the same improper prior, marginal posterior probabilities can be defined from their likelihood functions under finite integral conditions.

Appendix C. Truncated harmonic mean estimator

The supplied appendix passage introduces a positive joint density on R^d × R^m as the starting point for the truncated harmonic mean estimator.

  • The appendix begins with a positive probability density p defined on R^d × R^m.

C.1. Case of proper prior distributions.

For proper priors, the section develops Monte Carlo estimators for normalizing quantities and ratios using ergodic chains targeting the conditional distribution.

  • The truncated harmonic mean estimator uses an ergodic Markov chain targeting p(x|y) and converges almost surely to I_A(f, y).The estimator is defined for bounded Borel sets A.
  • For two positive distributions, the ratio p1(y)/p2(y) is estimated from their unnormalized versions through the corresponding normalizing-quantity identity.The construction introduces f1 and f2 and then estimates the ratio using the estimator in (34).
  • When the normalizing constant is unknown, MYULA or Px-MALA can compute it, and the computation may be performed offline when the relevant ratio is independent of y.This applies to the priors considered in the referenced experiment.

C.2. Case of improper prior distributions.

For improper priors, the section defines conditional distributions from an improper joint density and estimates normalizing quantities through truncated harmonic means.

  • With an improper prior, f acts as an improper joint density, while p(x|y) is defined by normalizing f over x.The resulting conditional distribution is used for the estimator construction.
  • The normalizing quantity satisfies ˜p(y) = Vol(A)/I_A(f, y), linking it to the truncated harmonic mean estimator.This identity holds for bounded Borel sets A under the stated construction.
Loading 1612.07471v1…