1
Tracing the Oracle: Improving Diffusion Timestep Scheduling for 3D CT Reconstruction Yujia Wu and Zhaoqiang Liu∗
arXiv:2606.06236v1 [cs.LG] 4 Jun 2026
Abstract Pretrained diffusion models demonstrate impressive potential in solving highly ill-posed 3D computed tomography (CT) inverse problems, while the inference process suffers from significant computational overhead. Furthermore, existing uniform timestep schedules fail to capture the non-uniform evolution of the reverse conditional diffusion stochastic differential equation, thereby introducing substantial truncation errors. To overcome this limitation, we propose Tracing the Oracle (TrO), a plug-and-play framework for improved timestep scheduling. Specifically, we treat densely sampled numerical integration trajectories on a few samples as the reference oracle. The optimized schedule is extracted by leveraging dynamic programming to globally minimize the cumulative error between the fewstep approximation and the oracle. This mechanism precisely allocates the limited sampling steps to critical evolution stages that are highly susceptible to truncation errors. Our extensive experiments on the AAPM dataset across multiple 3D CT reconstruction tasks demonstrate that, when combined with the state-of-the-art 3D CT reconstruction method DDS, our optimized timesteps significantly improve reconstruction fidelity and computational efficiency compared to existing heuristic schedules, especially under a strict budget of no more than 10 sampling steps.
I. I NTRODUCTION The reconstruction of 3D computed tomography (CT) is mathematically formulated as a linear inverse problem. This problem aims to reconstruct an unknown clean signal x∗0 ∈ Rd from a degraded observation y ∈ Rm . This degradation process is mathematically represented by the linear forward operator A ∈ Rm×d and additive Gaussian noise ξ ∈ Rm with ξ ∼ N (0, σy2 Im ), formalized as y = Ax∗0 + ξ.
(1)
Deep neural networks trained on paired datasets have been widely employed to directly learn mappings from degraded observations to clean signals [1], [2], [3]. However, acquiring large-scale paired datasets is often prohibitively expensive, and such data-driven models exhibit severe generalization vulnerabilities across different measurement operators [4]. Deep generative models including variational autoencoders [5] and generative adversarial networks [6] have achieved remarkable success in solving inverse problems, by exploiting data-driven priors learned from unpaired training datasets [7], [8], [9], [10], [11]. Recently, diffusion models [12], [13], [14] have emerged as powerful alternatives by modeling the data distribution through a stochastic differential equation (SDE). Leveraging ∗ Corresponding author.
The authors are with the School of Computer Science and Engineering, University of Electronic Science and Technology of China.
June 5, 2026
DRAFT
2
pre-trained unconditional diffusion models as generative priors enables highly effective training-free inference for inverse problems [15], [16], [17], [18], [19], [20], [21], [22], [23]. When addressing highly ill-posed 3D inverse problems, such as volumetric computed tomography (CT) reconstruction, directly training a 3D diffusion model imposes prohibitive computational costs [24]. Consequently, existing diffusion posterior sampling approaches for 3D inverse problems primarily utilize 2D diffusion priors augmented with inter-slice regularization constraints [24], [25]. For instance, Decomposed Diffusion Sampling (DDS) [25] improves reconstruction efficiency by using CGbased data-consistency updates initialized from the Tweedie denoised estimate. Under the local assumption that the tangent space around the denoised estimate can be represented or approximated by a Krylov subspace, these updates remain within the corresponding local subspace [26]. Nevertheless, existing approaches typically uniformly allocate the reverse sampling timesteps. Such heuristic timestep scheduling mechanisms fail to dynamically adapt to the complex interplay between the generative diffusion prior and the measurements of specific 3D inverse problems, thus yielding sub-optimal reconstruction quality under restricted sampling budgets. The sequential nature of the reverse diffusion process typically necessitates numerous neural function evaluations (NFEs), thereby imposing a substantial computational burden [27], [28], [29]. To maximize sampling efficiency, various sampler optimization strategies have been proposed [30], [31]. In addition, various timestep scheduling optimization strategies have been proposed [32], [33], [34] to accelerate the sampling process. For instance, the work [32] minimizes the explicit distance between numerical and analytic solutions, while GITS [33] leverages the geometric regularity of probability flow trajectories. However, solving inverse problems relies on measurementconditioned diffusion posterior sampling, where conditional measurements disrupt the unconditional generation manifold. Furthermore, achieving data consistency in inverse problems typically requires injecting stochastic noise during the reverse refinement process [16], [20], which invalidates the deterministic trajectory assumptions relied upon by unconditional scheduling algorithms, rendering them ineffective for complex conditional inverse solvers. To overcome these limitations, we propose Tracing the Oracle (TrO), which is a plug-and-play non-uniform timestep scheduling framework for 3D CT reconstruction (illustrated in Fig. 1). We initially construct a highfidelity reference trajectory, referred to as the oracle, using a densely sampled stochastic process on a subset of data. To accurately quantify the approximation error across these stochastic transitions, we introduce a noise reuse technique that isolates the discretization error between arbitrary timesteps. Subsequently, we formulate the scheduling task as a discrete optimization problem and deploy dynamic programming to minimize the cumulative trajectory error between the few-step approximation and the reference oracle. Given that diffusion trajectories for identically distributed data exhibit strong structural similarities [33], the optimized timestep schedule generalizes robustly across the entire dataset without incurring additional inference overhead. This optimization framework enables highly efficient and accurate 3D CT reconstructions. The primary contributions of this work are summarized as follows: •
We propose a plug-and-play framework that formulates timestep scheduling for measurement-conditioned diffusion as a discrete shortest-path problem, aiming to minimize cumulative trajectory errors in 3D inverse problems.
•
We introduce a noise reuse mechanism that isolates deterministic truncation errors from stochastic variance,
June 5, 2026
DRAFT
3
Noised Diffusion Direction
Conditional Direction
Final Direction
Gaussian Noise Ground Truth
Oracle Uniform-t TrO
Measurement
T
t
0
Fig. 1. Geometric interpretation of the proposed Tracing the Oracle (TrO) framework. The heuristic uniform timestep schedule diverges from the high-fidelity oracle trajectory during the conditional generation process. Our method improves the timestep allocation by minimizing the cumulative approximation error to closely track the oracle path. This optimized non-uniform scheduling strategy concentrates limited computational resources on critical evolution stages and effectively approximates the ground truth.
enabling precise quantification of discretization errors within SDE-based sampling trajectories. •
Extensive evaluations on 3D CT benchmarks demonstrate that our optimized schedules consistently improve reconstruction fidelity over predefined schedules without additional inference overhead. II. R ELATED W ORKS
A. Diffusion Models for 3D Inverse Problems Medical image reconstruction requires solving mathematically ill-posed linear inverse problems, under severe measurement constraints driven by patient safety protocols [35]. Diffusion priors have recently emerged as a prominent paradigm to regularize these underdetermined systems by enabling zero-shot inference over complex anatomical data manifolds [24]. While these zero-shot conditioning methods successfully steer initial unconditional noise vectors toward measurement consistency, they inherently rely on reversing SDEs, which requires thousands of sequential integrations [36]. When scaled to high-dimensional volumetric domains such as 3D CT, these numerical integrations incur severe memory consumption and substantial inference latency, making existing frameworks computationally demanding [37], [24]. To mitigate the memory limitations arising from native volumetric processing, recent literature explores patchbased and multi-view slice approximations. For instance, Song et al. [37] propose learning 3D patch image priors through position-aware diffusion score blending to ensure cross-slice consistency without relying on heuristic
June 5, 2026
DRAFT
4
regularizers. Further compressing the representational space, He et al. [38] introduce a triplane neural field integrated with diffusion priors to enable fast zero-shot reconstruction while managing volumetric complexities. To enforce measurement fidelity, Askari et al. [35] construct a bi-level guided diffusion model that approximates an inner-level conditional posterior mean before applying an outer-level proximal optimization. To circumvent the substantial computational burden of computing manifold-constrained gradients, Chung et al. [25] exploit the observation that, if the tangent space at a Tweedie denoised sample is locally represented by a Krylov subspace, then CG updates initialized from the denoised estimate remain in that local tangent subspace. By executing classical conjugate gradient (CG) operations strictly within this clean tangent space, the framework significantly reduces the required NFEs while preserving high-fidelity volumetric spatial features. However, DDS typically follows a fixed uniform timestep schedule, which can lead to suboptimal reconstruction fidelity when restricted to a limited computational budget (i.e., the few-step regime). Our proposed method can improve the timestep schedule to further enhance the sampling efficiency and reconstruction performance of the DDS framework. B. Sampling Acceleration for Unconditional Diffusion To improve inference speed, researchers have focused on refining numerical solvers for diffusion-based differential equations. DDIM [13] reformulates the generation pipeline from Markovian stochastic processes to deterministic numerical integration, allowing substantial reductions in temporal evaluations by skipping intermediate noise scales. DPM-Solver [27] leverages the inherent semi-linear nature of probability flow ODEs by employing exact exponential integrators to compute linear drifts analytically while approximating nonlinear network predictions via high-order Taylor expansions. UniPC [39] introduces a predictor-corrector framework that achieves high-order accuracy without requiring additional network evaluations per sampling step, significantly enhancing sample quality under limited sampling budgets. In addition to solver design, recent studies focus on optimizing discretized timestep schedules and sampler coefficients to enhance low-budget generation. Sabour et al. [40] optimize non-uniform schedules by minimizing a Kullback-Leibler upper bound between exact and approximate paths. Xue et al. [32] formulate a differentiable surrogate objective solved via trust-region methods to minimize the explicit distance between the numerical and analytical solutions. Zhou et al. [34] confirm that high-dimensional trajectories inherently reside within strict twodimensional subspaces and invoke the mean value theorem to distill exact intermediate scaling factors that eliminate polynomial truncation errors. The most relevant work to our approach is perhaps GITS [33], which focuses on the concept of trajectory geometry. GITS characterizes shape regularity within the implicit denoising trajectory and deploys dynamic programming to align the timestep schedule with the underlying structure of the probability flow ODE. While GITS demonstrates improved performance in unconditional generation, applying such deterministic trajectory geometry directly to conditional diffusion sampling corresponding to 3D inverse solvers poses significant challenges. We highlight two primary distinctions between our proposed framework and GITS designed to address these limitations. 1) GITS assumes a trajectory governed solely by the unconditional data manifold. In inverse problems, however, measurementinformed updates introduce external constraints that modify the trajectory geometry [31]. A schedule optimized for
June 5, 2026
DRAFT
5
unconditional priors may misalign with these conditional dynamics, especially when the measurement operator induces high curvature in the manifold transitions. Furthermore, our empirical results in Sec. V-C demonstrate that even distinct inverse problems exhibit varying preferences for timestep allocation, suggesting that the sampling dynamics are inherently coupled with the measurement process rather than being governed solely by the data prior. 2) GITS fundamentally relies on probability flow ODEs to compute deterministic state transitions. However, highfidelity inverse problem solvers typically require diffusion SDEs to prevent the generated samples from drifting off the data manifold [37], [16], [24]. Directly substituting SDE transitions into an ODE-based error metric, or computing naive geometric distances between stochastic SDE states, introduces severe random variance that biases the true trajectory error estimation. Our framework resolves this issue by proposing a noise reuse mechanism that accurately isolates the deterministic truncation error within the stochastic sampling procedure. III. P RELIMINARIES A. Diffusion Models Diffusion models establish a generative framework by progressively corrupting a data distribution into a tractable prior distribution through a continuous-time stochastic process. Let x0 ∈ Rd denote a clean data sample drawn from the underlying distribution pdata . The forward process over a continuous time variable t ∈ [0, T ] is governed by an SDE formulated as dxt = f (t)xt dt + g(t)dwt ,
(2)
where xt ∈ Rd represents the state vector at time t and wt ∈ Rd denotes the standard Wiener process. The function f (t) acts as the drift coefficient while g(t) functions as the diffusion coefficient. For each time t, the transition kernel from the initial data x0 to the latent state xt follows a Gaussian distribution defined by q(xt |x0 ) = N (xt ; αt x0 , σt2 Id ),
(3)
where αt and σt are differentiable functions representing the scale and noise schedule. Specifically, αt is strictly decreasing while σt is strictly increasing over t. To ensure that the SDE in Eq. (2) indeed induces the transition distribution q(xt |x0 ), the drift and diffusion coefficients are related to these functions according to [27]: f (t) =
d log αt , dt
g 2 (t) =
dσt2 d log αt − 2σt2 . dt dt
(4)
The boundary condition at the terminal time T ensures that the marginal distribution approximates a standard normal distribution characterized by xT ∼ N (0, Id ). The progression of the forward process is often measured by one half of the logarithm of signal-to-noise ratio (SNR), which is a strictly decreasing function of time denoted by λt = log(αt /σt ). To invert this forward transformation and recover the original data distribution, a neural network model ϵθ (xt , t) parameterized by θ is trained to predict the injected noise ϵ ∼ N (0, Id ). The network parameters are optimized by minimizing the following objective: Et∼U (0,T ),x0 ∼pdata ,ϵ∼N (0,Id ) ∥ϵθ (αt x0 + σt ϵ, t) − ϵ∥22 .
June 5, 2026
(5)
DRAFT
6
Song et al. [14] show that the reverse-time diffusion process evolving from T down to 0 can be formulated as an SDE: dxt = f (t)xt − g 2 (t)∇xt log pt (xt ) dt + g(t)dwt ,
(6)
where wt is a standard Wiener process in the reverse-time direction, pt denotes the marginal probability density of xt at time t, and the score function ∇xt log pt (xt ) can be approximated via the noise prediction network as ∇xt log pt (xt ) ≈ −ϵθ (xt , t)/σt . Furthermore, there exists a deterministic probability flow ODE that shares the same marginal densities as the SDE [14]: dxt =
1 f (t)xt − g 2 (t)∇xt log pt (xt ) dt. 2
(7)
In practice, Eq. (7) is integrated numerically over a user-defined fixed timestep schedule {ti }N i=0 , where ti > ti+1 for i = 0, 1, ..., N − 1 with t0 = T and tN = 0. At an arbitrary timestep ti , Tweedie’s formula [41] provides the expected posterior mean of the clean data, which is approximately derived as: x0|ti =
xti + σt2i ∇x log pti (xti ) αti
≈
xti − σti ϵθ (xti , ti ) αti
.
(8)
The state transition from the current step ti to the subsequent step ti+1 follows a deterministic numerical integration rule given by DDIM [13] as: xti+1 = αti+1 x0|ti + σti+1 ϵθ (xti , ti ).
(9)
B. Decomposed Diffusion Sampler (DDS) Training-free approaches solve the linear inverse problem in Eq. (1) by drawing samples from the conditional posterior distribution. However, achieving high reconstruction fidelity with such training-free solvers [36], [24], [42] necessitates multiple NFEs. DDS [25] mitigates these computational burdens by performing data-consistency updates on the Tweedie denoised representation. The key idea is that, under a local Krylov-subspace approximation of the tangent space around x0|ti , CG updates initialized from x0|ti provide an efficient local residual minimization strategy without explicitly computing manifold-constrained gradients. At timestep ti , to address highly ill-posed 3D inverse problems such as volumetric CT reconstruction, the algorithm seeks a refined estimate x̂ by solving a regularized optimization problem: 1 min ∥y − Ax∥22 + ζ∥Dz x∥1 , x 2
(10)
where Dz denotes the finite difference operator along the z-axis for Total Variation regularization [43], and ζ > 0 controls the regularization strength. To efficiently solve this while constraining the update direction strictly along the appropriate manifold tangent space, the algorithm utilizes the Alternating Direction Method of Multipliers (ADMM) [44]. By introducing an auxiliary variable z and a scaled dual variable w, the optimization is decomposed into tractable subproblems. Specifically, the data consistency update for x requires solving a linear system governed by ACG = A⊤ A+ρD⊤ z Dz and bCG = A⊤ y + ρD⊤ z (z − w), where ρ > 0 is the ADMM penalty parameter enforcing the constraint, while the
June 5, 2026
DRAFT
7
TV regularization strength is controlled by ζ. This system is efficiently resolved by applying M steps of the CG method [25] initialized at the unconditional Tweedie estimate x0|ti in Eq. (8): x̂0|ti = CG(ACG , bCG , x0|ti , M ).
(11)
Subsequently, the new auxiliary variable z′ is updated using a soft-thresholding operator Sζ/ρ (·) as follows: z′ = Sζ/ρ (Dz x̂0|ti + w),
(12)
where Sω (u) = sign(u) · max(|u| − ω, 0). In addition, the new scaled dual variable w′ is updated as follows: w′ = Dz x̂0|ti + w − z′ .
(13)
Implementation details of the ADMM and CG processes are provided in the supplementary material. Following this refinement, the unconditioned estimate in the discrete integration rule is substituted by the data-consistent estimate x̂0|ti . The algorithm transitions from the current timestep ti to the subsequent step ti+1 utilizing a stochastic generalized timestep schedule as DDIM [13], rather than the purely deterministic integration rule. The transition incorporates the measurement information defined by q xti+1 = αti+1 x̂0|ti + σt2i+1 − η 2 βt2i+1 ϵθ (xti , ti ) + ηβti+1 ϵi ,
(14)
where ϵi ∼ N (0, Id ) represents standard Gaussian noise. Here, η ∈ [0, 1] regulates the degree of stochasticity and βt represents the standard deviation of the noise at time t as defined in DDIM [13]. The deterministic trajectory is susceptible to accumulated approximation errors from the iterative data consistency projections. Introducing controlled stochastic noise mitigates these deviations and functions as a corrective mechanism to prevent samples from drifting off the clean data manifold. IV. M ETHODS To further enhance the sampling efficiency of diffusion-based 3D inverse solvers and enable high-fidelity reconstruction within the few-step regime, we propose Tracing the Oracle (TrO). Our framework formulates improved timestep scheduling as a discrete shortest-path optimization problem, identifying an optimized subset of evaluation points that minimizes the global truncation error for 3D inverse problems. A. Oracle Trajectory Continuous reverse diffusion processes trace trajectories along the data manifold from noise to real data. In the few-step regime, approximating this continuous evolution with sparse discrete transitions introduces truncation errors, causing the generated samples to deviate from the true data manifold and affecting reconstruction fidelity. Conversely, numerical integration with dense evaluations more closely tracks the continuous conditional distribution. Therefore, we utilize DDS [25] with a large number of evaluations N to generate a reference trajectory for a given measurement y. We define this densely sampled path as the oracle trajectory, denoted by Toracle = {xti }N i=0 . Our objective is to extract a target subset of L timesteps, {τj }L j=0 (L ≪ N ), such that the trajectory generated by this subset minimizes the cumulative approximation error relative to the oracle under our specified metric.
June 5, 2026
DRAFT
8
B. Discretization Error To extract the optimized subset, it is necessary to quantify the transition error incurred by skipping intermediate timesteps. Suppose we transition from a source timestep ti to a target timestep tj (ti > tj ), following the DDIM sampling rule in Eq. (14), the one-step approximation xtj |ti is derived as: q xtj |ti = αtj x̂0|ti + σt2j − η 2 βt2j ϵθ (xti , ti ) + ηβtj ϵ,
(15)
where x̂0|ti is the data-consistent estimate obtained through the CG optimization at time ti in Eq. (11), and ϵ ∼ N (0, Id ) is the injected standard Gaussian noise. The discretization error associated with this transition can be evaluated by measuring the geometric distance ∥xtj |ti − xtj ∥2 between the approximation xtj |ti and the corresponding intermediate state xtj from the oracle trajectory following GITS [33]. C. Noise Reuse However, the geometric distance between xtj |ti and xtj presents challenges in stochastic frameworks. Since both the dense oracle trajectory and the one-step approximation independently sample the stochastic noise term ϵ, the resulting difference xtj |ti − xtj is non-deterministic. Even if the deterministic components are fixed, independent target-step noise injections contribute Cov(ηβtj ϵ − ηβtj ϵ′ ) = 2η 2 βt2j Id , where ϵ and ϵ′ are independently sampled standard Gaussian noises. This covariance indicates that the geometric distance fluctuates, which obscures the deterministic truncation error. While setting the stochasticity parameter η to zero avoids this variance, executing purely deterministic sampling in conditional inverse solvers can trap the generation process in local optima and degrade reconstruction performance [20], [25], [16]. To address this, we employ a noise reuse mechanism. During the generation of the oracle trajectory Toracle , we explicitly record the stochastic noise vectors ϵi injected at each step ti . Specifically, in the error evaluation phase, to compute the approximation state xtj |ti targeting the oracle state xtj , we reuse the exact noise vector ϵj−1 that was originally utilized during the oracle transition from tj−1 to tj . This precise alignment neutralizes the random variance and ensures that the approximated state remains consistent with the marginal distribution of the actual SDE sampling trajectory, thereby isolating the deterministic truncation error within the stochastic sampling procedure. Consequently, we define the transition cost Ci,j from ti to tj as: 2 T Ci,j = ∥xtj |ti − xtj ∥2 + κ (ti − tj ) − . L
(16)
The first term represents the truncation error isolated via noise reuse. The second term is a regularization component scaled by a hyperparameter κ ≥ 0. This regularization penalizes large deviations from a uniform temporal stride (T /L) to prevent the neural network from operating outside its trained conditional manifold. The computational procedure for constructing the cost matrix is detailed in Algorithm 1.
June 5, 2026
DRAFT
9
D. Schedule Optimization Given the cost matrix C ∈ R(N +1)×(N +1) and {ti }N i=0 , the timestep scheduling task is formulated as a discrete shortest path problem to find the optimized timestep schedule {τj }L j=0 . Our objective is defined as: min
{τj }L j=0
s.t.
L−1 X
Ckj ,kj+1 ,
j=0
τj = tkj , kj ∈ {0, 1, . . . , N }, τj > τj+1 ,
j = 0, 1, ..., L
(17)
j = 0, 1, . . . , L − 1,
τ0 = t0 , τL = tN . To compute a schedule that minimizes this objective over the directed acyclic graph defined by C, we implement the dynamic programming algorithm (DP) [45] similar to GITS [33], presented in Algorithm 2. Specifically, the DP state records the minimum accumulated cost of reaching a timestep with a fixed number of transitions by enumerating all possible predecessor nodes. Therefore, given the pre-computed cost matrix C, the returned schedule is the global optimum of Eq. (17) rather than a greedy or locally optimized solution. Without dynamic programming, exact search −1 over all feasible schedules leads to a combinatorial time complexity of O L N . Given the precomputed cost L−1 matrix C, our DP solver instead obtains the global optimum in O(LN 2 ) time and O(LN ) memory. Data samples originating from the same underlying distribution typically exhibit consistent structural properties within their reverse diffusion trajectories [33], [30]. Thus, the timestep schedule extracted from a subset of calibration data generalizes across the test dataset. In our pipeline, the reference trajectory generation and the DP search are conducted in an offline preparatory phase. During the inference stage for 3D volume reconstruction, the pre-computed schedule {τj }L j=0 is applied as a target timestep schedule, ensuring that this scheduling approach introduces no additional computational overhead to the conditional generation process. V. E XPERIMENTAL R ESULTS A. Experimental Setups 1) 3D CT Reconstruction: We evaluate the proposed framework on sparse-view CT (SV-CT) and limited-angle CT (LA-CT) tasks. The forward operator A employs the discrete Radon transform. For SV-CT, we uniformly sample 2, 4, and 8 projection angles. For LA-CT, we restrict 120 projections to a 90◦ wedge. Experiments are conducted on the AAPM dataset [46] following standard protocols [25], [24]. We train the 2D diffusion prior on the training split, and designate one specific volume to compute the oracle trajectory (N = 200) and extract the optimized schedule. Subsequent evaluations are conducted on a distinct volume, which serves as the standard single-volume benchmark commonly adopted by prior works [25], [24]. Although community conventions restrict the evaluation to this single 3D volume, the reconstruction actually entails hundreds of independent 2D inferences across its high-resolution slices. Therefore, reporting metrics averaged over the axial, coronal, and sagittal planes provides a reliable quantitative assessment of our method.
June 5, 2026
DRAFT
10
Algorithm 1 Tracing the Oracle (TrO) N N N Require: Source timesteps {ti }N i=0 , {αti }i=0 , {σti }i=0 and {βti }i=0 , target number of steps L, regularization parameter κ, A, y, Dz , ρ, Sζ/ρ , M , η 1: Sample xt0 ∼ N (0, Id ) 2: z0 ← 0d , w0 ← 0d 3: for i ← 0 to N − 1 do 4:
ϵ̂i ← ϵθ (xti , ti )
5:
x0|ti ← (xti − σti ϵ̂i )/αti
6:
ACG ← A⊤ A + ρD⊤ z Dz
7:
bCG ← A⊤ y + ρD⊤ z (zi − wi )
8:
x̂0|ti ← CG(ACG , bCG , x0|ti , M )
9:
zi+1 ← Sζ/ρ (Dz x̂0|ti + wi )
10:
wi+1 ← wi + Dz x̂0|ti − zi+1
11:
Sample ϵi ∼ N (0, Id ) q xti+1 ← αti+1 x̂0|ti + σt2i+1 − η 2 βt2i+1 ϵ̂i + ηβti+1 ϵi
12:
▷ Tweedie denoising
▷ Data consistency
▷ DDIM sampling to ti+1
13: end for 14: C ← ∞(N +1)×(N +1) 15: for i ← 0 to N − 1 do 16:
for j ← i + 1 to N do
17:
xtj |ti ← αtj x̂0|ti +
q σt2j − η 2 βt2j ϵ̂i + ηβtj ϵj−1 ▷ DDIM sampling to tj with reusing noise ϵj−1
Ci,j ← ∥xtj |ti − xtj ∥2 + κ
18: 19:
2 (ti − tj ) − TL
end for
20: end for N 21: {τj }L j=0 ← DP(C, {ti }i=0 , N, L)
▷ Execute Dynamic Programming Algorithm to solve Eq. (17) 22: return {τj }L j=0
2) Implementation Details: We integrate our proposed timestep scheduling algorithm into DDS [25]. Since existing ODE-based optimized schedules are difficult to apply directly to stochastic conditional inverse solvers, we evaluate our approach against several widely adopted fixed schedules within the identical pipeline to ensure a fair comparison. These baselines include the uniform time schedule Uniform-t [13], the uniform half of log-SNR schedule Uniform-λ [27], the Quadratic schedule [13], the EDM schedule [47], and the Cosine schedule [48]. Detailed formulations for these schedules are provided in the supplementary material. To isolate the impact of the scheduling strategies, all methods utilize the identical pre-trained diffusion model, and the algorithmic hyperparameters are kept consistent with the original DDS implementation. Additionally, we compare our method against
June 5, 2026
DRAFT
11
Algorithm 2 Dynamic Programming Solver (DP) Require: Cost matrix C, Source time steps {ti }N i=1 , Source number of steps N , Target number of steps L 1: State ← ∞L×(N +1) , Path ← 0L×(N +1) 2: for j = 1 to N − 1 do 3:
State[0, j] ← C[0, j]
4:
Path[0, j] ← 0
5: end for 6: for i ← 1 to L − 1 do 7:
for j ← i + 1 to N do for k ← i − 1 to j − 1 do
8: 9:
if State[i − 1, k] + C[k, j] < State[i, j] then
10:
State[i, j] ← State[i − 1, k] + C[k, j]
11:
Path[i, j] ← k end if
12:
end for
13: 14:
end for
15: end for 16: Timesteps ← [tN ] 17: i ← N 18: for j ← L − 1 down to 0 do 19:
i ← Path[j, i]
20:
Append ti to Timesteps
21: end for 22: return reverse(Timesteps)
established diffusion solvers for 3D inverse problems, such as DiffMBIR [24], Score-Med [36], and MCG [42], which typically demand thousands of NFEs to achieve competitive reconstruction quality. Following approaches [24], [25], [42], [36], volumetric reconstruction quality is quantitatively assessed using the Peak Signal-to-Noise Ratio (PSNR) and the Structural Similarity Index Measure (SSIM), evaluated on Axial, Coronal, and Sagittal planes. B. Quantitative Results We comprehensively evaluate the volumetric reconstruction performance across SV-CT (e.g., 8-view and 4-view) and LA-CT tasks under strictly constrained computational budgets corresponding to 15, 10 and 8 NFEs. The quantitative comparisons, averaged across the axial, coronal, and sagittal planes, are summarized in Table I. The evaluation indicates that TrO yields competitive reconstruction fidelity across the tested configurations. Notably, the advantage of our optimized schedule over predefined counterparts becomes more evident as the NFE budget decreases. While this performance gap may tend to narrow as the number of NFEs increases, these results highlight
June 5, 2026
DRAFT
12
（
）
（
）
Measurement
Ground Truth
Uniform-t
Uniform-λ
Quadratic
EDM
Cosine
TrO Ours
Measurement
Ground Truth
Uniform-t
Uniform-λ
Quadratic
EDM
Cosine
TrO Ours
Fig. 2. Qualitative comparison of 3D CT reconstruction results using different timestep schedules. Red boxes highlight the detail recovery of our method.
Uniform-t Uniform-λ
1000
EDM Quadratic Cosine TrO (SV-CT) TrO (LA-CT)
Timestep Value (τi)
800
600
400
200
0 0
4
8
12
16
20
Sampling Step Index (i)
Fig. 3. Visualization of timestep allocations across different scheduling strategies. The graph illustrates the temporal evolution of continuous time τi over discrete evaluation steps i.
the utility of TrO in further enhancing the efficiency of the DDS framework with limited computational resources. Additionally, visual comparisons of the reconstructed tomographic slices are provided in Fig. 2, which further corroborate the enhanced artifact suppression and detail preservation capabilities of our adaptive schedule. C. Further Analysis Fig. 3 illustrates the actual temporal evolution of various scheduling strategies. The quantitative results in Table I indicate that the relative performance of static timestep schedules is highly dependent on the specific measurement operator and the available NFEs. For instance, while the Cosine schedule exhibits a noticeable performance drop in SV-CT tasks, it demonstrates relatively better stability in LA-CT settings compared to the EDM and Quadratic schedules.
June 5, 2026
DRAFT
13
TABLE I Q UANTITATIVE COMPARISON OF TIMESTEP SCHEDULES ON AAPM ACROSS DIFFERENT NFE S AND SV/LA-CT. B EST RESULTS ARE HIGHLIGHTED IN BOLD AND THE SECOND - BEST RESULTS ARE UNDERLINED .
8-view Axial
Coronal
4-view Sagittal
Axial
Coronal
Sagittal
Method (NFE)
Schedule
DDS (8)
Uniform-t Uniform-λ Quadratic EDM Cosine TrO (Ours)
37.08 36.68 37.58 35.59 36.60 38.53
0.940 0.937 0.941 0.928 0.933 0.952
37.89 38.03 38.59 36.95 37.64 39.78
0.935 0.931 0.931 0.922 0.926 0.964
36.74 36.94 37.87 36.09 36.70 38.62
0.936 0.932 0.934 0.923 0.929 0.948
32.09 31.99 32.66 30.42 31.77 33.24
0.894 0.894 0.901 0.830 0.890 0.917
33.05 33.14 33.60 30.27 33.02 34.66
0.887 0.889 0.890 0.815 0.885 0.910
31.13 31.16 31.57 29.79 31.24 32.78
0.889 0.891 0.891 0.821 0.887 0.914
DDS (10)
Uniform-t Uniform-λ Quadratic EDM Cosine TrO (Ours)
37.66 38.08 38.72 38.38 36.29 39.33
0.945 0.948 0.951 0.946 0.907 0.955
38.52 39.03 39.67 39.37 36.42 40.36
0.939 0.943 0.946 0.940 0.895 0.950
37.34 37.96 38.66 38.61 35.96 39.46
0.941 0.944 0.948 0.942 0.900 0.951
32.78 32.34 33.29 33.66 30.88 34.12
0.900 0.899 0.908 0.921 0.800 0.927
33.56 33.55 34.25 34.63 30.96 35.37
0.895 0.894 0.899 0.916 0.780 0.922
31.53 31.51 32.15 33.21 30.12 33.86
0.896 0.895 0.901 0.917 0.791 0.920
DDS (15)
Uniform-t Uniform-λ Quadratic EDM Cosine TrO (Ours)
38.54 39.01 39.63 39.82 36.63 40.02
0.950 0.953 0.955 0.956 0.905 0.957
39.47 40.04 40.61 40.90 36.69 40.96
0.945 0.948 0.951 0.951 0.894 0.952
38.36 39.01 39.69 40.06 36.31 40.14
0.947 0.949 0.952 0.953 0.900 0.954
33.82 33.54 34.84 34.84 30.68 34.99
0.911 0.907 0.923 0.929 0.763 0.932
34.67 34.48 35.68 36.50 30.58 36.48
0.906 0.904 0.917 0.924 0.744 0.923
32.68 32.36 33.66 34.49 30.03 34.61
0.907 0.905 0.918 0.925 0.758 0.927
DiffMBIR(4000) Score-Med(4000) MCG(4000)
/ / /
33.49 29.10 28.61
0.942 0.882 0.873
35.18 27.93 28.05
0.967 0.875 0.884
32.18 24.23 24.45
0.910 0.759 0.765
30.52 28.20 27.33
0.914 0.867 0.855
30.09 27.48 26.52
0.938 0.889 0.863
27.89 25.08 23.04
0.871 0.783 0.745
PSNR↑ SSIM↑ PSNR↑ SSIM↑ PSNR↑ SSIM↑ PSNR↑ SSIM↑ PSNR↑ SSIM↑ PSNR↑ SSIM↑
90◦ -Angle
2-view Axial
Coronal
Sagittal
Axial
Coronal
Sagittal
Method(NFE)
Schedule
DDS (8)
Uniform-t Uniform-λ Quadratic EDM Cosine TrO (Ours)
25.04 24.55 25.18 21.50 23.39 26.23
0.742 0.730 0.752 0.606 0.705 0.798
25.98 25.76 26.02 22.03 25.52 26.74
0.723 0.714 0.753 0.634 0.752 0.789
24.12 23.97 24.16 21.30 23.11 24.74
0.710 0.707 0.740 0.615 0.704 0.803
38.07 38.29 38.28 36.71 38.36 38.82
0.944 0.949 0.944 0.941 0.947 0.952
39.33 39.84 39.73 38.26 39.62 40.32
0.936 0.942 0.942 0.934 0.940 0.945
38.27 38.52 38.60 36.64 38.57 39.21
0.939 0.945 0.946 0.937 0.946 0.946
DDS (10)
Uniform-t Uniform-λ Quadratic EDM Cosine TrO (Ours)
25.52 25.25 26.05 23.27 24.06 26.70
0.731 0.741 0.781 0.728 0.558 0.811
26.30 25.85 27.03 23.68 24.34 27.60
0.712 0.724 0.798 0.762 0.525 0.819
24.55 24.26 24.73 22.50 23.45 25.36
0.720 0.729 0.796 0.739 0.546 0.818
38.52 38.61 38.54 37.76 38.68 39.17
0.948 0.950 0.947 0.946 0.950 0.954
39.85 40.14 40.17 39.54 40.07 40.51
0.941 0.944 0.945 0.940 0.943 0.947
38.77 38.91 38.87 38.06 38.93 39.47
0.943 0.946 0.945 0.943 0.946 0.949
DDS (15)
Uniform-t Uniform-λ Quadratic EDM Cosine TrO (Ours)
25.93 25.52 26.93 26.37 24.04 27.29
0.681 0.702 0.813 0.802 0.491 0.820
26.32 26.01 28.01 26.41 24.03 28.24
0.667 0.687 0.821 0.809 0.463 0.833
24.77 24.48 25.27 24.61 23.34 25.74
0.680 0.697 0.827 0.803 0.489 0.832
39.00 38.80 38.88 37.93 38.97 39.13
0.951 0.950 0.952 0.944 0.950 0.952
40.35 40.30 40.33 39.69 40.32 40.52
0.945 0.944 0.944 0.937 0.944 0.946
39.28 39.14 39.22 38.33 39.25 39.46
0.947 0.946 0.947 0.940 0.946 0.948
DiffMBIR(4000) Score-Med(4000) MCG(4000)
/ / /
24.11 24.07 24.69
0.810 0.808 0.821
23.15 23.70 23.52
0.841 0.822 0.806
21.72 20.95 20.71
0.766 0.720 0.685
34.92 27.80 26.01
0.956 0.852 0.838
32.48 25.69 24.55
0.947 0.869 0.823
28.82 22.03 21.59
0.832 0.735 0.706
PSNR↑ SSIM↑ PSNR↑ SSIM↑ PSNR↑ SSIM↑ PSNR↑ SSIM↑ PSNR↑ SSIM↑ PSNR↑ SSIM↑
June 5, 2026
DRAFT
14
TABLE II A BLATION STUDY OF κ FOR THE 8- VIEW SV-CT TASK ON AAPM. B EST RESULTS ARE HIGHLIGHTED IN BOLD AND THE SECOND - BEST RESULTS ARE UNDERLINED .
NFE 8
10
15
Metric
κ=0
κ=1
κ=5
κ = 10
Uniform-t
PSNR↑
38.42
38.98
38.65
37.64
37.24
SSIM↑
0.951
0.955
0.953
0.946
0.942
PSNR↑
39.32
39.64
39.72
38.96
37.84
SSIM↑
0.943
0.949
0.952
0.948
0.942
PSNR↑
39.63
39.77
40.37
39.83
38.80
SSIM↑
0.951
0.951
0.954
0.952
0.947
1) Ablation Study: We conduct an ablation study on the AAPM dataset to investigate the impact of the regularization hyperparameter κ in Eq. (16) for 8-view SV-CT reconstruction task. The reported PSNR and SSIM are averaged across the axial, coronal, and sagittal planes. As shown in Table II, the optimal κ varies depending on the NFEs. Generally, a larger NFE requires a correspondingly larger κ to enforce a similar level of constraint, which stems from the differences in N/L. Furthermore, the empirical results indicate that the introduction of the regularization term effectively prevents the generated schedule from degenerating into extreme timestep allocations, which enhances the overall reconstruction performance. 2) Timestep Analysis: We hypothesize that this divergence originates from the frequency reconstruction dynamics of diffusion models. The early stages of the reverse diffusion process primarily synthesize low-frequency global structures, while the later stages refine high-frequency details. In SV-CT, the uniformly distributed projections provide reliable global low-frequency constraints. Consequently, computational resources are better allocated to the later stages (e.g., as in EDM and Quadratic) for high-frequency artifact suppression. Conversely, LA-CT concentrates projections within a narrow region, introducing severe structural artifacts and unreliable low-frequency information in the unmeasured wedges. Allocating insufficient steps to the early phase implicitly over-trusts these distorted lowfrequency measurements, severely impeding late-stage optimization. The Cosine schedule mitigates this by densely populating the early stages, providing the generative prior with sufficient capacity to synthesize the missing reliable low-frequency structures. These observations reveal that a single fixed timestep allocation is sub-optimal across diverse inverse problems. Notably, our proposed TrO framework effectively captures this physical discrepancy. As clearly depicted in Fig. 3, TrO adaptively shifts the step allocation toward the earlier stages for the LA-CT task compared to the SV-CT task. By continuously evaluating the actual SDE truncation error, TrO dynamically extracts task-specific and budget-aware temporal configurations, thereby ensuring consistent reconstruction fidelity across varying degradation scenarios. VI. C ONCLUSION In this paper, we present Tracing the Oracle (TrO), a plug-and-play framework that formulates timestep scheduling for diffusion-based 3D inverse problems as a discrete shortest-path optimization problem. To quantify the trajectory discretization error, we propose a noise reuse mechanism that isolates the deterministic truncation error within the
June 5, 2026
DRAFT
15
SDE sampling procedure. By leveraging dynamic programming, TrO extracts an optimized timestep schedule from a densely sampled oracle trajectory. Evaluations on sparse-view and limited-angle CT reconstruction demonstrate that our proposed approach improves reconstruction fidelity, particularly under limited computational budgets. By reducing the performance gap between few-step sampling and high-NFE solvers without additional inference overhead, TrO facilitates the efficient and robust deployment of diffusion models in 3D medical imaging applications. A PPENDIX A T IMESTEP S CHEDULE To evaluate the performance of our proposed Tracing the Oracle (TrO) framework, we compare it against several widely adopted fixed timestep schedules. All schedules partition the continuous time interval [ϵ, T ] into N discrete steps {tn }N n=0 , where t0 = T and tN = ϵ ≈ 0. The formulations for these schedules are summarized as follows: A. Uniform-t Schedule The Uniform-t schedule splits the time interval [ϵ, T ] into N sub-intervals of equal length [13]. The n-th timestep is defined as: n (ϵ − T ), N
tn = T +
n = 0, 1, . . . , N.
(18)
B. Uniform-λ Schedule The Uniform-λ schedule is based on the half of log-signal-to-noise ratio (SNR), defined as λt = log(αt /σt ). This scheme uniformly discretizes the interval [λT , λϵ ] and maps the values back to the time domain [27]: n tn = tλ λT + (λϵ − λT ) , N
(19)
where tλ (·) denotes the inverse function of λt . C. EDM Schedule The EDM schedule employs a non-linear discretization strategy based on the variable substitution κt = σt /αt 1/ρ
[47]. It uniformly discretizes κt
for a positive integer ρ: ρ n 1/ρ 1/ρ tn = tκ κT + (κ1/ρ − κT ) , ϵ N
(20)
where tκ (·) is the inverse function of κt . Following the standard implementation, we set ρ = 7 in our experiments. D. Quadratic Schedule The Quadratic schedule allocates more steps toward the smaller noise levels (near t = 0) to refine high-frequency details [13]. The timesteps are determined by: tn =
June 5, 2026
√
T+
√ 2 n √ ( ϵ − T) . N
(21)
DRAFT
16
E. Cosine Schedule The Cosine schedule is designed to ensure a more gradual loss of information during the forward diffusion process [48]. It is typically defined through the cumulative noise schedule ᾱt : 2 f (t) t/T + s π , ᾱt = , f (t) = cos · f (0) 1+s 2
(22)
where s is a small offset (e.g., s = 0.008) to prevent singularities at t = 0. A PPENDIX B D ETAILED DDS A. ADMM To address the highly ill-posed 3D CT inverse problem, we incorporate total variation (TV) regularization along the z-axis to leverage cross-slice structural priors: 1 min ∥y − Ax∥22 + ζ∥Dz x∥1 , x 2
(23)
where A is the Radon transform operator, y denotes the projection measurements, Dz represents the finite difference operator along the axial direction, and ζ is the regularization parameter. 1) Constraint Reformulation: By introducing an auxiliary variable z ∈ Rd , we convert the unconstrained L1 problem into an equivalent constrained optimization framework: 1 min ∥y − Ax∥22 + ζ∥z∥1 , x,z 2
s.t. z = Dz x.
(24)
The corresponding Augmented Lagrangian function is defined as: Lρ (x, z, w) =
1 ρ ∥y − Ax∥22 + ζ∥z∥1 + ∥Dz x − z + w∥22 , 2 2
(25)
where w is the scaled dual variable and ρ > 0 denotes the ADMM penalty parameter. 2) Iterative Updates: The ADMM algorithm proceeds by alternating the minimization of Lρ with respect to each variable: a) Primal Variable x-update:: Fixing zj and wj , the update for x is obtained by solving: 1 ρ xj+1 = arg min ∥y − Ax∥22 + ∥Dz x − zj + wj ∥22 . x 2 2
(26)
Taking the gradient with respect to x and setting it to zero yields the following linear system: ⊤ ⊤ (A⊤ A + ρD⊤ z Dz )xj+1 = A y + ρDz (zj − wj ).
(27)
In our framework, this high-dimensional system is solved iteratively using the Conjugate Gradient (CG) method. b) Auxiliary Variable z-update:: Fixing xj+1 and wj , the update for z is given by: ρ zj+1 = arg min ζ∥z∥1 + ∥Dz xj+1 − z + wj ∥22 . z 2
(28)
This proximal subproblem has a closed-form solution via the soft-thresholding operator Sτ (·): zj+1 = Sζ/ρ (Dz xj+1 + wj ),
(29)
where Sτ (u) = sign(u) · max(|u| − τ, 0).
June 5, 2026
DRAFT
17
c) Dual Variable w-update:: The dual variable is updated to enforce the linear constraint: wj+1 = wj + Dz xj+1 − zj+1 .
(30)
To ensure computational efficiency during the reverse diffusion process, we perform only a single ADMM iteration per sampling step and share the variables {z, w} across timesteps, a technique known as fast variable sharing. B. CG The detailed implementation for the conjugate gradient (CG) is given in Algorithm 3. Algorithm 3 Conjugate Gradient (CG) Require: Symmetric positive definite matrix A, right-hand side b, initial estimate x0 , number of iterations M 1: r 0 ← b − Ax0 2: p0 ← r 0 3: for k = 0 to M − 1 do r⊤ r
4:
αk ← p⊤kApk
5:
xk+1 ← xk + αk pk
6:
r k+1 ← r k − αk Apk
7:
βk ←
8:
pk+1 ← r k+1 + βk pk
k
k
r⊤ k+1 r k+1 r⊤ k rk
9: end for 10: return xM
C. Full DDS in 3D CT Reconstruction After obtaining the optimized timesteps {τi }L i=0 , we utilize this schedule with DDS [25] to achieve 3D reconstruction. This step is analogous to computing the oracle trajectory in the TrO algorithm, with the complete procedure summarized in Algorithm 4.
June 5, 2026
DRAFT
18
Algorithm 4 DDS L L L Require: Target timesteps {τi }L i=0 , {ατi }i=0 , {στi }i=0 and {βτi }i=0 , A, y, Dz , ρ, Sζ/ρ , M , η 1: Sample xτ0 ∼ N (0, Id ) 2: z0 ← 0d , w0 ← 0d 3: for i ← 0 to L − 1 do 4:
ϵ̂i ← ϵθ (xτi , τi )
5:
x0|τi ← (xτi − στi ϵ̂i )/ατi ⊤
▷ Tweedie denoising
A + ρD⊤ z Dz
6:
ACG ← A
7:
bCG ← A⊤ y + ρD⊤ z (zi − wi )
8:
x̂0|τi ← CG(ACG , bCG , x0|τi , M )
9:
zi+1 ← Sζ/ρ (Dz x̂0|τi + wi )
10:
wi+1 ← wi + Dz x̂0|τi − zi+1
11:
Sample ϵi ∼ N (0, Id )
12:
xτi+1 ← ατi+1 x̂0|τi +
q
▷ Data consistency
στ2i+1 − η 2 βτ2i+1 ϵ̂i + ηβτi+1 ϵi
▷ DDIM sampling to τi+1
13: end for 14: return xτL
A PPENDIX C D ETAILS OF DYNAMIC P ROGRAMMING We provide additional details of the dynamic programming solver used in Algorithm 2. Let State[s, j] denote the minimum accumulated cost of reaching timestep tj from t0 using exactly s + 1 transitions. The corresponding predecessor index is stored in Path[s, j]. The initialization is given by State[0, j] = C0,j ,
Path[0, j] = 0,
j = 1, . . . , N.
(31)
For s = 1, . . . , L − 1, the recurrence is State[s, j] = min {State[s − 1, k] + Ck,j } , k<j
(32)
and Path[s, j] records the index k that attains the minimum. After computing State[L − 1, N ], the optimized timestep indices are obtained by backtracking from N through Path. This procedure enumerates all valid predecessor nodes at each transition count, and therefore returns the global optimum of Eq. (17) under the pre-computed cost matrix C. The time complexity is O(LN 2 ) and the memory complexity is O(LN ).
June 5, 2026
DRAFT
19
R EFERENCES [1] C. Dong, C. C. Loy, K. He, and X. Tang, “Image super-resolution using deep convolutional networks,” IEEE TPAMI, 2015. [2] Z. Tu, H. Talebi, H. Zhang, F. Yang, P. Milanfar, A. Bovik, and Y. Li, “MAXIM: Multi-axis MLP for image processing,” in CVPR, 2022. [3] S. W. Zamir, A. Arora, S. Khan, M. Hayat, F. S. Khan, M.-H. Yang, and L. Shao, “Multi-stage progressive image restoration,” in CVPR, 2021. [4] M. Delbracio and P. Milanfar, “Inversion by direct iteration: An alternative to denoising diffusion for image restoration,” TMLR, 2023. [5] D. P. Kingma, “Auto-encoding variational bayes,” in ICLR, 2014. [6] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial nets,” in NeurIPS, 2014. [7] A. Bora, A. Jalal, E. Price, and A. G. Dimakis, “Compressed sensing using generative models,” in ICML, 2017. [8] Z. Liu, X. Wang, and J. Liu, “Misspecified phase retrieval with generative priors,” in NeurIPS, 2022. [9] J. Chen, J. Scarlett, M. Ng, and Z. Liu, “A unified framework for uniform signal recovery in nonlinear generative compressed sensing,” in NeurIPS, 2023. [10] G. Daras, J. Dean, A. Jalal, and A. Dimakis, “Intermediate layer optimization for inverse problems using deep generative models,” in ICML, 2021. [11] J. Chen, M. K. Ng, and Z. Liu, “Solving quadratic systems with full-rank matrices using sparse or generative priors,” IEEE TSP, 2025. [12] J. Ho, A. Jain, and P. Abbeel, “Denoising diffusion probabilistic models,” in NeurIPS, 2020. [13] J. Song, C. Meng, and S. Ermon, “Denoising diffusion implicit models,” in ICLR, 2021. [14] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole, “Score-based generative modeling through stochastic differential equations,” in ICLR, 2021. [15] H. Chung, J. Kim, M. T. McCann, M. L. Klasky, and J. C. Ye, “Diffusion posterior sampling for general noisy inverse problems,” in ICLR, 2023. [16] Y. Zhu, K. Zhang, J. Liang, J. Cao, B. Wen, R. Timofte, and L. Van Gool, “Denoising diffusion models for plug-and-play image restoration,” in CVPRW, 2023. [17] X. Peng, Z. Zheng, W. Dai, N. Xiao, C. Li, J. Zou, and H. Xiong, “Improving diffusion models for inverse problems using optimal posterior covariance,” in ICML, 2024. [18] J. Zhang, J. Zhuang, C. Jin, G. Li, and Y. Gu, “Unleashing the denoising capability of diffusion prior for solving inverse problems,” in NeurIPS, 2024. [19] J. Chang, C. Duan, Y. Jiao, R. Li, J. Z. Yang, and C. Yuan, “Provable diffusion posterior sampling for Bayesian inversion,” arXiv preprint arXiv:2512.08022, 2025. [20] B. Zhang, W. Chu, J. Berner, C. Meng, A. Anandkumar, and Y. Song, “Improving diffusion inverse problem solving with decoupled noise annealing,” in CVPR, 2025. [21] Y. Zheng, W. Li, and Z. Liu, “Integrating intermediate layer optimization and projected gradient descent for solving inverse problems with diffusion models,” in ICML, 2025. [22] J. Li and C. Wang, “Image restoration via primal dual hybrid gradient and flow generative model,” in AAAI, 2026. [23] Y. Zheng, W. Li, and Z. Liu, “Image restoration via diffusion models with dynamic resolution,” in ICML, 2026. [24] H. Chung, D. Ryu, M. T. Mccann, M. L. Klasky, and J. C. Ye, “Solving 3D inverse problems using pre-trained 2D diffusion models,” in CVPR, 2023. [25] H. Chung, S. Lee, and J. C. Ye, “Decomposed diffusion sampler for accelerating large-scale inverse problems,” in ICLR, 2024. [26] G. Meurant and J. D. Tebbens, Krylov methods for nonsymmetric linear systems: From theory to computations.
Springer, 2020.
[27] C. Lu, F. Bao, J. Chen, C. Li, and J. Zhu, “DPM-Solver: A fast ODE solver for diffusion probabilistic model sampling in around 10 steps,” in NeurIPS, 2022. [28] F. Fu, T. Guo, and Z. Liu, “Learnable sampler distillation for discrete diffusion models,” in NeurIPS, 2025. [29] W. Sun, J. Li, and Z. Liu, “Just-in-time: Training-free spatial acceleration for diffusion Transformers,” in CVPR, 2026. [30] E. Frankel, S. Chen, J. Li, P. W. Koh, L. J. Ratliff, and S. Oh, “S4S: Solving for a diffusion model solver,” in ICML, 2025. [31] J. Zhang, Z. Liu, L. Yan, G. Li, and Y. Gu, “Improving diffusion-based inverse algorithms under few-step constraint via linear extrapolation,” in NeurIPS, 2025. [32] S. Xue, Z. Liu, F. Chen, S. Zhang, T. Hu, E. Xie, and Z. Li, “Accelerating diffusion sampling with optimized time steps,” in CVPR, 2024.
June 5, 2026
DRAFT
20
[33] D. Chen, Z. Zhou, C. Wang, C. Shen, and S. Lyu, “On the trajectory regularity of ODE-based diffusion sampling,” in ICML, 2024. [34] Z. Zhou, D. Chen, C. Wang, and C. Chen, “Fast ODE-based sampling for diffusion models in around 5 steps,” in CVPR, 2024. [35] H. Askari, F. Roosta, and H. Sun, “Training-free medical image inverses via bi-level guided diffusion models,” in WACV, 2025. [36] Y. Song, L. Shen, L. Xing, and S. Ermon, “Solving inverse problems in medical imaging with score-based generative models,” in ICLR, 2022. [37] B. Song, J. Hu, Z. Luo, J. A. Fessler, and L. Shen, “DiffusionBlend: learning 3D image prior through position-aware diffusion score blending for 3D computed tomography reconstruction,” in NeurIPS, 2024. [38] J. He, B. Li, G. Yang, and Z. Liu, “Blaze3DM: Integrating triplane representation with diffusion for solving 3D inverse problems in medical imaging,” in MICCAI, 2025. [39] W. Zhao, L. Bai, Y. Rao, J. Zhou, and J. Lu, “UniPC: A unified predictor-corrector framework for fast sampling of diffusion models,” in NeurIPS, 2023. [40] A. Sabour, S. Fidler, and K. Kreis, “Align your steps: Optimizing sampling schedules in diffusion models,” in ICML, 2024. [41] B. Efron, “Tweedie’s formula and selection bias,” JASA, 2011. [42] H. Chung, B. Sim, D. Ryu, and J. C. Ye, “Improving diffusion models for inverse problems using manifold constraints,” in NeurIPS, 2022. [43] E. Y. Sidky and X. Pan, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization,” PMB, 2008. [44] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” ML, 2010. [45] R. Bellman, “Dynamic programming,” Science, 1966. [46] American Association of Physicists in Medicine (AAPM) and Mayo Clinic, “2016 NIH-AAPM-Mayo clinic low dose CT grand challenge dataset,” https://www.aapm.org/GrandChallenge/LowDoseCT/, 2016. [47] T. Karras, M. Aittala, T. Aila, and S. Laine, “Elucidating the design space of diffusion-based generative models,” in NeurIPS, 2022. [48] A. Q. Nichol and P. Dhariwal, “Improved denoising diffusion probabilistic models,” in ICML, 2021.
June 5, 2026
DRAFT