Source-linked AI summary
Reducing Complexity of Pseudo-Marginal Sequential Monte Carlo via Hierarchical Tempering
Josef Martínek, Julian Hofstadler, Alexander M. G. Cox, Robert Scheichl
TL;DR
In inverse problems with intractable likelihoods, stochastic forward solvers can yield biased likelihood estimates, while repeated evaluations make posterior approximation costly. The paper develops hierarchical pseudo-marginal SMC using increasing solver accuracy across tempering and mutation stages, reporting a 14.3× speedup in photoacoustic imaging and a cost ratio that grows linearly with T under assumptions.
Problem
Inverse problems can have intractable likelihoods without an available unbiased likelihood estimator, making posterior sampling computationally challenging.
Source: §1 · p. 1Method
Hierarchical SMC uses stochastic solvers of increasing accuracy along the tempering sequence and combines cheap approximate with expensive target-invariant kernels within mutation stages.
Source: §1 · p. 3, §5.1 · p. 14, §7 · p. 25Results
14.3× speedup reduced photoacoustic imaging runtime from approximately 4.6 days to 7.7 hours, while posterior approximations remained close to the standard-SMC reference; theoretically, the cost ratio grows linearly in T under assumptions.
Source: §1 · p. 3, §7 · p. 26Takeaways & Limitations
In the nonlinear photoacoustic imaging problem, adaptive hierarchical SMC lowered runtime while retaining posterior approximations close to the standard-SMC reference.
Source: §7 · p. 26A joint complexity analysis of the algorithm using both acceleration mechanisms simultaneously remains open.
Source: §7 · p. 26
Abstract
from arXiv · showhide
The Bayesian approach to inverse problems provides a powerful framework for inferring unknown parameters from noisy observations. Among the computationally most challenging statistical inference problems are problems with intractable likelihoods, especially if no unbiased likelihood estimator is available. We focus on a class of such problems, where a biased likelihood estimator is obtained through a stochastic approximation of a forward map (stochastic solver). Obtaining an empirical posterior approximation in applications is very expensive due to the repeated stochastic solver evaluations required by standard sampling methods. To address this computational bottleneck, we develop a novel algorithm based on a density-tempered marginalized Sequential Monte Carlo (SMC) sampler. Our approach (hierarchical SMC) accelerates inference by leveraging a hierarchy of stochastic solvers with increasing accuracy both along the tempering sequence and within mutation stages of the SMC sampler. We establish, under additional assumptions, that the asymptotic speedup is proportional to the number of tempering stages $T$. We demonstrate computational gains on several examples, including a nonlinear, high-contrast inverse problem arising in photoacoustic medical imaging. On this problem, we demonstrate a speedup by a factor of 14.3, reducing the inference wall-clock time from about 4.6 days to 7.7 hours on an NVIDIA A40 GPU.
1 Introduction
The paper develops hierarchical pseudo-marginal SMC for inverse problems with intractable likelihoods, using stochastic solvers of increasing accuracy. It analyzes the resulting complexity and mutation-stage acceleration, and demonstrates a substantial runtime reduction on a photoacoustic imaging problem.
1 Introduction
- When no unbiased likelihood estimator is available, stochastic approximations of the forward map can still define an approximate posterior, but solver accuracy trades computational cost against posterior fidelity.Higher solver accuracy is more expensive, while smaller accuracy may insufficiently target the true posterior.Source: §1 · pp. 1–2
- The method uses increasingly accurate stochastic solvers across tempering stages, relying on simpler early targets to mitigate inexact early likelihoods.High precision is used later to support sufficiently accurate Markov transitions.Source: §1 · p. 3
- The paper provides a general density-tempered pseudo-marginal SMC framework and finite-sample error bounds matching those of standard SMC for a suitable accuracy sequence.The framework targets problems where stochastic forward solvers yield likelihood approximations.Source: §1 · p. 3
- The cost ratio of standard to hierarchical SMC grows linearly in T while both achieve the same error bound.This result holds under the paper’s assumptions for the selected accuracy schedule.Source: §1 · p. 3
- A hierarchy of Markov kernels can improve SMC mutation stages, with mixing bounds used to derive an optimal budget allocation.The analysis concerns hierarchical mutation stages.Source: §1 · p. 3
- For a Bayesian inverse problem in photoacoustic medical imaging, the algorithm reduced computational cost by roughly 14.3×, from 4.6 days to 7.7 hours, while numerically computed errors remained the same.Source: §1 · p. 3
2 Bayesian inverse problems with stochastic forward solvers
The section formulates Bayesian inverse problems with stochastic forward solvers and explains how their random outputs induce approximate posterior targets. It then introduces pseudo-marginal sampling and density-tempered targets that vary solver accuracy alongside the tempering sequence.
2 Bayesian inverse problems with stochastic forward solvers
- A Bayesian inverse problem seeks the conditional distribution of an unknown parameter given noisy data generated by a forward map.The section uses a prior on the parameter and Gaussian additive observation noise.Source: §2 · p. 4
2.1 Approximate posteriors using randomized likelihoods
- A stochastic solver approximates the forward map with latent randomness governed by an accuracy-indexed noise model, and its marginal posterior need not equal the true posterior.The accuracy parameter typically increases both solver accuracy and computational cost; the applications use unbiased additive solvers, though other estimators may be biased.Source: §2.1 · pp. 4–5
- Pseudo-marginal methods permit sampling from the marginal approximation without explicitly constructing the typically intractable noise model.The paper’s hierarchical schedule evolves through targets associated with different solver accuracies and ends at a selected marginal posterior approximation.Source: §2.1 · p. 5
2.2 Density-tempering for pseudo-marginal SMC
- Density-tempered pseudo-marginal SMC uses an extended-space sequence of targets, and the proposed extension varies solver accuracy across that sequence.The sequence bridges prior and posterior targets, while tempering can become expensive when both the solver accuracy and number of stages are large.Source: §2.2 · p. 6
3 Hierarchical tempering
Hierarchical density-tempered SMC increases solver accuracy along the tempering sequence and uses cheap and expensive kernels for mutation. The section reports lower costs across tested regimes and an asymptotic complexity gain under additional assumptions.
3 Hierarchical tempering
- Solver accuracies are nondecreasing across tempering stages; equal accuracies with no cheap-kernel transitions recover the standard SMC sampler.When accuracy increases, the update uses classical importance sampling if noise-kernel densities are known, and Generalized IS otherwise.Source: §3 · pp. 6–7
- After importance sampling and resampling, mutation applies cheap-kernel transitions based on lower accuracy, followed by expensive-kernel transitions at the current accuracy.The algorithm permits both transition types at a stage.Source: §3 · pp. 6–7
- Under a geometric solver hierarchy and additional assumptions, Theorem 3.1 states an asymptotic complexity gain over standard SMC.The reduction stems from using cheaper solver accuracies at earlier tempering stages.Source: §3 · p. 7
- Hierarchical tempering with Generalized IS consistently reduces cost across tested regimes, while additional two-stage mutation savings depend on the mutation budget and cheap-chain allocation.The mutation-budget allocation balances transitions using a cheap kernel and an expensive kernel.Source: §3 · p. 7
- The mutation analysis restricts the kernel hierarchy to one cheap and one expensive kernel, although further complexity reduction could be possible with more cheap kernels.This restriction yields a simple, interpretable budget-allocation problem.Source: §3 · pp. 7–8
4 Hierarchical acceleration along the tempering sequence
The section analyzes the cost and error guarantees of hierarchical tempering, then extends solver-accuracy changes to a general importance-sampling setting. Its complexity results compare hierarchical and standard SMC under stated assumptions.
4.1 Complexity analysis assuming known densities
- The initial complexity analysis assumes solver-noise kernels have known, positive densities relative to a common measure.Source: §4.1 · p. 8
- The finite-sample guarantee is derived under A1–A3 for a suitably increasing sequence of solver accuracies.Theorem 4.2 gives bounds for both standard and hierarchical SMC.Source: §4.1 · p. 9
- Hierarchical and standard SMC attain the same error estimate, while hierarchical tempering yields an asymptotic computational gain under the theorem’s assumptions.The result assumes a spectral-gap lower bound independent of the number of tempering stages.Source: §4.1 · p. 10
4.2 The general case
- When such densities are unavailable, generalized importance sampling handles accuracy changes; its practical approximations are heuristic and their added error is unanalyzed.The accuracy change is separated from the tempering change, with no resampling or mutation between them.Source: §4.2 · pp. 11–12
5 Hierarchical acceleration of mutation
This section studies hierarchical acceleration of mutation by using cheaper approximate kernels before more accurate, expensive ones, and analyzes when that allocation is effective. It verifies the mixing assumptions for PMMH and an extended MCWM kernel, while noting that the resulting budget regimes depend on controlling MCWM’s invariant-measure discrepancy.
5 Hierarchical acceleration of mutation
- Mutation likelihood evaluations are costly, and the method uses cheap and expensive mutation blocks to address particle-weight degeneracy without relying solely on more tempering stages.The section identifies likelihood evaluations during mutation as a major SMC cost and adopts the strategy of more effective particle mutation through cheap and expensive kernels.Source: §5 · p. 13
5.1 Accelerated mutation analysis
- Under geometric-ergodicity assumptions, Theorem 5.1 bounds the product-kernel mutation error when particles move first toward the cheap-kernel target and then toward the expensive-kernel target.The analysis assumes a common suitable function V for the cheap and expensive kernels; the staged strategy is motivated when the cheap target approximates the expensive target and cheap transitions cost less.Source: §5.1 · p. 14
- The optimal mutation-budget allocation depends on cost efficiency and relative bias: sufficiently cheap transitions can justify a more biased cheap chain, while the combined-chain regime keeps its optimal cheap-chain length fixed as budget grows.The cost-efficiency threshold is q < τ−(1−θ), equivalent to τc < τθ; in the combined-chain regime, additional budget goes to the expensive chain, while sufficiently reduced budgets can shift all spending to the cheap chain.Source: §5.1 · p. 15
5.2 Application to SMC for BIPs with randomized likelihoods
- PMMH is the expensive and MCWM the cheap kernel; under Assumption 5.3, both are uniformly ergodic, and PMMH’s spectral-gap lower bound is δD > 0, independent of S, t, T, and the tempering schedule.MCWM is cheaper than PMMH only if Sc < S/2; the paper uses values as small as Sc = S/20.Source: §5.2 · pp. 16–17
- The PMMH and extended MCWM kernels are uniformly ergodic and satisfy the kernel and invariant-measure assumptions of Theorem 5.1, with K = 2, τc = τ = 1 − δD, and V ≡ 1.The extension places MCWM on the same extended state space as PMMH.Source: §5.2 · p. 17
- The MCWM extension is theoretical; in implementation, the auxiliary variable at accuracy S need only be generated once after the cheap block and before the first PMMH transition.The cheap MCWM update does not use the auxiliary coordinate.Source: §5.2 · p. 18
- The mixing result does not control the discrepancy between the MCWM and PMMH invariant measures, so Theorem 5.2’s budget regimes require a perturbation bound or an empirical discrepancy estimate.The discrepancy is expressed as f WV (µc,N, µe,N) in the passage.Source: §5.2 · p. 18
- Assumption 5.3 is sufficient rather than necessary; for an unbounded parameter space with a random-walk proposal, applying Theorem 5.1 would require drift and small-set conditions, a common unbounded Lyapunov function, and uniform tail control.The passage says this extension is not pursued.Source: §5.2 · p. 18
6 Applications
The applications evaluate hierarchical SMC on controlled linear and nonlinear photoacoustic inverse problems, examining computational cost and posterior reconstruction accuracy. They also illustrate solver-accuracy adaptation during SMC and cost reductions on challenging nonlinear problems.
6.1 A linear problem
- For every tempering schedule and γ ≥ 1 in the linear setting, assumptions (A1)–(A3) hold, and the PMMH kernels have absolute spectral gaps uniformly bounded below.The Gaussian random-walk proposal conditioned to remain in U yields the stated reversibility and gap properties.Source: §6.1 · pp. 18–19
- Approximately 3-fold cost savings occur for d = 8, and approximately 5-fold savings occur for d = 64 when Monte Carlo error dominates; generalized IS retains matching accuracy in the higher-dimensional test.For d = 64, generalized IS has a maximum relative error discrepancy of about 10%, while classical IS breaks down at N = 100 due to importance-weight collapse.Source: §6.1 · p. 19, Figure 1 · p. 20
- At fixed total cost, hierarchical SMC reduces error across all three tempering schedules, while two-stage mutation adds substantial gains at small budgets for T = 24 and 37.For T = 50, a cheap-chain cap of 100 adds little improvement, whereas a cap of 1000 yields clear additional gains; benefits depend on mutation budget and cheap-chain allocation.Source: Figure 2 · p. 21, §6.1 · pp. 21–22
6.2 A nonlinear application in quantitative photoacoustic imaging
- The inferred log-absorption reconstruction is best near the domain boundaries where sources are located, and posterior marginals and pairwise coefficient dependence are compared with prior distributions.The displayed comparison uses the worst-case hierarchical SMC realization by mean-marginal W1 error among the ten runs.Source: Figure 5 · p. 24, §6.2 · pp. 24–25, Figure 6 · p. 25
- Reducing photon counts further in the cheap mutation kernel gives limited wall-clock savings when the GPU is underutilized, though selective two-stage mutation at later stages may help if it saves sufficient time.Later stages use higher photon counts, where cheap-kernel savings may be more substantial.Source: §6.2 · p. 24
- Four of ten adaptive hierarchical SMC runs reach Smax = 10^9 photons, while mean photon counts remain below that maximum over most of tempering and rise toward β = 1.Solver accuracy starts at 10^5 photons and is multiplied by 10 after a mutation stage when mean acceptance falls below 0.15, subject to Smax.Source: §6.2 · p. 24
- 0.1016 mean and 0.1413 maximum mean-marginal W1 errors accompanied 0.0177 mean and 0.0282 maximum relative root mean squared errors across ten adaptive realizations.The largest observed posterior-mean discrepancy was approximately 2.8% relative to the standard-SMC reference, and even the worst-case comparison remained visually close.Source: §6.2 · p. 25
- 14.3× speedup reduced mean inference wall-clock time from about 4.6 days to 7.7 hours versus standard SMC using 109 photons throughout.For the four hierarchical realizations reaching the maximum photon count Smax = 109, the reported speedup was 6.4×; simulations used an NVIDIA A40 GPU.Source: §6.2 · p. 25
- The experiment shows that hierarchical SMC can substantially reduce computational cost for a challenging nonlinear inference problem while preserving accuracy.Source: §6.2 · p. 25
7 Conclusions
The paper develops hierarchical SMC for Bayesian inverse problems with stochastic forward solvers, using increasing solver accuracy across tempering and combining kernels within mutation stages. In photoacoustic imaging, it achieves a 14.3-fold speedup while producing posterior approximations close to the standard-SMC reference.
7 Conclusions
- Increasing solver accuracy along tempering gives an asymptotic cost reduction proportional to T under the stated assumptions, while two-stage mutation can add savings depending on its budget and cheap-chain allocation.The framework combines cheap approximate and expensive target-invariant kernels within a fixed mutation budget.Source: §7 · p. 25
- 14.3-fold speedup reduced photoacoustic inference runtime from approximately 4.6 days to 7.7 hours while keeping posterior approximations close to the standard-SMC reference.The conclusion presents this as the result on the nonlinear photoacoustic imaging problem.Source: §7 · p. 26
- A joint complexity analysis of both acceleration mechanisms and broader multilevel solver constructions remain open directions.The authors also identify possible extensions to deterministic discretizations, learned surrogates, and other SMC settings.Source: §7 · p. 26
A Proofs and Auxiliary Results of Section 5
The appendix develops auxiliary convergence results for SMC particle kernels and supplies proofs for results on mutation-kernel budgets and ergodicity. It uses product-kernel constructions and coupling-based Wasserstein arguments.
A Proofs and Auxiliary Results of Section 5
- The product kernel preserves the product invariant distribution, and the auxiliary result establishes geometric ergodicity under its stated assumptions.The proof uses coupling and Wasserstein contraction arguments.Source: A Proofs and Auxiliary Results of… · p. 29
- Applying a geometrically ergodic kernel to individual particles yields a product kernel that converges to the product target measure.The appendix defines the product kernel and target as products of the individual kernel and invariant distribution.Source: A Proofs and Auxiliary Results of… · pp. 29–30
- A one-step Doeblin minorization gives the PMMH and MCWM kernels unique invariant distributions and uniform ergodicity.The proof identifies the associated contraction parameters from the minorization bound.Source: A Proofs and Auxiliary Results of… · p. 31
B Proofs of Section 6
These proofs verify conditions used in the paper’s complexity analysis for stochastic solvers, including bounds on likelihood approximation error. They also use compactness to bound the solver-noise covariance uniformly.
B Proofs of Section 6
- A Taylor expansion and the symmetry of the even noise density support the likelihood-approximation error bound used in verifying Assumption (A3).The proof also invokes a uniform moment bound for the noise.Source: B Proofs of Section 6 · pp. 32–33
- Compactness of the parameter set and continuity of the covariance map yield uniform positive lower and finite upper eigenvalue bounds.Source: B Proofs of Section 6 · p. 32
- For S ≥ R^-2, the proof establishes Assumption (A3) with C_B = 2 exp(KΦ) C_L.The preceding bound on the solver-noise moment is uniform in u and S.Source: B Proofs of Section 6 · p. 33