Weak neural variational inference for solving Bayesian inverse problems without forward models: applications in elastography

TL;DR

This work tackles high-dimensional Bayesian inverse problems governed by PDEs, notably elastography, by introducing Weak Neural Variational Inference (WNVI) that uses weighted residuals as virtual observations to inject physics without solving the forward model. The approach jointly infers latent state variables and material-property fields through an approximate posterior $q_{\boldsymbol{\xi}}(\boldsymbol{x},\boldsymbol{y})$, optimized via Stochastic Variational Inference and a Monte Carlo sampling of residuals to form an ELBO that combines virtual and real data likelihoods. Numerical examples show WNVI achieves accuracy comparable to forward-solver–based methods but with orders of magnitude fewer residual evaluations, can handle ill-posed forward problems (e.g., missing Dirichlet data), and extends naturally to nonlinear material behavior. The method offers a practitioner-friendly, probabilistic inversion tool that provides uncertainty quantification and avoids the complexities of forward solvers and adjoint computations, with potential extensions to adaptive residual weighting and tempering for real-time diagnostics.

Abstract

In this paper, we introduce a novel, data-driven approach for solving high-dimensional Bayesian inverse problems based on partial differential equations (PDEs), called Weak Neural Variational Inference (WNVI). The method complements real measurements with virtual observations derived from the physical model. In particular, weighted residuals are employed as probes to the governing PDE in order to formulate and solve a Bayesian inverse problem without ever formulating nor solving a forward model. The formulation treats the state variables of the physical model as latent variables, inferred using Stochastic Variational Inference (SVI), along with the usual unknowns. The approximate posterior employed uses neural networks to approximate the inverse mapping from state variables to the unknowns. We illustrate the proposed method in a biomedical setting where we infer spatially varying material properties from noisy tissue deformation data. We demonstrate that WNVI is not only as accurate and more efficient than traditional methods that rely on repeatedly solving the (non)linear forward problem as a black-box, but it can also handle ill-posed forward problems (e.g., with insufficient boundary conditions).

Figures (11)

• Figure 1: Schematic illustration of the proposed formulation. The actual data $\hat{\boldsymbol{u}}$ induces an actual likelihood $p(\hat{\boldsymbol{u}}| \boldsymbol{y})$ on the (representation of the) PDE-solution $\boldsymbol{y}$. This is combined with the virtual likelihood $p(\hat{\boldsymbol{R}}| \boldsymbol{y}, \boldsymbol{x})$ that assigns probabilities to $(\boldsymbol{x},\boldsymbol{y})$ pairs. The virtual likelihood is highest on the solution manifold, i.e., for all $(\boldsymbol{x},\boldsymbol{y})$ pairs that yield $0$ values for the weighted residuals $r(\boldsymbol{y},\boldsymbol{x})$ and decays at a rate determined by the precision $\lambda$ as one moves away from it. The combination of these two likelihoods (and the priors which are not depicted here) gives rise to the sought posterior $p(\boldsymbol{x} |\hat{\boldsymbol{u}}, \hat{\boldsymbol{R}})$ on (the representation of) the unknown material field $\boldsymbol{x}$.
• Figure 2: Illustration of 4 randomly selected weight functions $\boldsymbol{w}^{(j)}$ used in the algorithm. Note that for each weight function depicted are only the dimensions $w^{(j)}_1$ or $w^{(j)}_2$ which have non-zero values.
• Figure 3: (a): Problem configuration with Dirichlet BCs on boundaries $\Gamma_\mathrm{top}$ and $\Gamma_\mathrm{left}$ and Neumann BCs on boundaries $\Gamma_\mathrm{bot}$ and $\Gamma_\mathrm{right}$, respectively. The domain shows the ground truth material property field $m$ (e.g., the logarithm of Youngs modulus) with two circular inclusions. (b): Boundary notation convention.
• Figure 4: Convergence of the proposed method. Depicted is the change of the ELBO $\mathcal{L}$ (left) and the expected squared residual $\left< r^2 \right>$ (right) over the iterations. The expected residuals are calculated using $50$$\boldsymbol{x}$-$\boldsymbol{y}$ samples and all weight functions. The fluctuations are due to the Monte Carlo estimates employed at each iteration.
• Figure 5: The three columns contain posterior mean (first), posterior standard deviation (second) over the whole domain, and posterior estimates along the diagonal $s_1=s_2$ of the problem domain (i.e. posterior mean $\mu$ (blue line) and $95 \%$ credible intervals (blue shaded areas)). The first row pertains to the ground truth posterior obtained by black-box HMC run ($1.35 \times 10^{11}$ residual evaluations). The second row corresponds to the posterior obtained by the proposed method after $4 \times 10^8$ weighted residual evaluations. The third and fourth rows contain posterior estimates obtained with black-box HMC / SVI, respectively, and after the same number of residual evaluations as the proposed method. The results were obtained with the same synthetic data contaminated by noise with SNR$=30$ dB.