Bayesian imaging using Plug & Play priors: when Langevin meets Tweedie

Since the seminal work of Venkatakrishnan et al. [89] 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) PnPULA (Plug & Play Unadjusted Langevin Algorithm) for Monte Carlo sampling and MMSE inference; and 2) PnP-SGD (Plug & Play 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. ∗VDB was partially supported by EPSRC grant EP/R034710/1. RL was partially supported by grants from Région Ile-De-France. AD acknowledges support of the Lagrange Mathematical and Computing Research Center. MP was partially supported by EPSRC grant EP/T007346/1. JD and AA acknowledge support from the French Research Agency through the PostProdLEAP project (ANR-19-CE23-0027-01). Computer experiments for this work ran on a Titan Xp GPU donated by NVIDIA, as well as on HPC resources from GENCI-IDRIS (Grant 2020-AD011011641). †These authors contributed equally ‡Université de Paris, MAP5 UMR 8145, F-75006 Paris, France §Department of Statistics University of Oxford 24-29 St Giles OX1 3LB, Oxford United Kingdom ¶Institut Universitaire de France (IUF) ‖Centre Borelli, UMR 9010, École Normale Supérieure Paris-Saclay ∗∗School of Mathematical and Computer Sciences, Heriot-Watt University & Maxwell Institute for Mathematical Sciences, Edinburgh, United Kingdom 1 ar X iv :2 10 3. 04 71 5v 4 [ st at .M E ] 1 9 M ar 2 02 1