Sampling strategies in Bayesian inversion: A study of RTO and Langevin methods

Rémi Laumont; Yiqiu Dong; Martin Skovgaard Andersen

doi:10.3934/ipi.2024046

Article Contents

2025, Volume 19, Issue 3: 592-612. Doi: 10.3934/ipi.2024046

This issue Previous Article Subspace-based coupled tensor decomposition for hyperspectral blind fusion Next Article Diffusion posterior sampling for magnetic resonance imaging

Sampling strategies in Bayesian inversion: A study of RTO and Langevin methods

Department of Applied Mathematics and Computer Science, Technical University of Denmark, DK-2800 Kgs. Lyngby, Denmark

^*Corresponding author: Rémi Laumont
^*Corresponding author: Rémi Laumont

Received: June 2024

Revised: November 2024

Early access: November 2024

Published: June 2025

This work was funded by a Villum Investigator grant (no. 25893) from the Villum Foundation.

Abstract / Introduction Full Text(HTML) Figure(16) / Table(2) Related Papers Cited by

Abstract

This paper studies two classes of sampling methods for the solution of inverse problems, namely Randomize-Then-Optimize (RTO), which is rooted in sensitivity analysis, and Langevin methods, which are rooted in the Bayesian framework. The two classes of methods correspond to different assumptions and yield samples from different target distributions. We highlight the main conceptual and theoretical differences between the two approaches and compare them from a practical point of view by tackling two classical inverse problems in imaging: deblurring and inpainting. We show that the choice of the sampling method has a significant impact on the reconstruction and the proposed uncertainty.

Keywords:

Mathematics Subject Classification: Primary: 65K10, 62F15, 62C10, 68U10, 90C25, 65C05; Secondary: 65K05, 65D18, 68Q25.

Citation:

Full Text(HTML)

Figure 1. Truncated Gaussian posterior, a smooth approximation, and the RTO density for a Gaussian observation model. The vertical dotted line marks the target mean relative to the observation $ y $, which is the mode of the likelihood

Download: Full-size image PowerPoint slide

Figure 2. Original images used for the deblurring and inpainting experiments

Download: Full-size image PowerPoint slide

Figure 3. (Deblurring) degraded images

Download: Full-size image PowerPoint slide

Figure 4. (Deblurring) MMSE estimates computed respectively with RTO and MYULA together with MAP estimates

Download: Full-size image PowerPoint slide

Figure 5. (Deblurring) Marginal posterior standard deviation computed with the samples generated by RTO and MYULA

Download: Full-size image PowerPoint slide

Figure 6. (Deblurring) The comparison of ACFs of RTO and MYULA on both test images

Download: Full-size image PowerPoint slide

Figure 7. (Deblurring) Samples generated by RTO and MYULA

Download: Full-size image PowerPoint slide

Figure 8. (Deblurring) Relative errors of samples generated of RTO and MYULA with respect to MMSE estimates

Download: Full-size image PowerPoint slide

Figure 9. (Deblurring) Results generated by MYULA with $ n_{pd} = 10 $

Download: Full-size image PowerPoint slide

Figure 10. (Deblurring) Results generated by RTO with $ \texttt{tol} = 10^{-2} $

Download: Full-size image PowerPoint slide

Figure 11. (Inpainting) Observations for the inpainting problem

Download: Full-size image PowerPoint slide

Figure 12. (Inpainting) MMSE estimates respectively obtained of RTO and MYULA together with MAP estimates

Download: Full-size image PowerPoint slide

Figure 13. (Inpainting) The standard deviations computed with the samples generated by RTO and MYULA

Download: Full-size image PowerPoint slide

Figure 14. (Inpainting) ACFs from samples generated by RTO and MYULA, respectively

Download: Full-size image PowerPoint slide

Figure 15. (Hierarchical Gibbs sampler) MMSE and marginal posterior standard deviations computed using the hierarchical Gibbs sampler presented in Algorithm 1 for image deblurring problem defined in Section 3.1

Download: Full-size image PowerPoint slide

Figure 16. (Hierarchical Gibbs sampler) Analysis of the sample correlation of the augmented Gibbs sampler. Left: ACFs of the image samples. Middle and right: traces of the noise level and regularization parameters

Download: Full-size image PowerPoint slide

Table 1. Comparisons between RTO and MYULA from a theoretical point of view for a Gaussian likelihood

	RTO	MYULA
Log-prior	polyhedral hypograph	concave
Target density	implicit form	explicit form
Random perturbation	in data space	in image space
Samples	independent	correlated
Burn-in	no	required
Sample generation cost	high: computing (6b)	low: computing (10a)
Parameter selection	online	offline
	hierarchical approach	empirical approach
Parallelization	yes	possible but difficult

| Show Table

DownLoad: CSV

Table 2. (Hierarchical Gibbs sampler) Empirical mean and standard deviation of the noise level and regularization parameters $ \lambda $ and $ \gamma $

	Obtained from hierarchical Gibbs sampler	Used in RTO in Section 3.1
	mean (standard deviation)
$ \gamma $	3.97 (0.44)	5
$ \lambda $	1015.49(7.19)	1000

| Show Table

DownLoad: CSV

Related Papers

Cited by

References

[1]	C. Aguerrebere, A. Almansa, J. Delon, Y. Gousseau and P. Muse, A Bayesian Hyperprior Approach for Joint Image Denoising and Interpolation, With an Application to HDR Imaging, IEEE Transactions on Computational Imaging, 3 (2017), 633-646. doi: 10.1109/TCI.2017.2704439.
[2]	J. M. Bardsley and C. Fox., An MCMC method for uncertainty quantification in nonnegativity constrained inverse problems, Inverse Problems in Science and Engineering, 20 (2012), 477–498. doi: 10.1080/17415977.2011.637208.
[3]	J. M. Bardsley, Computational Uncertainty Quantification for Inverse Problems: An Introduction to Singular Integrals, Comput. Sci. Eng., 19 Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2018. doi: 10.1137/1.9781611975383.
[4]	J. M. Bardsley and P. C. Hansen, MCMC algorithms for computational UQ of nonnegativity constrained linear inverse problems, SIAM Journal on Scientific Computing, 42 (2020), A1269-A1288. doi: 10.1137/18M1234588.
[5]	J. M. Bardsley, A. Solonen, H. Haario and M. Laine, Randomize-then-optimize: A method for sampling from posterior distributions in nonlinear inverse problems, SIAM Journal on Scientific Computing, 36 (2014), A1895-A1910. doi: 10.1137/140964023.
[6]	H. H Bauschke, P. L. Combettes, H. H. Bauschke and P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, Second edition, With a foreword by Hédy Attouch CMS Books Math./Ouvrages Math. SMC Springer, Cham, 2017. doi: 10.1007/978-3-319-48311-5.
[7]	A. Blake, P. Kohli and C. Rother, Markov Random Fields for Vision and Image Processing, MIT Press, Cambridge, MA, 2011.
[8]	S. Boyd, N. Parikh, E. Chu, B. Peleato and J. Eckstein, Distributed optimization and statistical learning via the alternating direction method of multipliers, Foundations and Trends® in Machine Learning, 3 (2011), 1-122. doi: 10.1561/2200000016.
[9]	Z. Cai, J. Tang, S. Mukherjee, J. Li, C. B. Schönlieb and X. Zhang, NF-ULA: Langevin Monte carlo with normalizing flow prior for imaging inverse problems, SIAM J. Imaging Sci., 17 (2024), 820–860. doi: 10.1137/23M1581807.
[10]	A. Chambolle, An algorithm for total variation minimization and applications, Journal of Mathematical Imaging and Vision, 20 (2004), 89-97. doi: 10.1023/B:JMIV.0000011321.19549.88.
[11]	A. Chambolle and T. Pock, A first-order primal-dual algorithm for convex problems with applications to imaging, Journal of Mathematical Imaging and Vision, 40 (2011), 120-145. doi: 10.1007/s10851-010-0251-1.
[12]	T. Chen, E. Fox and C. Guestrin, Stochastic gradient Hamiltonian Monte Carlo, Proceedings of the 31st International Conference on Machine Learning, Proceedings of Machine Learning Research, (2014), 1683-1691.
[13]	V. De Bortoli, A. Durmus, M. Pereyra and A. F. Vidal, Maximum likelihood estimation of regularization parameters in high-dimensional inverse problems: An empirical Bayesian approach. Part II: Theoretical analysis, SIAM Journal on Imaging Sciences, 13 (2020), 1990-2028. doi: 10.1137/20M1339842.
[14]	A. Durmus and É. Moulines, Nonasymptotic convergence analysis for the unadjusted Langevin algorithm, Annals of Applied Probability, 27 (2017), 1551-1587. doi: 10.1214/16-AAP1238.
[15]	A. Durmus, E. Moulines and M. Pereyra, Efficient Bayesian computation by proximal Markov chain Monte Carlo: When Langevin meets Moreau, SIAM Journal on Imaging Sciences, 11 (2018), 473-506. doi: 10.1137/16M1108340.
[16]	J. M. Everink, Y. Dong and M. S. Andersen, Bayesian inference with projected densities, SIAM/ASA Journal on Uncertainty Quantification, 11 (2023), 1025-1043. doi: 10.1137/22M150695X.
[17]	J. M. Everink, Y. Dong and M. S. Andersen, Sparse Bayesian inference with regularized Gaussian distributions, Inverse Problems, 39 (2023), 115004, 28 pp. doi: 10.1088/1361-6420/acf9c5.
[18]	M. Girolami and B. Calderhead, Riemann manifold Langevin and Hamiltonian Monte Carlo methods, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73 (2011), 123-214. doi: 10.1111/j.1467-9868.2010.00765.x.
[19]	P. C. Hansen, J. Jørgensen and W. R. B. Lionheart, Computed Tomography: Algorithms, Insight, and Just Enough Theory, Fundam. Algorithms, 18. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2021. doi: 10.1137/1.9781611976670.
[20]	J. Ho, A. Jain and P. Abbeel, Denoising diffusion probabilistic models, Advances in Neural Information Processing Systems, 33 (2020), 6840-6851.
[21]	M. Holden, M. Pereyra and K. C. Zygalakis, Bayesian imaging with data-driven priors encoded by neural networks, SIAM Journal on Imaging Sciences, 15 (2022), 892-924. doi: 10.1137/21M1406313.
[22]	A. Houdard, C. Bouveyron and J. Delon, High-dimensional mixture models for unsupervised image denoising (HDMI), SIAM Journal on Imaging Sciences, 11 (2018), 2815-2846. doi: 10.1137/17M1135694.
[23]	S. Hurault, A. Leclaire and N. Papadakis, Proximal denoiser for convergent plug-and-play optimization with nonconvex regularization, International Conference on Machine Learning, (2022), 9483-9505.
[24]	P. E. Jacob, J. O'Leary and Y. F Atchadé, Unbiased markov chain monte carlo methods with couplings, Journal of the Royal Statistical Society Series B: Statistical Methodology, 82 (2020), 543-600. doi: 10.1111/rssb.12336.
[25]	D. P Kingma, M. Welling, et al., An introduction to variational autoencoders, Foundations and Trends® in Machine Learning, 12 (2019), 307-392. doi: 10.1561/2200000056.
[26]	R. Laumont, V. De Bortoli, A. Almansa, J. Delon, A. Durmus and M. Pereyra, Bayesian imaging using plug & play priors: When Langevin meets Tweedie, SIAM Journal on Imaging Sciences, 15 (2022), 701-737. doi: 10.1137/21M1406349.
[27]	C. Louchet and L. Moisan, Posterior expectation of the total variation model: Properties and experiments, SIAM Journal on Imaging Sciences, 6 (2013), 2640-2684. doi: 10.1137/120902276.
[28]	R. Molina, A. K. Katsaggelos and J. Mateos, Bayesian and regularization methods for hyperparameter estimation in image restoration, IEEE Transactions on Image Processing, 8 (1999), 231-246. doi: 10.1109/83.743857.
[29]	S. Mukherjee, S. Dittmer, Z. Shumaylov, S. Lunz, O. Öktem and C.-B. Schönlieb, Data-driven convex regularizers for inverse problems, ICASSP 2024-2024 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), (2024), 13386-13390. doi: 10.1109/ICASSP48485.2024.10447719.
[30]	G. Papandreou and A. L. Yuille, Gaussian sampling by local perturbations, Advances in Neural Information Processing Systems, 23 (2010).
[31]	M. Pereyra, Proximal Markov chain Monte Carlo algorithms, Statistics and Computing, 26 (2016), 745-760. doi: 10.1007/s11222-015-9567-4.
[32]	M. Pereyra, J. M. Bioucas-Dias and M. AT Figueiredo, Maximum-a-posteriori estimation with unknown regularisation parameters, 2015 23rd European Signal Processing Conference (EUSIPCO), (2015), 230-234. doi: 10.1109/EUSIPCO.2015.7362379.
[33]	M. Pereyra, L. V. Mieles and K. C. Zygalakis, Accelerating proximal Markov chain Monte Carlo by using an explicit stabilized method, SIAM Journal on Imaging Sciences, 13 (2020), 905-935. doi: 10.1137/19M1283719.
[34]	M. Pereyra, L. A. Vargas-Mieles and K. C. Zygalakis, The split Gibbs sampler revisited: Improvements to its algorithmic structure and augmented target distribution, SIAM Journal on Imaging Sciences, 16 (2023), 2040-2071. doi: 10.1137/22M1506122.
[35]	J.-C. Pesquet, A. Repetti, M. Terris and Y. Wiaux, Learning maximally monotone operators for image recovery, SIAM Journal on Imaging Sciences, 14 (2021), 1206-1237. doi: 10.1137/20M1387961.
[36]	A. Repetti, M. Pereyra and Y. Wiaux, Scalable Bayesian Uncertainty Quantification in Imaging Inverse Problems via Convex Optimization, SIAM Journal on Imaging Sciences, 12 (2019), 87-118. doi: 10.1137/18M1173629.
[37]	D. Rezende and S. Mohamed, Variational inference with normalizing flows, International Conference on Machine Learning, (2015), 1530-1538.
[38]	C. P. Robert and G. Casella, Monte Carlo Statistical Methods, Springer Verlag, 2004.
[39]	G. O. Roberts and R. L. Tweedie, Exponential convergence of Langevin distributions and their discrete approximations, Bernoulli, 2 (1996), 341-363. doi: 10.2307/3318418.
[40]	L. I. Rudin, S. Osher and E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D: Nonlinear Phenomena, 60 (1992), 259-268. doi: 10.1016/0167-2789(92)90242-F.
[41]	A. M. Teodoro, J. M. Bioucas-Dias and M. A. T. Figueiredo, Scene-adapted plug-and-play algorithm with guaranteed convergence: Applications to data fusion in imaging, (2018).
[42]	A. F. Vidal, V. D. Bortoli, M. Pereyra and A. Durmus, Maximum likelihood estimation of regularization parameters in high-dimensional inverse problems: An empirical Bayesian approach part I: Methodology and experiments, SIAM Journal on Imaging Sciences, 13 (2020), 1945-1989. doi: 10.1137/20M1339829.
[43]	M. Vono, N. Dobigeon and P. Chainais, Split-and-augmented Gibbs sampler—application to large-scale inference problems, IEEE Transactions on Signal Processing, 67 (2019), 1648-1661. doi: 10.1109/TSP.2019.2894825.
[44]	Z. Wang, J. M. Bardsley, A. Solonen, T. Cui and Y. M. Marzouk, Bayesian inverse problems with $l_1$ priors: A randomize-then-optimize approach, SIAM Journal on Scientific Computing, 39 (2017), S140-S166. doi: 10.1137/16M1080938.
[45]	Z. Wang and A. C. Bovik, Mean squared error: Love it or leave it? A new look at signal fidelity measures, IEEE Signal Processing Magazine, 26 (2009), 98-117. doi: 10.1109/MSP.2008.930649.
[46]	Z. Wang, A. C. Bovik, H. R. Sheikh and E. P. Simoncelli, Image quality assessment: From error visibility to structural similarity, IEEE Transactions on Image Processing, 13 (2004), 600-612. doi: 10.1109/TIP.2003.819861.
[47]	M. Welling and Y. W. Teh, Bayesian learning via stochastic gradient Langevin dynamics, Proceedings of the 28th International Conference on Machine learning (ICML-11), (2011), 681-688.
[48]	G. Yu, G. Sapiro and S. Mallat, Solving inverse problems with piecewise linear estimators: from Gaussian mixture models to structured sparsity, IEEE Transactions on Image Processing, 21 (2011), 2481-2499. doi: 10.1109/TIP.2011.2176743.
[49]	D. Zoran and Y. Weiss, From learning models of natural image patches to whole image restoration, 2011 International Conference on Computer Vision, (2011), 479-486. doi: 10.1109/ICCV.2011.6126278.