Dual-space posterior sampling for Bayesian inference in constrained inverse problems

This paper proposes a dual-space posterior sampling framework that integrates the Alternating Direction Method of Multipliers (ADMM) with Stein Variational Gradient Descent (SVGD) to perform Bayesian inference in constrained inverse problems, effectively translating hard physical constraints into relaxable penalties that enable well-calibrated uncertainty quantification even with noisy and incomplete data.

The Gist

Imagine you are a detective trying to solve a mystery: What does the inside of the Earth look like?

You can't dig a hole to the center of the planet, so instead, you listen to the echoes of earthquakes (seismic waves) bouncing off underground rocks. This is called an Inverse Problem. You have the sound (the data), and you need to figure out what made the sound (the underground structure).

The Problem: A Noisy, Tricky Puzzle

The problem is that the data is messy. It's like trying to hear a whisper in a crowded stadium. Plus, the math is "ill-conditioned," meaning there are thousands of different underground maps that could explain the same sound. If you just pick one map, you might be wrong, and you won't know how confident you should be.

Bayesian Inference is the solution. Instead of guessing one map, you generate a whole cloud of possible maps (a distribution). This cloud shows you not just the most likely answer, but also the "what-ifs." It tells you, "We are 95% sure the rock is here, but there's a small chance it's over there."

The Big Hurdle: The "Physics Rule"

Here's the catch: The Earth follows strict physical laws (the Wave Equation). Any map you guess must obey these laws. If your map says sound travels faster than light, it's physically impossible.

In traditional methods, you have to force every single guess to obey these laws perfectly before you even check if it fits the data. This is like trying to solve a maze while blindfolded, but you have to touch the walls perfectly at every step. It's slow, stiff, and often gets stuck.

Other methods try to "relax" the rules, letting the guesses break the laws a little bit to make the math easier. But then, you have to guess how much to relax the rules. If you relax them too much, your map is nonsense. If you don't relax them enough, the math breaks.

The Solution: ADMM-SVGD (The "Dual-Space" Dance)

This paper introduces a clever new method called ADMM-SVGD. Let's break it down with an analogy.

Imagine you are trying to find the best route through a city (the solution) while obeying traffic laws (the physics constraints).

1. The Team of Explorers (SVGD): Instead of sending one detective, you send a team of 1,000 explorers. They all start at different random spots. They talk to each other. If they all crowd into one dead-end, they push each other apart (to keep exploring different areas). If they find a good path, they all move toward it. This is Stein Variational Gradient Descent (SVGD). It's a smart way to explore all possibilities at once.

2. The Referee (ADMM): The problem is that the explorers keep breaking traffic laws. Enter the Referee (the Augmented Lagrangian method).

   • In the old way, the referee would stop the explorers immediately if they broke a rule, making the whole process grind to a halt.
   • In this new Dual-Space method, the referee is more flexible. He says, "Okay, you broke the rule, but I'm going to give you a penalty ticket (a multiplier)."
   • The explorers keep moving, but the penalty ticket gets bigger every time they break the rule.
   • Over time, the penalty gets so huge that breaking the rule becomes impossible. The explorers naturally learn to obey the laws without the referee stopping them at every single step.

Why is this a Game-Changer?

• It's Flexible: The explorers can wander a bit in the beginning (when the model is far from the truth) without getting stuck. They only get strictly forced to obey the laws at the very end.
• It's Honest: Because we have a whole team of explorers, we get a full picture of the uncertainty. We can see, "Oh, in this deep valley, the explorers are all over the place. We really don't know what's down there."
• It Works: The authors tested this on a fake mountain (Rosenbrock) and a real, complex geological map (Marmousi II).
  • They found that as they added more data (more "ears" listening), the cloud of possible maps shrank and focused on the truth.
  • They found that even with noisy data, the method didn't panic; it just said, "We are less sure, but here is the range of possibilities."

The Bottom Line

This paper gives scientists a new, smarter way to map the Earth's interior. Instead of forcing a single, rigid answer, it uses a team of explorers guided by a smart referee to find the truth while respecting the laws of physics. It tells us not just where the oil or water might be, but how sure we are about it, which is crucial for making safe and smart decisions.

Technical Summary

1. Problem Statement

Inverse problems constrained by Partial Differential Equations (PDEs), such as Full Waveform Inversion (FWI) in geophysics, are notoriously ill-conditioned and non-unique due to noisy, incomplete data and inherent physical limitations (e.g., shadow zones).

• The Challenge: Traditional deterministic approaches yield a single point estimate, ignoring solution non-uniqueness. Bayesian inference offers a solution by characterizing the posterior distribution, but sampling from this distribution in high-dimensional, PDE-constrained spaces is computationally prohibitive.
• The Bottleneck: Existing methods face a trade-off:
  • Reduced-space formulations: Enforce the PDE exactly at every iteration (e.g., $u = A(m)^{-1}b$). This creates a highly nonlinear, multimodal posterior landscape with spurious local modes, making sampling difficult.
  • Penalty-based methods: Relax the PDE constraint into a soft penalty. While this improves conditioning, it requires careful tuning of penalty parameters and often fails to satisfy the physics constraint exactly at convergence.
  • Gaussian approximations: Some recent works relax constraints and approximate the posterior as Gaussian. However, this fails to capture complex, multimodal, or non-Gaussian uncertainties common in geology.

2. Methodology: ADMM-SVGD

The authors propose ADMM-SVGD, a novel framework that combines the Alternating Direction Method of Multipliers (ADMM) with Stein Variational Gradient Descent (SVGD) to perform posterior sampling in a dual space.

Core Concept: Augmented Lagrangian Probabilistic Interpretation

Instead of enforcing the PDE constraint $A(m)u = b$ exactly at every step, the method treats the constraint via an Augmented Lagrangian formulation.

• Relaxation: The wave equation is relaxed, treating the model parameters ($m$) and wavefields ($u$) as independent variables.
• Probabilistic Bridge: The Augmented Lagrangian function is interpreted as the negative log-density of an unnormalized probability distribution. As the Lagrange multipliers (dual variables) are updated, this distribution evolves and converges to the true posterior distribution.
• Dual-Space Advantage: By working in the dual space, the method decouples the model and wavefield variables, significantly improving the conditioning of the subproblems compared to reduced-space approaches.

Algorithmic Workflow

The algorithm iterates through four main steps for an ensemble of $N_p$ particles:

1. Auxiliary Variable Update: For each particle, solve for auxiliary wavefields ($u$) and scaled multipliers ($\epsilon$) using the optimality conditions of the augmented Lagrangian. This involves solving linear systems involving the Helmholtz operator but avoids explicit adjoint solves required in reduced-space methods.
2. Gradient Computation: Compute the gradient of the log-posterior with respect to the model $m$. This gradient is derived directly from the dual iteration equations (combining data misfit, constraint penalty, and prior).
3. SVGD Update: Update the model particles using the SVGD perturbation direction. This involves:
   • A gradient term driving particles toward high-probability regions (fitting data and prior).
   • A repulsive kernel term (RBF kernel) that maintains diversity among particles, preventing mode collapse and capturing multimodal structures.
4. Multiplier Update: Update the scaled Lagrange multipliers by accumulating the PDE residuals. This progressively enforces the physics constraint, ensuring that at convergence, the samples satisfy the wave equation exactly.

Key Features:

• Parallelism: The auxiliary solves and particle updates are "embarrassingly parallel" across particles and sources.
• No Fixed Penalty: Unlike standard penalty methods, the multiplier updates progressively tighten the constraint, removing the need for manual calibration of a large penalty parameter.
• Non-parametric: SVGD represents the posterior as an ensemble of particles, capturing arbitrary geometries (multimodality, non-Gaussian tails) without assuming a specific distribution shape.

3. Key Contributions

1. Dual-Space Sampling Framework: The first systematic procedure to translate hard PDE constraints into a Bayesian sampling framework using an augmented Lagrangian formulation.
2. Integration of ADMM and SVGD: A seamless coupling of ADMM (for constraint handling) and SVGD (for uncertainty quantification), leveraging the favorable conditioning of dual-space solvers while enabling full posterior characterization.
3. Exact Constraint Satisfaction: The method ensures that the PDE constraint is satisfied exactly at convergence through multiplier updates, avoiding the approximation errors inherent in fixed-penalty approaches.
4. Handling Multimodality: By avoiding Gaussian approximations, the method successfully captures complex, multimodal posterior distributions that arise in geologically complex settings.

4. Results and Validation

The method was validated on three problems of increasing complexity:

A. Rosenbrock Conditional Inference (Stylized Test)

• Setup: Sampling from a Rosenbrock distribution (banana-shaped) under a nonlinear constraint.
• Findings: ADMM-SVGD produced posterior samples that adhered strictly to the nonlinear constraint manifold ($z \approx x_1^2$), whereas standard SVGD (without constraint splitting) exhibited broader spread and required smaller step sizes and more iterations to remain stable.

B. Gaussian Anomaly Model (FWI)

• Setup: Inverting for a low-velocity anomaly in a homogeneous background using frequency-domain FWI.
• Findings:
  • Uncertainty Calibration: The method produced well-calibrated credible intervals (68%, 95%, 99%) that contained the true velocity model.
  • Spatial Resolution: Pointwise standard deviations correctly identified high uncertainty at anomaly boundaries and deep regions (poor illumination) and low uncertainty in well-constrained areas.

C. Marmousi II Benchmark (Complex Geology)

• Setup: Inverting a complex geological model with faults, anticlines, and lateral velocity variations using varying numbers of sources ($N_s = 17, 34, 68$).
• Findings:
  • Posterior Contraction: As data coverage increased (more sources), the posterior distribution contracted, and the relative model error (RME) decreased.
  • Multimodality: At geologically complex locations, the posterior exhibited multimodal structures, revealing distinct plausible velocity values that a single-point estimate would miss.
  • Noise Robustness: Experiments with 10–20% noise showed that the method gracefully degrades, with uncertainty increasing linearly with noise levels, and the conditional mean remaining geologically coherent.
  • Convergence: The algorithm demonstrated stable convergence with posterior contraction, confirming that the multiplier updates successfully enforce constraints without trapping particles in local modes.

5. Significance and Impact

• Bridging Optimization and Sampling: The paper successfully bridges the gap between deterministic constrained optimization (ADMM) and probabilistic inference (SVGD), offering a robust solution for high-dimensional inverse problems.
• Practical Uncertainty Quantification: It provides a computationally feasible way to quantify uncertainty in complex geophysical problems (like FWI) where traditional MCMC is too slow and Gaussian approximations are insufficient.
• Physical Consistency: By ensuring exact satisfaction of the wave equation at convergence, the generated samples are physically meaningful, addressing a major limitation of previous relaxed-constraint methods.
• Scalability: The embarrassingly parallel nature of the algorithm makes it well-suited for modern high-performance computing architectures, paving the way for Bayesian inference in large-scale 3D seismic imaging.

In conclusion, ADMM-SVGD offers a principled, efficient, and robust framework for Bayesian inference in constrained inverse problems, enabling the extraction of reliable uncertainty estimates from complex physical data.