We consider the problem of sampling from a posterior distribution arising in Bayesian inverse problems in science, engineering, and imaging. Our method belongs to the family of independence Metropolis-Hastings (IMH) sampling algorithms, which are common in Bayesian inference. Relying on the existence of an approximate posterior distribution that is cheaper to sample from but may have significant bias, we introduce Proximal-IMH, a scheme that removes this bias by correcting samples from the approximate posterior through an auxiliary optimization problem. This yields a local adjustment that trades off adherence to the exact model against stability around the approximate reference point. For idealized settings, we prove that the proximal correction tightens the match between approximate and exact posteriors, thereby improving acceptance rates and mixing. The method applies to both linear and nonlinear input-output operators and is particularly suitable for inverse problems where exact posterior sampling is too expensive. We present numerical experiments including multimodal and data-driven priors with nonlinear input-output operators. The results show that Proximal-IMH reliably outperforms existing IMH variants.
Figures & tables
Figure 1: Sensitivity of the expected KL divergence of the Approx , Latent , and Proximal proposals (relative to the Exact posterior ) for different noise levels, operator errors, observation ratios, and dimensions. Details of the experimental setup are given in Table 2 .
Figure 2: Comparison of sampling performance in the bimodal test. Rows correspond to Tests I–III. From left to right, panels show Metropolis–Hastings acceptance rates for three IMH methods, convergence of the relative mean error for the IMH methods together with NUTS and MALA, and histograms of the projected samples w⊤x for two levels of operator error. In the histogram panels, Approx-NUTS denotes samples drawn from the Approx posterior πa(x∣y) . For Proximal-IMH , we set β=σ2 . Results in the left panel are averaged over 5 independent trials for each method.
Figure 3: Operator discrepancy in the bimodal test. The quantities ∥I−F−1F∥ and ∥I−K−1∥ are reported for Latent-IMH and Proximal-IMH , respectively, where K is defined in Equation 7 . For Proximal-IMH , the hyperparameter is set to β=σ2 .
Figure 4: Effect of β on Proximal-IMH in the bimodal test. Values in parentheses represent ∥A−A∥/∥A∥ . Blue curves show acceptance rates (left y-axis), and pink curves show relative mean errors after 105 MH steps (right y-axis).
Figure 5: Ground truth, posterior mean and MAP estimates for the MNIST inverse problem in Tests IV–VII . The MAP estimates are obtained from different random initial guesses for the optimization problem. We observe that in several cases the reconstructed MAP point is wrong. This is due to the nonconvexity of solving for a MAP point.
Figure 6: Sampling performance on MNIST imaging inverse problems: convergence of the relative mean error and relative second-moment error for different sampling methods. Dashed lines indicate the bias induced by the Approx posterior πa(x∣y) . For Proximal-IMH , we set β=σ2 . Results in the convergence plots are reported over 5 independent trials.
Figure 7: Effect of the hyperparameter β on Proximal-IMH in the inverse problem in MNIST. Values in parentheses represent ∥A−A∥/∥A∥ . Blue curves show acceptance rates (left y-axis), and pink curves show relative mean errors after 105 MH steps (right y-axis). The results are not sensitive to β .
Figure 8: Nonlinear Helmholtz inverse problem. Ground truth, MAP estimate, and posterior mean for the acoustic scattering test, together with convergence of the relative mean error.
Appendix figures & tables9 assets
Supplementary material from the paper’s appendix.
Appendix
Posterior
Distribution
Mean
Covariance
Pseudo-inverse
Exact
π(x∣y)=N(μ,Σ)
μ=A†y
Σ=I−A†A
A†=A⊤(AA⊤+σ2I)−1
Approx
πa(x∣y)=N(μa,Σa)
μa=Aa†y
Σa=I−Aa†A
Aa†=A⊤(AA⊤+σ2I)−1
Latent
πl(x∣y)=N(μl,Σl)
μl=Al†y
Σl=F−1FΣaF⊤F−⊤
Al†=F−1FAa†
Proximal
πp(x∣y)=N(μp,Σp)
μp=Ap†y
Σp=KΣaK⊤
Ap†=KAa†
Appendix
Table 1: Distributions of Exact , Approx , Latent and Proximal posteriors for Gaussian prior p(x)=N(0,I) and noise q(e)=N(0,σ2I) . For Proximal posterior , K=(A⊤A+βI)−1(A⊤A+βI) .
Effects
log10SNR
∥F∥∥F−F∥
dy/d
d
Noise level
0.5 —4.0
6%
0.2
500
Operator error
2.5
2% —21%
0.2
500
Observation ratio
2.5
6%
0.05 —0.5
500
Dimension
2.5
6%
0.2
100 —2000
Appendix
Table 2: Experimental setup for comparing Da , Dl and Dp .
Figure 9: Histogram of the projection of prior samples in the bimodal test (see Equation 21 ).
dz
dx
du
du
dy
128
784 ( =282 )
104 ( =1002 )
{352,202}
100
Appendix
Table 3: Dimensions of parameters in Tests VI and VII in Section 4.2 .
Figure 10: Matrix constructions of A and A in Tests VI and VII in Section 4.2 . Left: A=OFP . Right: A=OPuFPx . In our setting dy<dx<du , and du<du .
Figure 11: Jacobian log-determinant diagnostics for the Gauss–Newton correction. Each row corresponds to one test case. Left: Sorted values of the log determinant of the Jacobian of a single Gauss–Newton transformation, evaluated at 2000 randomly selected posterior samples from πa(x∣y) . Right: Order statistics of ratios of Jacobian determinants between 2000 randomly selected pairs of posterior samples from πa(x∣y) , i.e., detJGN(xi)/detJGN(xj) . The numbers shown in the legends indicate the relative operator approximation error ∥A−A∥/∥A∥ . In the right panels, dashed lines denote the 5%–95% quantile range.
dx
du
du
dy
1024 ( =322 )
16,384 ( =1282 )
4,096 ( =642 )
512
Appendix
Table 4: Dimensions of parameters in the nonlinear Helmholtz test in Section 4.3 .
Figure 12: Effect of preconditioning on GMRES convergence for the nonlinear Helmholtz test.
Figure 13: Wavefields {ui=L(x)−1fi}i=14 for the nonlinear Helmholtz test. Each column corresponds to a distinct source fi . Top: Ground-truth wavefields. Middle: Wavefields obtained from the MAP estimate. Bottom: Posterior mean wavefields computed from posterior samples.
Sampling from the posterior is the central task of computational Bayesian inverse problems. The standard workhorse in Bayesian inference - Markov chain Monte Carlo (MCMC) - is sequential, yields correlated samples, and must be rerun for each observation. Conditional flow matching offers an amortized alternative: a transport map, trained once on joint samples of parameter and data, that yields independent approximate posterior samples for any observation at negligible online cost, without new likelihood evaluations. We give a careful, MCMC-literate assessment of flow matching for PDE-based inverse problems with function-valued parameters. Exploiting the flow's tractable density, we derive computable accuracy estimates of the underlying approximate posterior in total-variation distance and Kullback-Leibler divergence and, moreover, propose a hybrid sampler that is asymptotically exact by Metropolization. We validate the accuracy estimates and demonstrate the amortization in several numerical examples, including electrical impedance tomography and a likelihood-free Lotka-Volterra model.
Jan Blechschmidt, Oliver G. Ernst, Moritz Poguntke +1
Diffusion models have excellent capacity to model complex distributions of natural data, which has made them a popular and effective choice for posterior sampling in imaging inverse problems. Existing methods can incorporate any measurement model at inference time but must use an inexact approximation for the likelihood at intermediate timesteps for computational tractability. Although these approximations can often work well empirically, their downstream effect on the sampled posterior is poorly understood and can result in unexplained failures. To understand when, why, and how these likelihood approximations propagate to erroneous posterior distributions, we introduce a finite-sample perspective on posterior sampling that approximates the posterior to arbitrary precision as training set size tends towards infinity, for any forward model and prior distribution. Using this finite-sample lens, we observe that popular posterior sampling approximations tend to under- or over-estimate the spread of the posterior at intermediate timesteps, causing downstream consequences including sensitivity to early stopping time, inaccurate relative weighting of posterior modes, and hallucination, both of prior modes that are not in the posterior and likelihood modes that are not supported by the prior. Moreover, we find that the cause of these posterior errors requires neither a nonlinear measurement model nor a multimodal posterior, but can arise solely due to a multimodal prior and inaccurate posterior spread at intermediate sampling times. Our finite-sample posterior sampling approach is agnostic to the type of likelihood approximation and the type of (linear or nonlinear) forward model, and can thus serve as a drop-in diagnostic to evaluate the accuracy and failure modes of existing and future posterior samplers.
This work introduces a new approximate proximal sampler that operates solely with zeroth-order information of the potential function. Prior theoretical analyses have revealed that proximal sampling corresponds to alternating forward and backward iterations of the heat flow. The backward step was originally implemented by rejection sampling, whereas we directly simulate the dynamics. Unlike diffusion-based sampling methods that estimate scores via learned models or by invoking auxiliary samplers, our method treats the intermediate particle distribution as a Gaussian mixture, thereby yielding a Monte Carlo score estimator from directly samplable distributions. Theoretically, when the score estimation error is sufficiently controlled, our method inherits the exponential convergence of proximal sampling under isoperimetric conditions on the target distribution. In practice, the algorithm avoids rejection sampling, permits flexible step sizes, and runs with a deterministic runtime budget. Numerical experiments demonstrate that our approach converges rapidly to the target distribution, driven by interactions among multiple particles and by exploiting parallel computation.
Hirohane Takagi, Atsushi Nitanda
Graduate School of Information Science and Technology, The University of Tokyo, Japan · Institute of High Performance Computing, Agency for Science, Technology and Research, Singapore · Centre for Frontier AI Research, Agency for Science, Technology and Research, Singapore +1