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

Hopkins Computational Imaging Group

Pre-print

Picard Proximal Monte Carlo for Parallel Bayesian Imaging with Score-Based Generative Priors Deliang Wei1 , Evan Bell1 , Wenhan Guo1 , Yifan Chen2 , and Yu Sun1, 1 Johns Hopkins University

2 University of California, Los Angeles

Corresponding author: ysun214@jh.edu

arXiv:2608.17666v1 [cs.LG] 18 Aug 2026

Abstract

1

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 timeparallel 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-logconcave 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 × 512 × 80 sparse-view computed tomography (CT) problem, annealed multi-block PiX-MC achieves up to a 50× runtime speedup over the standard Langevin sampler using eight GPUs.

Introduction

The recovery of an unknown image x ∈ Rd from a set of sparse and noisy measurements y ∈ Rc is a fundamental problem in computational imaging. It is often mathematically formulated as an inverse problem of the following system y = A(x) + e, (1) where the forward operator A : Rd → Rc models the imaging system response and e ∈ Rc denotes measurement noise. Due to the sparsity and noisiness of y, multiple plausible x can explain the same set of measurements, leading to significant uncertainty in the recovered image. Bayesian sampling methods offer a principled framework for addressing this ambiguity. Rather than seeking a single point estimate, they aim to characterize the full posterior distribution π(x|y) given an image prior p(x). This capability is particularly beneficial in high-stakes biomedical and scientific imaging applications, where different imaging results can lead to different scientific conclusions, clinical decisions, or downstream actions [7, 43]. Traditionally, the image prior is handcrafted, including total variation and sparsity-promoting penalties [72, 75, 32]; yet such priors often fall short in capturing the complexity of real-world image distributions. Recently, score-based diffusion models (SDMs) have emerged as a powerful class of generative models capable of modeling such complexity, inspiring a growing body of work on SDM-based posterior sampling for imaging inverse problems. One popular line of work incorporates 1

PiX-MC

Wei et al.

Picard Iteration

measurements as guidance to steer the generative process toward the posterior distribution [37, 16]; this approach is convenient to implement on top of pre-trained models and has been demonstrated across a wide range of imaging applications [17, 79, 35]. However, such guidance-based methods generally lack theoretical guarantees on whether the generated samples indeed follow the true posterior, significantly limiting their interpretability and reliability in high-stakes settings. In light of this, another class of methods integrates SDMs with classical Markov chain Monte Carlo (MCMC) frameworks, providing more rigorous posterior sampling guarantees [10]. Popular instantiations include combining Langevin dynamics with the score network of an SDM [40, 86], and incorporating a full SDM as an empirical proximal sampler in the split Gibbs process [19, 92]. Despite their algorithmic differences, these methods share a common sequential structure: samples must be evolved through long iterative trajectories, each requiring repeated evaluations of both the score model and the forward operator. As a result, the computational cost scales poorly with the problem dimension and the complexity of the forward model, making large-scale imaging applications particularly challenging. This limitation is increasingly pronounced as modern machine learning systems are routinely deployed on multiple graphics processing unit (multi-GPU) platforms that offer substantial parallel computational resources; yet Time axis existing samplers are not designed N to exploit such parallelism. While Langevin MC efforts have explored parallelism for (sequential processing) solvers applicable to inverse probPicard Proximal MC lems [53, 99, 85], they are largely re- (batch-time processing) stricted to deterministic optimization Time axis N frameworks, limiting the direct applicability of these techniques to posteParallel Computed rior sampling. In this work, we bridge this gap Picard Update by proposing Picard proXimal Monte Carlo (PiX-MC), a time-parallel posteK : GPU #1 : GPU #2 : GPU #2 rior sampling framework for scalable Bayesian imaging. PiX-MC is built upon a forward-backward splitting of Figure 1: Schematic illustration of the proposed Picard Proximal Langevin dynamics that decouples the Monte Carlo (PiX-MC) method. PiX-MC employs time-batch proprior and likelihood, enabling modu- cessing and exploits Picard iteration to enable parallel computation lar integration of a score-based genera- across time nodes. tive prior and a likelihood update via proximal mapping. To overcome the sequential bottleneck, PiX-MC further leverages Picard iteration to parallelize the sampling trajectory across time steps. A schematic illustration is provided in Figure 1. The key contributions of this work are as follows: • We formulate the sampling trajectory of PiX-MC as a fixed-point problem and solve it via Picard iteration. Rather than parallelizing across image pixels, PiX-MC introduces parallelism across time steps, enabling concurrent processing over multiple GPUs. Unlike existing Picardbased methods [77, 3], PiX-MC incorporates a forward–backward formulation that imposes measurement consistency through a proximal update. This design is particularly attractive for computational imaging, as it naturally exploits the structure of many imaging forward operators. We further develop annealed and multi-block variants of PiX-MC. The former employs 2

PiX-MC

Wei et al.

temperature annealing to accelerate the underlying sequential sampling trajectory, while the latter partitions the sampling trajectory into multiple blocks, allowing PiX-MC to still operate efficiently under limited GPU budgets. • We present a comprehensive theoretical analysis of the proposed suite of PiX-MC algorithms. Our analysis bounds the time-averaged relative Fisher information by an O(1/T ) transient term together with explicit approximation-error terms. The theory accommodates imperfect learned score networks and weakly convex imaging likelihoods, covering practical settings encountered in a wide range of imaging problems. Furthermore, our results provide a detailed error decomposition that explicitly characterizes the individual contributions of forward–backward splitting, Picard iterations, and time discretization, offering interpretable insight into the sources of sampling error. • We validate PiX-MC on a diverse set of computational imaging inverse problems, including both linear and nonlinear forward models. In addition, we consider large-scale problems ranging from 1024 × 1024 image deblurring to 512 × 512 × 80 sparse-view computed tomography (CT), covering both 2D and 3D imaging settings. These experiments are conducted on multi-GPU platforms to reflect realistic deployment conditions. The results demonstrate that PiX-MC generalizes across these settings while substantially reducing wall-clock time without sacrificing reconstruction quality. On the large-scale 3D CT task, APiX-MC achieves up to a 50× runtime speedup over the standard Langevin sampler using eight GPUs.

2

Background and Related Works

We begin by reviewing regularized image reconstruction from a Bayesian perspective. We then introduce SDMs and discuss posterior sampling methods for imaging inverse problems. Finally, we review related parallelization approaches for image reconstruction.

2.1

Bayesian Image Reconstruction

Classical approaches to imaging inverse problems are formulated as regularized optimization problems, seeking a solution that simultaneously fits the measurements and conforms to prior knowledge about the image. Within the Bayesian framework, this corresponds to maximum a posteriori (MAP) estimation  x̂ = arg max ℓ(y|x) p(x) = arg min L(x) + V (x) , (2) x∈Rd

x∈Rd

where ℓ(y|x) is the likelihood of the measurements, L(x) = −log ℓ(y|x) the data-fidelity term, and V (x) = −log p(x) is the regularizer. A common choice for the data-fidelity is the least-squares loss L(x) = 2β1 2 ∥y − A(x)∥22 , which arises from the Gaussian likelihood y|x ∼ N (A(x), β 2 I) with noise variance β 2 > 0. Popular regularizers include transform-domain sparsity penalties such as V (x) = τ ∥T x∥1 [11, 22], where T is a sparsifying transform and τ > 0 controls the regularization strength. Proximal methods [9] are widely adopted for solving (2) in computational imaging. Their popularity stems from the ability to handle nonsmooth regularizers, such as transform-domain sparsity [6], while exploiting the composite structure of (2). Moreover, proximal algorithms naturally accommodate the structure of imaging forward operators, often enabling efficient implementations when A 3

PiX-MC

Wei et al.

admits fast transforms or closed-form inverses [1]. The fundamental building block is the proximal operator   1 2 proxηF (z) ≜ arg min ∥x − z∥2 + ηF (x) , (3) 2 x∈Rd where the quadratic term encourages the output to remain close to the input z, and η > 0 controls the influence of F . The proximal operator admits an interpretation as a backward or implicit gradient step, as its optimality condition yields x ∈ z − η∂F (x)

x = (I + η∂F )−1 (z),

⇔

where ∂F (x) denotes the subdifferential of F at x, and (I + η∂F )−1 is known as the resolvent of ∂F . A common proximal algorithm is forward-backward splitting (FBS), which combines an explicit gradient step with an implicit proximal step. For instance, given a smooth V and well-structured A, the update of FBS reads as
xn+1 = proxγL xn − γ∇V (xn ) = (I + γ∂L)−1 ◦ (I − γ∇V )(xn ).

(4)

where ◦ denotes function composition, and γ > 0 denotes the stepsize. Note that the prior and likelihood information are incorporated separately through the gradient and proximal steps, respectively. Classic convex optimization theory has established convergence of FBS to a minimizer of (2) under various conditions [61, 5].

2.2

Score-Based Diffusion Models as Learned Image Priors

SDMs offer a powerful framework for learning image priors from data. Rather than modeling the image distribution p directly, SDMs learn the score of a Gaussian-smoothed version of p. Specifically, let x ∼ p be a clean image and let z = x + σζ with ζ ∼ N (0, I) be a noisy observation at noise level σ > 0. Denote by pσ = p ∗ N (0, σ 2 I) the distribution of z. By Tweedie’s formula [26], the score of pσ admits the closed-form expression ∇log pσ (z) = (E[x|z] − z) /σ 2 ,

(5)

where E[x|z] is the minimum mean-squared error (MMSE) estimator of the clean image given the noisy observation. Equation (5) establishes a direct connection between score estimation and Gaussian denoising. This observation motivates denoising score matching (DSM) [88], which learns a neural network approximation of the smoothed score, i.e. Sθ ( · ; σ) ≈ −∇log pσ where θ denotes the weights of the network. A common DSM training objective is h i L(θ) = E ∥Sθ (z ; σ) − ζ/σ∥2 , with z = x + σζ, (6) which can be optimized efficiently using a dataset of clean images. As σ → 0, the smoothed density pσ approaches p, yielding Sθ ( · ; σ) ≈ −∇log p = ∇V . As a result, the learned score network provides a data-driven approximation of the gradient of the prior potential. With the learned score network Sθ (·; σ), one can construct sampling algorithms that generate samples from the learned prior. A popular approach is based on Langevin dynamics [63, 80], which evolves the sample according to the following stochastic differential equation (SDE) √ dxt = −Sθ (xt , σt ) dt + 2 dBt , (7) 4

PiX-MC

Wei et al.

where σt denotes a time-dependent noise level and Bt ∈ Rd is standard Brownian motion. While one could choose σt to be a small constant, empirical studies have shown that gradually decreasing σt toward zero during sampling leads to substantially improved performance [60, 80]. This strategy, commonly known as annealing, allows the sampler to first explore the global structure of the distribution at high noise levels before progressively refining samples at lower noise levels. An alternative family of SDMs formulates sampling through a reverse-time diffusion process [36, 82, 44]. Let pt denote the density of a forward noising process governed by dxt = u(t) xt dt + v(t) dBt , where u(t) and v(t) denote the drift and diffusion coefficients, respectively [44]. The corresponding reverse-time dynamics is given by   dxt = u(t)xt + v(t)2 Sθ (xt , σt ) dt + v(t) dBt , (8) where σt corresponds to the noise level associated with pt , which is predetermined by the forward process. Starting from Gaussian noise, simulating (8) progressively transforms the sample toward the target image distribution. Both Langevin-based and reverse-diffusion formulations have demonstrated remarkable empirical success in generating high-quality images across a wide range of domains, including natural images [81, 20], medical imaging [47, 23], and scientific imaging applications [28, 86]. This strong generative capability makes score-based diffusion models particularly attractive as image priors for solving inverse problems, as discussed in the next section.

2.3

Posterior Sampling Methods for Imaging

Existing posterior sampling methods can largely be viewed as constructing stochastic processes whose stationary or terminal distribution approximates the target posterior π(x|y). These methods can be roughly categorized into three families according to the type of their stochastic process: Langevinbased methods, split Gibbs samplers, and reverse-time diffusion approaches. Langevin-Based Methods The underlying idea of this family is to construct a Langevin diffusion process whose invariant distribution coincides with the posterior. Following (7), the drift is decomposed into a likelihood gradient ∇L and a learned score Sθ , and a forward Euler discretization yields the Langevin Monte Carlo (L-MC) iteration
p xn+1 = xn − γ ∇L(xn ) + Sθ (xn ; σn ) + 2γ Wn , (9) √ where γ > 0 denotes the step size and where Wn = (B(n+1)γ − Bnγ )/ γ ∼ N (0, I) is a d-dimensional standard Gaussian vector. Numerous variants of L-MC have been developed for imaging inverse problems. Recent works have incorporated customized annealing strategies for the likelihood gradient ∇L and the score network Sθ to accelerate empirical convergence [40, 86], explored score-based priors in latent spaces [38, 18, 83], and extended Langevin sampling to non-Gaussian likelihoods with highly nonconvex landscapes [58]. These methods have been successfully applied across a broad range of imaging tasks [40, 18, 86]. On the theoretical side, substantial effort has been devoted to establishing quantitative convergence guarantees for both continuous-time and discretized Langevin dynamics; see, e.g., [31, 24]. While classical analyses often rely on strong log-concavity assumptions, recent advances have significantly broadened the theoretical understanding beyond this regime [4]. Many of these results are formulated through functional inequalities satisfied by the target distribution, including logarithmic Sobolev [87] and Poincaré inequalities [14]. One closely related line of work exploits the composite structure of the posterior potential through proximal updates. When the prior potential V is nonsmooth, [25] proposed to replace the prior score 5

PiX-MC

Wei et al.

with the gradient of the Moreau-Yosida envelope, yielding the following update   p
1 xn+1 = xn − γ ∇L(xn ) + xn − proxλV (xn ) + 2γ Wn , (10) λ
where λ1 xn −proxλV (xn ) is the gradient of the Moreau-Yosida envelope of V . Subsequent works have replaced the proximal operator with pretrained deep denoisers, giving rise to plug-and-play MoreauYosida Langevin samplers [50]. Another line of work considers an asymmetric discretization of the Langevin dynamics, applying a backward step to one term and a forward step to the other [24, 91, 73]. This yields a proximal Langevin update of the form
p xn+1 = xn − γ ∇L(xn+1 ) + ∇V (xn ) + 2γ Wn ,   p = proxγL xn − γ∇V (xn ) + 2γ Wn . (11) which applies a proximal step on the negative log-likelihood L and a gradient step on the negative log-prior V . This forward-backward structure draws a natural connection to FBS in (4). PiX-MC complements this line of research by integrating Picard iteration into the proximal Langevin framework, enabling parallel trajectory evaluation for large-scale imaging inverse problems. Split Gibbs Samplers An alternative approach to posterior sampling is based on variable splitting and Gibbs sampling [89, 65, 8, 27]. These methods introduce an auxiliary variable z and couple it 1 with x through a quadratic penalty 2ρ ∥x − z∥2 , where ρ > 0 controls the coupling strength. Sampling is then performed by alternating between the conditional distributions   1 2 xn+1 ∼ exp −L(x) − ∥x − zn ∥ (12a) 2ρ   1 2 zn+1 ∼ exp −V (z) − ∥z − xn+1 ∥ . (12b) 2ρ The resulting decoupling allows the likelihood and prior terms to be handled independently, making it possible to exploit the structure of the forward model while incorporating sophisticated learned priors in a modular fashion. Recent works have employed diffusion models to approximate the prior conditional distribution in (12b), establishing convergence to the augmented posterior under various conditions [92, 19, 93], with extensions to blind inverse problem settings [39]. Despite these advances, the conditional updates remain computationally demanding, as each Gibbs iteration typically requires solving an additional sampling subproblem. Guided Reverse Diffusion Rather than constructing a Markov chain that directly targets the posterior distribution, guided reverse-diffusion methods modify the reverse diffusion trajectory using measurement information. Starting from a learned generative diffusion process, the reverse dynamics are guided toward measurement-consistent samples through likelihood-based corrections [16] or approximations of the posterior score [37]. This viewpoint underlies a broad class of diffusion posterior sampling and guided diffusion methods, which have demonstrated strong empirical performance across a wide range of imaging tasks [17, 15, 46, 71, 55, 79, 95]. Compared with Langevin-based approaches, however, these methods typically rely on approximations of the time-dependent posterior score and therefore do not always admit a straightforward characterization in terms of a target invariant distribution. 6

PiX-MC

2.4

Wei et al.

Related Parallel Imaging Methods & Picard Iteration

Most existing parallel methods consider a general optimization formulation applicable to the MAP problem in (2). They accelerate computation by exploiting parallelism within individual iterations, including stochastic gradient methods [67, 54], block-coordinate updates [29, 70], and asynchronous procedures [64, 84, 34, 85]. Some of these techniques have also been extended to sampling [73, 74, 21]. However, they do not reduce the intrinsic sequential structure of iterative algorithms, which remains the fundamental computational bottleneck for large-scale posterior sampling. A more recent line of work seeks to parallelize the sampling trajectory itself. In image generation, Picard iteration has been used to evaluate multiple denoising steps simultaneously [77, 78]. Consider a dynamical system dxt = ft (xt ) dt, t ∈ [0, T ], with initial condition x0 . Picard iteration generates a sequence of trajectories according to (0)

xt

≡ x0 , (k+1)

Note that each point xt

(k+1)

xt

Z t = x0 +

fs (x(k) s ) ds,

for k = 0, . . . , K − 1.

(13)

(k)

(k+1)

0

depends only on the previous trajectory (xs )s∈[0,t] , rather than on xs (k+1)

for s < t. This formulation breaks the sequential dependency across time and allows xt to be computed in parallel across time nodes after discretization. On the theoretical side, [3] established parallelization guarantees for Langevin dynamics under functional inequalities such as logarithmic Sobolev, with subsequent works extending these ideas to diffusion models [12], diagonal trajectory updates [98], and higher-order Langevin schemes [56]. Despite these advances, existing trajectoryparallel methods are designed primarily for either unconditional generative models or posterior distributions with explicit unnormalized densities. They do not directly accommodate the structure of Bayesian imaging posteriors, which often combine proximal likelihood updates, learned image priors, and annealing schedules. PiX-MC significantly extends this line of work to posterior sampling for imaging inverse problems, and provides a comprehensive theoretical analysis that accommodates all these components.

3

Method: Picard Proximal Monte Carlo (PiX-MC)

This section presents the proposed PiX-MC framework. We first construct a proximal-likelihood drift motivated by the FBS-style update, which induces a sequential proximal Langevin dynamics. We then reformulate the resulting dynamics using Picard iteration, enabling parallel computation across sequential time steps. We further establish finite-horizon bounds on the time-averaged relative Fisher information of PiX-MC under transparent assumptions. Finally, we introduce multi-block and annealed variants of PiX-MC to further improve scalability and imaging performance. The analysis extends to both variants under their stated variant-specific assumptions.

3.1

Proximal-Likelihood Drift via Forward-Backward Splitting

We consider the FBS update in (4) to construct our drift term. To distinguish the FBS splitting parameter from the Langevin stepsize γ introduced later, we denote the former by η > 0 throughout this section. Replacing ∇V (x) with a pretrained score network Sθ (x; σ) yields
xn+1 = proxηL xn − ηSθ (xn ; σ) , 7

(14)

PiX-MC

Wei et al.

or equivalently,

 xn+1 − xn 1 = − xn − proxηL xn − ηSθ (xn ; σ) . η η

(15)

The incremental representation motivates the following ordinary differential equation (ODE) dxt = −Tη,σ (xt )dt,

where

Tη,σ ≜

 1 I − proxηL ◦ I − ηSθ (·; σ) . η

(16)

When the noise level σ is fixed throughout the iteration, we omit its dependence and simply write Sθ and Tη . The drift Tη combines the learned prior score model Sθ with a proximal likelihood update, enabling efficient incorporation of likelihood terms through their problem-specific proximal operators. We next characterize the discrepancy between Tη and the ideal posterior drift −∇log π(x|y) ≜ ∇L(x) + ∇V (x). To proceed, we impose the following transparent assumptions. Assumption 1. The negative log-likelihood L is αL -weakly convex and βL -smooth. R The negative log-prior V is αV -weakly convex and βV -smooth. In addition, assume 0 < Z := Rd exp[−L(x) − V (x)] dx < ∞ and define π(dx) := Z −1 exp[−L(x) − V (x)] dx. The conditions on the likelihood are standard and hold for many linear and nonlinear imaging models; see [69] for examples. Moreover, if L is αL -weakly convex and 0 < η < 1/αL , then the proximal operator proxηL = (I + η∇L)−1 is well defined, single valued, and 1/(1 − ηαL )-Lipschitz; see Proposition 2 of [30]. We note that the smoothness of the likelihood is commonly required in the analysis of proximal Langevin algorithms [91, 73]. Assumption 2. Given an arbitrary σ > 0, the score network Sθ (x; σ) satisfies: (i) There exists δ > 0 such that, for all x ∈ Rd , we have ∥Sθ (x; σ) − ∇V (x)∥ ≤ δ; (ii) There exists LS > 0 such that for any x1 , x2 ∈ Rd , ∥Sθ (x1 ; σ) − Sθ (x2 ; σ)∥ ≤ LS ∥x1 − x2 ∥.

(17)

Assumption 2(i) allows for imperfectly learned score networks with δ-accuracy to ∇V , while Assumption 2(ii) is widely used for convergence analysis with score networks [94, 51, 13, 86, 69, 68]. Note that, under Assumptions 1–2, the proximal drift Tη is naturally Λη -Lipschitz Λη ≜ LS + βL

1 + ηLS . 1 − ηαL

(18)

See Lemma 4 in Appendix A for details. The following lemma summarizes the key approximation property of the drift Tη used throughout the subsequent analysis. Lemma 1. Let Assumptions 1 and 2(i) hold. Let η be chosen to satisfy 1 − ηαL > 0 and 1 − 3η 2 βL2 > 0. Then, for any fixed σ > 0 and any x ∈ Rd , ∥Tη (x) − ∇L(x) + ∇V (x))∥2 ≤

6η 2 βL2 2δ 2 + ∥∇L(x) + ∇V (x)∥2 . 2 1 − 3η 2 βL 1 − 3η 2 βL2 8

(19)

PiX-MC

Wei et al.

Proof. See proof in Appendix A. Lemma 1 shows that Tη (x) approximates the ideal posterior drift ∇L(x) + ∇V (x), with the approximation error controlled by the score estimation error δ and the splitting error induced by the FBS-type treatment. Adding Brownian motion Bt to the ODE in (16) yields the following continuoustime SDE √ dxt = −Tη (xt )dt + 2dBt . (20) By applying the forward Euler–Maruyama discretization, we obtain the proXimal MC (X-MC) p xn+1 = xn − γTη (xn ) + 2γWn ,

(21)

which serves as the sequential baseline for the proposed PiX-MC, introduced in the next section. It is also worth noting that, in the exact-score case Sθ = ∇V , a change of variables [49] rewrites X-MC in a form similar to the proximal Langevin algorithm [91]; see also (11). Let zn = proxηL (xn − η∇V (xn )). We have √ xn+1 = xn − γTη (xn ) + 2γWn   p γ γ −1 ⇔ (I − η∇V ) ◦ (I + η∇L)(zn+1 ) = 1 − xn + zn + 2γWn (22) η η    p γ γ ⇔ (I + η∇L)(zn+1 ) = (I − η∇V ) ◦ 1 − xn + zn + 2γWn , η η where (I − η∇V )−1 is well-defined by setting η < 1/βV . Letting γ = η, we arrive at   p
p zn+1 = proxγL zn − γ∇V zn + 2γWn + 2γWn .

(23)

This formulation shows that X-MC effectively applies the proximal Langevin update to zn , up to the √ residual noise term 2γWn appearing inside the evaluation of ∇V . From this perspective, X-MC recovers the proximal Langevin algorithm in this special case, apart from this additional perturbation in the gradient.

3.2

PiX-MC: Algorithm Formulation and Convergence Analysis

Starting in this section, we introduce the family of PiX-MC algorithms. Our presentation follows a two-part structure: we first introduce the formulation of each algorithm, and then establish its companion convergence analysis, with additional assumptions and lemmas introduced as needed. This structure will be maintained in subsequent sections. Algorithm formulation. We apply Picard iteration to parallelize the computation of the X-MC trajectory across time nodes. Specifically, the Picard iteration for the continuous-time SDE in (20) is given by Z t √ (k+1) xt = x0 − Tη (x(k) (24) s )ds + 2Bt , 0

where we note that all trajectory iterates share the same realization of the Brownian process (Bt )t∈[0,T ] . Let tn = nγ for n = 0, . . . , N , where N γ = T . Discretizing (24) using left-point Riemann sums yields

9

PiX-MC

Wei et al.

Algorithm 1 PiX-MC (single block) Input: Initial distribution µ0 , drift Tη , stepsize γ > 0, number of time steps N , and number of Picard iterations K. Output: xN i.i.d. 1: Sample x0 ∼ µ0 , and independently of x0 , sample W1 , . . . , WN ∼ N (0, I) n P√ 2: Set B0 ← 0, Bnγ ← γWj , for n = 1, . . . , N j=1 (0) (k) 3: Set xn ← x0 , for n = 0, . . . , N , and x0 ← x0 for k = 0, . . . , K

4: for k = 0, . . . , K − 1 do 5: 6: 7:

(k)

(k)

(k)

(k)

Tη (x0 ), . . . , Tη (xN −1 ) ← ParallelCompute(x0 , . . . , xN −1 ) for n = 1, . . . , N do n−1 √ P (k) (k+1) Tη (xi ) + 2Bnγ xn ← x0 − γ i=0

8:

end for 9: end for (K) 10: Return xN = xN the PiX-MC update x(0) n = x0 ,

(k)

x0 = x0 ,

x(k+1) = x0 − γ n

n−1 X

(k)

Tη (xi ) +

√ 2Bnγ

(25)

i=0

with n = 1, . . . , N, k = 0, . . . , K. Pn √ Here, Bnγ = i=1 γWi denotes the cumulative sum of the discrete Brownian increments. Unlike (k)

the sequential update rule of X-MC, PiX-MC has the entire previous trajectory {xn }N n=0 available (k+1) N (k) when computing {xn }n=0 . Hence, the drift evaluations Tη (x0 ), . . . , (k) Tη (xN −1 ) can be computed in parallel across time nodes and accumulated for the next trajectory refinement. Suppose sufficient parallel computational resources are available such that each time node can be processed concurrently; the computational cost of PiX-MC then no longer scales with N but with the number of Picard iterations K. This is equivalent to viewing the entire trajectory as a single time block that is fully parallelized. Algorithm 1 summarizes the algorithmic details of PiX-MC under this single-block condition. In §3.3, we will introduce a multi-block version that better balances computational resource constraints against the degree of parallelization. Convergence analysis. The goal of the presented analysis is to quantify the distance between the (k) evolving probability distribution of the PiX-MC iterates xn and the target posterior distribution. To this end, we use the relative Fisher information to measure the distance between two probability distributions. Definition 1. Let µ and π denote two probability distributions such that µ ≪ π. The Fisher information (FI) of µ relative to π is defined as Z dµ 2 FI(µ∥π) = ∇log dµ. (26) dπ 10

PiX-MC

Wei et al.

If µ and π have positive and smooth densities, FI(µ ∥ π) = 0 if and only if µ = π, motivating its use as a criterion for measuring convergence between distributions. From an optimization perspective, Langevin dynamics can be viewed as a gradient flow in the space of probability distributions equipped with the Wasserstein metric, with FI playing the role of the squared gradient norm. Thus, relative FI provides a stationarity criterion for the algorithm. Lemma 2. Let µ0 be the initial distribution of x0 and suppose Assumptions 1-2 hold. Set the parameters η, γ > 0 sufficiently small such that 1 − ηαL > 0,

1 − 3η 2 βL2 > 0,

qsing ≜ γN Λη < 1.

Define the error across two successive Picard refinements at the k-th iteration as h i (k−1) 2 εk ≜ max E ∥x(k) ∥ . n − xn n=1,2,...,N

(27) (28)

Then, the Picard iteration error for PiX-MC satisfies (29)

2K−2 εK ≲ qsing 1 + FI(µ0 ∥π) . Proof. See Appendix B for the detailed proof and range of the hyperparameters.

Lemma 2 shows that the Picard iteration error εK decays geometrically with the total number of Picard iterations. As shown by the numerical results in Fig. 3 and §5, this geometric rate leads to rapid stabilization of the trajectories within only a few iterations. We next present the main theorem for PiX-MC. Theorem 1. Let (µt )t≥0 denote the law of the continuous interpolation of {xn (K) }N n=0 generated by PiX-MC, where N > 0 is the total number of time steps. Assume Assumptions 1-2 hold, and the parameters η, γ > 0 are chosen sufficiently small. Then, the following inequality holds for PiX-MC:   1Z T 2 2 0.75 − O(η ) − O(γ ) FI(µt ∥π)dt T 0
1 2K−2 (30) KL(µ0 ∥π) + O(η 2 ) + O(δ 2 ) + O(γ) + O(qsing ) 1 + FI(µ0 ∥π) . ≤ | {z } | {z } | {z } T | {z } Splitting discrepancy

Score mismatch error

Discretization error

Finite Picard refinement error

Proof. See Appendix B for the detailed proof and range of the hyperparameters. Theorem 1 establishes a finite-time stationarity guarantee for PiX-MC. The theorem also explicitly quantifies the errors stemming from the splitting discrepancy O(η 2 ), score mismatch O(δ 2 ), and discretization error O(γ), which are inherited from X-MC. The additional finite Picard-refinement
2K−2 error O(qsing ) 1 + FI(µ0 ∥π) is specific to PiX-MC. When η, γ are chosen sufficiently small such that
the coefficient 0.75 − O(η 2 ) − O(γ 2 ) is positive and qsing < 1, the bound is nontrivial, and the finite Picard refinement error decays geometrically with K according to Lemma 2. Additionally, as FI(·∥π) is convex, Jensen’s inequality yields Z 1 T FI(µ̄T ∥π) ≤ FI(µt ∥π)dt, (31) T 0 RT where µ̄T = (1/T ) 0 µt dt is the time-averaged marginal law. A possible approach for drawing samples from this distribution has been discussed in [4]. Theorem 1 therefore also bounds the stationarity gap of µ̄T relative to π. 11

PiX-MC

3.3

Wei et al.

Multi-Block PiX-MC: Partitioning the Trajectory into Multiple Blocks

As discussed in the previous section, the vanilla PiX-MC treats the entire trajectory as a single block to be fully parallelized, offering the largest possible degree of parallelism. Yet this design comes with two practical trade-offs. First, the degree of parallelism interacts with the choice of hyperparameters. Recall from Lemma 2 and Theorem 1 that controlling the Picard error requires qsing < 1. Since qsing = T Λη , this contraction condition directly restricts the time horizon that can be handled in a single block. Second, refining the full trajectory at once is more demanding on computational resources, as it requires storing and updating all time nodes simultaneously, which can be memory-intensive for large-scale imaging problems. To better balance parallelism with these practical considerations, we introduce multi-block PiX-MC that partitions the trajectory into multiple time blocks. Algorithm formulation. We partition the full interval [0, T ] into M successive blocks, each spanning N time steps, so that T = M N γ, τ = N γ, (32) where τ denotes the length of each block. Starting from the endpoint of the preceding block, multiblock PiX-MC then performs K Picard refinements only within the current block, rather than over the full trajectory. As a result, both the trajectory storage and the concurrent drift evaluations scale with the block size N rather than the total number of time steps M N , allowing N to be chosen according to the available GPU memory and parallel computing resources. This block-wise refinement also relaxes the condition required for convergence: the local Picard contraction factor becomes qmult ≜ τ Λη . Note that when M > 1, qmult = (τ /T )qsing < qsing , the multi-block scheme admits a larger stepsize than its single-block counterpart. Algorithm 2 summarizes the complete procedure of multi-block PiX-MC. As shown, parallelization occurs within each time block, while the algorithm proceeds sequentially across blocks. Convergence analysis. The convergence analysis of multi-block PiX-MC requires an additional stability control compared with the single-block setting. In single-block PiX-MC, every Picard trajectory is initialized at the same point x0 , and the initial Picard residual can be controlled directly through the initial Fisher information. In contrast, the starting point xmN of the m-th block is itself a random variable, determined by the output of all preceding blocks. Therefore, we impose the following stability assumption on the block starting points, which allows the initial Picard residual to be controlled uniformly across blocks. Assumption 3. The block starting points have uniformly bounded second moments. Specifically, there exists a constant Cx > 0, independent of the number of blocks M and the total time horizon T , such that (0) 2 sup E xmN ≤ Cx . (33) 0≤m≤M −1

Assumption 3 is a standard condition for Langevin-based methods, and can be implied by coercivity of the posterior potential or dissipativity of the drift; see e.g., [66, 68]. We also provide an example condition for Assumption 3 to hold in Appendix C. This assumption is consistent with practical imaging settings, where the image intensity typically lies in a bounded range such as [0, 1]d . Under Assumption 3, the block-wise Picard initialization remains uniformly controlled over the sampling horizon. We next establish a result analogous to Lemma 2 for multi-block PiX-MC. 12

PiX-MC

Wei et al.

Algorithm 2 multi-block PiX-MC Input: Initial distribution µ0 , drift Tη , stepsize γ > 0, number of blocks M , time step number N inside each block, and Picard iteration number K. i.i.d. 1: Sample x0 ∼ µ0 and, independently of x0 , sample W1 , . . . , WM N ∼ N (0, I) i √ P γWj for i = 1, . . . , M N . 2: Set B0 ← 0, and Biγ ← j=1

3: for m = 0, . . . , M − 1 do 4: 5: 6: 7: 8:

(0)

(k)

Set xmN +i ← xmN for i = 1, . . . , N , and xmN ← xmN for k = 0, . . . , K for k = 0, . . . , K − 1 do (k) (k) (k) (k) Tη (xmN ), . . . , Tη (x(m+1)N −1 ) ← ParallelCompute(xmN , . . . , x(m+1)N −1 ) for i = 1, . . . , N do i−1 √
P (k+1) (k) xmN +i = xmN − γ Tη (xmN +j ) + 2 B(mN +i)γ − BmN γ j=0

9:

end for 10: end for (K) 11: x(m+1)N ← x(m+1)N 12: end for (K) 13: Return xM N = xM N Lemma 3. Let Assumptions 1–3 hold, and set the parameters η, γ > 0 sufficiently small such that 1 − ηαL > 0,

1 − 3η 2 βL2 > 0,

qmult = γN Λη < 1.

For each block m = 0, . . . , M − 1, define the in-block Picard error as i h (k) (k−1) εm,k ≜ max E ∥xmN +i − xmN +i ∥2 . i=1,...,N

(34)

(35)

2 εm,k . Consequently, Then, for multi-block PiX-MC, the in-block Picard error satisfies εm,k+1 ≤ qmult for each block m = 0, . . . , M − 1,   3τ 2 (δ 2 + G2 ) 2K−2 2K−2 εm,K ≤ qmult 2τ d + = qmult R. (36) 1 − 3η 2 βL2

Here G > 0 is a constant independent of the total time T ; see Appendix D for details. Proof. See Appendix D for the detailed proof and explicit expressions for the constants. Lemma 3 establishes geometric Picard contraction within each local block. We are now ready to present the convergence result for multi-block PiX-MC. N Theorem 2. Let (µt )t≥0 denote the law of the continuous interpolation of {xn (K) }M n=0 generated by multi-block PiX-MC, where N > 0 denotes the number of time steps in each time block and M > 0 is the number of blocks. Assume Assumptions 1-3 hold and choose sufficiently small parameters

13

PiX-MC

Wei et al.

η, γ > 0. Then, the following inequality holds for multi-block PiX-MC:  ≤

 1Z T

0.75 − O(η 2 ) − O(γ 2 )

KL(µ0 ∥π) + O(η 2 ) + | {z } T Splitting discrepancy

T

FI(µt ∥π)dt

0

O(δ 2 ) | {z }

Score mismatch error

+

O(γ) | {z }

Discretization error

2K−2 + O(qmult )R . {z } |

(37)

Finite Picard refinement error

where R is the constant from Lemma 3. Proof. See Appendix D for the detailed proof and range of the hyperparameters. Theorem 2 extends the single-block guarantee to the multi-block case by replacing the global contraction coefficient with the in-block qmult < 1. Note that qmult is determined by the local time horizon τ rather than the total sampling time T , and the finite Picard error still decays geometrically within each block. Compared with the result in Theorem 1, the associated constant is now controlled by R established in Lemma 3, rather than by the initial Fisher information.

3.4

Multi-Block Annealed PiX-MC: Improving Sample Quality via Annealing

PiX-MC has so far used a fixed noise level in the score network Sθ . In this section, we introduce multi-block annealed PiX-MC (multi-block APiX-MC) that employs a decreasing noise schedule to first capture coarse image structures and then progressively refine details. Algorithm formulation. Let n = mN +i denote the global time node index, where m = 0, . . . , M −1 is the block index and i = 0, . . . , N − 1 is the local node index. We use the geometric schedule σn = max{σmax ξ n , σmin },

(38)

N −1 where 0 < ξ < 1 is the annealing parameter. We further associate {σn }M with a non-increasing n=0 M N −1 score-weight schedule {αn }n=0 satisfying αn ≥ 1 and αn → 1 when n → ∞. At the n-th time node, the corresponding weighted proximal-likelihood drift is defined as αn Tη,σ ≜ n

 1 I − proxηL ◦ (I − ηαn Sθ (·; σn )) . η

(39)

Rather than restarting the annealing schedule at each block boundary, the global indexing in (38) allows the noise level and score weight to decrease continuously over the entire sampling trajectory. Within the m-th block, the local Picard iteration thus uses the time-dependent drifts αmN , . . . , T α(m+1)N −1 , while retaining the same block-wise trajectory structure as multi-block PiXTη,σ η,σ(m+1)N −1 mN MC. The resulting algorithm, multi-block APiX-MC, is summarized in Algorithm 3. Convergence analysis. The convergence analysis follows the same overall argument as that of multiblock PiX-MC, with modifications accounting for the time-varying denoising-strength schedule σn α now varies across time nodes, we replace and score-weight schedule αn . Since the proximal drift Tη,σ the single score-mismatch error δ with the node-dependent error and adjust the uniform Lipschitz bound of the proximal drift accordingly. To control these effects, we replace Assumption 2 with the following annealed counterpart. 14

PiX-MC

Wei et al.

Algorithm 3 multi-block APiX-MC Input: Initial distribution µ0 , stepsize γ > 0, number of blocks M , time step number N inside each α with denoising strength schedule {σ }M N −1 and score-weight schedule {α }M N −1 , block, drift Tη,σ i i=0 i i=0 and Picard iteration number K. i.i.d. 1: Sample x0 ∼ µ0 and, independently of x0 , sample W1 , . . . , WM N ∼ N (0, I) i √ P γWj for i = 1, . . . , M N . 2: Set B0 ← 0, and Biγ ← j=1

3: for m = 0, . . . , M − 1 do 4: 5: 6: 7: 8:

(0)

(k)

Set xmN +i ← xmN for i = 1, . . . , N , and xmN ← xmN for k = 0, . . . , K − 1 for k = 0, . . . , K − 1 do (k) (k) αmN (x(k) ), . . . , T α(m+1)N −1 (x(k) Tη,σ η,σ(m+1)N −1 mN (m+1)N −1 ) ← ParallelCompute(xmN , . . . , x(m+1)N −1 ) mN for i = 1, . . . , N do i−1 √
P αmN +j (k) (k+1) xmN +i = xmN − γ Tη,σmN +j (xmN +j ) + 2 B(mN +i)γ − BmN γ j=0

9:

end for 10: end for (K) 11: x(m+1)N ← x(m+1)N 12: end for (K) 13: Return xM N = xM N M N −1 N −1 N −1 and {αn }M Assumption 4. Given the schedules {σn }M n=0 , there exist {δn }n=0 , LS > 0, and n=0 RS > 0 such that the score network Sθ satisfies:

(i) for every n = 0, . . . , M N − 1 and x ∈ Rd , ∥Sθ (x; σn ) − ∇V (x)∥ ≤ δn ;

(40)

(ii) for every n = 0, . . . , M N − 1 and x, z ∈ Rd , ∥Sθ (x; σn ) − Sθ (z; σn )∥ ≤ LS ∥x − z∥;

(41)

(iii) for every n = 0, . . . , M N − 1 and x ∈ Rd , (42)

∥Sθ (x; σn )∥ ≤ RS . Assumption 4(iii) controls the additional discrepancy induced by the score weights. Indeed, ∥αn Sθ (x; σn ) − ∇V (x)∥ ≤ δn + (αn − 1) RS .

(43)

Accordingly, define the maximal and time-averaged annealed-score errors by δann ≜

max

0≤n≤M N −1

{δn + (αn − 1) RS } ,

2 δ̄ann ≜

MX N −1 1 [δn + (αn − 1) RS ]2 . MN n=0

15

(44)

PiX-MC

Wei et al.

2 The maximal score error δann controls the in-block Picard error, whereas δ̄ann appears in the timeM N −1 averaged convergence bound. Since {αn }n=0 is non-increasing, the proximal drifts are uniformly Λη,ann -Lipschitz, where

Λη,ann ≜ α0 LS + βL

1 + ηα0 LS , 1 − ηαL

qann ≜ τ Λη,ann .

(45)

Here, qann is the uniform in-block Picard contraction factor for multi-block APiX-MC. The following theorem states the convergence result for multi-block APiX-MC. N Theorem 3. Let (µt )t≥0 denote the law of the continuous interpolation of {xn (K) }M n=0 generated by multi-block APiX-MC, where N > 0 denotes the number of time steps in each time block and M > 0 is the number of blocks. Assume Assumptions 1, 3, and 4 hold, and choose sufficiently small parameters η, γ > 0. Then, the following inequality holds for multi-block APiX-MC:

 1Z T 0.75 − O(η 2 ) − O(γ 2 ) FI(µt ∥π)dt T 0 KL(µ0 ∥π) 2 + O(η 2 ) + O(δ̄ann ) + O(γ) | {z } | {z } | {z } T

 ≤

Splitting discrepancy

where Rann = 2τ d +

Annealed score mismatch error

Discretization error

2K−2 + O(qann ) Rann . | {z }

(46)

Finite Picard refinement error

2 +G2 3τ 2 (δann ) . 2 1−3η 2 βL

Proof. See Appendix E for the detailed proof and full expression of the terms. Theorem 3 extends the multi-block guarantee in Theorem 2 to the weighted annealing setting. Compared with multi-block PiX-MC, the score-mismatch term is replaced by the time-averaged 2 , while the finite Picard-refinement error is governed by the uniform annealed-score discrepancy δ̄ann contraction factor qann . In particular, 2 δ̄ann ≤

MX N −1 M N −1 2RS2 X 2 δn2 + (αn − 1)2 . MN MN n=0

(47)

n=0

The first term represents the approximation error of score networks, whereas the second term is induced by the score weights. If αn = 1 after a finite annealing stage, which is true in our experiments, the latter contribution decreases as M N increases. Meanwhile, with qann < 1, the finite Picardrefinement error decays geometrically with K. More broadly, annealing and Picard refinement act on complementary aspects of the computation: annealing improves the sequential trajectory itself, reducing the effective sequential effort required to attain a given reconstruction quality, whereas Picard refinement reduces the wall-clock cost of evaluating that trajectory through time parallelism. Their gains can therefore compound in practice. As demonstrated in §5, APiX-MC yields additional wall-clock acceleration beyond trajectory parallelization.

4

Numerical Validation of the Theory

We consider a compressed sensing problem whose target posterior is available in closed form. This setting allows us to directly compare the empirical distributions generated by PiX-MC algorithms 16

PiX-MC

Wei et al.

1.0

0.05

std dev

0.0

0.00

1.0

0.05

0.0

0.00

PiX-MC+Pretrained Score

PiX-MC+Analytical Score 1.0

0.05

1.4 APiX-MC +Pretrained Score PiX-MC +Pretrained Score PiX-MC +Analytical Score

0.0

0.00

1.0

0.05

0.0

0.00

APiX-MC+Pretrained Score

MMD2

True Posterior mean

0.0

Time Step

0

1000

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 MMD2 , while APiX-MC reduces MMD2 more rapidly than PiX-MC with the learned score.

with the true posterior, providing direct empirical support for the convergence results established in our theorems. As the single-block implementation of PiX-MC is computationally demanding, we implement all variants of PiX-MC using their multi-block schemes, namely Algorithms 2–3. For brevity, we omit the qualifier "multi-block" hereafter. All PiX-MC algorithms use eight NVIDIA RTX PRO 6000 Blackwell GPUs. Experimental Setup. We construct a Gaussian image prior (48)

p(x) = N x; µprior , Σprior ,

where µprior and Σprior are estimated from 70,000 FFHQ images [45], after converting the images to grayscale and x∗ ∼
resizing them to 32 × 32. A ground-truth image is∗ drawn according to 100×1024 N µprior , Σprior . We then generate the linear measurements y = Ax + e, where A ∈ R
is a Gaussian random matrix whose rows are normalized to unit norm, and e ∼ N 0, σy2 I , with σy = 0.001. Under this Gaussian prior and likelihood, the posterior distribution is also Gaussian,
π(x|y) = N x; µpost , Σpost , (49) where

 Σpost =

−1

1 ⊤ Σ−1 prior + σ 2 A A y

 ,

µpost = Σpost

1 ⊤ Σ−1 prior µprior + σ 2 A y y

 .

(50)

We train a score network following the EDM framework [44] by using samples from the empirical Gaussian prior. We evaluate three instantiations of PiX-MC: PiX-MC and APiX-MC equipped with the
learned score network, and PiX-MC with the analytical score ∇V (x) = Σ−1 x − µ . Comparing prior prior PiX-MC and APiX-MC with the learned score network additionally illustrates the effect of annealing on the evolution of the sampled distribution. Distributional Evaluation. For each algorithm, we run 150 independent chains and record the samples every 10 time steps. We also draw 150 independent samples directly from the closed-form 17

PiX-MC

Wei et al.

posterior. At the final time step, we compute the empirical posterior mean and pixelwise standard deviation from the samples. To assess the full distribution beyond these first two moments, we compute the squared maximum mean discrepancy (MMD) [41] at every recorded time step. Given e = {e em } and exact posterior samples X = {x1 , . . . , xn }, we use the empirical samples X x1 , . . . , x estimator m

e X) = MMD2 (X,

n

m

m

n

n

1 XX 1 XX 2 XX ej ) + 2 k(e xi , xj ). k(e xi , x k(xi , xj ) − 2 m n mn i=1 j=1

i=1 j=1

(51)

i=1 j=1

  ∥x1 −x2 ∥22 with σMMD = 3.0. We use the Gaussian kernel k(x1 , x2 ) = exp − 2σ 2 MMD

Distribution-Level Results Figure 2 reports the empirical posterior mean and pixelwise standard deviation at the final time step, together with the evolution of MMD2 over the sampling trajectory. All three PiX-MC instantiations accurately recover the posterior mean and uncertainty maps. As expected, PiX-MC with the analytical score most closely matches the true posterior, owing to the absence of score approximation error. The MMD2 curves show that the empirical distributions produced by all three methods move toward the target posterior, providing empirical support for the finite-horizon stationarity bounds established in Theorems 2–3. At the final time step, PiX-MC with the analytical score achieves the lowest MMD2 of 0.0036, while APiX-MC and PiX-MC with the learned score attain 0.1361 and 0.1480, respectively. The remaining gap between the learned-score variants and the analytical variant is consistent with the score approximation error analyzed in Section 3. Compared with PiX-MC, APiX-MC reduces MMD2 more rapidly and reaches a lower final value, corroborating the discussion in Section 3.4 that annealing improves the underlying sampling trajectory while Picard refinement enables its parallel evaluation.

5

Real-World Imaging Experiments

We now present our experimental results, which are organized around three goals. First, we verify the behavior of Picard refinement on MRI and Rician noise removal, covering both linear and nonlinear inverse problem setups. Second, we evaluate practical acceleration on large-scale 1024 × 1024 image deblurring, where we vary the GPU budget and the block size to examine how the multi-block implementation affects wall-clock runtime. Third, we apply PiX-MC to a large-scale 512 × 512 × 80 sparse-view CT problem to demonstrate its scalability on high-dimensional imaging tasks. All hyperparameters are tuned to maximize the peak signal-to-noise ratio (PSNR). Additional implementation details and hyper-parameter settings are provided in Appendix F. Experimental Setup. All experiments are conducted on a single compute node equipped with eight NVIDIA RTX PRO 6000 Blackwell GPUs and two AMD EPYC 9334 32-Core processors. All Picard-based methods are implemented using the multi-block scheme and executed on all eight GPUs, whereas sequential methods are executed on a single GPU. We use DPS [16] and DAPS [97] as external baselines. For the deblurring and 3D CT experiments, we additionally include sequential and ablation variants to isolate the effects of Picard time parallelization, the proximal-likelihood drift, and the annealing strategy. We use the following nomenclature throughout this section. L-MC denotes the standard sequential Langevin sampler with score network Sθ (Eq. (9)), and X-MC denotes its proximallikelihood counterpart (Eq. (21)). PiL-MC is obtained by applying Picard time parallelization to L-MC 18

PiX-MC

105

Wei et al.

MRI × 8

Non-annealed PiX-MC PiL-MC

103

Annealed APiX-MC APiL-MC

𝜀𝜀5,𝑘𝑘

10−5 0

𝜀𝜀3,𝑘𝑘

Picard iteration 𝑘𝑘

75 0

Picard iteration 𝑘𝑘

75

10−3 0

Rician Denoising Annealed APiX-MC APiL-MC

Non-annealed PiX-MC PiL-MC

Picard iteration 𝑘𝑘

50 0

Picard iteration 𝑘𝑘

50

Figure 3: Visualization of the in-block Picard error ε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.

by replacing the drift Tη in Algorithm 2 with Sθ + ∇L. The prefix “A” denotes the corresponding annealed variants, yielding AL-MC, AX-MC, and APiL-MC. Unless otherwise specified, all methods use the same score network for each imaging task. Task-specific score network choices and multiblock configurations are described in the corresponding subsections, while complete hyperparameter settings are provided in Appendix F.1.

5.1

Behavior of Picard Refinement

We first present results on linear MRI reconstruction and nonlinear Rician denoising. These experiments empirically verify the geometric decay of the Picard error predicted by Lemmas 2 and 3, demonstrating that the predicted convergence behavior is agnostic to the linearity of the inverse problem. Implementation Details. For MRI, we use the pretrained score network provided by [86]1 . We use ten brain images with size 320 × 320 from the FastMRI dataset [96] as the test set. We use a radial subsampling mask corresponding to 8× acceleration. We additionally add additive white Gaussian noise (AWGN) with standard deviation σy = 0.05. For Rician denoising, we use the pretrained 256 × 256 unconditioned diffusion model2 . The tested dataset (CBSD10) is composed of a subset of 10 images of the CBSD68 dataset [57]. We add Rician noise with σy = 0.1. We set M = 10, N = 75, K = 100 for MRI, and M = 5, N = 100, K = 50 for Rician denoising. We refer to Appendix F.2 for more technical details on these two tasks. Geometric Decay of the Picard Error. Lemma 3 establishes geometric contraction of the in-block Picard error εm,k defined in (35). Figure 3 reports the empirical evolution of εm,k on a logarithmic scale for the MRI reconstruction and Rician denoising tasks. We report the block indexed by m = 5 for MRI reconstruction and the block indexed by m = 3 for Rician denoising, with εm,k averaged over the test dataset. After a short initial transient, the Picard error decreases approximately linearly on the logarithmic scale for all methods, including both proximal and non-proximal variants as well 1 2

The pretrained score network is available at https://github.com/sunyumark/PnP-MonteCarlo The pretrained diffusion model is available at https://github.com/openai/guided-diffusion

19

PiX-MC

Wei et al.

MRI × 8

Rician Denoising X-M

MC

APi

PiX-

C

PiX-

APi

MSE

1

MSE

1

C X-M

MC

10-3 0

Pic

ard

0

Ite rat

ion

100

750

e Tim

p Ste

0

Pic

ard

10-3 0

0

Ite

rat

ion

100

750

e Tim

p Ste

Pic

ard

0

Ite

rat

ion

50

500

e Tim

p Ste

0

Pic

ard

0

Ite

rat

ion

50

500

eS Tim

tep

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. Rician Denoising (σy = 0.1) Method PSNR ↑ SSIM ↑ TV 27.64 0.7678 DAPS 29.66 0.7598 MRI (Accel.= 8× with σy = 0.05) PiL-MC 30.21 0.8201 APiL-MC 31.30 0.8353 PiX-MC 31.41 0.8439 0.8472 APiX-MC 31.69

Method PSNR ↑ LPIPS ↓ SSIM ↑ DPS 26.02 0.1429 0.7865 DAPS 27.98 0.1723 0.7533 PiL-MC 30.51 0.1139 0.8657 APiL-MC 30.62 0.1025 0.8710 29.75 0.1632 0.8363 PiX-MC 0.0826 0.8677 APiX-MC 30.43

Table 1: Average quantitative results of PiX-MC and its variants for Rician noise removal (σy = 0.1) on the CBSD10 dataset and MRI reconstruction (8×, σy = 0.05) on 10 FastMRI brain images. The best values are shown in bold.

as their annealed counterparts. This behavior indicates rapid convergence of the Picard iteration in practice and is empirically consistent with the geometric contraction established in Lemma 3. Trajectory Evolution. We further visualize the evolution of the trajectories during the Picard iteration. Given an image x and the ground-truth image xGT , we compute the mean-squared error (MSE) MSE := d1 ∥x − xGT ∥22 , and Figure 4 visualizes the resulting MSE curves. At the initial Picard iteration (k = 0), all time nodes within each block are initialized to the same point, yielding piecewise-constant MSE curves. As the Picard iteration proceeds, the trajectories are progressively refined, and the piecewise-constant curves evolve into smooth trajectories across successive blocks. This evolution provides an intuitive illustration of the Picard refinement mechanism and is consistent with the convergence guarantees established in Lemma 3 and Theorems 2–3. Quantitative and Qualitative Comparison. Now we evaluate the final reconstruction quality of the proposed methods. Table 1 reports the average quantitative results over ten test images for each task. For Rician denoising, APiL-MC achieves the highest PSNR and SSIM, while APiX-MC attains the lowest LPIPS. For accelerated MRI, APiX-MC achieves the best PSNR and SSIM among all compared methods. Overall, the proposed methods are competitive with diffusion-based baselines, and annealing consistently improves reconstruction quality over the corresponding non-annealed variants. Figures 5–6 provide visual comparisons. For accelerated MRI, TV exhibits staircase artifacts, 20

PiX-MC

Wei et al.

PiL-MC

TV

26.57 dB

APiL-MC

29.67 dB

31.02 dB

29.30 dB

31.47 dB

31.15 dB

23.16 dB

GT

APiX-MC

PiX-MC

DAPS

ZF

PSNR

Figure 5: Visual comparison for 8× accelerated MRI reconstruction with measurement noise level σ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.

while DAPS leaves noticeable residual noise. Both PiL-MC and APiL-MC remove most blur and noise but lose fine image details. In contrast, PiX-MC and APiX-MC recover sharper edges with fewer residual artifacts, as highlighted in the zoomed-in regions of Figure 5. A similar trend is observed for Rician denoising. DPS exhibits noticeable color shifts, and DAPS leaves residual noise in the regions highlighted in purple. While PiL-MC and APiL-MC effectively suppress the Rician noise, PiX-MC and APiX-MC preserve finer structures. In particular, APiX-MC successfully reconstructs the letter “A” together with the surrounding dot-like markings (red box), whereas these details are blurred into a thin horizontal line by PiX-MC and are absent from the remaining reconstructions. These visual observations are consistent with the quantitative results in Table 1.

5.2

Practical Acceleration on Large-Scale Deblurring

In this section, we investigate whether Picard time parallelization can accelerate posterior sampling without sacrificing reconstruction quality, and how this acceleration scales with the GPU budget and the number of time nodes per block.

21

PiX-MC

Wei et al.

DPS

38.98 dB

37.84 dB

28.81 dB

PiX-MC

DAPS

32.74 dB

Noisy

APiL-MC

PiL-MC

20.03 dB

GT

APiX-MC

39.19 dB

37.95 dB

PSNR

Figure 6: Visual comparison for Rician denoising with measurement noise level σ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.

Implementation Details. For the score network, we use the same diffusion model as in the Rician noise removal experiment. We use ten natural images from the DIV2K validation set from [2]3 . We center-crop each image to size 1024 × 1024. We apply a camera-shake blur from [52] to the test images, and add AWGN noise with standard deviation σy = 0.05. For the proximal methods, the proximal operator on likelihood can be solved efficiently in the Fourier domain according to [62]. For each block of the Picard-based methods, we adopt an adaptive stopping criterion: the Picard refinement is terminated when either the maximum number of Picard iterations K = 20 is reached, or the normalized Picard error satisfies εm,k /d < 3 × 10−4 . More implementation details are provided in Appendix F. Acceleration under Different Parallelization Configurations. We first investigate the practical acceleration under different numbers of GPUs and time nodes N within each block. Table 2 reports the runtime for reaching an average PSNR of 26 dB and the corresponding speedup over each sequential counterpart under different configurations. It can be seen that increasing the number of GPUs improves the speedup across all tested methods and block configurations. In particular, using eight GPUs with N = 8 yields speedups between 2.87× and 2.99× across the four Picard-based methods, as all eight time nodes can be evaluated concurrently. Moderate oversubscription with N = 16 remains effective, but the benefit decreases as N increasingly exceeds the number of available GPUs. For example, with four GPUs and N = 32, all four Picard-based methods are slower than their sequential counterparts. These results demonstrate the practical trade-off between the block length and the available computational resources. Figure 7 plots the average PSNR against runtime in seconds under the configuration of eight GPUs and N = 8. As shown, Picard parallelization consistently accelerates convergence. In particular, 3

The DIV2K dataset is available at https://data.vision.ee.ethz.ch/cvl/DIV2K/

22

PiX-MC

Wei et al.

26.0

26.0

PSNR (dB)

28.5

PSNR (dB)

28.5

Non-annealed PiX-MC PiL-MC X-MC L-MC

7.5 0

127

365

171

491

Runtime (s)

718

7.5 0

Annealed APiX-MC APiL-MC AX-MC AL-MC 54

63

157

190

Runtime (s)

718

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 4 GPUs 8 GPUs Methods 4 GPUs 8 GPUs

PiL-MC vs. L-MC (491s) N =8 264s / 1.86× 171s / 2.87×

N = 16 328s / 1.50× 175s / 2.81×

APiL-MC vs. AL-MC (190s)

N = 32 550s / 0.89× 295s / 1.66×

PiX-MC vs. X-MC (365s) N =8 195s / 1.87× 127s / 2.87×

N = 16 230s / 1.59× 138s / 2.64×

N =8 99s / 1.92× 63s / 2.99×

N = 16 139s /1.37× 81s / 2.35×

N = 32 250s / 0.76× 151s / 1.26×

APiX-MC vs. AX-MC (157s)

N = 32 440s / 0.83× 188s / 1.94×

N =8 84s / 1.87× 54s / 2.91×

N = 16 120s / 1.31× 70s / 2.24×

N = 32 194s / 0.81× 94s / 1.67×

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.

APiX-MC reaches an average PSNR of 26 dB in only 54 seconds, compared with 157 seconds for AX-MC, 365 seconds for X-MC, and 491 seconds for L-MC. Beyond the pairwise speedups, this comparison demonstrates that Picard time parallelization can compound with mechanisms that accelerate the underlying sequential trajectory. For example, the proximal update and annealing together reduce the runtime from 491 seconds for L-MC to 157 seconds for AX-MC, corresponding to a 3.13× sequential-side speedup. Picard parallelization further reduces this runtime to 54 seconds, providing an additional 2.91× parallelization-side speedup. Together, they yield an overall 9× speedup over L-MC. These results empirically support the discussion at the end of §3.4 that Picard parallelization provides a complementary direction that can be combined with sequential strategies. This is particularly valuable for large-scale imaging tasks, where sequential trajectory evaluation is computationally expensive. Quantitative and Qualitative Comparison. We next show that PiX-MC methods accelerate posterior sampling without compromising reconstruction quality. Table 3 reports the average final PSNR, SSIM, and LPIPS over the test set. The L-MC family is also included as references. PiX-MC and APiX-MC closely match the reconstruction quality of their corresponding sequential counterparts, with PSNR differences of at most 0.01 dB, while remaining competitive with the external baselines DPS and DAPS. 23

PiX-MC

Wei et al.

DPS

L-MC

AL-MC

X-MC

AX-MC

Blurry

DAPS

PiL-MC

APiL-MC

PiX-MC

APiX-MC

GT

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 × 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) ↑

LPIPS ↓

SSIM ↑

Figure 8 compares visual results under a DPS 26.43 0.2498 0.7142 fixed runtime budget of 200 seconds. While DAPS 28.17 0.1632 0.7444 DPS and DAPS require 334 and 741 seconds, reL-MC 26.82 0.2605 0.7273 spectively, to converge, PiX-MC and APiX-MC AL-MC 27.92 0.2042 0.7737 already recover sharp edges and fine structures X-MC 27.24 0.2686 0.7342 AX-MC 27.98 0.2138 0.7719 with substantially reduced blur and noise within the 200-second runtime budget. Compared with PiL-MC 26.82 0.2604 0.7268 APiL-MC 27.92 0.2048 0.7738 their sequential counterparts, Picard time parPiX-MC 27.23 0.2726 0.7343 allelization consistently produces higher-quality APiX-MC 27.98 0.2120 0.7718 reconstructions under the same wall-clock budget by advancing further along the sampling tra- Table 3: Average quantitative results for large-scale dejectory. These observations, together with the blurring on ten 1024 × 1024 test images from DIV2K quantitative results, show that Picard time par- validatation set. allelization achieves substantial acceleration without sacrificing the final reconstruction quality.

5.3

Scaling to 3D CT Reconstruction

Finally, we evaluate PiX-MC on a large-scale sparse-view CT reconstruction problem with volume size 512 × 512 × 80, demonstrating its practical scalability to high-dimensional 3D imaging.

24

PiX-MC

Wei et al.

35.0

35.0 32.5

18.0 0

PSNR (dB)

PSNR (dB)

30.5 Non-annealed PiX-MC PiL-MC X-MC L-MC 37

116

154

Runtime (min)

474

526

18.0 0

Annealed APiX-MC APiL-MC AX-MC AL-MC 9 28

103

305

Runtime (min)

526

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× wall-clock speedup over L-MC.

Implementation Details. For the score network, we trained a diffusion model based on the EDM framework [44], using CT images from the AMOS [42] and LDCT [59] datasets. The image values in Hounsfield units were clamped to the range [−1000, 1000], then linearly rescaled to [−1, 1]. We consider a sparse-view parallel-beam CT reconstruction problem with 30 uniformly spaced projection views over 180◦ . The target volume consists of 80 axial slices, each of size 512 × 512. Measurements are generated by applying the parallel-beam projection operator and adding i.i.d. Gaussian noise with standard deviation σy = 0.5 in the sinogram domain. Filtered back-projection (FBP) is used for initialization. Forward projection, back projection, and FBP are implemented using LEAP [48]4 . To encourage consistency across adjacent slices, we additionally incorporate a one-dimensional inter-slice Huber-TV (HTV) regularizer H(x) into the proximal drift. The detailed formulation of the HTV regularizer and the resulting proximal drift are provided in Appendix F.2. Accordingly, the proximal drift Tη,σ is given by
 1 Tη,σ ≜ I − proxηL ◦ I − ηSθ (·; σ) − ηβH ∇H , (52) η where βH > 0 is a balancing parameter for H. All Picard-based methods use eight GPUs with N = 8 time nodes per block. Runtime Acceleration. Figure 9 plots the mean slice-wise PSNR over the 80 axial slices against wall-clock time. The curves distinguish the effects of the three algorithmic components: the proximallikelihood drift, Picard time parallelization, and annealing. Compared with their likelihood-gradient counterparts, X-MC and AX-MC converge substantially faster and attain higher final PSNR, demonstrating that the proximal-likelihood drift provides a more effective reconstruction trajectory for this CT task. Building upon these improved trajectories, PiX-MC and APiX-MC further reduce the wall-clock time while preserving the attainable reconstruction quality. Specifically, PiX-MC reaches 4

LEAP is available at https://github.com/LLNL/LEAP

25

PiX-MC

L-MC

30.71 dB (526 mins)

PiL-MC

31.95 dB (380 mins)

Wei et al.

AL-MC

32.16 dB (526 mins)

APiL-MC

32.39 dB (200 mins)

X-MC

AX-MC

32.59 dB (300 mins)

PiX-MC

33.60 dB (100 mins)

APiX-MC

32.70 dB (100 mins)

FBP

33.64 dB (35 mins)

19.16 dB

GT

PSNR (Runtime)

Residual

Figure 10: Visual comparison of PiX-MC and baseline methods for large-scale sparse-view 3D CT reconstruction. The 20th slice of a 512 × 512 × 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.

its reference PSNR in only 37 minutes, compared with 116 minutes for X-MC, corresponding to a 3.14× speedup. Likewise, APiX-MC requires only 9 minutes, compared with 28 minutes for AXMC, yielding a 3.11× speedup. The same pattern is observed for the likelihood-gradient methods, where PiL-MC and APiL-MC accelerate L-MC and AL-MC by 3.08× and 2.96×, respectively. Annealing further improves the attainable reconstruction quality while preserving the approximately 3× speedup provided by Picard time parallelization. Overall, APiX-MC achieves 32.5 dB in only 9 minutes, demonstrating how the proximal-likelihood drift, annealing, and Picard time parallelization compound to deliver approximately a 50× wall-clock speedup over the standard L-MC. Axial-Slice and Volumetric Reconstruction. Figures 10 and 11 compare the reconstructed axial slice after convergence and the reconstructed volume under a fixed wall-clock budget of 35 minutes, respectively. Across all four sequential/Picard pairs, Picard time parallelization achieves comparable or better reconstruction quality in substantially less wall-clock time. The visual comparisons further separate the effects of the remaining algorithmic components: compared with the likelihood-gradient methods, the proximal-likelihood variants more effectively suppress sparse-view streak artifacts and recover finer anatomical structures, while annealing further reduces residual artifacts and improves visual fidelity. These improvements are consistently observed in both the selected axial slice and the full 3D volume. Overall, APiX-MC achieves the best visual quality, and the two figures together 26

PiX-MC

L-MC

Wei et al.

𝑥𝑥

𝑧𝑧

𝑦𝑦

AL-MC

22.72 dB PiL-MC

X-MC

AX-MC

28.57 dB APiL-MC

25.64 dB

26.89 dB PiX-MC

30.87 dB

FBP

32.89 dB APiX-MC

30.42 dB

19.50 dB GT

34.14 dB

PSNR

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.

demonstrate that the proximal-likelihood drift, annealing, and Picard time parallelization complement each other to enable high-quality large-scale 3D CT reconstruction.

6

Conclusion

This paper introduces PiX-MC, a principled time-parallel posterior sampling framework for Bayesian imaging. Motivated by forward-backward splitting, we construct a proximal drift that exploits efficient proximal operators available for imaging likelihoods. Instead of advancing the samples sequentially along a long trajectory, PiX-MC parallelizes the trajectory across time nodes, enabling multiple time nodes to be processed concurrently across GPUs. To improve practical scalability under limited GPU budgets, we further develop annealed and multi-block variants. The annealed variant empirically accelerates the convergence to reach a target reconstruction quality, while the multi-block variant partitions the trajectory into blocks, allowing PiX-MC to operate efficiently under limited GPU budgets. We establish finite-time stationarity guarantees for all variants of PiX-MC under mild assumptions. Extensive experiments across linear, nonlinear, and large-scale 2D and 3D imaging problems demonstrate practical scalability of PiX-MC and substantial wall-clock acceleration while maintaining reconstruction quality.

7

Acknowledgements

Y. Sun acknowledges support from the U.S. National Science Foundation (NSF) under Grant CCF2542022. E. Bell is supported by the U.S. Department of Energy Computational Science Graduate Fellowship under Award No. DE-SC0026073. Y. Chen acknowledges support from the AMS-Simons Travel Grant and National Science Foundation award DMS-2608264. 27

PiX-MC

A

Wei et al.

Tool Lemmas

Lemma 4 (Lipschitz continuity of the proximal drifts). Suppose Assumption 1 holds and 1 − ηαL > 0. (i) If Assumption 2 (ii) holds, the proximal drift Tη,σ is Λη -Lipschitz, where Λη ≜ L S + β L

1 + ηLS . 1 − ηαL

(53)

N −1 (ii) If Assumption 4 (ii) holds and {αn }M is non-increasing, then every weighted proximal drift n=0 α n Tη,σn is uniformly Λη,ann -Lipschitz, where

Λη,ann ≜ α0 LS + βL

1 + ηα0 LS . 1 − ηαL

(54)

Proof. We prove (i) first. For x1 , x2 ∈ Rd , define z1 = proxηL (x1 − ηSθ (x1 ; σ)),

z2 = proxηL (x2 − ηSθ (x2 ; σ)).

(55)

Since L is αL -weakly convex and 1 − ηαL > 0, proxηL is (1 − ηαL )−1 -Lipschitz. Thus, ∥z1 − z2 ∥ ≤

1 + ηLS 1 ∥x1 − x2 − η (Sθ (x1 ; σ) − Sθ (x2 ; σ))∥ ≤ ∥x1 − x2 ∥. 1 − ηαL 1 − ηαL

(56)

The proximal optimal condition yields Tη (x1 ) = ∇L(z1 ) + Sθ (x1 ; σ),

Tη (x2 ) = ∇L(z2 ) + Sθ (x2 ; σ).

(57)

Hence, ∥Tη,σ (x1 ) − Tη,σ (x2 )∥ ≤ ∥∇L (z1 ) − ∇L (z2 )∥ + ∥Sθ (x1 ; σ) − Sθ (x2 ; σ)∥ ≤ βL ∥z1 − z2 ∥ + LS ∥x  1 − x2 ∥ 1 + ηLS ≤ LS + β L ∥x1 − x2 ∥. 1 − ηαL

(58)

For (ii), replace Sθ (·; σ) by αn Sθ (·; σn ). Since αn ≤ α0 , the weighted score is uniformly α0 LS Lipschitz. Repeating the argument yields (ii). The following lemma from [14, Lemma 16] is essential for the proof. Lemma 5. Under Assumption 1, ∇L(x) + ∇V (x) is (βL + βV )-Lipschitz. For any probability measure µ ≪ π and x ∼ µ, we have # " 2   dµ + 2d(βL + βV ) = FI(µ∥π) + 2d(βL + βV ). (59) E ∥∇L(x) + ∇V (x)∥2 ≤ E ∇ ln dπ We introduce the following technical lemma. Lemma 6. Assume that Assumptions 1 and 2 (i) are satisfied, and η > 0 is chosen sufficiently small such that 1 − ηαL , 1 − 3η 2 βL2 > 0. Then, for any x ∈ Rd , ∥Tη (x)∥2 ≤

3 3δ 2 + ∥∇L(x) + ∇V (x)∥2 . 1 − 3η 2 βL2 1 − 3η 2 βL2 28

(60)

PiX-MC

Wei et al.

Proof. Let

z = (I + η∇L)−1 (x − ηSθ (x; σ)).

(61)

z + η∇L(z) = x − ηSθ (x; σ).

(62)

z − x = −η∇L(z) − ηSθ (x).

(63)

∥z − x∥2 = η 2 ∥∇L(z) + Sθ (x)∥2 = η 2 ∥∇L(z) − ∇L(x) + Sθ (x; σ) − ∇V (x) + ∇L(x) + ∇V (x)∥2 ≤ 3η 2 ∥∇L(z) − ∇L(x)∥2 + 3η 2 ∥Sθ (x; σ) − ∇V (x)∥2 + 3η 2 ∥∇L(x) + ∇V (x)∥2 ≤ 3η 2 βL2 ∥z − x∥2 + 3η 2 δ 2 + 3η 2 ∥∇L(x) + ∇V (x)∥2 .

(64)

Then Tη (x) = η1 (x − z). We have

That is, Thus,

This gives us (1 − 3η 2 βL2 )∥z − x∥2 ≤ 3η 2 δ 2 + 3η 2 ∥∇L(x) + ∇V (x)∥2 3 3δ 2 + ∥∇L(x) + ∇V (x)∥2 . ∥Tη (x)∥2 ≤ 2 2 2 2 1 − 3η βL 1 − 3η βL

⇐⇒

(65)

We now prove Lemma 1. Proof. Let z = proxηL (x − ηSθ (x; σ)). Thus z + η∇L(z) = x − ηSθ (x; σ). Therefore Tη (x) = ∇L(z) + Sθ (x; σ). It follows that
Tη (x) − ∇L(x) + ∇V (x) = ∇L(z) − ∇L(x) + Sθ (x; σ) − ∇V (x). (66) Then we have

∥Tη (x) − ∇L(x) + ∇V (x) ∥2 ≤ 2∥∇L(z) − ∇L(x)∥2 + 2∥Sθ (x; σ) − ∇V (x)∥2 ≤ 2βL2 ∥z − x∥2 + 2δ 2 .

(67)

Since z − x = −ηTη (x), this gives
∥Tη (x) − ∇L(x) + ∇V (x) ∥2 ≤ 2η 2 βL2 ∥Tη (x)∥2 + 2δ 2 .

(68)

Applying Lemma 6 yields
6η 2 βL2 δ 2 6η 2 βL2 ∥Tη (x) − ∇L(x) + ∇V (x) ∥2 ≤ 2δ 2 + + ∥∇L(x) + ∇V (x)∥2 . 1 − 3η 2 βL2 1 − 3η 2 βL2 Finally, note that 2δ 2 +

(69)

6η 2 βL2 δ 2 2δ 2 = . Hence, 1 − 3η 2 βL2 1 − 3η 2 βL2

∥Tη (x) − ∇L(x) + ∇V (x) ∥2 ≤

6η 2 βL2 2δ 2 + ∥∇L(x) + ∇V (x)∥2 . 1 − 3η 2 βL2 1 − 3η 2 βL2

29

(70)

PiX-MC

B

Wei et al.

Proof of single-block PiX-MC

We present the convergence result of single-block PiX-MC. Before that, we define some essential constants used throughout the proof. Define the constants cscore , csplit : 2δ 2 , 1 − 3η 2 βL2

cscore =

csplit =

6η 2 βL2 . 1 − 3η 2 βL2

Define the constant ρ:  ρ=

18γ 2 (βL + βV )2 1− 1 − 3η 2 βL2

−1

(71)

.

(72)

2(βL + βV )2 (1 + 9η 2 βL2 ) . 1 − 3η 2 βL2

(73)

Define the constant A1 : A1 = 2(βL + βV )2 (1 + 2csplit ) = Define the constants A4 , A5 : 18ρA1 γ 2 A4 = 4csplit + , 1 − 3η 2 βL2

 A5 = 2cscore + A1 ρ

 9γ 2 δ 2 + 6γd . 1 − 3η 2 βL2

(74)

Define the constants B2 , B3 :  2K−2 · 2T d + B2 = qsing

  2 2  1 2 3T δ + 6T d(β + β ) , L V 1 − 3η 2 βL2

and 2K−2 · B3 = qsing

3T 2 . 1 − 3η 2 βL2

(75)

(76)

Finally, we define the constant A: A := 3ρA1 B2 + A5 + 2A4 d(βL + βV ).

(77)

Recall that qsing = T Λη is the contraction factor introduced in §3.2. We now state a detailed version of Theorem 1 with explicit constants. Theorem 4 (Detailed version of Theorem 1). Let (µt )t≥0 denote the law for the continuous interpolation of {xn (K) }N n=0 generated by PiX-MC, where N > 0 is the total number of time steps. Assume Assumptions 1-2 hold, and the parameters η, γ > 0 are chosen sufficiently small such that 1 − ηαL > 0, 1 − 3η 2 βL2 > 0, 1 − Then, 1 (0.75 − A4 ) T

18γ 2 (βL + βV )2 > 0, qsing < 1, A4 < 0.75. 1 − 3η 2 βL2

(78)

1 KL(µ0 ∥π) + A + 3ρA1 B3 · FI(µ0 ∥π). T

(79)

Z T FI(µt ∥π)dt ≤ 0

30

PiX-MC

Note that

Wei et al.

36γ 2 (βL + βV )2 1 + 9η 2 βL2 24η 2 βL2 . A4 = +
2  2 (β +β )2 L V 1 − 3η 2 βL2 1 − 3η 2 βL2 1 − 18γ1−3η 2β2

(80)

L

Hence, A4 → 0 as (η, γ) → (0, 0), and therefore A4 < 0.75 for sufficiently small η and γ. More explicitly, direct simplification gives 0.75 − A4 = 0.75 ·

1 − 35η 2 βL2 − 66γ 2 (βL + βV )2 1 − 3η 2 βL2 − 18γ 2 (βL + βV )2

.

(81)

Consequently, the sufficient conditions η 2 βL2 <

1 , 35

1 − 35η 2 βL2 66

γ 2 (βL + βV )2 <

(82)

ensure that A4 < 0.75. Before the proof of Theorem 4, we first introduce a technical lemma, which shows that the error decays exponentially fast in the parallel refinement, if qsing < 1. Lemma 7 (Detailed version of Lemma 2). Let µ0 be the initial distribution of x0 and assume Assumptions 1-2 hold. Choose the parameters η, γ > 0 sufficiently small such that 1 − ηαL > 0,

1 − 3η 2 βL2 > 0,

(83)

qsing < 1.

Define the error across two successive Picard refinements at the k-th iteration as i h (k−1) 2 εk ≜ max E ∥x(k) ∥ . n − xn n=1,2,...,N

(84)

Then, the Picard iteration error for PiX-MC satisfies (85)

εK ≤ B2 + B3 · FI(µ0 ∥π). Proof. For each n = 1, 2, . . . , N , (k+1) (k) E∥xn − x n ∥2 =

E γ

≤ nγ 2 ≤ ≤

n−1 Ph

i=0 n−1 P

i 2

(k) (k−1) Tη (xi ) − Tη (xi ) (k)

(k−1)

E∥Tη (xi ) − Tη (xi

)∥2

i=0 n−1 P (k) (k−1) 2 nγ 2 Λ2η E∥xi − xi ∥ i=0 2 ε . n2 γ 2 Λ2η εk ≤ (nγΛη )2 εk ≤ qsing k

(86)

(0)

2 ε . Note that according to Algorithm 1, x Thus, εk+1 ≤ qsing n ≡ x0 . In particular, we have k (1) (0) E∥xn − xn ∥2 = E

γ

n−1 P i=0

=

Tη (x0 ) +

√

2

2Bnγ

n2 γ 2 E∥Tη (x0 )∥2 + 2nγd ≤ T 2 E∥Tη (x0 )∥2 + 2T d. 31

(87)

PiX-MC

Wei et al.

Note that since x0 is independent of Bnγ and E[Bnγ ] = 0, the cross term vanishes. Thus ε1 ≤ 2 ε and ε ≤ T 2 E∥T (x )∥2 + 2T d, by induction we have T 2 E∥Tη (x0 )∥2 + 2T d. Since εk+1 ≤ qsing 1 η 0 k (88)

2K−2 εK ≤ qsing T 2 E∥Tη (x0 )∥2 + 2T d . We bound E∥Tη (x0 )∥2 using Lemma 5. By Lemma 6, we have that E∥Tη (x0 )∥2 ≤

3 3δ 2 + E∥∇L(x0 ) + ∇V (x0 )∥2 . 2 2 1 − 3η βL 1 − 3η 2 βL2

(89)

Since ∇L(x) + ∇V (x) is (βL + βV )-smooth, by Lemma 5, we have (90)

E∥∇L(x0 ) + ∇V (x0 )∥2 ≤ 2d(βL + βV ) + FI(µ0 ∥π). Therefore, E∥Tη (x0 )∥2 ≤

3δ 2 3 + [2d(βL + βV ) + FI(µ0 ∥π)] . 2 2 1 − 3η βL 1 − 3η 2 βL2

(91)

Substituting (91) into (88) completes the proof. Now we are ready to prove Theorem 4. Proof. We will use the interpolation method. Let Xt∈[0,T ] denote the continuous-time interpolation (K)

(K)

of the final Picard trajectory {xn }N n=0 . Let Xnγ = xn . For any t ∈ [nγ, (n + 1)γ], define √ Xt := x(K) − (t − nγ)Tη (x(K−1) ) + 2(Bt − Bnγ ). n n

(92)

(K−1)

Now, let µt be the law of Xt , µt ≡ law(Xt ). With respect to the natural filtration, Tη (xn ) is Fnγ measurable and constant in time on this interval. The marginal process therefore has conditional drift h i (K−1)

−E Tη (xn )|Xt . Applying the entropy-dissipation identity and conditional Jensen’s inequality, as in [87, proof of Lemma 3], gives:   µt (Xt ) (K−1) ∂t KL(µt ∥π) = −FI(µt ∥π) + E ∇L(Xt ) + ∇V (Xt ) − Tη (xn ), ∇log . (93) π(Xt )

Thus we have that 2 3 ∂t KL(µt ∥π) ≤ − FI(µt ∥π) + E ∇L(Xt ) + ∇V (Xt ) − Tη (x(K−1) ) . n 4

(94)

We need to bound the last term. (K−1)

E∥∇L(Xt ) + ∇V (Xt ) − Tη (xn )∥2 (K−1) (K−1) 2 ≤ 2E∥∇L(Xt ) + ∇V (Xt ) − ∇L(xn ) − ∇V (xn )∥ (K−1) (K−1) (K−1) 2 +2E∥∇L(xn ) + ∇V (xn ) − Tη (xn )∥ .

(95)

We bound each of the two terms. For the first term, E∥∇L(Xt ) + ∇V (Xt ) − ∇L(x(K−1) ) − ∇V (x(K−1) )∥2 ≤ (βL + βV )2 E∥Xt − x(K−1) ∥2 . n n n 32

(96)

PiX-MC

Wei et al.

For the second term, note that (97)

∇L(x) + ∇V (x) − Tη (x) = ∇V (x) − Sθ (x; σ) + ∇L(x) − ∇L(z),
where z = (I + η∇L)−1 x − ηSθ (x; σ) . Thus, ∥∇L(x) + ∇V (x) − Tη (x)∥2 ≤ 2∥∇V (x) − Sθ (x; σ)∥2 + 2∥∇L(x) − ∇L(z)∥2 ≤ 2δ 2 + 2βL2 ∥x − z∥2 . By Lemma 6,

3δ 2 + 3∥∇L(x) + ∇V (x)∥2 . 1 − 3η 2 βL2

(99)

6η 2 βL2 2δ 2 + ∥∇L(x) + ∇V (x)∥2 1 − 3η 2 βL2 1 − 3η 2 βL2 = cscore + csplit ∥∇L(x) + ∇V (x)∥2 .

(100)

∥Tη (x)∥2 ≤ Thus,

(98)

∥∇L(x) + ∇V (x) − Tη (x)∥2 ≤ We have

(K−1)

E∥∇L(Xt ) + ∇V (Xt ) − Tη (xn

(K−1) 2 ∥ + 2c

≤ 2(βL + βV )2 E∥Xt − xn (K−1)

Now we bound E∥∇L(xn (K−1)

(K−1)

score + 2csplit E∥∇L(xn

(K−1)

) + ∇V (xn

)∥2 (K−1)

) + ∇V (xn

)∥2 .

)∥2 :

(K−1)

E∥∇L(xn ) + ∇V (xn )∥2 (K−1) (K−1) ≤ 2E∥∇L(xn ) + ∇V (xn ) − ∇L(Xt ) − ∇V (Xt )∥2 + 2E∥∇L(Xt ) + ∇V (Xt )∥2 (K−1) ≤ 2(βL + βV )2 E∥xn − Xt ∥2 + 2E∥∇L(Xt ) + ∇V (Xt )∥2 . and

(101)

(102)

(K−1)

E∥∇L(Xt ) + ∇V (Xt ) − Tη (xn )∥2   (K−1) 2 ≤ 2(βL + βV )2 + 4csplit (βL + βV )2 E∥Xt − xn ∥ 2 +4csplit E∥∇L(Xt ) + ∇V (Xt )∥ + 2cscore .

(103)

By the definition of A1 , we have (K−1)

E∥∇L(Xt ) + ∇V (Xt ) − Tη (xn )∥2 (K−1) 2 ≤ A1 E∥Xt − xn ∥ + 4csplit E∥∇L(Xt ) + ∇V (Xt )∥2 + 2cscore . (K−1) 2 ∥ .

Now we consider E∥Xt − xn

(104)

√ 2(Bt − Bnγ ).

(105)

E∥Xt − x(K−1) ∥2 ≤ 3E∥x(K) − x(K−1) ∥2 + 3γ 2 E∥Tη (x(K−1) )∥2 + 6γd. n n n n

(106)

Xt − x(K−1) = x(K) − x(K−1) − (t − nγ)Tη (x(K−1) )+ n n n n By Lemma 6, (K−1)

(K−1) 2 E∥Tη (xn )∥ ≤

(K−1)

3δ 2 + 3E∥∇L(xn ) + ∇V (xn 2 1 − 3η βL2 (K−1)

≤

3δ 2 + 6(βL + βV )2 E∥xn

)∥2

− Xt ∥2 + 6E∥∇L(Xt ) + ∇V (Xt )∥2 ; 1 − 3η 2 βL2

33

(107)

PiX-MC

Wei et al.

(K−1)

E∥Xt − xn ∥2 (K) (K−1) 2 ≤ 3E∥xn − xn ∥ + 6γd (K−1) 2 2 2 9γ δ + 18γ (βL + βV )2 E∥xn − Xt ∥2 + 18γ 2 E∥∇L(Xt ) + ∇V (Xt )∥2 + . 1 − 3η 2 βL2

(108)

Rearranging this, we obtain (K−1) 2 ∥

E∥Xt − xn (K)

≤ 3ρE∥xn

(K−1) 2 ∥ +ρ

− xn





9γ 2 δ 2 2 + 6γd 1−3η 2 βL

2

18ργ 2 + 1−3η 2 β 2 E∥∇L(Xt ) + ∇V (Xt )∥ .

(109)

L

Thus, (K−1)

E∥∇L(Xt ) + ∇V (Xt ) − Tη (xn ≤

3ρA1 E∥x(K) − x(K−1) ∥2 + A 1 ρ n n



)∥2

 9γ 2 δ 2 + 6γd 1 − 3η 2 βL2

(110)

18ρA1 γ 2 E∥∇L(Xt ) + ∇V (Xt )∥2 + 4csplit E∥∇L(Xt ) + ∇V (Xt )∥2 + 2cscore . + 1 − 3η 2 βL2 The coefficient of E∥∇L(Xt ) + ∇V (Xt )∥2 is A4 , the constant term is A5 . That is, (K−1)

By Lemma 5,

E∥∇L(Xt ) + ∇V (Xt ) − Tη (xn )∥2 (K) (K−1) 2 ≤ 3ρA1 E∥xn − xn ∥ + A4 E∥∇L(Xt ) + ∇V (Xt )∥2 + A5 .

(111)

E∥∇L(Xt ) + ∇V (Xt )∥2 ≤ FI(µt ∥π) + 2d(βL + βV ).

(112)

Thus we have (K−1)

E∥∇L(Xt ) + ∇V (Xt ) − Tη (xn )∥2 (K) (K−1) 2 ≤ 3ρA1 E∥xn − xn ∥ + A4 FI(µt ∥π) + A5 + 2A4 d(βL + βV ).

(113)

Back to (94), (K)

(K−1)

∂t KL(µt ∥π) ≤ − 34 FI(µt ∥π) + 3ρA1 E∥xn − xn ∥2 + A4 FI(µt ∥π) + A5 + 2A4 d(βL + βV )
(K) (K−1) 2 = A4 − 43 FI(µt ∥π) + 3ρA1 E∥xn − xn ∥ + A5 + 2A4 d(βL + βV ).

(114)

By Lemma 7, E∥x(K) − x(K−1) ∥2 ≤ B2 + B3 · FI(µ0 ∥π). n n

(115)

∂t KL(µt ∥π) ≤ (A4 − 0.75)FI(µt ∥π) + A + 3ρA1 B3 · FI(µ0 ∥π),

(116)

Thus, where A = 3ρA1 B2 + A5 + 2A4 d(βL + βV ). Now we integrate over t ∈ [0, T ]: Z T KL(µT ∥π) − KL(µ0 ∥π) ≤ (A4 − 0.75) FI(µt ∥π)dt + AT + 3ρT A1 B3 · FI(µ0 ∥π).

(117)

0

Note that KL(µT ∥π) ≥ 0, thus (0.75 − A4 )

1 T

Z T FI(µt ∥π)dt ≤ 0

1 KL(µ0 ∥π) + A + 3ρA1 B3 · FI(µ0 ∥π). T

This completes the proof. 34

(118)

PiX-MC

C

Wei et al.

Sufficient conditions for Assumption 3

We present Lemma 8 as an example for Assumption 3 to hold. Lemma 8 (Sufficient condition for Assumption 3). Suppose Assumptions 1-2 hold. Also suppose that E∥x0 ∥2 < ∞. Assume additionally that there exist constants κ > αL + αV and b ≥ 0 such that L(x) + V (x) ≥

κ ∥x∥2 − b, 2

(119)

∀x ∈ Rd .

Then, there exist η0 , γ0 > 0 such that, for any 0 < η ≤ η0 and 0 < γ ≤ γ0 satisfying qmult = τ Λη < 1, there exists K0 < ∞ such that, for every K ≥ K0 , the block starting points generated by multi-block PiX-MC satisfy (0) sup E∥xmN ∥2 ≤ Cx (120) m≥0

for some finite constant Cx > 0 independent of the number of blocks M and the total time horizon T . Hence, Assumption 3 holds. Proof. We first show that the quadratic coercivity condition implies dissipativity of the ideal posterior drift. Since L + V is (αL + αV )-weakly convex, for any x ∈ Rd , αL + α V ∥x∥2 . 2

(121)

αL + α V ⟨x, ∇L(x) + ∇V (x)⟩ ≥ L(x) + V (x) − L(0) − V (0) − ∥x∥2 2 κ − α L − αV ≥ ∥x∥2 − b − |L(0) + V (0)|. 2

(122)

L(0) + V (0) ≥ L(x) + V (x) − ⟨∇L(x) + ∇V (x), x⟩ − Combining this inequality with (119) gives

Then, ⟨x, ∇L(x) + ∇V (x)⟩ ≥ a∥x∥2 − b0 , a ≜

κ − α L − αV > 0, b0 ≜ b + |L(0) + V (0)|. 2

(123)

We next transfer this dissipativity property to the proximal drift Tη . By Lemma 1,   2 Tη (x) − ∇L(x) + ∇V (x) ≤ cscore + csplit ∥∇L(x) + ∇V (x)∥2 .

(124)

Since ∇L + ∇V is (βL + βV )-Lipschitz,

Therefore, where

∥∇L(x) + ∇V (x)∥2 ≤ 2∥∇L(0) + ∇V (0)∥2 + 2(βL + βV )2 ∥x∥2 .

(125)

  2 ≤ Aη + Bη ∥x∥2 , Tη (x) − ∇L(x) + ∇V (x)

(126)

Aη ≜ cscore + 2csplit ∥∇L(0) + ∇V (0)∥2 ,

35

Bη ≜ 2csplit (βL + βV )2 .

(127)

PiX-MC

Wei et al.

Since csplit = O(η 2 ), we can choose η0 > 0 sufficiently small such that Bη < a2 whenever 0 < η ≤ η0 . Using (123), (126), and Young’s inequality, we obtain   ⟨x, Tη (x)⟩ = ⟨x, ∇L(x) + ∇V (x)⟩ + x, Tη (x) − ∇L(x) + ∇V (x)   ≥ a∥x∥2 − b0 − ∥x∥ Tη (x) − ∇L(x) + ∇V (x)   2 1 a (128) ≥ a∥x∥2 − b0 − ∥x∥2 − Tη (x) − ∇L(x) + ∇V (x) 2 2a   a2 − Bη Aη ≥ . ∥x∥2 − b0 + 2a 2a Hence, we have ⟨x, Tη (x)⟩ ≥ aη ∥x∥2 − bη , aη ≜

a2 − B η Aη > 0, bη ≜ b0 + . 2a 2a

(129)

We also note that the Lipschitz continuity of Tη directly implies a linear-growth bound. Indeed, by Lemma 4, ∥Tη (x)∥2 ≤ 2∥Tη (0)∥2 + 2Λ2η ∥x∥2 . (130) We now consider one block of multi-block PiX-MC. Fix m ≥ 0, and define a sequential reference em,0 = xmN = trajectory {e xm,i }N i=0 using the same Brownian increments as the m-th Picard block. Set x (0) xmN , and, for i = 0, . . . , N − 1, √
em,i+1 = x em,i − γTη (e x xm,i ) + 2 B(mN +i+1)γ − B(mN +i)γ . (131) This is precisely the sequential X-MC trajectory over the current block. em,i and has Since the Brownian increment over [(mN + i)γ, (mN + i + 1)γ] is independent of x zero mean, by (129) and (130),

E∥e xm,i+1 ∥2 ≤ 1 − 2aη γ + 2γ 2 Λ2η E∥e xm,i ∥2 + γ 2bη + 2γ∥Tη (0)∥2 + 2d . (132) Choose γ0 > 0 sufficiently small such that, for 0 < γ ≤ γ0 , define 2

rη,γ ≜ 1 − 2aη γ + 2γ Λ2η ∈ (0, 1), Cseq ≜

γ 2bη + 2γ∥Tη (0)∥2 + 2d . 1 − rη,γ

(133)

Iterating the preceding inequality gives, for every i = 0, . . . , N , i E∥e xm,i ∥2 ≤ rη,γ E∥xmN ∥2 + Cseq ≤ E∥xmN ∥2 + Cseq .

(134)

It remains to control the finite-Picard approximation of this sequential block. For convenience, set (k) xmN = xmN for all k = 0, . . . , K, and define (k)

em,i Dm,k ≜ max E xmN +i − x 0≤i≤N

2

.

(135)

The sequential reference trajectory satisfies the fixed-point representation em,i = xmN − γ x

i−1 X

Tη (e xm,j ) +

j=0

36

√

2 B(mN +i)γ − BmN γ .

(136)

PiX-MC

Wei et al.

Comparing this identity with the Picard update in Algorithm 2, we obtain, for i = 1, . . . , N , (k+1)

em,i E xmN +i − x

2

≤ iγ 2 ≤

i−1 X

(k)

E Tη (xmN +j ) − Tη (e xm,j )

j=0 i−1 X iγ 2 Λ2η E j=0

2

2 (k) 2 em,j ≤ τ 2 Λ2η Dm,k = qmult Dm,k . xmN +j − x

(137)

Taking the maximum over i yields the following contraction that comes directly from the Lipschitz continuity of Tη : 2 Dm,k+1 ≤ qmult Dm,k . (138) (0)

Since the zeroth Picard trajectory satisfies xmN +i = xmN , (134) gives Dm,0 ≤ 2E∥xmN ∥2 + 2 max E∥e xm,i ∥2 ≤ 4E∥xmN ∥2 + 2Cseq .

(139)

2K 4E∥xmN ∥2 + 2Cseq . Dm,K ≤ qmult

(140)

0≤i≤N

Therefore, by (138),

(K)

Finally, Algorithm 2 sets x(m+1)N = x(m+1)N . For any υ > 0, Young’s inequality, together with (134) for i = N and (140), gives (K)

em,N E∥x(m+1)N ∥2 ≤ (1 + υ)E∥e xm,N ∥2 + (1 + υ −1 )E x(m+1)N − x   N −1 2K 2 ≤ (1 + υ)rη,γ + 4(1 + υ )qmult  E∥xmN ∥ 2K C . + (1 + υ) + 2(1 + υ −1 )qmult seq

2

(141)

N < 1, choose υ > 0 sufficiently small such that (1 + υ)r N < 1. Since q Since rη,γ mult < 1, there exists a η,γ 2K0 N −1 finite K0 such that (1 + υ)rη,γ + 4(1 + υ )qmult < 1. Hence, for every K ≥ K0 , there exist constants χ ∈ (0, 1) and C < ∞, independent of the block index m, such that

(142)

E∥x(m+1)N ∥2 ≤ χE∥xmN ∥2 + C. Iterating (142) over the block index gives E∥xmN ∥2 ≤ χm E∥x0 ∥2 + C

m−1 X j=0

Therefore, sup E∥xmN ∥2 ≤ E∥x0 ∥2 + m≥0

χj ≤ E∥x0 ∥2 +

C . 1−χ

C ≜ Cx < ∞. 1−χ

The constant Cx does not depend on M or T = M τ . This proves Assumption 3.

37

(143)

(144)

PiX-MC

D

Wei et al.

Proof of multi-block PiX-MC

We present the convergence result of multi-block PiX-MC. Before that, we define some essential constants. Recall the constants ρ, A1 , A4 , and A5 defined in Appendix B. Under Assumption 3, define G2 ≜ 2∥∇L(0) + ∇V (0)∥2 + 2(βL + βV )2 Cx .

(145)

Since ∇L + ∇V is (βL + βV )-Lipschitz, Assumption 3 yields h i sup E ∥∇L (xmN ) + ∇V (xmN )∥2 ≤ G2 .

(146)

0≤m≤M −1

We define the constant D. 2K−2 D = 3ρA1 qmult R + A5 + 2A4 d(βL + βV ).

(147)

Recall that qmult = τ Λη is the in-block contraction factor introduced in Section 3.3. Theorem 5 (Detailed version of Theorem 2). Let (µt )t≥0 denote the law of the continuous interpolation of {xn (K) }N n=0 generated by multi-block PiX-MC. Suppose Assumptions 1-3 hold. Let the parameters η, γ be chosen sufficiently small such that 1 − ηαL > 0, 1 − 3η 2 βL2 > 0, 1 − Then, 1 (0.75 − A4 ) T

18γ 2 (βL + βV )2 > 0, qmult < 1, A4 < 0.75. 1 − 3η 2 βL2

(148)

KL(µ0 ∥π) + D. T

(149)

Z T FI(µt ∥π)dt ≤ 0

The conditions involving ρ and A4 are identical to those in Theorem 4. Before proving Theorem 5, we establish the following detailed version of Lemma 3, which shows that the error decays geometrically fast in the parallel refinement. Lemma 9 (Detailed version of Lemma 3). Let Assumptions 1-3 hold. Choose the parameter η, γ > 0 sufficiently small such that 1 − ηαL > 0,

1 − 3η 2 βL2 > 0,

qmult = τ Λη < 1.

(150)

For each block m = 0, . . . , M − 1 and each Picard refinement k = 1, . . . , K, define (k)

(k−1)

εm,k ≜ max E xmN +i − xmN +i 1≤i≤N

2

.

2 Then for k = 1, . . . , K − 1, εm,k+1 ≤ qmult εm,k . Moreover "
# 2 δ 2 + G2 3τ 2K−2 2K−2 εm,K ≤ qmult 2τ d + = qmult R. 1 − 3η 2 βL2

Consequently, max

0≤m≤M −1 1≤i≤N

(K−1) 2

(K)

E xmN +i − xmN +i 38

2K−2 ≤ qmult R.

(151)

(152)

(153)

PiX-MC

Wei et al.

Proof. Fix a block index m ∈ {0, 1, . . . , M − 1}. For any k = 1, . . . , K − 1 and i = 1, . . . , N , multi-block PiX-MC gives (k+1) (k) E∥xmN +i − xmN +i ∥2 =

E γ

≤ iγ 2 ≤ ≤

i−1 Ph

j=0 i−1 P

i 2

(k) (k−1) Tη (xmN +j ) − Tη (xmN +j ) (k)

(k−1)

E∥Tη (xmN +j ) − Tη (xmN +j )∥2

(154)

j=0 i−1 P (k) (k−1) iγ 2 Λ2η E∥xmN +j − xmN +j ∥2 j=0 2 τ 2 Λ2η εm,k = qmult εm,k .

2 Thus, εm,k+1 ≤ qmult εm,k . It remains to bound the initial Picard error εm,1 . Since the zeroth Picard (0) trajectory is initialized at the block starting point, xmN +i = xmN for i = 0, . . . , N , we have, for i = 1, . . . , N , i−1 X √
(1) (0) (155) xmN +i − xmN +i = −γ Tη (xmN ) + 2 B(mN +i)γ − BmN γ . j=0

Note that (1)

(0)

E xmN +i − xmN +i

2

= i2 γ 2 E ∥Tη (xmN )∥2 + 2iγd ≤ τ 2 E ∥Tη (xmN )∥2 + 2τ d.

(156)

By Lemma 6, we have E ∥Tη (xmN )∥2 ≤

3 3δ 2 + E ∥∇L (xmN ) + ∇V (xmN )∥2 , 1 − 3η 2 βL2 1 − 3η 2 βL2

(157)

and that E ∥∇L (xmN ) + ∇V (xmN )∥2 ≤ G2 . As a result,
3 δ 2 + G2 E ∥Tη (xmN )∥ ≤ . 1 − 3η 2 βL2 2

2

2

(158)

2

(δ +G ) Thus εm,1 ≤ 2τ d + 3τ1−3η = R. We finally obtain 2β2 L

(159)

2K−2 2K−2 εm,K ≤ qmult εm,1 ≤ qmult R.

Since R is independent of the total time horizon T and the bound is uniform over m = 0, . . . , M − 1, the global estimate follows. This completes the proof. Now we are ready to prove Theorem 5. Proof. Let m ∈ {0, . . . , M − 1} and i ∈ {0, . . . , N − 1}. For t ∈ [(mN + i)γ, (mN + i + 1)γ], define the continuous-time interpolation by   √ (K) (K−1) Xt = xmN +i − [t − (mN + i)γ] Tη xmN +i + 2(Bt − B(mN +i)γ ). (160) (K)

Let µt = Law(Xt ). By multi-block PiX-MC, this interpolation satisfies X(mN +i)γ = xmN +i and (K)

X(mN +i+1)γ = xmN +i+1 . 39

PiX-MC

Wei et al.

Applying the same conditional-drift argument as in the proof of Theorem 4, we obtain, for almost every t ∈ [(mN + i)γ, (mN + i + 1)γ], 3 (K−1) 2 ∂t KL(µt ∥π) ≤ − FI(µt ∥π) + E ∇L(Xt ) + ∇V (Xt ) − Tη (xmN +i ) . 4

(161)

The local interpolation estimate established in the proof of Theorem 4 also applies here without modification, and yields   2 (K−1) E ∇L(Xt ) + ∇V (Xt ) − Tη xmN +i (K−1) 2

(K)

≤ 3ρA1 E xmN +i − xmN +i

(162)

+ A4 E∥∇L(Xt ) + ∇V (Xt )∥2 + A5 .

By Lemma 5, we have E∥∇L(Xt ) + ∇V (Xt )∥2 ≤ FI(µt ∥π) + 2d(βL + βV ). (K)

(K−1)

For i = 0, we have xmN = xmN obtain

(163)

= xmN . Hence for i = 0, . . . , N − 1, together with Lemma 9, we (K−1) 2

(K)

E xmN +i − xmN +i

2K−2 R. ≤ qmult

(164)

Substituting this back, we arrive at 2K−2 ∂t KL(µt ∥π) ≤ (A4 − 0.75) FI(µt ∥π) + 3ρA1 qmult R + A5 + 2A4 d(βL + βV ).

By the definition of D, ∂t KL(µt ∥π) ≤ (A4 − 0.75) FI(µt ∥π) + D.

(165) (166)

Integrating (166) over each interval [(mN + i)γ, (mN + i + 1)γ] and summing over m = 0, . . . , M − 1 and i = 0, . . . , N − 1, we obtain Z T KL(µT ∥π) − KL(µ0 ∥π) ≤ (A4 − 0.75) FI(µt ∥π)dt + DT. (167) 0

Since KL(µT ∥π) ≥ 0, rearranging and dividing by T gives (0.75 − A4 )

1 T

Z T FI(µt ∥π)dt ≤ 0

KL(µ0 ∥π) + D. T

(168)

This completes the proof.

E

Proof of multi-block APiX-MC

We present the detailed convergence result for multi-block APiX-MC with weighted annealing. We first introduce several quantities specific to the annealed analysis. Recall the constants ρ, A1 , and A4 defined in Appendix B. Define the maximal and averaged score error by δann ≜

max

{δn + (αn − 1)RS },

0≤n≤M N −1

2 δ̄ann ≜

MX N −1 1 [δn + (αn − 1)RS ]2 . MN n=0

40

(169)

PiX-MC

Wei et al.

Assumption 4 gives ∥αn Sθ (x; σn ) − ∇V (x)∥ ≤ ∥Sθ (x; σn ) − ∇V (x)∥ + (αn − 1)∥Sθ (x; σn )∥ ≤ δn + (αn − 1)RS . Since ∇L + ∇V is (βL + βV )-Lipschitz, Assumption 3 gives h i sup E ∥∇L (xmN ) + ∇V (xmN )∥2 ≤ G2 .

(170)

(171)

0≤m≤M −1 N −1 Since {αn }M is non-increasing, define n=0

Λη,ann ≜ α0 LS + βL

1 + ηα0 LS , 1 − ηαL

qann ≜ τ Λη,ann ,

Rann ≜ 2τ d +

2 + G2 ) 3τ 2 (δann , 1 − 3η 2 βL2

(172)

where G is defined above through the uniform moment bound in Assumption 3. Finally define the constant   4 + 9ρA1 γ 2 2 2K−2 δ̄ann + 6ρA1 γd + 2A4 d(βL + βV ). (173) Dann = 3ρA1 qann Rann + 1 − 3η 2 βL2 Theorem 6 (Detailed version of Theorem 3). Let (µt )t≥0 denote the law of the continuous interpolation of {xn (K) }N n=0 generated by multi-block APiX-MC. Suppose Assumptions 1, 3, and 4 hold. Choose parameters η, γ > 0 sufficiently small such that 1 − ηαL > 0,

1 − 3η 2 βL2 > 0,

Then, 1 (0.75 − A4 ) T

1−

18γ 2 (βL + βV )2 > 0, 1 − 3η 2 βL2

Z T FI(µt ∥π)dt ≤ 0

qann < 1,

A4 < 0.75.

KL(µ0 ∥π) + Dann . T

(174)

(175)

The effective discrepancy satisfies 2 δ̄ann ≤

MX N −1 M N −1 2RS2 X 2 δn2 + (αn − 1)2 . MN MN n=0

(176)

n=0

Hence, the first term on the right captures the score approximation error, while the second summation represents the transient error introduced by weighted annealing. In particular, this annealing error vanishes as the post-annealing sampling horizon increases, provided that αn = 1 after a finite number of iterations. Before proving Theorem 6, we establish the following annealed counterpart of Lemma 9. It shows that the in-block Picard refinement error continues to decay geometrically. Lemma 10. Let Assumptions 1, 3, and 4 hold. Choose η, γ small enough such that 1 − ηαL > 0,

1 − 3η 2 βL2 > 0,

qann = τ Λη,ann < 1.

(177)

For each block m = 0, . . . , M − 1 and each Picard refinement k = 1, . . . , K, define (k)

(k−1)

εm,k := max E xmN +i − xmN +i 1≤i≤N

41

2

.

(178)

PiX-MC

Wei et al.

2 ε Then, for k = 1, . . . , K − 1, εm,k+1 ≤ qann m,k . Moreover, "
# 2 + G2 3τ 2 δann 2K−2 2K−2 = qann Rann . εm,K ≤ qann 2τ d + 1 − 3η 2 βL2

Consequently,

(K−1) 2

(K)

max

E xmN +i − xmN +i

(K)

(K−1)

0≤m≤M −1 1≤i≤N

At each block starting point, xmN = xmN nodes n = 0, . . . , M N − 1. max

0≤n≤M N −1

(179)

(180)

2K−2 ≤ qann Rann .

= xmN . Hence, the same bound also holds over all global

E x(K) − x(K−1) n n

2

(181)

2K−2 ≤ qann Rann .

Proof. Fix a block index m ∈ {0, . . . , M − 1}. For any k = 1, . . . , K − 1 and i = 1, . . . , N , multi-block APiX-MC gives (k+1)

(k)

xmN +i − xmN +i = −γ

i−1 h    i X α αmN +j (k) (k−1) +j Tη,σmN x − T x . η,σmN +j mN +j mN +j mN +j

(182)

j=0

By the Cauchy-Schwarz inequality and the uniform Lipschitz bound of the weighted drifts, (k+1)

(k)

E xmN +i − xmN +i

2

≤ iγ 2

i−1 P j=0

≤

  2   α α (k−1) (k) +j +j xmN +j xmN +j − Tη,σmN E Tη,σmN mN +j mN +j

iγ 2 Λ2η,ann

i−1 P

E

j=0

(k) (k−1) 2 xmN +j − xmN +j .

(183)

(k)

Since xmN = xmN holds for all k, (k+1)

(k)

E xmN +i − xmN +i

2

2 ≤ i2 γ 2 Λ2η,ann εm,k ≤ τ 2 Λ2η,ann εm,k = qann εm,k .

(184)

2 ε Taking the maximum over i = 1, . . . , N yields εm,k+1 ≤ qann m,k . It remains to bound εm,1 . Since the (0) zeroth Picard trajectory is initialized by xmN +i = xmN for i = 0, . . . , N , we have (1)

(0)

xmN +i − xmN +i = −γ

i−1 X

α

+j Tη,σmN mN +j (xmN ) +

√

2 B(mN +i)γ − BmN γ .

(185)

j=0

Since xmN is independent of the Brownian increments within the current block, the cross term vanishes. Hence, E

2 (1) (0) xmN +i − xmN +i = γ 2 E

i−1 X

2 α +j Tη,σmN mN +j (xmN )

+ 2iγd.

(186)

j=0

Applying Cauchy-Schwarz inequality gives γ2E

i−1 X

2 α +j Tη,σmN mN +j (xmN )

≤ iγ 2

j=0

i−1 X j=0

42

α

+j E Tη,σmN mN +j (xmN )

2

.

(187)

PiX-MC

Wei et al.

Applying the argument of Lemma 6 to the weighted score αmN +j Sθ (·; σmN +j ), we obtain α

+j E Tη,σmN mN +j (xmN )

By (171)

2

≤

2 3δann 3 + E ∥∇L(xmN ) + ∇V (xmN )∥2 . 2 2 2 2 1 − 3η βL 1 − 3η βL

(188)

E ∥∇L (xmN ) + ∇V (xmN )∥2 ≤ G2 .

(189)

2 + G2 δann . 1 − 3η 2 βL2

(190)

2 + G2 3τ 2 δann εm,1 ≤ 2τ d + = Rann . 1 − 3η 2 βL2

(191)

Thus E

3 2 α +j Tη,σmN ≤ mN +j (xmN )

Consequently, since iγ ≤ τ ,

Combining this estimate with the geometric recursion gives (192)

2K−2 2K−2 εm,K ≤ qann εm,1 ≤ qann Rann .

Since the bound is uniform over m = 0, . . . , M − 1, the global estimate follows. This completes the proof. Now we are ready to prove Theorem 6. Proof. For each global time node n = 0, . . . , M N − 1 and t ∈ [nγ, (n + 1)γ], define the continuous-time interpolation by   √ αn (K−1) Xt = x(K) − (t − nγ)T (193) x + 2 (Bt − Bnγ ) . n η,σn n (K)

(K)

Let µt = Law(Xt ). By multi-block APiX-MC, Xnγ = xn , X(n+1)γ = xn+1 . Applying the same conditional-drift argument as in the proofs of Theorems 4 and 5, we obtain, for almost every t ∈ [nγ, (n + 1)γ], 2 3 αn (x(K−1) ) . ∂t KL(µt ∥π) ≤ − FI(µt ∥π) + E ∇L(Xt ) + ∇V (Xt ) − Tη,σ n n 4

(194)

At the n-th node, Assumption 4 gives ∥αn Sθ (x; σn ) − ∇V (x)∥ ≤ δn + (αn − 1) RS .

(195)

Therefore, repeating the local interpolation estimate in the proof of Theorem 4, with δ replaced by δn + (αn − 1) RS gives us   2 (K−1) αn E ∇L (Xt ) + ∇V (Xt ) − Tη,σ x n n (K)

(K−1) 2

≤ 3ρA1 E xn − xn + A4 E ∥∇L (Xt ) + ∇V (Xt )∥2   4 + 9ρA1 γ 2 + [δn + (αn − 1) RS ]2 + 6ρA1 γd. 1 − 3η 2 βL2

43

(196)

PiX-MC

Wei et al.

By Lemma 5 and Lemma 10, ∂t KL(µt ∥π) ≤

2K−2 A4 − 34 FI(µt ∥π)  + 3ρA1 qann Rann 2 4 + 9ρA1 γ + [δn + (αn − 1) RS ]2 + 6ρA1 γd + 2A4 d (βL + βV ) 1 − 3η 2 βL2

(197)

Integrating over [nγ, (n + 1)γ], and summing over n = 0, . . . , M N − 1 yields KL (µT ∥π) − KL (µ0 ∥π)
RT 2K−2 R ≤ A4 − 43 0 FI (µt ∥π) dt + 3ρA1 qann ann T   MX N −1 2 4 + 9ρA1 γ + γ [δn + (αn − 1) RS ]2 + [6ρA1 γd + 2A4 d (βL + βV )] T. 1 − 3η 2 βL2 n=0 2 , we have γ By the definition of δ̄ann

Dann ,

MP N −1 n=0

(198)

2 . Further, by the definition of [δn + (αn − 1) RS ]2 = T δ̄ann

 KL (µT ∥π) − KL (µ0 ∥π) ≤

3 A4 − 4

Z T FI (µt ∥π) dt + T Dann .

(199)

0

Finally by the non-negativity of KL, the proof is complete.

F

Additional experimental details

F.1

Common Implementation Details and Hyperparameter Selections

In all experiments, the images are scaled to [0, 1]. Before evaluating a score network, we rescale the input image to [−1, 1] and map the network output back to [0, 1] accordingly. The pretrained models used for natural images, MRI, and CT adopt different output parameterizations, whose model-specific implementations are described below. Score-network implementation. For the natural-image experiments, we employ a pretrained diffusion model as in DDPM [36]. Let ᾱt denote the cumulative diffusion coefficient. The connection between the score network Sθ and noise predictor εθ in DDPM can be formulated approximately as   2 2x − 1 Sθ (x; σ) = εθ √ ;σ , (200) σ 1 + σ2 r 1 − ᾱt where σ = . For MRI reconstruction, we employ the pretrained score network from PMC ᾱt [86]. Denote its output in the network coordinates by SθMRI (2x − 1; σ). Since the network directly estimates the positive score, the negative-score function used in our formulation is Sθ (x; σ) = −2SθMRI (2x − 1; σ).

(201)

For CT images, we train an MMSE denoiser Dθ using the standard EDM framework [44]. The relation between Sθ and Dθ is 2 Sθ (x; σ) = 2 [2x − 1 − Dθ (2x − 1; σ)] . (202) σ 44

PiX-MC

Wei et al.

Balancing parameter for likelihood. To account for the different numerical scales of the likelihood and score networks, we introduce a task-dependent balancing parameter λ > 0 for the likelihood. For example, MRI, deblurring, and CT, the likelihood is L(x) =

λ ∥Ax − y∥2 , 2

(203)

where A is a task-specific forward operator. Annealing Schedule. Let n = mN + i denote the global time node index, where m = 0, . . . , M − 1 is the block index and i = 0, . . . , N − 1 is the local node index. Recall that we use the geometric schedule σn = max{σmax ξ n , σmin }, (204) where 0 < ξ < 1 is the annealing parameter. For the score network weight schedule αn , we also use a corresponding geometric sequence αn = (σn /σmin )a (205) This schedule is non-increasing and becomes one once σn reaches σmin . For the MRI and CT experiments, we set a = 2; for the deblurring and Rician denoising experiments, we set a = 5/2. Picard stopping criteria. For deblurring and CT experiments, we employ adaptive Picard refinement. We compute the empirical in-block residual εm,k = max

1≤i≤N

(k)

(k−1)

xmN +i − xmN +i .

(206)

The refinement within a block is terminated when either k reaches K = 20 or the normalized residual εm,k /d < 3 × 10−4 . Hyperparameter selection. The hyperparameters are tuned separately for each method and imaging task to obtain the best reconstruction results. The gradient-based and proximal variants use different drift parameterizations and are therefore assigned method-specific stepsizes. In particular, the implicit treatment of the likelihood in the proximal variants may improve numerical stability and permit larger stepsizes in practice. For each paired sequential and Picard-based method, we use the same Langevin stepsize, score schedule, likelihood weight, and total number of time nodes. The complete hyperparameter settings are reported in Table 4.

F.2

Likelihood and Proximal Implementations

We next provide the task-specific implementations of the likelihood gradient and proximal operator. The datasets, score networks, and hyperparameters are reported in $ 5 and Appendix F.1. MRI reconstruction. The forward operator A = M F consists of the discrete Fourier transform F followed by a radial subsampling mask. The likelihood and its gradient are L(x) =

λ ∥M F x − y∥2 , 2

∇L(x) = λF H M T (M F x − y) .

45

(207)

PiX-MC

Wei et al.

Inverse problem Method

λ

γ

η

ξ

σmin

σmax

M

N

K

MRI

PiL-MC – APiL-MC 1.6 × 103 6.25 × 10−5 – PiX-MC 3.5 × 103 8.75 × 10−5 4.375 × 10−4 APiX-MC 3.5 × 103 8.75 × 10−5 4.375 × 10−4

1 0.98 1 0.98

0.05 0.05 0.05 0.05

0.05 0.2 0.05 0.2

10 10 10 10

75 75 75 75

100 100 100 100

Rician

PiL-MC APiL-MC PiX-MC APiX-MC

2.2 2.2 2.2 2.2

1 × 10−4 1 × 10−4 1 × 10−4 1 × 10−4

– – 1 × 10−3 1 × 10−3

1 0.98 1 0.98

0.1 0.1 0.1 0.1

0.1 0.4 0.1 0.4

5 5 5 5

100 100 100 100

50 50 50 50

Deblurring

PiL-MC APiL-MC PiX-MC APiX-MC

800 1100 900 1150

1 × 10−4 1 × 10−4 1 × 10−4 1 × 10−4

– – 1 × 10−3 1 × 10−3

1 0.98 1 0.98

0.1 0.1 0.1 0.1

0.1 2.0 0.1 2.0

125 125 125 125

8 8 8 8

20 20 20 20

3D CT

PiL-MC APiL-MC PiX-MC APiX-MC

1.11 1.11 2.22 2.22

1 × 10−5 1 × 10−5 4.5 × 10−5 4.5 × 10−5

– – 9 × 10−4 9 × 10−4

1 0.03 0.03 3875 0.995 0.03 0.5 3875 1 0.085 0.085 3875 0.999 0.085 0.5 3875

8 8 8 8

20 20 20 20

1.6 × 103

6.25 × 10−5

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

The proximal operator has the closed form proxηL (z) = ηλF H M T M F + I

−1

z + ηλF H M T y .

(208)

Since F is unitary and M T M is diagonal in the Fourier domain, the inverse above can be evaluated elementwise in k-space, without the need to solve a large-scale linear system. Rician denoising. Let ereal and eimag be two i.i.d. Gaussian noise vectors with zero mean and standard deviation σy . The Rician-noise-corrupted observation y is obtained by: q y = (x + ereal )2 + e2imag , (209) where theoperations are elementwise. According to [33], the conditional probability density function (PDF) of y is     d Y yr x2r + yr2 x r yr p(y|x) = exp − I , (210) 0 σy2 2σy2 σy2 r=1

which is known as the Rice or Rician distribution. Here I0 is the modified Bessel function of the first kind of order zero. Consider the modified Bessel’s differential equations: x2 y ′′ + xy ′ − (x2 + n2 )y = 0,

n ≥ 0.

The solutions of the previous equation In (x) = y(x) continuous in zero are called the modified Bessel function of the first kind of order n. After omitting terms independent of x, the negative log-likelihood L for Rician noise removal is    x⊙y x2 . (211) L(x) = λ 1, 2 − log I0 2σy σy2 46

PiX-MC

Wei et al.

Note that according to Lemma 12 in [69], L is Lipschitz smooth, thus satisfies Assumption 1. Following [90], its gradient is evaluated using the ratio B(x) = I1 (x)/I0 (x), and the pixel-wise proximal subproblem is solved using the iterative reweighted ℓ1 (IRL1) method. The IRL1 algorithm for solving proxηL (z) is given as follows:   t y x ⊙y t+ 21 x = 2B 2 σy σ y  1 x2 1 x⊙y t+1 x = arg min ηλ 1, 2 − + ∥x − z∥2 − ηλ⟨xt+ 2 , x⟩ 2 x 2σy σy 2 (212) ηλ t+ 21 z + 2 y + ηλx σy . = ηλ 1+ 2 σy In experiments, we set the maximum inner iteration number for IRL1 to be 10 because it typically solves the proximal sub-problem within 5 iterations. Image deblurring. Following [69], we model the forward operator as a circular convolution, which is diagonalized by the discrete Fourier transform. The likelihood is L(x) = λ2 ∥Kx − y∥2 , where K is the circular convolution operator associated with a blur kernel k. The likelihood gradient and proximal operator are
−1 (z + ηλK T y). (213) ∇L(x) = λK T (Kx − y), proxηL (z) = ηλK T K + I The proximal operator can be evaluated efficiently in the Fourier domain through elementwise multiplication and division:   F(z + ηλK T y) proxηL (z) = F H . (214) F(ηλK T K + 1) 3D CT reconstruction. The quadratic likelihood potential and its gradient are given by L(x) =

λ ∥Ax − y∥2 , ∇L(x) = λAT (Ax − y). 2 The proximal operator is evaluated by solving the following linear system:

(215)

x = proxηL (z) ⇐⇒ (I + ηλAT A)x = z + ηλAT y.

(216)

We solve this linear system using five conjugate gradient iterations [76]. To promote consistency between adjacent slices, we additionally use a one-dimensional inter-slice Huber-TV (HTV) regularizer. For x ∈ R512×512×80 , HTV regularizer H(x) with smoothing parameter δH > 0 is H(x) =

512 X 512 X 79 X

HδH (x(i, j, s + 1) − x(i, j, s)),

(217)

i=1 j=1 s=1

where x(i, j, s) denotes the element located at i-th row, j-th column in s-th slice. For any x ∈ R, HδH (x) is defined as  x2   , if |x| ≤ δH ; 2δH HδH (x) := (218)   |x| − δH , if |x| > δ . H 2 47

PiX-MC

Wei et al.

We set δH = 10−8 in the experiments. We view HTV as an inter-slice prior for CT volumes characterizing correlations between adjacent slices. Accordingly, the proximal drift Tη,σ is Tη,σ ≜

 1 I − proxηL ◦ I − ηSθ (·; σ) − ηβH ∇H , η

(219)

where βH > 0 is a balancing parameter for H. We set βH = 0.28 for all methods.

References [1] M. V. Afonso, J. M.Bioucas-Dias, and M. A. T. Figueiredo. Fast image recovery using variable splitting and constrained optimization. IEEE Trans. Image Process., 19(9):2345–2356, Sep. 2010. [2] Eirikur Agustsson and Radu Timofte. Ntire 2017 challenge on single image super-resolution: Dataset and study. In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR) Workshops, July 2017. [3] Nima Anari, Sinho Chewi, and Thuy-Duong Vuong. Fast parallel sampling under isoperimetry. In The Thirty Seventh Annual Conference on Learning Theory, pages 161–185. PMLR, 2024. [4] Krishna Balasubramanian, Sinho Chewi, Murat A Erdogdu, Adil Salim, and Shunshi Zhang. Towards a theory of non-log-concave sampling: first-order stationarity guarantees for langevin monte carlo. In Conference on Learning Theory, pages 2896–2923. PMLR, 2022. [5] A. Beck and M. Teboulle. Convex Optimization in Signal Processing and Communications, chapter Gradient-Based Algorithms with Applications to Signal Recovery Problems, pages 42–88. Cambridge, 2009. [6] A. Beck and M. Teboulle. Fast gradient-based algorithm for constrained total variation image denoising and deblurring problems. IEEE Trans. Image Process., 18(11):2419–2434, November 2009. [7] Edmon Begoli, Tanmoy Bhattacharya, and Dimitri Kusnezov. The need for uncertainty quantification in machine-assisted medical decision making. Nature Machine Intelligence, 1(1):20–23, 2019. [8] Charles A Bouman and Gregery T Buzzard. Generative plug and play: Posterior sampling for inverse problems. In 2023 59th annual Allerton conference on communication, control, and computing (Allerton), pages 1–7. IEEE, 2023. [9] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributed optimization and statistical learning via the alternating direction method of multipliers. Found. Trends Mach. Learn., 3(1):1– 122, 2011. [10] Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng. Handbook of markov chain monte carlo. CRC press, 2011. [11] E. J. Candès, J. Romberg, and T. Tao. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. IEEE Trans. Inf. Theory, 52(2):489–509, February 2006. 48

PiX-MC

Wei et al.

[12] Haoxuan Chen, Yinuo Ren, Lexing Ying, and Grant M Rotskoff. Accelerating diffusion models with parallel sampling: Inference at sub-linear time complexity. Advances in Neural Information Processing Systems, 37:133661–133709, 2024. [13] Sitan Chen, Sinho Chewi, Holden Lee, Yuanzhi Li, Jianfeng Lu, and Adil Salim. The probability flow ode is provably fast. Advances in Neural Information Processing Systems, 36:68552–68575, 2023. [14] Sinho Chewi, Murat A Erdogdu, Mufan Li, Ruoqi Shen, and Matthew S Zhang. Analysis of langevin monte carlo from poincare to log-sobolev. Foundations of Computational Mathematics, 25(4):1345–1395, 2025. [15] Hyungjin Chung, Jeongsol Kim, Sehui Kim, and Jong Chul Ye. Parallel diffusion models of operator and image for blind inverse problems. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 6059–6069, 2023. [16] Hyungjin Chung, Jeongsol Kim, Michael T Mccann, Marc L Klasky, and Jong Chul Ye. Diffusion posterior sampling for general noisy inverse problems. arXiv preprint arXiv:2209.14687, 2022. [17] Hyungjin Chung and Jong Chul Ye. Score-based diffusion models for accelerated mri. Medical image analysis, 80:102479, 2022. [18] Florentin Coeurdoux, Nicolas Dobigeon, and Pierre Chainais. Normalizing flow sampling with langevin dynamics in the latent space. Machine Learning, 113(11):8301–8326, 2024. [19] Florentin Coeurdoux, Nicolas Dobigeon, and Pierre Chainais. Plug-and-play split gibbs sampler: embedding deep generative priors in bayesian inference. IEEE Transactions on Image Processing, 33:3496–3507, 2024. [20] Prafulla Dhariwal and Alexander Quinn Nichol. Diffusion models beat GANs on image synthesis. In Advances in Neural Information Processing Systems, 2021. [21] Zhiyan Ding, Qin Li, Jianfeng Lu, and Stephen J Wright. Random coordinate langevin monte carlo. In Conference on learning theory, pages 1683–1710. PMLR, 2021. [22] D. L. Donoho. Compressed sensing. IEEE Trans. Inf. Theory, 52(4):1289–1306, April 2006. [23] Zolnamar Dorjsembe, Hsing-Kuo Pao, Sodtavilan Odonchimed, and Furen Xiao. Conditional diffusion models for semantic 3d brain mri synthesis. IEEE Journal of Biomedical and Health Informatics, 28(7):4084–4093, 2024. [24] Alain Durmus, Szymon Majewski, and Błażej Miasojedow. Analysis of langevin monte carlo via convex optimization. Journal of Machine Learning Research, 20(73):1–46, 2019. [25] Alain Durmus, Eric Moulines, and Marcelo Pereyra. Efficient bayesian computation by proximal markov chain monte carlo: when langevin meets moreau. SIAM Journal on Imaging Sciences, 11(1):473–506, 2018. [26] Bradley Efron. Tweedie’s formula and selection bias. Journal of the American Statistical Association, 106(496):1602–1614, 2011.

49

PiX-MC

Wei et al.

[27] Elhadji C Faye, Mame Diarra Fall, and Nicolas Dobigeon. Regularization by denoising: Bayesian model and langevin-within-split gibbs sampling. IEEE Transactions on Image Processing, 34:221– 234, 2024. [28] Berthy T Feng, Jamie Smith, Michael Rubinstein, Huiwen Chang, Katherine L Bouman, and William T Freeman. Score-based diffusion models as principled priors for inverse imaging. In Proceedings of the IEEE/CVF international conference on computer vision, pages 10520–10531, 2023. [29] Olivier Fercoq and Peter Richtárik. Accelerated, parallel, and proximal coordinate descent. SIAM Journal on Optimization, 25(4):1997–2023, 2015. [30] Rémi Gribonval and Mila Nikolova. A characterization of proximity operators. Journal of Mathematical Imaging and Vision, 62(6):773–789, 2020. [31] Leonard Gross. Logarithmic sobolev inequalities. American Journal of Mathematics, 97(4):1061– 1083, 1975. [32] Shuhang Gu, Lei Zhang, Wangmeng Zuo, and Xiangchu Feng. Weighted nuclear norm minimization with application to image denoising. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2862–2869, 2014. [33] Hákon Gudbjartsson and Samuel Patz. The rician distribution of noisy mri data. Magnetic resonance in medicine, 34(6):910–914, 1995. [34] Robert Hannah, Fei Feng, and Wotao Yin. A2bcd: Asynchronous acceleration with optimal complexity. In International Conference on Learning Representations, 2019. [35] Chunming He, Yuqi Shen, Chengyu Fang, Fengyang Xiao, Longxiang Tang, Yulun Zhang, Wangmeng Zuo, Zhenhua Guo, and Xiu Li. Diffusion models in low-level vision: A survey. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2025. [36] Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in neural information processing systems, 33:6840–6851, 2020. [37] Jonathan Ho and Tim Salimans. Classifier-free diffusion guidance. In NeurIPS 2021 Workshop on Deep Generative Models and Downstream Applications, 2021. [38] Matthew Holden, Marcelo Pereyra, and Konstantinos C Zygalakis. Bayesian imaging with datadriven priors encoded by neural networks. SIAM Journal on Imaging Sciences, 15(2):892–924, 2022. [39] Yuanyun Hu, Evan Bell, Guijin Wang, and Yu Sun. Prism: Probabilistic and robust inverse solver with measurement-conditioned diffusion prior for blind inverse problems. In ICASSP 2026-2026 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 11432–11436. IEEE, 2026. [40] Ajil Jalal, Marius Arvinte, Giannis Daras, Eric Price, Alexandros G Dimakis, and Jon Tamir. Robust compressed sensing mri with deep generative priors. Advances in neural information processing systems, 34:14938–14954, 2021.

50

PiX-MC

Wei et al.

[41] Sadeep Jayasumana, Srikumar Ramalingam, Andreas Veit, Daniel Glasner, Ayan Chakrabarti, and Sanjiv Kumar. Rethinking fid: Towards a better evaluation metric for image generation. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 9307–9315, 2024. [42] Yuanfeng Ji, Haotian Bai, Chongjian Ge, Jie Yang, Ye Zhu, Ruimao Zhang, Zhen Li, Lingyan Zhanng, Wanling Ma, Xiang Wan, et al. AMOS: A large-scale abdominal multi-organ benchmark for versatile medical image segmentation. Advances in neural information processing systems, 35:36722–36732, 2022. [43] Kostas Kampourakis and Kevin McCain. How uncertainty makes science advance. In Uncertainty: How It Makes Science Advance. Oxford University Press, 12 2019. [44] Tero Karras, Miika Aittala, Timo Aila, and Samuli Laine. Elucidating the design space of diffusion-based generative models. Advances in neural information processing systems, 35:26565– 26577, 2022. [45] Tero Karras, Samuli Laine, and Timo Aila. A style-based generator architecture for generative adversarial networks. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 4401–4410, 2019. [46] Bahjat Kawar, Michael Elad, Stefano Ermon, and Jiaming Song. Denoising diffusion restoration models. In Advances in Neural Information Processing Systems, 2022. [47] Amirhossein Kazerouni, Ehsan Khodapanah Aghdam, Moein Heidari, Reza Azad, Mohsen Fayyaz, Ilker Hacihaliloglu, and Dorit Merhof. Diffusion models in medical imaging: A comprehensive survey. Medical image analysis, 88:102846, 2023. [48] Hyojin Kim and Kyle Champley. Differentiable forward projector for x-ray computed tomography. In ICML 2023 Workshop on Differentiable Almost Everything: Differentiable Relaxations, Algorithms, Operators, and Simulators. [49] Tim Tsz-Kit Lau, Han Liu, and Thomas Pock. Non-log-concave and nonsmooth sampling via langevin monte carlo algorithms. In INdAM Workshop: Advanced Techniques in Optimization for Machine learning and Imaging, pages 83–149. Springer, 2022. [50] Rémi Laumont, Valentin De Bortoli, Andrés Almansa, Julie Delon, Alain Durmus, and Marcelo Pereyra. Bayesian imaging using plug & play priors: when langevin meets tweedie. SIAM Journal on Imaging Sciences, 15(2):701–737, 2022. [51] Holden Lee, Jianfeng Lu, and Yixin Tan. Convergence for score-based generative modeling with polynomial complexity. Advances in Neural Information Processing Systems, 35:22870–22882, 2022. [52] Anat Levin, Yair Weiss, Fredo Durand, and William T Freeman. Understanding and evaluating blind deconvolution algorithms. In 2009 IEEE conference on computer vision and pattern recognition, pages 1964–1971. IEEE, 2009. [53] Xiangru Lian, Yijun Huang, Yuncheng Li, and Ji Liu. Asynchronous parallel stochastic gradient for nonconvex optimization. Advances in neural information processing systems, 28, 2015.

51

PiX-MC

Wei et al.

[54] Xiangru Lian, Ce Zhang, Huan Zhang, Cho-Jui Hsieh, Wei Zhang, and Ji Liu. Can decentralized algorithms outperform centralized algorithms? a case study for decentralized parallel stochastic gradient descent. Advances in neural information processing systems, 30, 2017. [55] Jiaming Liu, Rushil Anirudh, Jayaraman J Thiagarajan, Stewart He, K Aditya Mohan, Ulugbek S Kamilov, and Hyojin Kim. DOLCE: A model-based probabilistic diffusion framework for limitedangle ct reconstruction. In Proceedings of the IEEE/CVF international conference on computer vision, pages 10498–10508, 2023. [56] Jaideep Mahajan, Kaihong Zhang, Feng Liang, and Jingbo Liu. Fast and efficient parallel sampling using higher order langevin dynamics, 2025. [57] David Martin, Charless Fowlkes, Doron Tal, and Jitendra Malik. A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics. In Proceedings eighth IEEE international conference on computer vision. ICCV 2001, volume 2, pages 416–423. Ieee, 2001. [58] Savvas Melidonis, Paul Dobson, Yoann Altmann, Marcelo Pereyra, and Konstantinos Zygalakis. Efficient bayesian computation for low-photon imaging problems. SIAM Journal on Imaging Sciences, 16(3):1195–1234, 2023. [59] Taylor R Moen, Baiyu Chen, David R Holmes III, Xinhui Duan, Zhicong Yu, Lifeng Yu, Shuai Leng, Joel G Fletcher, and Cynthia H McCollough. Low-dose CT image and projection dataset. Medical physics, 48(2):902–911, 2021. [60] Radford M. Neal. Annealed importance sampling. Statistics and Computing, 11(2):125–139, 2001. [61] Y. Nesterov. Introductory Lectures on Convex Optimization: A Basic Course. Kluwer Academic Publishers, 2004. [62] Jinshan Pan, Zhe Hu, Zhixun Su, and Ming-Hsuan Yang. l_0-regularized intensity and gradient prior for deblurring text images and beyond. IEEE transactions on pattern analysis and machine intelligence, 39(2):342–355, 2016. [63] Giorgio Parisi. Correlation functions and computer simulations. Nuclear Physics B, 180(3):378– 384, 1981. [64] Zhimin Peng, Yangyang Xu, Ming Yan, and Wotao Yin. Arock: an algorithmic framework for asynchronous parallel coordinate updates. SIAM Journal on Scientific Computing, 38(5):A2851– A2879, 2016. [65] Marcelo Pereyra, Luis A Vargas-Mieles, and Konstantinos C Zygalakis. The split gibbs sampler revisited: improvements to its algorithmic structure and augmented target distribution. SIAM Journal on Imaging Sciences, 16(4):2040–2071, 2023. [66] Maxim Raginsky, Alexander Rakhlin, and Matus Telgarsky. Non-convex learning via stochastic gradient langevin dynamics: a nonasymptotic analysis. In Conference on Learning Theory, pages 1674–1703. PMLR, 2017.

52

PiX-MC

Wei et al.

[67] Benjamin Recht, Christopher Re, Stephen Wright, and Feng Niu. Hogwild!: A lock-free approach to parallelizing stochastic gradient descent. Advances in neural information processing systems, 24, 2011. [68] Marien Renaud, Valentin De Bortoli, Arthur Leclaire, and Nicolas Papadakis. From stability of langevin diffusion to convergence of proximal mcmc for non-log-concave sampling. Advances in Neural Information Processing Systems, 38:115709–115773, 2026. [69] Marien Renaud, Julien Hermant, Deliang Wei, and Yu Sun. Provably accelerated imaging with restarted inertia and score-based image priors. In The Fourteenth International Conference on Learning Representations, 2026. [70] Peter Richtárik and Martin Takáč. Parallel coordinate descent methods for big data optimization. Mathematical Programming, 156(1):433–484, 2016. [71] Litu Rout, Negin Raoof, Giannis Daras, Constantine Caramanis, Alex Dimakis, and Sanjay Shakkottai. Solving linear inverse problems provably via posterior sampling with latent diffusion models. In Thirty-seventh Conference on Neural Information Processing Systems, 2023. [72] Leonid I Rudin, Stanley Osher, and Emad Fatemi. Nonlinear total variation based noise removal algorithms. Physica D: nonlinear phenomena, 60(1-4):259–268, 1992. [73] Adil Salim, Dmitry Kovalev, and Peter Richtárik. Stochastic proximal langevin algorithm: Potential splitting and nonasymptotic rates. Advances in Neural Information Processing Systems, 32, 2019. [74] Adil Salim and Peter Richtarik. Primal dual interpretation of the proximal stochastic gradient langevin algorithm. Advances in Neural Information Processing Systems, 33:3786–3796, 2020. [75] Ken Sauer and Charles Bouman. Bayesian estimation of transmission tomograms using segmentation based optimization. IEEE Transactions on Nuclear Science, 39(4):1144–1152, 1992. [76] Jonathan Richard Shewchuk et al. An introduction to the conjugate gradient method without the agonizing pain, 1994. [77] Andy Shih, Suneel Belkhale, Stefano Ermon, Dorsa Sadigh, and Nima Anari. Parallel sampling of diffusion models. Advances in Neural Information Processing Systems, 36:4263–4276, 2023. [78] Junhyuk So, Jiwoong Shin, Chaeyeon Jang, and Eunhyeok Park. Pcm: Picard consistency model for fast parallel sampling of diffusion models. In Proceedings of the Computer Vision and Pattern Recognition Conference, pages 23313–23322, 2025. [79] Bowen Song, Soo Min Kwon, Zecheng Zhang, Xinyu Hu, Qing Qu, and Liyue Shen. Solving inverse problems with latent diffusion models via hard data consistency. In International Conference on Learning Representations, volume 2024, pages 7624–7654, 2024. [80] Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. In Advances in Neural Information Processing Systems, volume 32, 2019. [81] Yang Song and Stefano Ermon. Improved techniques for training score-based generative models. In Advances in Neural Information Processing Systems, 2020. 53

PiX-MC

Wei et al.

[82] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:2011.13456, 2020. [83] Alessio Spagnoletti, Jean Prost, Andrés Almansa, Nicolas Papadakis, and Marcelo Pereyra. Latino-pro: Latent consistency inverse solver with prompt optimization. In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 19597–19607, 2025. [84] Tao Sun, Robert Hannah, and Wotao Yin. Asynchronous coordinate descent under more realistic assumptions. Advances in Neural Information Processing Systems, 30, 2017. [85] Yu Sun, Jiaming Liu, Yiran Sun, Brendt Wohlberg, and Ulugbek Kamilov. Async-RED: A provably convergent asynchronous block parallel stochastic method using deep denoising priors. In International Conference on Learning Representations, 2021. [86] Yu Sun, Zihui Wu, Yifan Chen, Berthy T Feng, and Katherine L Bouman. Provable probabilistic imaging using score-based generative priors. IEEE Transactions on Computational Imaging, 10:1290–1305, 2024. [87] Santosh Vempala and Andre Wibisono. Rapid convergence of the unadjusted langevin algorithm: Isoperimetry suffices. Advances in neural information processing systems, 32, 2019. [88] Pascal Vincent. A connection between score matching and denoising autoencoders. Neural Comput., 23(7):1661–1674, 2011. [89] Maxime Vono, Nicolas Dobigeon, and Pierre Chainais. Split-and-augmented gibbs sampler—application to large-scale inference problems. IEEE Transactions on Signal Processing, 67(6):1648–1661, 2019. [90] Deliang Wei, Shiyang Weng, and Fang Li. Nonconvex rician noise removal via convergent plug-and-play framework. Applied Mathematical Modelling, 123:197–212, 2023. [91] Andre Wibisono. Proximal langevin algorithm: Rapid convergence under isoperimetry. arXiv preprint arXiv:1911.01469, 2019. [92] Zihui Wu, Yu Sun, Yifan Chen, Bingliang Zhang, Yisong Yue, and Katherine L Bouman. Principled probabilistic imaging using diffusion models as plug-and-play priors. Advances in Neural Information Processing Systems, 37:118389–118427, 2024. [93] Xingyu Xu and Yuejie Chi. Provably robust score-based diffusion posterior sampling for plugand-play image reconstruction. Advances in Neural Information Processing Systems, 37:36148–36184, 2024. [94] Kaylee Yingxi Yang and Andre Wibisono. Convergence of the inexact langevin algorithm and score-based generative models in kl divergence. arXiv preprint arXiv:2211.01512, 2022. [95] Liutao Yang, Jiahao Huang, Guang Yang, and Daoqiang Zhang. Ct-sdm: A sampling diffusion model for sparse-view ct reconstruction across various sampling rates. IEEE Transactions on Medical Imaging, 44(6):2581–2593, 2025.

54

PiX-MC

Wei et al.

[96] Jure Zbontar, Florian Knoll, Anuroop Sriram, Tullie Murrell, Zhengnan Huang, Matthew J Muckley, Aaron Defazio, Ruben Stern, Patricia Johnson, Mary Bruno, et al. fastmri: An open dataset and benchmarks for accelerated mri. arXiv preprint arXiv:1811.08839, 2018. [97] Bingliang Zhang, Wenda Chu, Julius Berner, Chenlin Meng, Anima Anandkumar, and Yang Song. Improving diffusion inverse problem solving with decoupled noise annealing. In Proceedings of the Computer Vision and Pattern Recognition Conference, pages 20895–20905, 2025. [98] Huanjian Zhou and Masashi Sugiyama. Parallel simulation for sampling under isoperimetry and score-based diffusion models, 2025. [99] Zhengyuan Zhou, Panayotis Mertikopoulos, Nicholas Bambos, Peter Glynn, Yinyu Ye, Li-Jia Li, and Li Fei-Fei. Distributed asynchronous optimization with unbounded delays: How slow can you go? In International Conference on Machine Learning, pages 5970–5979. PMLR, 2018.

55

Related documents

• GSNR: Graph Smooth Null-Space Representation for Inverse Problems#986825
• Bayesian Analysis of Simple Random Densities#612660
• Bayesian Analysis of Simple Random Densities#611183
• FairBED: A Bayesian Experimental Design Approach to Gathering Fairer Data#299949
• Non-Gaussian Parameter Bias: Formalism and Validation#1157452