Picard Proximal Monte Carlo for Parallel Bayesian Imaging with Score-Based Generative Priors

Deliang Wei; Evan Bell; Wenhan Guo; Yifan Chen; Yu Sun

cs.LGAug 18, 2026

Picard Proximal Monte Carlo for Parallel Bayesian Imaging with Score-Based Generative Priors

Authors: Deliang Wei, Evan Bell, Wenhan Guo, Yifan Chen, Yu Sun

Organizations: Johns Hopkins University · University of California, Los Angeles

Abstract

Bayesian imaging inverse problems often require sampling from high-dimensional posterior distributions. While recent score-based and diffusion models provide expressive Bayesian priors, their sampling procedures remain inherently sequential and computationally expensive for large-scale imaging applications. We propose PiX-MC, a time-parallel posterior sampling framework based on proximal Langevin dynamics and Picard iteration. The proximal-likelihood formulation exploits the fact that many imaging likelihoods admit efficient, problem-specific proximal operators, while Picard refinement exposes parallelism across discretization nodes and naturally supports multi-GPU implementation. To further improve practical scalability and sampling performance, we develop multi-block and annealed variants of the proposed framework. We establish convergence guarantees under transparent assumptions, accommodating non-log-concave posteriors, imperfect learned score models, multi-block implementations, and annealing schedules. Experiments on a diverse collection of imaging inverse problems demonstrate that PiX-MC substantially reduces wall-clock time while preserving reconstruction quality. On a $512\times512\times80$ sparse-view computed tomography (CT) problem, annealed multi-block PiX-MC achieves up to a $50\times$ runtime speedup over the standard Langevin sampler using eight GPUs.

Figures & tables

Figure 1: Schematic illustration of the proposed Picard Proximal Monte Carlo (PiX-MC) method. PiX-MC employs time-batch processing and exploits Picard iteration to enable parallel computation across time nodes.

Figure 2: Distribution-level comparison of PiX-MC and APiX-MC against the closed-form Gaussian posterior. Each method runs 150 independent chains, with samples recorded every 10 time steps. The left panels show the posterior mean and pixel-wise standard deviation at the final time step, with the exact posterior shown for reference. The right panels plot the squared MMD to the exact posterior over the sampling trajectory. Note how the analytical-score variant achieves the lowest $\mathsf{MMD}^{2}$ , while APiX-MC reduces $\mathsf{MMD}^{2}$ more rapidly than PiX-MC with the learned score.

Figure 3: Visualization of the in-block Picard error $\varepsilon_{m,k}$ for PiX-MC methods on accelerated MRI reconstruction and Rician denoising. Langevin-based methods are included for reference. The error is plotted against the Picard iteration $k$ for representative blocks, using the block indexed by $m=5$ for MRI reconstruction and the block indexed by $m=3$ for Rician denoising. After a short initial transient, all PiX-MC methods exhibit geometric decay, consistent with Lemma 3 .

Figure 4: MSE evolution of PiX-MC and APiX-MC over time steps and Picard iterations for accelerated MRI reconstruction ( left ) and Rician denoising ( right ). At

k=0

, all time nodes within each block share the same initialization, resulting in a piecewise-constant MSE profile. As the Picard iteration proceeds, the block-wise trajectories are progressively refined into smooth curves.

Figure 5: Visual comparison for $8\times$ accelerated MRI reconstruction with measurement noise level $\sigma_{{\bm{y}}}=0.05$ . Green boxes highlight the zoomed-in regions, and blue panels show the corresponding residual maps. Note the sharper local structures and fewer residual artifacts recovered by PiX-MC and APiX-MC compared with the baseline methods.

Figure 6: Visual comparison for Rician denoising with measurement noise level $\sigma_{{\bm{y}}}=0.1$ . Purple and red boxes highlight the zoomed-in regions. Note the stronger Rician noise suppression and better preservation of fine markings achieved by PiX-MC and APiX-MC compared with the baseline methods.

Figure 7: Average PSNR against wall-clock time for large-scale deblurring. The left and right panels show the non-annealed and annealed methods, respectively. The horizontal dashed line marks the reference PSNR of $26$ dB, and the vertical dotted lines indicate the corresponding threshold-crossing time in seconds. Note how PiX-MC and APiX-MC substantially accelerate convergence without sacrificing reconstruction quality.

Methods	PiL-MC vs. L-MC (491s)			APiL-MC vs. AL-MC (190s)
	$N=8$	$N=16$	$N=32$	$N=8$	$N=16$	$N=32$
4 GPUs	264s / $1.86\times$	328s / $1.50\times$	550s / $0.89\times$	99s / $1.92\times$	139s / $1.37\times$	250s / $0.76\times$
8 GPUs	171s / $2.87\times$	175s / $2.81\times$	295s / $1.66\times$	63s / $2.99\times$	81s / $2.35\times$	151s / $1.26\times$
Methods	PiX-MC vs. X-MC (365s)			APiX-MC vs. AX-MC (157s)
	$N=8$	$N=16$	$N=32$	$N=8$	$N=16$	$N=32$
4 GPUs	195s / $1.87\times$	230s / $1.59\times$	440s / $0.83\times$	84s / $1.87\times$	120s / $1.31\times$	194s / $0.81\times$

Table 2: Runtime and speedup required to reach an average PSNR of

26

dB under different GPU budgets and numbers of time nodes

N

per block. Each speedup is computed relative to the corresponding sequential method shown in the column header.

Figure 8: Visual comparison of the reconstructions by PiX-MC and baselines for the large-scale deblurring task, where the underlying image is cropped into a size of $1024\times 1024$ pixels. DPS (334 seconds) and DAPS (741 seconds) methods are run until convergence. The other methods are run for 200 seconds. Within this limited runtime budget, the proposed methods better preserve local structures with fewer noise residuals.

Method	PSNR (dB) $\uparrow$	LPIPS $\downarrow$	SSIM $\uparrow$
DPS	26.43	0.2498	0.7142
DAPS	28.17	0.1632	0.7444
L-MC	26.82	0.2605	0.7273
AL-MC	27.92	0.2042	0.7737
X-MC	27.24	0.2686	0.7342
AX-MC	27.98	0.2138	0.7719

Table 3: Average quantitative results for large-scale deblurring on ten

1024\times 1024

test images from DIV2K validatation set.

Figure 9: Average PSNR against wall-clock time for 3D CT reconstruction. The left and right panels show the non-annealed and annealed methods, respectively. The horizontal dashed lines mark the reference PSNR levels of $30.5$ dB and $32.5$ dB, and the vertical dotted lines indicate the corresponding threshold-crossing times. Note how PiX-MC and APiX-MC attain even higher reconstruction quality in substantially less wall-clock time than their sequential counterparts. In particular, APiX-MC achieves the 32.5 dB reference PSNR in only $9$ minutes, delivering approximately a $50\times$ wall-clock speedup over L-MC.

Figure 10: Visual comparison of PiX-MC and baseline methods for large-scale sparse-view 3D CT reconstruction. The 20th slice of a $512\times 512\times 80$ CT volume is shown. Green boxes indicate the zoomed-in regions, together with the corresponding residual maps. The PSNR values and runtime (minutes) after convergence are reported. Note how APiX-MC achieves the best visual quality and the highest PSNR of $33.64$ dB in only $35$ minutes, by combining Picard iteration, annealing, and the proximal-likelihood drift.

Figure 11: Visual comparison of PiX-MC and baseline methods for large-scale sparse-view 3D CT reconstruction under a 35-minute time budget. Note how APiX-MC substantially suppresses the streak artifacts and noise while preserving clear structures.

Appendix figures & tables1 asset

Supplementary material from the paper’s appendix.

Appendix

Inverse problem	Method	$\lambda$	$\gamma$	$\eta$	$\xi$	$\sigma_{\mathop{\mathsf{min}}}$	$\sigma_{\mathop{\mathsf{max}}}$	$M$	$N$	$K$
MRI	PiL-MC	$1.6\times 10^{3}$	$6.25\times 10^{-5}$	–	$1$	$0.05$	$0.05$	$10$	$75$	$100$
	APiL-MC	$1.6\times 10^{3}$	$6.25\times 10^{-5}$	–	$0.98$	$0.05$	$0.2$	$10$	$75$	$100$
	PiX-MC	$3.5\times 10^{3}$	$8.75\times 10^{-5}$	$4.375\times 10^{-4}$	$1$	$0.05$	$0.05$	$10$	$75$	$100$
	APiX-MC	$3.5\times 10^{3}$	$8.75\times 10^{-5}$	$4.375\times 10^{-4}$	$0.98$	$0.05$	$0.2$	$10$	$75$	$100$
Rician	PiL-MC	$2.2$	$1\times 10^{-4}$	–	$1$	$0.1$	$0.1$	$5$	$100$	$50$
	APiL-MC	$2.2$	$1\times 10^{-4}$	–	$0.98$	$0.1$	$0.4$	$5$	$100$	$50$

Appendix

Table 4: Algorithmic hyperparameter settings for the real-world imaging experiments.

Explore similar work

Kinetic Langevin Meets Split Gibbs: Accelerated Posterior Sampling for Imaging Inverse Problems with Diffusion Priors
Oct 7, 2026Dai Hai Nguyen, Duc Dung NguyenBayesian Inverse Problems Markov Chain Monte Carlo
Hokkaido University · Institute of Information Technology, VAST
Deep unfolding of MCMC kernels: scalable, modular & explainable GANs for high-dimensional posterior sampling
Feb 24, 2026Jonathan Spence, Tobías I. Liaudat, Konstantinos Zygalakis +1Posterior Sampling Generative Adversarial Networks
Maxwell Institute for Mathematical Sciences, Bayes Centre, 47 Potterrow, Edinburgh, Scotland, UK · School of Mathematical and Computer Sciences, Heriot-Watt University, Edinburgh, Scotland, UK · IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France +1
Proximal-Based Generative Modeling for Bayesian Inverse Problems
May 13, 2026Boyang Zhang, Zhiguo Wang, Ya-Feng LiuScore Matching Bayesian Inverse Problems
School of Advanced Interdisciplinary Sciences, University of Chinese Academy of Sciences, Beijing, China · Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing, China · School of Mathematics, Sichuan University, Chengdu, China +2

Oct 7, 2026J/K move · Enter open · S save

Kinetic Langevin Meets Split Gibbs: Accelerated Posterior Sampling for Imaging Inverse Problems with Diffusion Priors

Dai Hai Nguyen, Duc Dung Nguyen

Hokkaido University · Institute of Information Technology, VAST

Bayesian Inverse Problems Markov Chain Monte Carlo Diffusion Posterior Sampling Langevin Dynamics Image Reconstruction

Open paper AskPDF arXiv

Split Gibbs sampling (SGS) is a popular framework for posterior sampling in Bayesian imaging inverse problems. It decouples a Gaussian data-fidelity term from a complex prior through an auxiliary variable, so the data variable is updated exactly and only the prior-side conditional is hard to sample. Existing samplers treat this conditional in one of two ways. Plug-and-play SGS runs a multi-step diffusion denoiser at every iteration, which is expensive and lacks non-asymptotic guarantees. Langevin-within-SGS takes cheap overdamped Langevin steps but needs many iterations. We propose RED-KLwSGS, which keeps the exact Gaussian update for the data variable and updates the auxiliary variable with underdamped (kinetic) Langevin diffusions driven by a one-shot denoising score, at the same per-iteration cost as Langevin-within-SGS. We prove non-asymptotic Wasserstein-2 convergence in continuous and discrete time for strongly log-concave priors. We also introduce Joint-RED-KLwSGS, which applies kinetic Langevin diffusions to both variables. Experiments with Denoising diffusion probabilistic models as diffusion priors on FFHQ and ImageNet datasets show faster convergence and high-quality image reconstruction.

Abstract

Explore similar work

Kinetic Langevin Meets Split Gibbs: Accelerated Posterior Sampling for Imaging Inverse Problems with Diffusion Priors

Deep unfolding of MCMC kernels: scalable, modular & explainable GANs for high-dimensional posterior sampling

Proximal-Based Generative Modeling for Bayesian Inverse Problems