WorksMathematics

Inverse domain uncertainty quantification in hemodynamics: sampling virtual patient cohorts under constraints

Francesco Romor, Jose Rodrigo Rojo-Garcia, André-Alexander Zepernick, Ana Djurdjevac, Alfonso Caiazzo

Published Oct 20, 2026Fresh — added to the Hub within the last 30 daysv149 min read

AI contributions

Claude Code

Abstract

We consider the following inverse domain uncertainty quantification problem: given a prescribed value or measurement of a hemodynamic biomarker, characterize the distribution of vascular geometries compatible with it. The problem is formulated as a Bayesian inverse problem in which a high-dimensional latent parameter encodes the cardiovascular shape and accounts for inter-patient variability. The parameter-to-observable map composes a conditional LDDMM stochastic interpolant generating the aortic surface with a high-performance numerical finite element solver for the unsteady incompressible Navier–Stokes equations, from which biomarkers of interest are evaluated. The map is nonlinear, high-dimensional and computationally expensive. To address the problem, we compare three ensemble strategies: ensemble adjusted Kalman inversion (EAKI), adaptive Bayesian updating with structural reliability (aBUS), and a Wasserstein–Fisher–Rao (WFR) sequential Monte Carlo scheme. The first two approaches are derivative-free. For the third method, we introduce a novel differentiable flow-matching surrogate that reduces the cost of a single forward map evaluation at the price of a distribution shift depending on the training distribution used for the surrogate. We apply the framework to the generation of virtual patient cohorts matching prescribed pressure drops along aortic geometries affected by coarctation, and we validate and compare the three methods on synthetic data.

Introduction

Computational blood flow models, in combination with medical imaging, can provide relevant insights into available medical data and support clinical decision making. Modern numerical methods, together with the constant increase in computational resources and the increasing quality of medical data, contributed to making in silico models more and more relevant for several clinical and preclinical applications, including non-invasive diagnosis and staging of cardiovascular pathologies, simulation of pathological conditions and medical devices, as well as treatment planning, optimization, and monitoring. Moreover, the possibility of generating patient-specific and disease-specific synthetic data makes computational models a valuable tool for the development and validation of scientific machine learning algorithms, especially in a field where availability of real data remains a relevant challenge.

However, the quality and reliability of these approaches depend heavily on the model parameters, on the quality of the images used to generate the computational domain, and on measurement accuracy. Uncertainty quantification (UQ) has therefore been increasingly studied for forward simulations, i.e., to assess the robustness of hemodynamic predictions with respect to uncertain model inputs, such as boundary conditions, physical, or geometric parameters (see, e.g., [1, 2, 3]). In particular, several recent works focused on studying the impact of geometric uncertainty, such as segmentation errors, on the simulation output, using machine learning methods to generate geometries from image data [4], perturbing a patient-specific aortic surface using Gaussian random fields [5], and considering a generative model, trained on a synthetic patient cohort, for generating patient-specific three-dimensional aortic geometries [6].

From the mathematical perspective, the propagation of geometrical uncertainties has been studied more generally in the context of partial differential equations (PDEs) posed on random domains, attracting growing interest in recent years. In the case of small boundary perturbations, methods from shape calculus and perturbation techniques have been exploited to derive approximations of statistical moments of the PDE solution [7, 8]. Alternatively, domain-mapping approaches describe each realization as a deformation of a fixed reference domain and transfer the geometric uncertainty to the coefficients and operators of the transformed PDE [9, 10]. The latter approach has provided a basis for rigorous regularity and convergence analyses for high-dimensional approximation techniques with respect to the parameters describing the geometry [11].

In this work, we consider an inverse domain uncertainty quantification problem, in which observations of the PDE solution are used to infer the underlying domain geometry [12, 13, 14]. In other words, rather than asking how uncertainty in the vascular geometry affects a given biomarker, the goal is to characterize the distribution of patient geometries compatible with a measured quantity of interest, taking also into account the possible uncertainty in the measurements. We formulate this task as a Bayesian inverse problem for uncertain computational domains.

Focusing on blood flow in the thoracic aorta, we introduce and validate a novel efficient framework for creating virtual patient cohorts with prescribed hemodynamic conditions, to support in silico trials, e.g., for the design and testing of medical devices, and to enable cardiovascular shape optimization under geometric uncertainty.

The Bayesian approach tackles both the quantification of uncertainty and the ill-posedness of the inverse problem. For the latter, the regularization of the problem is introduced by modeling the parameter of interest (i.e., the parameters describing the domain) as a random variable distributed according to a prior distribution (probability measure). Then, in conjunction with a so-called likelihood function, which quantifies the error between the results of the forward model and the available measurements, Bayes’ theorem is applied to determine a posterior distribution of the parameter [15, 16]. A rigorous analysis of well-posedness via Bayesian inversion can be found in [18] [17], where suitable conditions on the mathematical model and prior measure ensure the well-posedness of the posterior measure defined through Bayes’ theorem. The uncertainty quantification process is based on realizations drawn from the posterior measure. From these samples, moments and statistical estimators of the parameter and quantities of interest derived from the forward model can be inferred. Additionally, the samples provide information about the posterior distribution, such as multimodalities and correlations. For linear forward models, Gaussian prior distributions and Gaussian additive noise, the posterior distribution is also Gaussian with a closed-form expression for its hyperparameters [16]. In this setting, the sampling process becomes straightforward. This is no longer the case for more general nonlinear forward operators, for which classical approaches are based on sampling with Markov chain Monte Carlo methods (MCMC) (e.g., [19]). However, these techniques may require a high number of evaluations to adequately explore the target distribution, and the cost can be prohibitively high. Several adaptive methods have been proposed to tackle these issues (see, e.g., [20, 21, 22, 23]).

The inverse UQ problem is thus particularly challenging in the present setting. Firstly, one has to handle the high dimensionality of the parameter space (hundreds of degrees of freedom to describe the features of the aortic surface). Secondly, the map from latent variables to observations is highly nonlinear and involves shape generation, volumetric mesh deformation, and the numerical solution of the unsteady incompressible Navier–Stokes equations. Thirdly, a given observation of a biomarker may be compatible with multiple geometries, potentially leading to multimodal posterior distributions.

The high-dimensionality of the problem drastically increases the complexity of the sampling process. Two widely used approaches for sampling in high-dimensional spaces are Hamiltonian Monte Carlo (HMC), and Preconditioned Crank-Nicolson (pCN) (see, e.g., [24]). HMC requires gradient information of the forward model. In contrast, pCN is a derivative-free method, but it may still require a large number of forward model evaluations to adequately explore the posterior distribution. In all cases, the sequential nature of the MCMC methods limits their direct parallelization, and as a consequence, makes simulations within a reasonable time difficult for computationally costly forward models. A popular alternative to the methods described above is the Laplace approximation, which consists of a local Gaussian approximation around the Maximum a Posteriori point (MAP) [25, 26] [15]. By construction, the approximation is expected to perform well for unimodal posterior distributions with a local geometrical structure similar to a Gaussian density. On the other hand, this method requires sufficient differentiability of the likelihood function, which may not always be the case, and it may become computationally expensive in high-dimensional spaces.

Our framework builds on two recent works handling arbitrary non-parametric registration of aortic shapes [27] and generative modeling of geometries based on LDDMM stochastic interpolants [6], which allow generating aortic shape distributions from a latent representation of the geometry. Combined with a finite element solver for the blood flow in the resulting domain, we can then construct a forward map that links the latent geometrical parameter to selected observations on the fluid solution. In this setting, we investigate three complementary state-of-the-art approaches for high-dimensional Bayesian inversion to solve the inverse domain uncertainty quantification problem (see Figure 1): the Ensemble Adjusted Kalman Inversion (EAKI), the Adaptive Bayesian Updating with Structural reliability (aBUS), and the Wasserstein–Fisher–Rao sequential Monte Carlo (WFR).

Overview of the forward model and Bayesian inference framework, including latent parameters, aortic geometry, blood flow solution, predicted biomarker, prior distribution, Bayesian inference, posterior distribution, compatible aortic shapes, and target observation and uncertainty — Figure 1. Overview of the inverse domain uncertainty quantification framework.

Kalman filters are based on ensembles of particles that iteratively construct Gaussian approximations of the target distribution. The iterative process is derivative-free and typically requires a relatively small number of iterations and model evaluations. These methods were initially considered for Bayesian data assimilation of dynamical systems with uncertainty [28], and their use has been extended to inverse problems by imposing artificial dynamics, leading to Ensemble Kalman Inversion (EKI) methods (see, e.g., [29, 30, 31]). Inspired by EKI, the EAKI method [32] introduces suitable artificial dynamics specifically designed for Bayesian inference. Given a Gaussian prior and an ensemble generated from it, the iterative Gaussian densities converge to a Gaussian approximation of the posterior distribution, which, in the case of a linear forward model, coincides with the true posterior density.

An alternative approach to Bayesian inference is the Bayesian Updating with Structural reliability (BUS) methods [33], which recast the target posterior distribution as a rare event under a standard Gaussian density. This formulation requires rare-event simulation methods such as Subset Simulation (SuS) [34], which internally relies on MCMC methods. The adaptive BUS (aBUS) [35] relies on an adaptive strategy for the hyperparameters of BUS to improve the speed of convergence. Unlike the Kalman filter, aBUS does not rely on the assumption of an approximate Gaussian geometry or unimodality of the posterior distribution. Therefore, this approach can efficiently sample complex and multimodal distributions. Moreover, it is derivative-free and, although it internally uses MCMC methods, no burn-in truncation is required. In addition, this method can be efficiently parallelized [36].

Recently, particle-based sampling and optimization methods based on gradient flows have gained significant traction in the context of high-dimensional uncertainty quantification. While classical Langevin dynamics approximate Wasserstein gradient flows to transport particles via diffusion, they often suffer from severe metastability in non-convex or multimodal loss landscapes. Conversely, Fisher–Rao gradient flows enable non-local mass redistribution among different modes. Building on recent studies [37, 38], we consider a Wasserstein–Fisher–Rao (WFR) sequential Monte Carlo (SMC) method that combines the two gradient flows (metrics). With respect to MCMC methods (e.g., Hamiltonian Monte Carlo), the proposed approach has a better performance in cases of severe multimodality, such as the one considered in this work (see section 4.2.2). Unlike the EAKI and the aBUS, the Wasserstein step of the WFR requires the evaluation of the gradients of the forward map. In the current application, where the forward map combines several complex steps, including generative modeling for surface deformation, nonlinear extension of the bulk mesh, and solution of the incompressible Navier–Stokes equations, we circumvent this issue by training a surrogate model that approximates the map and provides cheap batched gradient evaluations on GPUs. The machine learning model chosen for this surrogate is an autoregressive flow-matching (stochastic interpolant) model with a transformer-based U-Net, similar to [39], but extended to the case of variable geometries. To this purpose, we utilize an LDDMM registration pre-processing step [27], which provides meshes with the same topology on different shapes, a natural geometrical encoding that speeds up the training of such heavy architectures.

We test the proposed methods in two different settings: using observations of the pressure drop along the whole vessel centerline and specifying only the maximum value of the pressure drop along the vessel centerline. The first problem is solved with EAKI, as the observed/specified values of the pressure drop along the vessel centerline are informative enough to reconstruct a single geometry, while aBUS and WFR are used for the second case for sampling multiple geometries associated with the same observed/prescribed value of the maximum pressure drop along the vessel centerline. The results and main contributions of this work can be summarized as follows.

We apply three approaches for domain UQ in high-dimensional parameter spaces that do not require gradients of the high-fidelity forward map, comparing them in terms of efficiency and sampling quality.
We show that, even when high-performance numerical solvers are used on relatively coarse meshes, a surrogate forward map provides a remarkable speedup at the price of a small loss of accuracy for the inverse problem under study.
We show that, for topologically equivalent shapes, tailored generative priors and surrogate models that exploit the geometrical encoding provided by registration maps represent a novel and efficient paradigm.

The rest of the paper is organized as follows. Section 2 formulates the inverse domain UQ problem and introduces the three considered methods in detail. Section 3 describes the forward map which defines the aortic shape inversion problem, as well as the surrogate version used for WFR. Section 4 presents the numerical results considering a given pressure drop along the aorta centerline and a given maximum pressure drop. Finally, section 5 draws the conclusions of our study.

Inverse domain uncertainty quantification

Problem setting

In what follows, let $\Omega_\theta\subset\mathbb{R}^{d}$ ( $d = 2, 3$ ) denote a computational domain, depending on a parameter vector $\theta\in\mathcal{X} \subset\mathbb{R}^{n_\theta}$ , $n_\theta\in\mathbb{N}$ , where $\mathcal{X}$ stands for a suitable parameter space. The relation between $\theta$ and $\Omega_\theta$ may correspond to a standard parametrization, such as a free-form deformation or a perturbation induced by a truncated Karhunen–Loève expansion of a covariance kernel specified on the domain, or, as in the case considered in this work, to a generative model (generative prior) [40], in which $\theta$ plays the role of a latent variable.

On the domain $\Omega_\theta$ , we consider a PDE in the general form

\mathcal{R}(\Omega_\theta; s_\theta) = 0, \tag*{(1)}

where the solution $s_\theta\in V(\Omega_\theta)$ belongs to an abstract function space (e.g., $s_\theta= (u_\theta, p_\theta) \in H^1(\Omega_\theta) \times L^2(\Omega_\theta)$ for the velocity and pressure fields of the Navier–Stokes equations). We also define an observation operator $\mathcal{O}(\theta; \cdot) : V(\Omega_\theta) \to\mathbb{R}^{n_y}$ , $s_\theta\mapsto\mathcal{O}(s_\theta)$ , for $n_y \in\mathbb{N}$ , which extracts selected quantities of interest from the PDE solution. The observation operator $\mathcal{O}(\theta; \cdot)$ is assumed to be bounded, not necessarily linear. For the sake of simplicity, we will omit the dependency on $\theta$ in the observation operator in the following description.

The domain uncertainty quantification problem consists in inferring the posterior distribution of $\theta$ given a vector $y \in\mathbb{R}^{n_y}$ of noisy observations or prescribed values of a quantity of interest. The solution corresponds to a probability distribution on the set of admissible domains $\{\Omega_\theta\}_{\theta\in\mathcal{X}}$ . Formally, we introduce a parameter-to-observable map

\mathcal{G} : \mathcal{X} \subset\mathbb{R}^{n_\theta} \to\mathbb{R}^{n_y}, \qquad\mathcal{G}(\theta) := \mathcal{O}(s_\theta), \tag*{(2)}

which associates with a particular realization of the parameter $\theta\in\mathcal{X}$ the observation of the corresponding PDE solution $s_\theta$ on the domain $\Omega_\theta$ . We model the observations as

y = \mathcal{G}(\theta) + \eta, \qquad\eta\sim\mathcal{N}(0,\Sigma_\eta), \tag*{(3)}

where $\Sigma_\eta\in\mathbb{R}^{n_y \times n_y}$ is a symmetric positive-definite noise covariance matrix.

Assuming a Gaussian prior $\theta\sim\mathcal{N}(m_0,C_0)$ on $\mathcal{X}$ , with corresponding measure $\mu_0$ , Bayes’ theorem yields the posterior measure $\mu^y$ , which is absolutely continuous with respect to $\mu_0$ and satisfies

d\mu^y = \frac{d\mu^y}{d\mu_0}d\mu_0, \tag*{(4)}

where

\frac{d\mu^y}{d\mu_0}(\theta) \propto\exp\left(-\frac{1}{2}\lVert y-\mathcal{G}(\theta)\rVert_{\Sigma_\eta^{-1}}^2\right), \tag*{(5)}

and the weighted norm is defined as $\lVert\cdot\rVert_{\Sigma_\eta^{-1}}^2 := \langle\cdot,\Sigma_\eta^{-1}\cdot\rangle$ .

In the application considered in this work, $\mathcal{G}$ is a nonlinear map acting on a high-dimensional parameter space, which requires the numerical solution of a computationally expensive PDE (see Figure 1). In such situations, the evaluation of the gradient is sometimes computationally prohibitive or undefined if the map is not guaranteed to be differentiable with respect to $\theta$ . We therefore focus on derivative-free methods tailored to high-dimensional parameter spaces: an ensemble Kalman inversion method (EAKI) that constructs a sequence of Gaussian approximations (§ 2.2), an adaptive Bayesian updating method with structural reliability (aBUS) (§ 2.3), and a Wasserstein–Fisher–Rao sequential Monte Carlo method (WFR, § 2.4) that leverages a geometry-conditioned LDDMM stochastic interpolant surrogate for efficient gradient evaluation (§ 3.7).

For ease of notation, we introduce the potential functions $\Phi_{\mathrm{like}},\Phi_{\mathrm{post}} : \mathcal{X} \times\mathbb{R}^{n_y} \to[0,\infty]$ , defined as

\Phi_{\mathrm{like}}(\theta,y) := \frac{1}{2}\lVert y-\mathcal{G}(\theta)\rVert_{\Sigma_\eta^{-1}}^2, \tag*{(6)}

and

\Phi_{\mathrm{post}}(\theta,y) := \Phi_{\mathrm{like}}(\theta,y) + \frac{1}{2}\lVert\theta-m_0\rVert_{C_0^{-1}}^2. \tag*{(7)}

Remark 2.1. These methods are not specifically designed for domain UQ problems: for example, they can be employed also when the parameter to be inferred $\theta$ modifies the PDE model, but the domain $\Omega_\theta$ remains fixed. Domain UQ problems are characterized by additional challenges, since the parameter-to-observable map requires, besides the resolution of a PDE, additional geometry and mesh processing steps, and a careful design of the shape space taking into account possible geometrical constraints. In this work, we consider only geometries with the same topology and for which the underlying PDE can be defined with the same type of boundary conditions. Domain UQ and shape optimization problems with varying topologies pose further challenges [41].

Ensemble Adjusted Kalman Inversion (EAKI)

Rather than sampling the posterior $\mu^y$ directly, the ensemble Kalman inversion (EKI) [29, 32, 42] recasts the inverse problem as a filtering problem for an artificial dynamical system, in which the single available observation $y$ is assimilated repeatedly. The conditional distributions of the state given the data assimilated up to iteration $t$ are approximated by a sequence of Gaussians $\{\mathcal{N}(m_t,C_t)\}_{t\in\mathbb{N}}$ , which converges exponentially to the exact posterior mean and covariance if the forward map is linear [32 Theorem 1], and provides a Gaussian approximation of $\mu^y$ otherwise, as in our setting. Only the first two moments of these distributions enter the iteration, and they are estimated empirically from an ensemble of $N$ so-called particles ( $N \ge2$ ). The method is therefore derivative-free and requires, at each iteration, $N$ independent evaluations of the forward map, which can be run in parallel. At the last iteration $T$ of the artificial dynamics, the empirical mean and covariance of the ensemble provide the estimates $m_T$ and $C_T$ of the sought Gaussian approximation. Usually, $T$ is specified by looking at convergence metrics, e.g. data misfit or $\Phi_{\mathrm{post}}$ .

For the formulation of the method, one considers the augmented operator

\mathcal{F} : \mathcal{X} \to\mathbb{R}^{n_y+n_\theta}, \qquad\mathcal{F}(\theta) = \begin{pmatrix} \mathcal{G}(\theta) \\ \theta\end{pmatrix},

and the artificial augmented discrete dynamics

\begin{aligned} \theta_{t+1} &= \theta_t + \varepsilon_{t+1},\\ z_{t+1} &= \mathcal{F}(\theta_{t+1}) + \nu_{t+1}. \tag*{(8)} \end{aligned}

for $t=0,1,\ldots,T$ , with an initial Gaussian guess on the parameter $\theta_0 \sim\mathcal{N}(m_0,C_0)$ . The second block of the augmented operator $\mathcal{F}$ consists of the parameter itself, so that each iteration assimilates not only the datum $y$ , but also the prior mean $m_0$ . This definition turns the potential $\Phi_{\mathrm{like}}$ (6) into the regularized potential $\Phi_{\mathrm{post}}$ (7) and prevents the repeated assimilation of the same datum from destroying the prior information. The noise terms in (8) are sampled from $\nu_{t+1} \sim\mathcal{N}(0,\Sigma_{\nu,t+1})$ and $\varepsilon_{t+1} \sim\mathcal{N}(0,\Sigma_{\varepsilon,t+1})$ , respectively, where

\Sigma_{\nu,t+1} := \frac{1}{\Delta\tau} \begin{pmatrix} \Sigma_\eta& 0\\ 0 & C_0 \end{pmatrix} \tag*{(9)}

depends on the measurement noise and on the prior on the parameter, and

\Sigma_{\varepsilon,t+1} := \frac{\Delta\tau}{1-\Delta\tau} C_t, \tag*{(10)}

depends on the current estimate of the parameter covariance $C_t$ . The parameter $\Delta\tau\in(0,1)$ in (9) and (10) acts as a time step for the augmented dynamics (8): it controls the fraction of the datum and of the prior that is assimilated at each iteration, see also remark 2.2.

Let us introduce the augmented observation $z$ defined as

z_t^\dagger:= \begin{pmatrix} y\\m_0 \end{pmatrix}, \quad\forall t=0,1,\ldots

and the observation set at time $t$

Y_t := \{z_0^\dagger,\ldots,z_t^\dagger\}.

The method is based on the assumption that the random variables $(\theta,z_{t+1})$ conditioned on $Y_t$ are distributed, for each iteration $t$ , according to

\begin{aligned} \theta_t \mid Y_t &\sim\rho(\theta_t \mid Y_t) \approx\mathcal{N}(m_t,C_t),\\ z_{t+1} \mid\theta_{t+1},Y_t &\sim\rho(z_{t+1}\mid\theta_{t+1},Y_t) \propto\rho(\theta_t\mid Y_t)e^{-\Delta\tau\Phi_{\mathrm{post}}(\theta_{t+1},y)}. \tag*{(11)} \end{aligned}

where $\Phi_{\mathrm{post}}(\theta_{t+1},y)$ is the regularized potential defined in (7).

Each iteration is composed of a prediction and an analysis step. In the prediction step, the mean is left unchanged, $\widehat{m}_{t+1}=m_t$ , while the covariance is inflated to $\widehat{C}_{t+1}=C_t/(1-\Delta\tau)$ , as prescribed by the artificial noise $\varepsilon_{t+1}$ in (8). In the analysis step, the current joint distribution $\rho(\theta_{t+1},z_{t+1}\mid Y_t)$ , which is in general not Gaussian, is projected onto the Gaussian distribution that matches its mean and covariance, and a standard Kalman update conditions the latter on the augmented observation $z_t^\dagger$ . In this way, the conditional distribution $\theta_{t+1}\mid Y_{t+1}$ is approximated by $\mathcal{N}(m_{t+1},C_{t+1})$ .

Different versions of the method have been proposed, depending on how the ensemble is updated to reproduce the analysis mean and covariance. In the stochastic EKI, each particle is updated using a perturbed observation, and the empirical covariance of the updated ensemble matches $C_{t+1}$ only in a statistical sense, see [32 Remark 5]. The EAKI employed here belongs instead to the family of square-root filters. Namely, the particles $\{\theta_{t+1}^j\}_{j=1}^{N}$ are updated deterministically by applying a linear operator $A:\mathbb{R}^{n_\theta}\to\mathbb{R}^{n_\theta}$ to the ensemble anomalies, i.e., at each iteration $t$ ,

\theta_{t+1}^j = m_{t+1} + A(\widehat{\theta}_{t+1}^j-\widehat{m}_{t+1}), \quad j=1,\ldots,N.

The operator $A$ is assembled from two singular value decompositions of the ensemble anomalies and guarantees that the empirical covariance of the updated ensemble matches $C_{t+1}$ exactly. Moreover, to prevent the collapse of the particles, we apply a jittering step at the end of each iteration, in which a perturbation drawn from a centered Gaussian distribution is added to each component of the particles whose standard deviation is below a given threshold. The EAKI method implemented in this work is described in detail in algorithm 1.

Remark 2.2 (Relation with the Fisher–Rao gradient flow). The splitting into a prediction step and an analysis step of the ensemble Kalman inversion method approximates a discrete-time fractional step of the Fisher–Rao gradient flow [42] minimizing the Kullback–Leibler divergence $\operatorname{KL}(\mu^y\Vert\mu^y)$ , where $\mu^y$ is the posterior distribution defined in (4) and $\mu_t^y$ is the one associated with $\theta_t|Y_t$ . The fractional step method is approximated as follows:

\hat{\rho}(\theta_{t+1}|Y_t) \propto\rho(\theta_t|Y_t)^{1-\Delta\tau}, \qquad\text{(prediction step)}

\rho(\theta_{t+1}|Y_{t+1}) \propto\hat{\rho}(\theta_{t+1}|Y_t)e^{-\Delta\tau\Phi_{\mathrm{post}}(\theta_{t+1},y)}, \qquad\text{(analysis step)}

where $\Delta\tau$ is the time step of the fractional step method.

Algorithm 1 Ensemble adjusted Kalman inversion (EAKI)

Input: $T_{\max}$ (EAKI iterations), $N$ (ensemble size), $\Delta\tau$ (time step), $\sigma_{\mathrm{thr}}$ (jitter threshold), $m_0$ , $C_0$ (prior mean and covariance)

Output: $m_{T_{\max}}$ , $C_{T_{\max}}$ (mean and covariance of the approximate posterior $\theta_{T_{\max}}\sim\mathcal{N}(m_{T_{\max}},C_{T_{\max}})$ ).

Initialization: sample $\theta_0^j\sim\mathcal{N}(m_0,C_0)=\mathcal{N}(0,I_d)$ , $j=1,\ldots,N$
for $t=0,1,\ldots,T_{\max}-1$ do Prediction step
$\hat{m}_{t+1}=m_t$
$\hat{\theta}_{t+1}^j=m_t+\frac{1}{\sqrt{1-\Delta\tau}}(\theta_t^j-m_t)$ , $j=1,\ldots,N$
$\hat{z}_{t+1}^j=F_k(\hat{\theta}_{t+1}^j)$ , $j=1,\ldots,N$ ▷ parallel forward map
Ensemble generation
$\hat{z}_{t+1}=\frac{1}{N}\sum_{j=1}^{N}\hat{z}_{t+1}^j$
$\hat{Z}_{t+1}=\frac{1}{\sqrt{N-1}}\left(\hat{\theta}_{t+1}^1-\hat{m}_{t+1}\ \hat{\theta}_{t+1}^2-\hat{m}_{t+1}\ \cdots\ \hat{\theta}_{t+1}^N-\hat{m}_{t+1}\right)$
$\hat{Y}_{t+1}=\frac{1}{\sqrt{N-1}}\left(\hat{z}_{t+1}^1-\hat{z}_{t+1}\ \hat{z}_{t+1}^2-\hat{z}_{t+1}\ \cdots\ \hat{z}_{t+1}^N-\hat{z}_{t+1}\right)$
apply SVD to evaluate $P$ , $\hat{D}$ , $V$ : $\hat{Z}_{t+1}=P\hat{D}^{1/2}V^T$
apply SVD to evaluate $U$ , $D$ : $V^T(I_d+(\hat{Y}_{t+1})^T(\Sigma_{\nu,t+1})^{-1}\hat{Y}_{t+1})^{-1}V=UDU^T$
$A=P\hat{D}^{1/2}UD^{1/2}\hat{D}^{-1/2}P^T$
Analysis step
$\hat{C}_{t+1}^{\theta,z}=\hat{Z}_{t+1}(\hat{Y}_{t+1})^T$
$\hat{C}_{t+1}^{z,z}=\hat{Y}_{t+1}(\hat{Y}_{t+1})^T+\Sigma_{\nu,t+1}$
mean: $m_{t+1}=\hat{m}_{t+1}+\hat{C}_{t+1}^{\theta,z}(\hat{C}_{t+1}^{z,z})^{-1}(\bar{z}_0-\hat{z}_{t+1})$
variance: $C_{t+1}=\hat{C}_{t+1}-\hat{C}_{t+1}^{\theta,z}(\hat{C}_{t+1}^{z,z})^{-1}(\hat{C}_{t+1}^{\theta,z})^T$
$\theta_{t+1}^j=m_{t+1}+A(\hat{\theta}_{t+1}^j-\hat{m}_{t+1})$
Apply jitter
for $j=1,\ldots,N$ do
for $l=1,\ldots,n_\theta$ (components of parameter vector) do
if $\operatorname{std}((\theta_{t+1})_l)<\sigma_{\mathrm{thr}}$ then
$(\theta_{t+1}^j)_l=(\theta_{t+1}^j)_l+\xi,\quad\xi\sim\mathcal{N}(0,1)$
end if
end for
end for
Convergence check
end for

Adaptive Bayesian Updating with Structural reliability (aBUS)

Bayesian Updating with Structural reliability (BUS) methods recast the inference of the posterior into a rare-event estimation problem, as studied in the context of structural reliability [33, 35]. In structural reliability analysis, the state of a system is described by a limit-state function, which is negative on the so-called failure domain, and the quantity of interest is the failure probability $\mathbb{P}(h\leq0)$ computed with respect to a standard normal reference density.

In this context, the difficulty is that failure is typically a rare event. BUS builds on the observation that the posterior can be characterized in the very same form: after augmenting the parameter with one auxiliary variable, the acceptance criterion of rejection sampling defines a domain whose conditional distribution is exactly $\mu^y$ . Unlike EAKI, no Gaussianity assumption is made on the posterior.

To describe the method, let us introduce a bijective map $T:\mathbb{R}^{n_\theta}\to\mathbb{R}^{n_\theta}$ that maps the standard normal variables $u_{1:n_\theta}$ back into the physical parameter space, i.e., such that $\theta=T(u_{1:n_\theta})$ for $u_{1:n_\theta}\sim\mathcal{N}(0,\mathrm{Id})$ . BUS considers an augmented parameter vector with an independent standard-normal variable $u_{n_\theta+1}\sim\mathcal{N}(0,1)$ and encodes the observation through the so-called limit-state function [33]:

\begin{aligned} h(u,\ell)&:=\log\Phi_{\mathrm{CDF}}(u_{n_\theta+1})-\log\left(e^{-\ell}\mathcal{L}(y\mid T(u_{1:n_\theta}))\right)\\ &=\log\Phi_{\mathrm{CDF}}(u_{n_\theta+1})+\ell-\log\mathcal{L}(y\mid T(u_{1:n_\theta})). \tag*{(12)} \end{aligned}

where $\mathcal{L}(y\mid\theta)$ is the likelihood function induced by the forward map (2)

\mathcal{L}(y\mid\theta):=\exp\left(-\Phi_{\mathrm{like}}(\theta,y)\right), \tag*{(13)}

$\ell$ is a scaling constant such that $e^{-\ell}\leq1/\sup_\theta\mathcal{L}(y\mid\theta)$ , and $\Phi_{\mathrm{CDF}}$ is the standard-normal CDF. Based on (12), we define the failure domain as the event

\mathcal{Y}_u:=\{u\in\mathbb{R}^{n_\theta+1}\mid h(u,\ell)\leq0\}. \tag*{(14)}

The corresponding conditional probability is given by

\mathbb{P}(u\in\mathcal{Y}_u\mid u_{1:n_\theta})=\mathbb{P}\left(\Phi_{\mathrm{CDF}}(u_{n_\theta+1})\leq e^{-\ell}\mathcal{L}(y\mid T(u_{1:n_\theta}))\right). \tag*{(15)}

Since $u_{n_\theta+1}\sim\mathcal{N}(0,1)$ , we have $\Phi_{\mathrm{CDF}}(u_{n_\theta+1})\sim\mathcal{U}(0,1)$ (standard uniform). As a consequence,

\mathbb{P}(u\in\mathcal{Y}_u\mid u_{1:n_\theta})=e^{-\ell}\mathcal{L}(y\mid T(u_{1:n_\theta})), \tag*{(16)}

which is equivalent to a standard rejection sampling criterion: a candidate $\theta=T(u_{1:n_\theta})$ drawn from the prior is accepted whenever the augmented sample $u$ falls in $\mathcal{Y}_u$ , that is, with probability $e^{-\ell}\mathcal{L}(y\mid\theta)$ . The condition on $\ell$ ensures that this quantity does not exceed one.

Then, we introduce the augmented density given by

\rho_{\mathcal{Y}_u}(u)\propto\mathbf{1}_{\mathcal{Y}_u}(u)\phi_{n_\theta+1}(u), \tag*{(17)}

where $\phi_{n_\theta+1}$ is the $(n_\theta+1)$ -dimensional standard normal density. Marginalizing the augmented density with respect to $u_{n_\theta+1}$ gives

\begin{aligned} \rho_{\mathcal{Y}_u}^{1:n_\theta}(u_{1:n_\theta})&\propto\phi_{n_\theta}(u_{1:n_\theta})\int_{\mathbb{R}}\mathbf{1}_{\mathcal{Y}_u}(u)\phi_1(u_{n_\theta+1})\,\mathrm{d}u_{n_\theta+1}\\ &\propto\phi_{n_\theta}(u_{1:n_\theta})\mathbb{P}(u\in\mathcal{Y}_u\mid u_{1:n_\theta})\\ &\propto\phi_{n_\theta}(u_{1:n_\theta})\mathcal{L}(y\mid T(u_{1:n_\theta})). \tag*{(18)} \end{aligned}

The marginal density thus corresponds to the sought density of $\mu^y$ under a change of coordinates. Hence, to obtain samples of $\mu^y$ , the first $n_\theta$ components of the samples drawn from the augmented density are pushed through the map $T$ . This formulation in the augmented space is particularly convenient, since the resulting density is expressed with respect to a standard normal reference density, while the prior information is preserved through the map $T$ .

Sampling $\mu^y$ by plain rejection requires, on average, one accepted sample every $\left(e^{-\ell}\int_{\mathcal{X}}\mathcal{L}(y\mid\theta)\,\mathrm{d}\mu_0\right)^{-1}$ evaluations of the forward map. For informative data, the likelihood is concentrated with respect to the prior, the domain $\mathcal{Y}_u$ becomes a rare event under the standard normal reference, and this cost, like that of standard MCMC methods in high-dimensional spaces, becomes prohibitive. We therefore consider the combination of BUS with Subset Simulation (SuS) [34] proposed in [33], in which the rare event is reached through a sequence of nested and progressively less probable events.

Namely, we consider a series of intermediate threshold levels

+\infty=:h_0>h_1>\cdots>h_{n_L}:=0,

and define the nested events

F_i:=\{u\mid h(u,\ell)\leq h_i\},\quad i=0,1,\ldots,n_L,

as well as the densities

\rho_i(u)\propto\mathbf{1}_{F_i}(u)\phi_{n_\theta+1}(u),\quad i=0,1,\ldots,n_L. \tag*{(19)}

By definition, the level sets $F_i$ form a nested chain, i.e.,

\mathbb{R}^{n_\theta+1}=F_0\supseteq F_1\supseteq\cdots\supseteq F_{n_L}=\mathcal{Y}_u.

Hence, the acceptance probability telescopes into the product of the level-conditional probabilities, i.e.,

\mathbb{P}(h\leq0)=\mathbb{P}(F_{n_L})=\prod_{k=1}^{n_L}\mathbb{P}(F_k\mid F_{k-1})=\prod_{k=1}^{n_L}P_k. \tag*{(20)}

In practice, at each level $i=0,1,\ldots,n_L-1$ , we adaptively select $h_{i+1}$ as the $p_0$ -quantile of the ordered samples at level $i$ , so that a fraction $p_0$ of them satisfies the event $F_{i+1}$ . These samples are retained as seeds

Algorithm 2 Parallel adaptive BUS with subset simulation (aBUS–SuS). The forward map is detailed in § 3 and algorithm 6.

Input: samples per level $N$ , conditional probability $p_0 = P(F_i \mid F_{i-1})$ , with $i = 1,\ldots,n_L$ , prior transform $T(\tilde{\theta}) = m_c + L_c\tilde{\theta}$ , log-likelihood $\log\mathcal{L}(y \mid\theta) = -\Phi_{\mathrm{like}}(\theta,y)$

Output: thresholds $\{h_i\}$ , posterior samples $\{\theta_j\}$

$u = (u_{1:n_\theta},u_{n_\theta+1}) \in\mathbb{R}^{n_\theta+1}$ , with $u_{n_\theta+1} \sim\mathcal{N}(0,1)$ ▷ augment the standard space
$h(u,\hat{\ell}) = \log\Phi_{\mathrm{CDF}}(u_{n_\theta+1}) + \hat{\ell} - \log\mathcal{L}(y \mid T(u_{1:n_\theta}))$ ▷ limit-state function (12)

Initial Monte Carlo level

draw $u_j \sim\mathcal{N}(0,\mathbf{I}_d)$ , $j = 1,\ldots,N$
$\{\log\mathcal{L}_j\} = \operatorname{ForwardMap}(\{u_j\})$ ; ▷ forward map, in parallel
$\hat{\ell} = \max_j \log\mathcal{L}_j$ ; $i = 0$ ; $h_0 = +\infty$

Subset simulation

while $h_i > 0$ do
$i = i + 1$ ; $g_j = h(u_j,\hat{\ell})$ for all $j = 1,\ldots,N$
$h_i = p_0$ -quantile of $\{g_j\}$ ; $n_F = \#\{j:g_j \le\max(h_i,0)\} \approx p_0 \cdot N$
if $h_i < 0$ then $h_i = 0$ , $P_i = n_F/N$ else $P_i = p_0$
$\{\omega_k\}_{k=1}^{n_F}$ = the $n_F$ samples with smallest $g_j$ ▷ set up the seeds
$\{u_j\},\{\log\mathcal{L}_j\},\lambda= \operatorname{aCS}(N,\lambda,h_i,\{\omega_k\}_{k=1}^{n_F},\hat{\ell})$ ▷ adaptive MCMC, algorithm 3
$\ell' = \max(\hat{\ell},\max_j \log\mathcal{L}_j)$ ; $h_i = h_i - \hat{\ell} + \ell'$ ; $\hat{\ell} = \ell'$ ▷ update threshold
$u_{j,n_\theta+1} = \Phi_{\mathrm{CDF}}^{-1}(p_j)$ , $p_j \sim\mathcal{U}(0,\min(1,e^{\log\mathcal{L}_j-\hat{\ell}+h_i}))$ ▷ refresh BUS variable
end while
$n_L = i$ ; $\theta_j = T(u_{j,1:n_\theta})$

and used to initialize multiple MCMC chains with target distribution $\rho_{i+1}$ . By construction, the seeds already belong to $F_{i+1}$ and are distributed according to $\rho_{i+1}$ : the chains therefore start in stationarity and no burn-in truncation period is required.

The efficiency of the MCMC sampling process depends on the proposal distribution. In the context of SuS, some specialized methods have been proposed (see, e.g., [43, 44, 45]), in particular, we use adaptive conditional sampling (aCS) as introduced in [46].

The overall efficiency of BUS in conjunction with SuS depends on the choice of the proportionality constant $e^{-\ell}$ : the smaller it is, the rarer the event $\mathcal{Y}_u$ and the larger the number of levels, while its admissibility condition $e^{-\ell} \le1/\sup_\theta\mathcal{L}(y \mid\theta)$ involves the maximum of the likelihood, which is not known a priori. The adaptive BUS (aBUS) method addresses this issue by estimating the constant at runtime [35]: $\ell$ is replaced by the largest log-likelihood value $\hat{\ell}$ observed so far and, whenever a larger value is found, the current threshold is shifted accordingly and the auxiliary variable $u_{n_\theta+1}$ is resampled from its conditional distribution (last two steps of the loop in algorithm 2), so that the samples remain distributed according to the current intermediate density.

Each evaluation of the limit-state function (12) requires one evaluation of the forward map $\mathcal{G}$ , which includes the numerical solution of the Navier–Stokes equations on the different domains (§ 3). Since the candidates of a level, or of a batch of MCMC steps, are mutually independent, they can be evaluated concurrently: following [36], each batch of forward evaluations $\theta\mapsto\mathcal{G}(\theta)$ is dispatched at once, so that the $N$ solutions of the PDE underlying the forward map are run in parallel. The implemented aBUS algorithm is summarized in algorithm 2.

Wasserstein–Fisher–Rao sequential Monte Carlo (WFR)

2.4. Wasserstein–Fisher–Rao sequential Monte Carlo (WFR). As pointed out in remark 2.2, EAKI approximates a fractional step scheme of the Fisher–Rao gradient flow of the Kullback–Leibler divergence $\operatorname{KL}(\mu_t^y\Vert\mu^y)$ . In this section, we introduce the Wasserstein–Fisher–Rao sequential Monte Carlo (WFR) method, which can be seen as a further extension [38, 47] of EAKI or of the Fisher–Rao gradient flow (see remark 2.2) that approximates the full Wasserstein–Fisher–Rao gradient flow of $\operatorname{KL}(\rho\Vert\mu^y)$ , for a probability measure $\rho$ .

The idea behind the method is to evolve an empirical measure $\rho_t$ towards the posterior $\mu^y$ of (4) along the gradient flow of the Kullback–Leibler divergence $\operatorname{KL}(\rho\Vert\mu^y)$ under the Wasserstein–Fisher–Rao metric [38, 47], combining the Wasserstein (optimal-transport) and Fisher–Rao (birth–death) geometries,

Algorithm 3 Adaptive conditional sampling (aCS) called by algorithm 2. The forward map is detailed in § 3 and algorithm 6.

Input: target acceptance value ( $\hat{a}_{\mathrm{target}} = 0.44$ ), initial proposal std ( $\sigma_{\mathrm{aCS},0} = 1$ ), batch size ( $n_B = 100$ ).

function $\mathrm{aCS}(N,\text{ proposal scaling parameter }\lambda,\tau,\{\omega_k\}_{k=1}^{n_F},\hat{\ell})$
distribute $N$ samples over $n_F$ chains; seed each chain head from $\{\omega_k\}_{k=1}^{n_F}$
$\sigma_{\mathrm{aCS}} = \min(\lambda\sigma_{\mathrm{aCS},0},1)$ ; $\varrho= \sqrt{1-\sigma_{\mathrm{aCS}}^2}$
split the $N$ chain steps into $n_A = \lceil N/n_B\rceil$ batches of (at most) $n_B$ proposals each
for $k=1,\ldots,n_A$ do
propose $v^\star\sim\mathcal{N}(\varrho u_{t-1},\sigma_{\mathrm{aCS}}^2)$ for every active chain ▷ conditional proposal
$\{\log\mathcal{L}^\star\}=\operatorname{ForwardMap}(\{v^\star\})$ ▷ forward map, in parallel
accept $u_t=v^\star$ if $h(v^\star,\hat{\ell})\le\tau$ , else $u_t=u_{t-1}$ (carry $\log\mathcal{L}$ )
evaluate batch acceptance rate $\hat{a}$
$\lambda=\exp\left(\log\lambda+\frac{1}{\sqrt{k}}(\hat{a}-\hat{a}_{\mathrm{target}})\right)$ ; update $\sigma_{\mathrm{aCS}},\varrho$
end for
return $\{u_j\},\{\log\mathcal{L}_j\},\lambda$
end function

i.e.,

\partial_t\rho_t=\nabla\cdot\left(\rho_t\nabla\log\left(\frac{\rho_t}{\mu^y}\right)\right)-\rho_t\left(\log\frac{\rho_t}{\mu^y}-\mathbb{E}_{\rho_t}\left[\log\frac{\rho_t}{\mu^y}\right]\right). \tag*{(21)}

The measure is approximated by an ensemble $\{\theta_j\}_{j=1}^N$ of $N$ particles, i.e., $\rho_t\approx\frac{1}{N}\sum_j\delta_{\theta_j}$ , and the evolution is realized via an operator splitting of the two half-flows over a step $\gamma_t$ [37]. Unlike EAKI and aBUS, WFR requires the evaluation of the likelihood score, i.e., the gradient of the log-likelihood with respect to the parameters, which in turn requires the computation of the Jacobian $\partial\mathcal{G}/\partial\theta$ of the forward map. To reduce the computational cost, we replace $\mathcal{G}$ with a differentiable neural network surrogate $\widetilde{\mathcal{G}}$ (§ 3.7) that approximates the PDE solution on the deforming domain.

Let

\widetilde{\Phi}_{\mathrm{like}}(\theta,y):=\frac{1}{2}\|y-\widetilde{\mathcal{G}}(\theta)\|_{\Sigma_\eta^{-1}}^2

denote the potential (6) evaluated through the differentiable surrogate $\widetilde{\mathcal{G}}$ . Writing the log-posterior, up to an additive constant, as $\log\mu^y(\theta)=-\widetilde{\Phi}_{\mathrm{like}}(\theta,y)+\log\pi_0(\theta)$ , and the prior $\pi_0=\mathcal{N}(0,\mathrm{Id})$ , and approximating the ensemble density $\rho_t$ by the Gaussian-mixture kernel estimate $q(\theta)=\frac{1}{N}\sum_i\mathcal{N}(\theta;\theta_i,2\gamma_t\mathrm{Id})$ , where $\mathcal{N}(\theta;\theta_i,2\gamma_t\mathrm{Id})$ denotes the density of the Gaussian distribution with mean $\theta_i$ and covariance $2\gamma_t\mathrm{Id}$ evaluated at $\theta$ , at one analysis step applies the following updates in sequence:

Fisher–Rao (birth–death). Each particle is scored by the log density-ratio

\ell_j=\log\mu^y(\theta_j)-\log q_j=-\widetilde{\Phi}_{\mathrm{like},j}+\log\pi_0(\theta_j)-\log q_j

between the posterior and the ensemble, positive where $\theta_j$ is under-represented relative to $\mu^y$ and negative where it is over-represented. The particles are reweighted and normalized,

w_j=\frac{\exp(\alpha_{\mathrm{FR}}\ell_j)}{\sum_{k=1}^{N}\exp(\alpha_{\mathrm{FR}}\ell_k)}, \tag*{(22)}

then resampled (systematic resampling), duplicating under-represented particles and pruning over-represented ones. The exponent $\alpha_{\mathrm{FR}}$ is the Fisher–Rao step size: its nominal value $1-e^{-\gamma_t}$ discretizes the birth–death half-flow over $\gamma_t$ , but to prevent weight degeneracy it is decreased to the largest $\alpha_{\mathrm{FR}}\le1-e^{-\gamma_t}$ for which the effective sample size

\operatorname{ESS}(\alpha_{\mathrm{FR}})=\left(\sum_{j=1}^{N}w_j^2\right)^{-1}\in[1,N], \tag*{(23)}

stays above $N/2$ ; ESS $=N$ for uniform weights and ESS $=1$ when a single particle dominates.

Wasserstein (Langevin). Each resampled particle is transported by one Langevin step

\theta_{t+1}=\theta_t+\gamma_t\left(-\nabla_\theta\widetilde{\Phi}_{\mathrm{like}}+\nabla_\theta\log\pi_0\right)+\sqrt{2\gamma_t}\,\xi_t,\qquad\xi_t\sim\mathcal{N}(0,\mathrm{Id}), \tag*{(24)}

with likelihood score

-\nabla_\theta\widetilde{\Phi}_{\mathrm{like}}(\theta)=\left(\partial\widetilde{\mathcal{G}}/\partial\theta\right)^\mathsf{T}\Sigma_\eta^{-1}\left(y-\widetilde{\mathcal{G}}(\theta)\right),

Algorithm 4 Parallel Wasserstein–Fisher–Rao (WFR) sequential Monte Carlo sampling. The surrogate forward map is defined in § 3.7, see also algorithm 7. The analysis step is detailed in algorithm 5.

Input: ensemble size $N$ , iterations $T_{\max}$ , noise scale $\sigma_{\eta}$ ( $\sigma_{\eta}^{2}=\Sigma_{\eta}$ for the scalar observable $O_{2}$ ), observed datum $y$ , parameter space maps $(L_c,m_c,U)$ of (25), differentiable surrogate $\widetilde{G}$ of $G$

Output: posterior ensemble $\{\widetilde{\theta}_{j}\}_{j=1}^{N}$

$\widetilde{\theta}_{j}\sim\mathcal{N}(0,\mathrm{Id})$ ▷ prior init $\pi_{0}=\mathcal{N}(0,\mathrm{Id})$
for $t=0,\ldots,T_{\max}-1$ do
$\gamma_{t}=\gamma_{\min}+\frac{1}{2}(\gamma_{\max}-\gamma_{\min})(1+\cos(\pi t/(T_{\max}-1)))$ ▷ cosine step scheduler
$\{\widetilde{G}(\widetilde{\theta}_{j})\},\{\nabla\widetilde{G}(\widetilde{\theta}_{j})\}=\operatorname{SurrogateForwardMap}(\{\widetilde{\theta}_{j}\})$ ▷ batched differentiable GPU job
$\{\widetilde{\theta}_{j}\}=\operatorname{WFRanalysis}(\{\widetilde{\theta}_{j}\},\{\widetilde{G}(\widetilde{\theta}_{j})\},\{\nabla\widetilde{G}(\widetilde{\theta}_{j})\},y,\sigma_{\eta},\gamma_{t})$
end for
return $\{\widetilde{\theta}_{j}\}$

Algorithm 5 One WFR sequential MC step called by algorithm 4.

function $\operatorname{WFRanalysis}(\{\widetilde{\theta}_{j}\},\{\widetilde{G}(\widetilde{\theta}_{j})\},\{\nabla\widetilde{G}(\widetilde{\theta}_{j})\},y,\sigma_{\eta},\gamma)$

(1) Likelihood potential and its score

$r_{j}=\widetilde{G}(\widetilde{\theta}_{j})-y$ ; $\widetilde{\Phi}_{\mathrm{like},j}=\frac{1}{2}(r_{j}/\sigma_{\eta})^{2}$ ▷ Gaussian likelihood potential (6)
$s_{j}=-\nabla_{\widetilde{\theta}}\widetilde{\Phi}_{\mathrm{like},j}=-\frac{r_{j}}{\sigma_{\eta}^{2}}\nabla\widetilde{G}(\widetilde{\theta}_{j})$ ▷ likelihood score

(2) Fisher–Rao flow: birth–death by reweighting and resampling

the FR velocity is $\log(\mu^{y}/\rho)$ , with target $\mu^{y}\propto\pi_{0}e^{-\widetilde{\Phi}_{\mathrm{like}}}$ and $\rho$ the current particle density

$\log q_{j}=-\log N+\operatorname{logsumexp}_{1\leq i\leq N}\left[-\frac{n_{\theta}}{2}\log(4\pi\gamma)-\frac{1}{4\gamma}\lVert\widetilde{\theta}_{j}-\widetilde{\theta}_{i}\rVert^{2}\right]$ ▷ $q\approx\rho$

a Gaussian-mixture density estimate over the ensemble, kernel $\mathcal{N}(\cdot;\widetilde{\theta}_{i},2\gamma\mathrm{Id})$

$\log\pi_{0}(\widetilde{\theta}_{j})=-\frac{1}{2}\lVert\widetilde{\theta}_{j}\rVert^{2}$ ▷ $\widetilde{\theta}\sim\mathcal{N}(0,\mathrm{Id})$ prior, up to an additive constant
$\ell_{j}=-\widetilde{\Phi}_{\mathrm{like},j}+\log\pi_{0}(\widetilde{\theta}_{j})-\log q_{j}$ ▷ $\log\mu^{y}-\log q$ ; constants cancel in $w$ below
$\overline{w}_{j}(\alpha_{\mathrm{FR}})=\frac{\exp(\alpha_{\mathrm{FR}}\ell_{j})}{\sum_{i}\exp(\alpha_{\mathrm{FR}}\ell_{i})}$ , $\operatorname{ESS}(\alpha_{\mathrm{FR}})=\left(\sum_{j}\overline{w}_{j}(\alpha_{\mathrm{FR}})^{2}\right)^{-1}\in[1,N]$
if $\operatorname{ESS}(1)\geq N/2$ then
$\alpha_{\mathrm{FR}}=1$
else
$\alpha_{\mathrm{FR}}=$ largest $\alpha_{\mathrm{FR}}\in[0,1]$ with $\operatorname{ESS}(\alpha_{\mathrm{FR}})\geq N/2$ , by 40 bisection steps ▷ ESS decreases in $\alpha_{\mathrm{FR}}$
end if
$w_{j}=\overline{w}_{j}(\alpha_{\mathrm{FR}})$ ▷ $w_{j}\propto(\mu^{y}/q)^{\alpha_{\mathrm{FR}}}$ at $\widetilde{\theta}_{j}$
draw $u\sim U(0,1)$ ; $\operatorname{idx}_{j}=\min\{i:\sum_{l\leq i}w_{l}\geq(j+u)/N\}$ , $j=0,\ldots,N-1$ ▷ systematic resampling
$\widetilde{\theta}_{j}=\widetilde{\theta}_{\operatorname{idx}_{j}}$ , $s_{j}=s_{\operatorname{idx}_{j}}$ ▷ scores travel with their particles

(3) Wasserstein flow: one Langevin step on the resampled particles

$v_{j}=s_{j}-\widetilde{\theta}_{j}$ ▷ $\nabla_{\widetilde{\theta}}\log\mu^{y}$ : likelihood score + prior score $\nabla_{\widetilde{\theta}}\log\pi_{0}=-\widetilde{\theta}$
$\xi_{j}\sim\mathcal{N}(0,\mathrm{Id})$ ; $\widetilde{\theta}_{j}=\widetilde{\theta}_{j}+\gamma v_{j}+\sqrt{2\gamma}\,\xi_{j}$ ▷ Euler–Maruyama, unadjusted
return $\{\widetilde{\theta}_{j}\}_{j=1}^{N}$
end function

evaluated from the Jacobian $\partial\widetilde{G}/\partial\theta$ of the surrogate by automatic differentiation (see remark 3.5), and prior score $\nabla_{\theta}\log\pi_{0}=-\theta$ , so that the drift equals $\nabla_{\theta}\log\mu^{y}$ and the invariant law of the continuous dynamics is the posterior $\mu^{y}$ .

Unlike the ensembles of § 2.2–§ 2.3, the approximate forward map $\widetilde{G}$ and its Jacobian $\partial\widetilde{G}/\partial\theta$ are here evaluated in a single batched, differentiable forward pass of the PDE surrogate introduced in § 3.7. Consequently, the flow does not require a separate high-fidelity PDE solution for each ensemble member. An additional advantage of a neural network surrogate is the use of batched evaluations of observables and their gradients, which is more efficient than the sequential evaluation of the high-fidelity PDE solver. Algorithms 4 and 5 sketch the steps of the WFR method.

Forward problem

This work is motivated by the problem of finding a posterior distribution compatible with a measured hemodynamic biomarker in aortic blood flows, such as pressure drops. In the formalism introduced in section 2, the quantities of interest are represented by the map $\mathcal{G}(\theta)$ , while the parameter vector $\theta\in\mathcal{X}$ defines a particular aortic shape. The map is built on two recent works concerning diffeomorphic registration of aortic shapes [27] and generative modeling of shapes based on stochastic interpolants [6]. In this section, we first introduce the geometry dataset used for training the registration algorithm and the generative model (§ 3.1), as well as for defining a latent parameter space encoding the main geometrical features (§ 3.2). In § 3.3 we detail the different steps of the forward map, entailing a generative model for the construction of a surface mesh, a smoothing step for the registration of a reference hexahedral mesh, the numerical solution of the Navier–Stokes equations, and the evaluation of an observation functional (e.g., a pressure drop between selected locations). The WFR method of § 2.4 requires a differentiable surrogate map for the evaluation of the likelihood score. The steps to define this surrogate and the differences with the full forward map are discussed in § 3.7. A sketch of the different steps defining the forward map is provided in Figure 2.

Figure 2. Schematic representation of the forward map, from the latent representation of the shape to the observation, using either a high-fidelity numerical solver (§ 3.3) or a surrogate model (§ 3.7).

Geometry dataset and shape registration

The building blocks of the forward map rely on a dataset of synthetic aortic shapes created using statistical shape modeling (SSM) from an initial cohort of 228 subjects (106 aortic coarctation patients, 37 post-treatment patients, and 85 healthy individuals). The procedure resulted in a dataset of 1257 aortic shapes, described via a surface triangular mesh, a set of centerline control points (same amount for all shapes), and for each centerline point, the radius of the corresponding inscribed sphere. For the detailed methodology, we refer the reader, e.g., to [48, 49, 50, 51, 52].

All shapes describe the same region of interest, comprising the vessel surface of the aortic arch up to the thoracic aorta (TA), including three main branches and the corresponding boundary surfaces (brachiocephalic artery, BCA, left common carotid artery, LCCA, left subclavian artery, LSA) (Figure 3, left). The dataset of 1257 geometries is then split into 1205 training and 52 test shapes, as in [6]. One training shape is selected as template and indexed with $i = 0$ . A sample of the test dataset is shown in Figure 4.

Considered region of the aorta and dimension reduction for the parameter space — Figure 3. Left. Considered region of the aorta, highlighting the four outlet boundaries (BCA: brachiocephalic artery, LCCA: left common carotid artery, LSA: left subclavian artery, TA: thoracic aorta). Right. Dimension reduction for the definition of the parameter space, encoding 94 4-dimensional centerline coordinates into a latent variable ( $\tilde{\theta}$ ) of dimension 250.

Figure 4. Selection of 30 of the 52 test geometries.

Based on the SSM construction, each synthetic geometry is parametrized via a centerline encoding $c \in\mathbb{R}^{d_c}$ , with $d_c := 94 \cdot4 = 376$ , gathering the three coordinates and the radii of the corresponding inscribed spheres of 94 selected control points along the centerline (main vessel and branches, Figure 3, right). For each shape in the database, we compute a bi-Lipschitz homeomorphic registration map on the template surface using the multi-level ResNet-LDDMM algorithm proposed in [27]. The template mesh is registered onto the point cloud describing other shapes, allowing us to define push-forward and pull-back maps between the different computational domains and finite element solutions. This step is crucial to define the generative model introduced in § 3.3.1.

Remark 3.1. The formulation of the method does not depend on the particular choice of the template. However, since the template shape defines the reference space onto which the other shapes and numerical solutions are mapped, one should avoid selecting shapes with extremely marked features. For a more detailed discussion on the choice of the template we refer to [27].

Parameter space

The parameter space $\mathcal{X}$ is defined in two steps. First, we reduce the dimension of the centerline encoding using Principal Component Analysis (PCA) on the training set and retaining the first $n_c = 250$ modes. Next, to further increase the variability of the considered population, we include an additional parameter describing the degree of a stenosis. The stenosis is included via a one-dimensional parametric deformation, gradually reducing the radius of the inscribed sphere around a fixed location along the aorta centerline (see Figure 5).

Figure 5. Parametric stenosis deformation for increasing severity levels $\alpha\in\{0\%,25\%,50\%,75\% \}$ (from left to right).

The resulting parameter space $\mathcal{X} = \mathbb{R}^{n_\theta}$ has dimension $n_\theta= n_c + 1$ . For a parameter $\theta= (\tilde{\theta},\tilde{\alpha}) \in\mathcal{X}$ , the component $\tilde{\theta} \in\mathbb{R}^{n_c}$ will be referred to as the latent coordinate, whilst $\tilde{\alpha} \in[0,1]$ refers to the stenosis degree.

Formally, let us denote by $\{c_i \in\mathbb{R}^{d_c}\}_{i=1}^{n_{\mathrm{train}}}$ the centerline encoding of the training set, and with $U \in\mathbb{R}^{d_c \times n_c}$ the resulting PCA basis. Moreover, we denote by $m_c \in\mathbb{R}^{n_c}$ and $L_c \in\mathbb{R}^{n_c \times n_c}$ the mean and the square root of the covariance of the reduced representations of the training distribution $\{U^\mathsf{T}c_i\}_{i=1}^{n_{\mathrm{train}}} \subset\mathbb{R}^{n_c}$ , respectively. The prior distribution for the parameter $\theta= (\tilde{\theta},\tilde{\alpha})$ is then defined as a standard normal distribution on the space $\mathcal{X}$ , after centering and scaling the latent coordinate with $m_c \in\mathbb{R}^{n_c}$ and $L_c \in\mathbb{R}^{n_c \times n_c}$ .

Remark 3.2. The UQ methods introduced in section 2 require sampling the parameter space. The dimension reduction, from 376 to 250 coordinates, has been implemented with the purpose of retaining the most relevant geometrical features and limiting the dimension of the parameter space. To this purpose, the truncation at 250 was chosen by computing the PCA reconstruction error for a varying number of modes 50, 100, $\ldots$ , 350 and selecting the value at which the decay of the singular values shows the largest spectral gap.

Remark 3.3. The location and the form of the additional stenosis are arbitrary. The purpose of including the stenosis in the parameter space is to demonstrate that any additional parametric deformation can be combined with the LDDMM stochastic interpolant used as generative prior for the surface.

Parameter-to-observable map

Generative model for surface mesh

The first step concerns the generation of a surface geometry. The approach is based on the LDDMM stochastic interpolant framework introduced in [6] in the context of aortic blood flow simulations.

Given $\theta= (\tilde{\theta},\tilde{\alpha})$ , we first recover the corresponding centerline encoding $c(\theta) \in\mathbb{R}^{d_c}$ via

c(\theta) = U\left(L_c\tilde{\theta} + m_c\right). \tag*{(25)}

The vector $c(\theta)$ is used as a conditioning variable to define the LDDMM stochastic interpolants [6], which are obtained by solving an SDE of the form

\begin{aligned} \mathrm{d}I_t &= b_t^{w}(I_t,c(\theta))\,\mathrm{d}t + \sigma_t\,\mathrm{d}W_t,\\ I_{t=0} &= X_0. \tag*{(26)} \end{aligned}

In (26), $X_0$ denotes the surface of the template shape and $b_t^w(I_t,c(\theta))$ is a drift coefficient trained using surface registration maps between the shapes of the training dataset and the template [6]. Moreover, $\sigma_t=\sqrt{t(1-t)}$ and $W_t$ is a standard Brownian motion. Equation (26) is solved using the forward Euler method with 10 time steps up to $t=1$ and the resulting surface displacement is given by $\Psi_s=\mathbb{E}[I_{t=1}]$ .

Next, we apply the deformation related to the stenosis parameter. The stenosis degree is obtained from the variable $\tilde{\alpha}\in\mathbb{R}$ through the standard-normal CDF, rescaled to a physical range, i.e.,

\alpha=\alpha_{\min}+(\alpha_{\max}-\alpha_{\min})\Phi_{\mathrm{CDF}}(\tilde{\alpha}).

Based on the value of $\alpha$ , the radii of the inscribed spheres are reduced gradually, up to a factor $\alpha$ at the selected location (see Figure 5), directly on the generated surface mesh.

Mesh

The next step handles the extension of the surface deformation to a suitable hexahedral mesh of the bulk domain. In particular, since the mesh is required for the numerical solution of the Navier–Stokes equations, it is necessary to guarantee an acceptable mesh regularity and avoid element degradation, especially in the highly deformed regions. The mesh generation process consists of two main steps.

First, the surface is preprocessed, applying a Taubin smoothing with the graph Laplacian and correcting the inlet and outlet faces by projecting the mesh nodes on planar surfaces.

The hexahedral mesh in the bulk domain is defined by mapping the template volume mesh onto the new surface. We solve a multi-step linear elasticity problem for an internal displacement field, starting from the template mesh and gradually deforming it towards the generated surface $\Psi_s$ :

\begin{aligned} \nabla\cdot\left[\frac{E(\mathbf{x})}{2(1+\nu_P)}\left(\nabla\Psi_n+\nabla\Psi_n^{\mathsf{T}}\right)+\frac{E(\mathbf{x})\nu_P}{(1+\nu_P)(1-2\nu_P)}(\nabla\cdot\Psi_n)\mathbf{I}\right]&=0 &&\text{in }\Psi_{n-1}(\Omega_0),\\ \Psi_n&=\frac{e^{-n/N_{\mathrm{steps}}-1}}{e^{-1}-1}\Psi_s &&\text{on }\partial\Psi_{n-1}(\Omega_0). \tag*{(27)} \end{aligned}

In (27), $\Omega_0$ denotes the template domain, $N_{\mathrm{steps}}=8$ is the number of pseudo-time steps, and $\Psi_s$ is the surface generated from (26), including, for ease of notation, also the smoothing and planar correction steps mentioned above. The Poisson ratio is defined as $\nu_P=0.3$ , and the stiffness $E(\mathbf{x})$ is set cell-wise to the largest ratio between the cell diameter and the length of one of its edges, so that already distorted cells are stiffer and deform less than well-shaped ones. The procedure is designed to gradually deform the template volume mesh onto a mesh compatible with the new surface mesh, preserving the original mesh topology. However, it does not guarantee that elements do not degenerate after each iteration. Therefore, after each pseudo-time step $n$ , the cells with non-positive Jacobian determinant or with aspect ratio above 100 are repaired by at most 40 gradient-descent steps on a cell aspect-ratio functional. We refer to [6] for the details of the mesh generation procedure. In what follows, we will denote by $\Omega_\theta$ the computational domain generated for the parameter $\theta\in\mathcal{X}$ .

Remark 3.4 (Handling of invalid geometries). The forward map $\mathcal{G}$ is defined on the whole parameter space $\mathbb{R}^{n_\theta}$ . However, not every parameter vector corresponds to a valid geometry, since the prior distribution is a standard normal on the PCA space of the geometric conditioning vector $\mathbf{c}\in\mathbb{R}^{d_c}$ (section 3.2), which is only an approximation of the true distribution of admissible geometries. As a result, the map can produce invalid geometries, i.e., shapes with self-intersecting surfaces, with negative inscribed radii, or surfaces that cannot be meshed with the multi-step linear elasticity extension (27). In our study, when the multi-step linear elasticity extension fails to produce a valid mesh at a step $n<N_{\mathrm{steps}}$ , we stop the iteration at that stage. If the extension fails at the first step, the forward map returns the template shape. Invalid parameters are thus mapped to an interpolation between the template and the generated surface, so that the forward map always returns an admissible domain and no prior mass is mapped to invalid geometries.

Notice that the numerical solver (section 3.3.3) requires a hierarchy of meshes for the geometric multigrid preconditioner. Hence, the mesh extension (27) is applied recursively to the entire hierarchy generated from the template mesh.

Numerical solution

On the computational domain $\Omega_\theta$ , we model the blood as an incompressible Newtonian fluid, described by a velocity field $\mathbf{u}:\Omega_\theta\times[0,T]\to\mathbb{R}^3$ and a pressure field $p:\Omega_\theta\times[0,T]\to\mathbb{R}$ , which satisfy the Navier–Stokes equations

\begin{aligned} \rho_f\partial_t\mathbf{u}+\rho_f(\mathbf{u}\cdot\nabla)\mathbf{u}-2\nabla\cdot(\mu_f\mathbb{D}(\mathbf{u}))+\nabla p&=0, &&\text{in }\Omega_\theta,\\ \nabla\cdot\mathbf{u}&=0, &&\text{in }\Omega_\theta. \tag*{(28)} \end{aligned}

In (28), $\mathbb{D}(\mathbf{u})$ is the symmetric part of the deformation tensor, $\rho_f$ is the blood density, and $\mu_f$ the dynamic viscosity, which includes a turbulent viscosity contribution based on a $\sigma$ -model [53]. We consider a cardiac period $T=1$ s. Notice that, for ease of notation, we omitted the dependence on $\theta$ of the solution $(\mathbf{u},p)$ .

Algorithm 6 Forward map

1: function ForwardMap(θ = (θ̃, α̃)) ▷ latent coordinate and stenosis degree
2: latent reconstruction c = U(L_c θ̃ + m_c)
3: generate surface mesh via LDDMM stochastic interpolants (26)
4: include stenosis defined by α̃
5: generate volume mesh (smoothing and multi-step elasticity (27))
6: compute (u, p) solution to (28)-(29)
7: evaluate observations 𝒢(θ) = 𝒪(p)
8: return 𝒢(θ)
9: end function

The system (28) is completed by the boundary conditions

\begin{cases} \mathbf{u}=\mathbf{u}_{\mathrm{in}} & \text{on }\Gamma_{\theta,\mathrm{in}},\\ \mathbf{u}=0 & \text{on }\Gamma_{\theta,\mathrm{wall}},\\ \left(-p\operatorname{Id}+\mu_f(\nabla\mathbf{u}+\nabla\mathbf{u}^{T})\right)\mathbf{n}=-P_i\mathbf{n} & \text{on }\Gamma_{\theta,i},\ i\in\{1,2,3,4\}. \tag*{(29)} \end{cases}

At the inlet $\Gamma_{\theta,\mathrm{in}}$ , we impose a given profile in the form $\mathbf{u}_{\mathrm{in}}(t,\mathbf{x})=a(t)\boldsymbol{\beta}(\mathbf{x})$ , obtained from a reference spatial profile $\boldsymbol{\beta}(\mathbf{x})$ , registered to the current inlet shape, and multiplied by a time-periodic function $a(t)$ [27]. On the vessel wall $\Gamma_{\theta,\mathrm{wall}}$ we impose no-slip Dirichlet boundary conditions, i.e., neglecting the interaction with the vessel wall. On the four outlet surfaces $\Gamma_{\theta,i}$ we impose an average pressure $P_i$ obtained via a so-called three-element Windkessel model, i.e., as the solution of the ODE

C_{d,i}\frac{\mathrm{d}}{\mathrm{d}t}(P_i-R_{p,i}Q_i)+\frac{P_i-R_{p,i}Q_i}{R_{d,i}}=Q_i. \tag*{(30)}

Equation (30) depends on the outgoing flow $Q_i$ on the outlet $\Gamma_{\theta,i}$ and on three parameters, $R_{p,i}$ , $R_{d,i}$ , and $C_{d,i}$ , which model the resistance and capacitance of the vascular system downstream. The Windkessel parameters are calibrated across the different shapes as described in [27], depending on a given flow split among the four outlets.

The system of equations (28)–(29) is discretized using a dual-splitting scheme in time and a high-order discontinuous Galerkin (DG2-DG1) discretization in space, with matrix-free geometric multigrid preconditioners in ExaDG [54] [6] with AVX SIMD vectorization and MPI parallelism on 32 cores. The geometric multigrid preconditioner requires a hierarchy of meshes. The finest mesh is composed of 2271 nodes and 1792 hexahedra, resulting in 159488 spatial degrees of freedom.

Observations

The last step in evaluating the forward map $\mathcal{G}(\theta)$ is defined by evaluating an observation operator $\mathcal{O}(\mathbf{u},p)$ on the numerical solution of (28)–(29).

For the examples considered in this work, we define a subset of 42 centerline points (shown in Figure 2 (observation operator) on the aortic arch of a sampled geometry) with corresponding cross-sections $\{S_i\}_{i=1}^{42}$ , and consider the following observations.

The pressure drop along the centerline at $t=0.25\,\mathrm{s}$ with respect to the average inlet pressure, defined as

\mathcal{O}_1(\mathbf{u},p)=\mathcal{O}_1(p)=\left(\frac{1}{|\Gamma_{\theta,\mathrm{in}}|}\int_{\Gamma_{\theta,\mathrm{in}}}\left(p(t=0.25\,\mathrm{s},\mathbf{x})-\frac{1}{|S_i|}\int_{S_i}p(t=0.25\,\mathrm{s},\mathbf{x})\,\mathrm{d}x\right)\,\mathrm{d}x\right)_{i=1}^{42}. \tag*{(31)}

The maximum pressure drop across all 42 centerline sections:

\mathcal{O}_2(\mathbf{u},p)=\mathcal{O}_2(p)=\max_{i=1,\ldots,42}\left|\mathcal{O}_1(p)_i\right|. \tag*{(32)}

Notice that we select the time instant $t=0.25\,\mathrm{s}$ as it lies right after the systolic peak.

We focus on observations only dependent on the pressure field. However, the framework can be analogously applied to other quantities of interest (centerline velocities, secondary-flow indices, wall shear stresses, etc.).

The steps of the forward map described in the previous sections are summarized in algorithm 6 and Figure 2.

Geometry-conditioned autoregressive LDDMM stochastic interpolant

The WFR method (§ 2.4) requires a differentiable surrogate forward map $\widetilde{G}$ to evaluate the likelihood score $-\nabla_{\theta}\Phi_{\mathrm{like}}(\theta)$ in (24).

To this end, we replace the high-fidelity PDE solver (§ 3.5) with a machine learning surrogate that emulates the pressure field over the cardiac cycle, conditioned on the geometry given by the centerline encoding vector.

The pressure surrogate is a flow-matching (stochastic-interpolant) model [55] of the transient lumen pressure field of a patient-specific aorta, and it relies heavily on the multi-step LDDMM registration of [27] to produce the training dataset of pressure fields supported on the same template.

The surrogate is trained to predict the pressure field $p$ on the whole domain. To increase the resolution of the predictions, the training pressure fields are interpolated onto a DG2 finite element space on the same mesh. Hence, conditioned on the geometry, the surrogate produces the scalar pressure $p$ at each of the $n_{\mathrm{nodes}} = 48\,384$ support points of the DG2 space for the physical time instants $t \in\{0, 0.05, 0.1, 0.15, 0.2, 0.25\}\,\mathrm{s}$ within the cardiac period. The SDE describing the model reads

\mathrm{d}I_t^p = b_t^p(I_t^p, x, c)\,\mathrm{d}t + \sigma_t^p\,\mathrm{d}W_t, \qquad I_{t=0}^p = p_0, \tag*{(33)}

where $x \in\mathbb{R}^{n_{\mathrm{nodes}}\times3}$ are the coordinates of the mesh vertices and $c$ the centerline encoding.

Rather than predicting each snapshot directly, the network learns the time-dependent drift $b^p$ that transports one snapshot into the next. For a source/target pair of snapshots at times $t$ and $t+\Delta$ , the interpolant is defined as

I^p(s) = (1-s)p_t + s p_{t+\Delta} + \sigma^p(s)\,\xi,

at interpolation time $s \in[0,1]$ , with Karras scheduler $\sigma^p(s) = \sigma_{\max} \left(\sigma_{\min}/\sigma_{\max}\right)^{s\xi}$ [56] and $\xi\sim\mathcal{N}(0,\mathrm{Id})$ . The neural network regresses the target drift $(p_{t+\Delta}-p_t) + \dot{\sigma}^p(s)\,\xi$ . At inference, the learned drift is integrated by the explicit Euler method with 10 substeps per interval to roll the pressure field forward across the instants $t \in\{0, 0.05, 0.1, 0.15, 0.2, 0.25\}\,\mathrm{s}$ .

At each evaluation, the drift network receives, for each node,

the current pressure state (scalar),
the node coordinates $x \in\mathbb{R}^{n_{\mathrm{nodes}}\times3}$ , lifted to Fourier features, and
a single global conditioning token built from the geometry conditioning vector $c(\theta) \in\mathbb{R}^{d_c}$ (the centerline control points and inscribed radii) together with two sinusoidal time embeddings (global cycle time and local interpolation time),

and outputs one scalar per node, i.e., the pressure drift field.

The node features are reordered along a Morton (Z-order) space-filling curve and packed into a dense 3D voxel grid. The ordering does not reproduce the connectivity of the mesh, but it makes the fields available on a regular grid, on which convolutional and transformer blocks are computationally more efficient than on the unstructured mesh, while keeping spatially close nodes close in memory. A depth-4 U-Net of conditioned transformer blocks processes the volume, with channel widths $64 \to128 \to256 \to512$ and a 1024-channel bottleneck. The output is un-permuted back to the nodes and projected to the scalar drift. The reported architecture uses 6 transformer blocks per stage and 6 bottleneck blocks ( $\approx112$ M parameters, see Table 1).

Component	Parameters	Share
Fourier PE (position + pressure)	2,128	< 0.1%
Conditioning embedding (geometry + 2×time)	98,624	0.1%
U-Net input embed + conditioning proj.	106,880	0.1%
Encoder 4 stages (ch. 64→512)	19,453,888	17.4%
Bottleneck (ch. 1024; conv + attention)	63,507,456	56.8%
Decoder 4 stages (ch. 512→64)	28,683,200	25.6%
Output (norm + linear → 1)	193	< 0.1%
Total	111,852,369	100%

Table 1. Parameter count of the stochastic interpolant drift network.

For the surrogate forward map, we restrict to the case $\tilde{\alpha}=0$ (without additional stenosis) degree, considering only the parameters of the latent space describing the geometric encoding (centerline and inscribed radii), i.e., $\theta= (\tilde{\theta}, 0)$ . This restriction is necessary, since the surrogate model is trained on numerical solutions obtained with the original training dataset, which contains the shapes without the stenosis. An extension with an additional parametric deformation would require the generation of a new training dataset. The strong dependence on the training distribution is one important limitation of this surrogate approach, and retraining or fine-tuning is needed whenever the training distribution must be updated. A further difference with respect to the forward map introduced in § 3.3 is that the volume mesh is generated using radial basis function (RBF) interpolation to extend the registered surface displacement into the interior of the mesh, instead of using the multi-step elasticity (27). RBF interpolation is needed to make the forward map amenable to batched execution on GPUs. Notice that the surrogate map does not require smoothing steps to improve the mesh quality, since the numerical solution of the Navier–Stokes equation, which might be sensitive to the computational mesh, is no longer needed.

Figure 6 reports the relative $l^2$ error on the pressure field and the relative error on the maximum pressure drop for the 52 test cases. Predicted and high-fidelity pressure snapshots are shown in Figure 7 for the cases with the maximum and the median $l^2$ -relative error.

Surrogate accuracy vs. time instant — Figure 6. LDDMM stochastic interpolant pressure surrogate: $l^2$ -relative error (left) and relative error on the maximum pressure drop (right), i.e., the output of the observation operator (32).

Snapshots of the maximum and median test cases — Figure 7. LDDMM stochastic interpolant pressure surrogate: snapshots of the cases with the maximum $l^2$ -relative error (left, case 18) and the median $l^2$ -relative error (right, case 43), averaged over all time instants.

The distribution of the pressure drop obtained from the surrogate for the training and test datasets (without stenosis) is shown in Figure 8.

Pressure-drop distribution of the surrogate dataset — Figure 8. Pressure drop distribution of the training and test datasets.

Notice that, as a result of the RBF interpolation, the surrogate is evaluated on volumetric meshes that differ from the ones it has been trained on. The errors in Figure 6 and Figure 7 are measured on the meshes produced by the multi-step elasticity extension. To quantify the effect of this substitution on the accuracy of the surrogate, we compare the two options case by case against the same high-fidelity solution. Figure 9 shows that the RBF interpolation affects the relative errors on the pressure drop and on the observed quantities $O_1$ and $O_2$ by at most 1%. Thus the surrogate is only mildly sensitive to the choice of the mesh extension approach. The surrogate forward map $\widetilde{G}$ used by the WFR method is summarized in algorithm 7 and Figure 2.

Effect of the RBF mesh extension on the surrogate error — Figure 9. Effect of the RBF interpolation on the accuracy of the surrogate, measured against the same high-fidelity solution. Left: difference between the $l^2$ relative errors on the pressure field, at each time instant, with mean, $\pm$ standard deviation band and min-max range over the 52 test cases. Right: difference between the relative errors on the centerline pressure drop $O_1$ and on the maximum pressure drop $O_2$ . Each marker corresponds to one test case, with median and interquartile range. Both left and right charts are in percentage points (pp): on the left, values of 1 correspond to 1% variation of the relative $l^2$ error.

Remark 3.5 (Approximation of the gradient of the surrogate output with respect to the parameters). At each discrete time instant in which the pressure is inferred, the LDDMM stochastic interpolant takes as input the coordinates of the nodes of the generated mesh, $x(\theta)$ , and the centerline encoding $c(\theta)$ , which is used as input to the conditioning branch. To approximate the gradient of the surrogate $\mathcal{T}_{t=0.25}^{p}(x(\theta),c(\theta))$ , we neglect the partial derivative with respect to the mesh nodes and consider a recursive gradient evaluation:

\frac{\partial\mathcal{T}_{t_k}^{p}}{\partial\theta} = \frac{\partial\mathcal{T}_{t_{k-1}}^{p}}{\partial\theta} + \Delta s\left( \frac{\partial b_{t_k}^{p}(\mathcal{T}_{t_{k-1}}^{p},x,c)}{\partial\mathcal{T}_{t_{k-1}}^{p}} \frac{\partial\mathcal{T}_{t_{k-1}}^{p}}{\partial\theta} + \frac{\cancel{\partial b_{t_k}^{p}(\mathcal{T}_{t_{k-1}}^{p},x,c)}}{\partial x} \frac{\partial x}{\partial\theta} + \frac{\partial b_{t_k}^{p}(\mathcal{T}_{t_{k-1}}^{p},x,c)}{\partial c} \frac{\partial c}{\partial\theta} \right),

:::`

Algorithm 7 Surrogate forward map

function $\operatorname{SurrogateForwardMap}(\theta= (\widetilde{\theta}, \widetilde{\alpha} = 0))$ ▷ latent coordinate only
latent reconstruction $c = U(L_c\widetilde{\theta} + m_c)$
generate surface mesh via LDDMM stochastic interpolants (26)
RBF interpolation to extend the surface mesh to the volume
$p =$ flow-matching pressure surrogate (33)
evaluate observations $\widetilde{G}(\widetilde{\theta}) = \mathcal{O}(p)$
batched evaluation of $\nabla_{\widetilde{\theta}}\widetilde{G}(\widetilde{\theta})$ by backward automatic differentiation
return $\{\widetilde{G}(\widetilde{\theta})\}, \{\nabla\widetilde{G}(\widetilde{\theta})\}$
end function

so the approximate gradient reduces to:

\frac{\partial}{\partial\theta}I^{p}_{t=0.25}\,s(x(\theta),c(\theta)) \approx\Delta s \sum_{k=1}^{\lfloor1/\Delta s\rfloor} \left[ \prod_{m=k+1}^{\lfloor1/\Delta s\rfloor} \left( \operatorname{Id}+\Delta s \frac{\partial b^{p}_{t_m}\left(I^{p}_{t_{m-1}},\mathbf{x},c\right)} {\partial I^{p}_{t_{m-1}}} \right) \right] \frac{\partial b^{p}_{t_k}\left(I^{p}_{t_{k-1}},\mathbf{x},c\right)}{\partial c} \frac{\partial c}{\partial\theta}.

Numerical results

This section presents the numerical results of the aortic shape inversion problem for the two selected types of observables (centerline and scalar pressure drops). We compare the chosen methods in terms of convergence, posterior distributions of the parameters and observables, and representative aortic geometries from the posterior parameter distributions.

The EAKI (section 2.2) is used for the case of centerline pressure drop (operator defined in (31)), as the observed/specified values of the pressure drop along the vessel centerline are informative enough to reconstruct a single geometry, while aBUS (section 2.3) and WFR (section 2.4) are used for the scalar pressure drop (32), for sampling multiple geometries associated with the same observed/prescribed value of the maximum pressure drop along the vessel centerline.

For EAKI and aBUS, the forward map is based on the high-fidelity PDE solver and on the mesh extended from the template (section 3). For WFR, which additionally requires the gradient $\frac{\partial G}{\partial\theta}$ , we use the surrogate map described in section 3.4. The high-fidelity solver for the forward map is run in parallel for each particle of the ensembles on an HTCondor queue with 4 nodes having 72, 128, 64, 128 cores, with each run requiring 32 CPU cores. The surrogate model runs in batches of 16 particles (forward and gradient evaluation) on a single NVIDIA H200 NVL4 141GB GPU.

In all cases, the prior distribution is the standard normal distribution in the latent parameter space, and the observation noise is relative, $\Sigma_{\eta} = \operatorname{diag}((0.1\,y)^2)$ , i.e., 10% of the assimilated value. The stenosis degree ( $\alpha$ ) is not considered in the surrogate model, as described in section 3.4. The setup of the three methods and the key differences are summarized in Table 2.

Method	Observation	Stenosis	Surrogate	Gradient	Gaussian posterior
EAKI	centerline pressure drop (31)	yes	no	no	yes
aBUS	pressure drop (32)	yes	no	no	no
WFR	pressure drop (32)	no	yes	yes	no

Table 2. Key properties of the three considered inversion methods.

Centerline pressure drop

We study the convergence behavior and the performance of the EAKI considering pressure data along the centerline.

To this end, we focus on two particular cases, i.e., the worst and the median test shapes based on the error of the pressure surrogate (cases 18 and 43, respectively, see Figure 7), with an additional stenosis with degree $\alpha= 0.7$ . With a little abuse of notation, we will refer to test case 18 and 43 in the following part also to indicate the geometries modified with the stenosis. For these two cases, we solve the Navier–Stokes equations (28) and evaluate the centerline pressure drop, defining the synthetic measurements as $y = \mathcal{O}_1(p)$ at $t = 0.25\,\mathrm{s}$ as defined in (31). The vector $y$ is then used as the observation for the EAKI method.

It is important to observe that the goal is not to reconstruct the original shapes (and, in general, this problem does not have a unique solution), but rather to obtain an ensemble of shapes – in the EAKI case, approximating a Gaussian measure – whose observation is close to the target one.

The most relevant parameters for the EAKI method are

the maximum number of iterations ( $T_{\max}$ ),
the ensemble size ( $N$ ),
the time step ( $\Delta\tau$ ), and
the jitter threshold ( $\sigma_{\mathrm{thr}}$ ).

In the numerical tests, we vary the ensemble size and the jitter threshold. At the first iteration, the EAKI ensemble particles are sampled from Latin hypercube sampling inside $[-3,3]^{n_\theta}$ . The number of iterations needed is usually small (below 10), and we set $\Delta\tau=1/2$ . For parameter studies related to $\Delta\tau$ and its interpretation, we refer to the original reference [32].

We show the convergence of EAKI for test case 18 in Figure 10 for $\sigma_{\mathrm{thr}}=0.1$ . The EAKI delivers, as a final result, an average shape given by the ensemble mean. The variance of the ensemble can be arbitrarily controlled through the jitter threshold $\sigma_{\mathrm{thr}}$ .

Centerline pressure drops for the EAKI ensemble over 10 iterations — Figure 10. Test case 18. Centerline pressure drops for the EAKI ensemble over $T_{\max}=10$ iterations. Jittering threshold $\sigma_{\mathrm{thr}}=0.1$ .

The results of jittering with $\sigma_{\mathrm{thr}}=0.1, 0.05, 0.01$ are shown in Figure 11 and Figure 12. The jittering can be used to avoid ensemble collapse, and one can see that the ensemble variance reduces considerably when decreasing the threshold from $\sigma_{\mathrm{thr}}=0.1$ to $\sigma_{\mathrm{thr}}=0.01$ . A higher ensemble variance is useful in practice to better study or model the uncertainty of the shape inversion problem: since we are evaluating the pressure drop with respect to the inlet pressure, the uncertainty is mostly distributed in the descending aortic arch rather than in the ascending arch, which is closer to the inlet, as can be seen in Figure 12.

Centerline pressure drops for the EAKI ensemble at the final iteration, varying the jittering threshold — Figure 11. Test case 18. Centerline pressure drops for the EAKI ensemble at the final iteration ( $T_{\max}=10$ ), varying the jittering threshold $\sigma_{\mathrm{thr}}$ .

Jittering results for test case 18 at three values of $\sigma_{\mathrm{thr}}$ — Figure 12. Test case 18. Jittering results corresponding to Figure 11. For each value of $\sigma_{\mathrm{thr}}$ , the picture shows the true (blue), the mean (red), and the farthest (green) EAKI ensemble members overlaid (left), the true geometry (center), and the EAKI mean geometry (right). The ensemble collapses to its mean from $\sigma_{\mathrm{thr}} = 0.1$ to $\sigma_{\mathrm{thr}} = 0.01$ .

The convergence of EAKI with respect to the relative misfit $\lVert y-G(\theta)\rVert_2/\lVert y\rVert_2$ is shown for the considered cases in Figure 13. It can be seen that the pressure drop profile at time $t = 0.25$ s is sufficient to identify a shape close to the true one after only 6 iterations.

Relative misfit over EAKI iterations for test cases 18 and 43 — Figure 13. Misfit $\lVert y - G(\theta) \rVert_2/\lVert y \rVert_2$ over the iterations. Left: test case 18. Right: test case 43.

The computational costs for the EAKI are reported in Table 3, detailing the ensemble size, the number of evaluations of the forward map, the wall-clock time (total and per evaluation), and the CPU time (in 32-core node hours).

Table 3. Computational cost of the aortic Bayesian inversion studies. $N_{\mathrm{ens}}$ is the ensemble size ( $N$ for EAKI, samples per level for aBUS, number of particles for WFR), $N_{\mathrm{fwd}}$ the total number of forward map evaluations, $T_{\mathrm{wall}}$ the wall-clock time of the whole study, $\bar{t}_{\mathrm{fwd}}$ the mean wall-clock time of a single forward map evaluation, and $T_{\mathrm{wall}}/N_{\mathrm{fwd}}$ the same time amortized over the study, which measures how effectively the evaluations are executed concurrently. The last two columns report the consumed resources: CPU time normalized to 32-core nodes and GPU time. For EAKI, a few forward evaluations did not terminate correctly, so that $N_{\mathrm{fwd}}$ is lower than $10\,N_{\mathrm{ens}}$ .

Run	$N_{\mathrm{ens}}$	$N_{\mathrm{fwd}}$	$T_{\mathrm{wall}}$ [h:mm]	$\bar{t}_{\mathrm{fwd}}$ [s]	$T_{\mathrm{wall}}/N_{\mathrm{fwd}}$ [s]	CPU [32-core node h]	GPU [h]
EAKI $O_1$ , $T_{\max}=10$ iterations, one high-fidelity solve per ensemble member
test 18, $\sigma_{\mathrm{thr}}=0.1$	100	1000	3:42	114.0	13.3	31.7	—
test 18, $\sigma_{\mathrm{thr}}=0.01$	252	2520	6:55	115.8	9.9	81.1	—
test 18, $\sigma_{\mathrm{thr}}=0.05$	252	2520	14:29	113.1	20.7	79.2	—
test 18, $\sigma_{\mathrm{thr}}=0.1$	252	2520	7:20	118.4	10.5	82.9	—
test 18, $\sigma_{\mathrm{thr}}=0.1$	500	4963	13:23	120.2	9.7	166.2	—
test 43, $\sigma_{\mathrm{thr}}=0.1$	100	1000	2:52	110.1	10.3	30.6	—
test 43, $\sigma_{\mathrm{thr}}=0.1$	252	2520	6:35	111.8	9.4	78.2	—
test 43, $\sigma_{\mathrm{thr}}=0.1$	500	4955	11:34	103.7	8.4	143.4	—
aBUS–SuS $O_2$ , one high-fidelity solve per sample and level
$O_2=5$ mmHg	500	1500	4:09	103.5	10.0	43.1	—
$O_2=10$ mmHg	500	2000	6:36	121.9	11.9	67.7	—
$O_2=20$ mmHg	250	1000	4:03	147.6	14.6	41.0	—
$O_2=20$ mmHg	500	2000	16:21	155.4	29.4	86.3	—
$O_2=20$ mmHg	1000	3000	10:08	129.3	12.2	107.8	—
$O_2=20$ mmHg	2000	6000	18:18	117.7	11.0	196.1	—
WFR $O_2$ , $T_{\max}=10$ iterations, one batched GPU job per iteration
$O_2=5$ mmHg	250	2500	4:03	5.7	5.8	0.5	4.0
$O_2=5$ mmHg	500	5000	7:57	5.7	5.7	1.0	7.9
$O_2=5$ mmHg	1000	10000	15:43	5.6	5.7	2.0	15.6
$O_2=10$ mmHg	250	2500	4:03	5.7	5.8	0.5	4.0
$O_2=10$ mmHg	500	5000	7:58	5.7	5.7	1.0	7.9
$O_2=10$ mmHg	1000	10000	15:50	5.7	5.7	2.0	15.7

Table 3.

A few forward evaluations did not terminate correctly in the two runs with $N_{\mathrm{ens}} = 500$ (only 4963 and 4955 valid evaluations instead of 5000). In fact, on a few meshes with particularly pronounced geometrical features, which require a very small adaptive time step, the numerical solver did not converge. These failures occurred only at the first EAKI iteration, since only in that case the particles are generated by Latin hypercube sampling, and not for the other considered methods.

Scalar pressure drop

In the next test, we consider the shape inversion problem based on a single observed value corresponding to the scalar pressure drop $O_2$ : the values are chosen by hand from $\{5\ \mathrm{mmHg}, 10\ \mathrm{mmHg}, 20\ \mathrm{mmHg}\}$ and corresponding to cases ranging from healthy to severe conditions. For this setting, we study the performance of the aBUS and the WFR methods.

aBUS

The main parameter that affects the convergence of aBUS is the number of samples ( $N$ ). For the other parameters, we use the values suggested in the literature [36], fixing the target acceptance rate $\hat{a}_{\mathrm{target}} = 0.44$ , the conditional probability between consecutive levels $p_0 = 0.1$ , and the batch size $n_B = 100$ .

We consider a pressure drop of 20 mmHg, for which the 10% observation noise gives $\sigma_{\eta}=2$ mmHg. The convergence with respect to the number of samples $N=250,500,1000,2000$ is shown in Figure 14. This pressure drop is considerably larger than physiological values for healthy subjects, and it can be obtained only in the presence of an additional stenosis. Figure 14 shows also the values of the observations associated to samples from the prior distribution (Figure 8).

Convergence plots of observed pressure drop for sample sizes N = 250, 500, 1000, and 2000 — Figure 14. Convergence of the aBUS method with respect to the sample size $N=250,500,1000,2000$ , with $\mathcal{O}_2=20$ mmHg.

The posterior distribution, as well as selected samples are shown in Figure 15. In particular, Figure 15, top-left, depicts the samples embedded in $\mathbb{R}^2$ with Uniform Manifold Approximation and Projection (UMAP). Unlike EAKI, the samples from the aBUS posterior distribution comprise multiple qualitatively different geometries rather than a single ensemble-averaged one.

UMAP embedding, posterior pressure-drop distribution, and selected aBUS shape-inversion samples — Figure 15. aBUS shape inversion with maximum pressure drop $O_2 = 20$ mmHg. Posterior distribution from aBUS and selected samples.

4.2.2. WFR. The WFR uses the surrogate forward map (§ 3.7), which does not include the additional stenosis. The method is therefore tested considering two moderate values of pressure drop, $\mathcal{O}_2 \in\{5,10\}$ mmHg, which do not lie too far from the prior/training distribution on which the surrogate model is trained (Figure 8). Without the stenosis, higher values of the pressure drop would correspond to a worse accuracy of the surrogate and, therefore, of the posterior distribution.

The main parameter considered for the convergence study is the number of samples $N=250,500,1000$ . Other relevant parameters are the number of steps of the gradient flow dynamics (fixed to $T_{\max}=10$ ), the minimum and maximum values of the cosine scheduler (fixed as $\gamma_{\min}=0.01$ , $\gamma_{\max}=0.1$ ), and the threshold on the effective sample size (ESS) set to $N/2$ .

Figure 16 and Figure 17 show the results for pressure drops of 5 mmHg and 10 mmHg, respectively. After sampling the parameters $\tilde{\theta}_j \in\mathbb{R}^{n_c}$ from the WFR posterior, the pressure values are computed by generating the corresponding computational domain and computing the pressure with the high-fidelity solver.

WFR scalar pressure drop distributions for ensemble sizes N=250, 500, and 1000, showing surrogate and high-fidelity posteriors, prior, means, target, and tolerance band — Figure 16. WFR, scalar pressure drop $O_2 = 5$ mmHg. Convergence of the WFR posterior distribution with respect to the ensemble size $N = 250, 500, 1000$ (from top to bottom).

Notice that values of the pressure drop farther from those of the training distribution correspond to a lower accuracy of the surrogate. As a result, when calculating the actual pressure drop corresponding to the posterior distribution, we observe a distribution shift. The case $\mathcal{O}_2=5$ mmHg is less affected than $\mathcal{O}_2=10$ mmHg. As observed in § 3.7 (see Figure 9), this shift is not due to the different approach used to generate the volume mesh, but should rather be attributed to the sensitivity of the surrogate to the training dataset.

These results can be improved by increasing the amount of training data associated with the prescribed values of the pressure drop, or by preventing the WFR method from exploring regions of the parameter space too far from the prior distribution. Nonetheless, the predicted posterior distributions match the prescribed values with good accuracy, which indicates that the distribution shift observed with the high-fidelity values is mainly due to the low accuracy of the surrogate model.

As for aBUS, the WFR posterior distributions can be visualized in $\mathbb{R}^2$ with a UMAP embedding, see Figure 18 and Figure 19.

UMAP 2d embedding, posterior density predicted and true plots, and representative sampled geometries labeled id=565, 4.64 mmHg; id=374, 2.10 mmHg; id=40, 3.63 mmHg; id=170, 3.28 mmHg; id=28, 4.00 mmHg; id=49, 2.94 mmHg; id=624, 3.92 mmHg; id=333, 4.57 mmHg; id=949, 3.60 mmHg; id=151, 9.82 mmHg — Figure 18. Representative sampled geometries from the WFR posterior distribution with $O_2 = 5$ mmHg.

UMAP 2d embedding, posterior density predicted and true plots, and representative sampled geometries labeled id=999, 11.71 mmHg; id=289, 7.85 mmHg; id=949, 13.57 mmHg; id=320, 5.77 mmHg; id=45, 9.52 mmHg; id=920, 10.51 mmHg; id=516, 7.07 mmHg; id=948, 9.97 mmHg; id=724, 8.01 mmHg; id=444, 27.58 mmHg — Figure 19. Representative sampled geometries from the WFR posterior distribution with pressure drop of 10 mmHg.

Comparisons: aBUS vs. WFR

We can compare the aBUS and WFR posterior distributions by restricting to the parameter space without the stenosis ( $\alpha= 0$ ). The resulting posterior distributions for pressure drops of 5 mmHg and 10 mmHg are shown in Figure 20, considering the same number of samples ( $N = 500$ ). In terms of misfit, aBUS is more accurate. WFR is accurate on the predicted posterior values but not on the posterior values computed with the high-fidelity solver. Hence, the higher misfit of the WFR results is mainly due to the limited accuracy of the surrogate outside the training distribution. This can be seen in the central plots in Figure 20, where the aBUS, the WFR, and the prior distributions are shown embedded in the first two principal components of a joint PCA, showing that aBUS remains much closer to the prior distribution in the parameter space.

Comparison of the posterior distributions of WFR and aBUS for pressure drops of 5 mmHg and 10 mmHg — Figure 20. Comparison of the posterior distributions of WFR and aBUS for pressure drops of 5 mmHg (top) and 10 mmHg (bottom).

In terms of computational costs, Table 3 provides a detailed comparison of the methods, considering the ensemble size, the number of evaluations of the forward map, the wall-clock time (total and per evaluation), and the CPU time (in 32-core node hours).

Using the surrogate forward map considerably reduces the computational cost. The HTCondor queue has 4 x64/openSUSE15 nodes with 72, 128, 64, 128 cores, respectively: each high-fidelity solver run requires 32 CPU cores, and a single CPU node can run multiple forward maps simultaneously. The surrogate runs in batches on a single NVIDIA H200 NVL4 141GB GPU.

Although this imbalance of resources results in a shorter wall-clock time for aBUS ( $\sim$ 4–6 hours against 8 hours for WFR), the compute time, measured in 32-core node hours and GPU hours, is higher for aBUS than for WFR, with an average forward evaluation time of 5.7 s for WFR against around 120 s for aBUS. Notice that the WFR steps are performed on CPU, while only the forward and gradient evaluations of the surrogate are performed on GPU. For example, compare the 67.7 node hours required by aBUS for the setting with pressure drop 10 mmHg and ensemble size $N_{\mathrm{ens}} = 500$ (2000 forward map evaluations), against 0.5 node hours plus 4 GPU hours spent by WFR for pressure drop 10 mmHg with $N_{\mathrm{ens}} = 250$ , which requires more forward map evaluations (2500), including also the gradients of the observation operator.

Conclusions

We addressed an inverse domain uncertainty quantification problem to determine virtual patient cohorts from target biomarkers in aortic coarctation. This approach has potential applications for data augmentation in the context of in silico testing of medical devices and treatments on patients with prescribed pathological conditions. Moreover, framed as a shape optimization problem, it can be applied to the estimation of the effect of geometric uncertainty introduced by the acquisition process of MRI and CT scans on numerical simulations. Using an initial cohort of synthetic patients, we define a latent parameter space (250-dimensional) that implicitly encodes centerlines and inscribed radii of the given population. The forward map for the generation of shapes is based on a ResNet-LDDMM aortic shape registration method [27], combined with LDDMM stochastic interpolants generative modeling [6], trained on the available set of shapes and corresponding pairwise registration maps. We considered three of the most promising UQ methods for large-scale and high-dimensional applications: EAKI, which relies on Gaussian approximations of the sought posterior distribution, and aBUS and WFR, which sample from multimodal and high-dimensional posterior distributions. We studied two different scenarios, considering observations of the pressure drop along the centerline (EAKI) or only its maximum value (aBUS and WFR). The results show that the posterior distributions reproduce the target measurements in all considered cases. Moreover, we observed that EAKI is more suitable for quickly extracting a single geometry with prescribed biomarkers from the generative model, whereas aBUS and WFR can be used to sample an entire cohort of aortic shapes associated with a prescribed biomarker of interest.

The present studies were conducted considering only the pressure drop. However, the general framework can be extended to other quantities, such as the velocity field and the wall shear stress, among others. The proposed methodology is also not limited to the considered application: it can be applied to any family of topologically equivalent shapes for which a training dataset is available. In this setting, registration maps onto a common template provide a natural geometrical encoding, shared by the generative prior and the surrogate model, and remove the need for a hand-crafted shape parametrization.

The restriction to a single topology is, at the same time, one of the main limitations of the framework: extending the generative prior to shapes with different topologies, e.g., aortas with a varying number of supra-aortic branches, requires different geometrical encodings. A further limitation concerns geometrical constraints that the generative prior does not enforce by construction, such as the absence of self-intersections or a minimum vessel radius. These constraints are currently satisfied only to the extent that they are represented in the training data, while invalid geometries are handled a posteriori by the forward map (remark (3.4)). Incorporating them explicitly, for instance through constrained generative models or additional terms in the likelihood, is a further direction of future research.

From a broader perspective, this work fits into the emerging field of generative design, also referred to as AI-assisted design, in which a generative model of admissible shapes is coupled with a physics-based forward model to find, or sample, geometries that attain prescribed performance targets. Further research directions include the improvement of the Fisher–Rao step of the WFR method through a more accurate particle density estimator, and the validation of the surrogate model with more training samples, which are needed to improve its accuracy.

Acknowledgements

Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – The Berlin Mathematics Research Center MATH+ (EXC-2046/1, EXC-2046/2, project ID: 390685689).

The authors used Claude Code (Anthropic) to polish the text and to check its grammar, and in part to support the implementation of the uncertainty quantification methods. The authors have reviewed and edited the output and take full responsibility for the content of this work.

References

[18]M. Dashti and A. M. Stuart. The Bayesian Approach to Inverse Problems. Ed. by R. Ghanem, D. Higdon, and H. Owhadi. Handbook of Uncertainty Quantification. Springer, 2017, pp. 311–428. DOI: 10.1007/978-3-319-12385-1_7.DOI
[19]C. P. Robert and G. Casella. Monte Carlo statistical methods. Vol. 2. Springer, 1999. DOI: 10.1007/978-1-4757-3071-5.DOI
[20]J. A. Christen and C. Fox. “A general purpose sampling algorithm for continuous distributions (the t-walk)”. In: Bayesian Analysis 5.2 (2010), pp. 263–281. DOI: 10.1214/10-BA603.DOI
[21]H. Haario, E. Saksman, and J. Tamminen. “Adaptive proposal distribution for random walk Metropolis algorithm”. In: Computational statistics 14.3 (1999), pp. 375–395. DOI: 10.1007/s001800050022.DOI
[22]H. Haario, E. Saksman, and J. Tamminen. “An adaptive Metropolis algorithm”. In: Bernoulli 7.2 (2001), pp. 223–242. DOI: 10.2307/3318737.DOI
[23]H. Haario et al. “DRAM: efficient adaptive MCMC”. In: Statistics and computing 16.4 (2006), pp. 339–354. DOI: 10.1007/s11222-006-9438-0.DOI
[24]S. L. Cotter et al. “MCMC methods for functions: modifying old algorithms to make them faster”. In: Statistical Science (2013), pp. 424–446.DOI
[25]T. Helin and R. Kretschmann. “Non-asymptotic error estimates for the Laplace approximation in Bayesian inverse problems”. In: Numerische Mathematik 150.2 (2022), pp. 521–549. DOI: 10.1007/s00211-021-01266-9.DOI
[26]P. Wacker. “Laplace’s method in Bayesian inverse problems”. In: arXiv preprint arXiv:1701.07989 (2017). DOI: 10.48550/arXiv.1701.07989.DOI
[27]F. Romor et al. “Shape-Informed Graph Neural Networks and Data Assimilation: Application to Velocity and Pressure Reconstruction in Aortic Blood Flow”. In: SIAM Journal on Imaging Sciences 19.2 (2026), pp. 710–751. DOI: 10.1137/25M1741820.DOI
[28]S. Särkkä and L. Svensson. Bayesian filtering and smoothing. Vol. 17. Cambridge university press, 2023. DOI: 10.1017/9781108917407.DOI
[29]M. A. Iglesias, K. J. Law, and A. M. Stuart. “Ensemble Kalman methods for inverse problems”. In: Inverse Problems 29.4 (2013), p. 045001. DOI: 10.1088/0266-5611/29/4/045001.DOI
[30]M. A. Iglesias. “A regularizing iterative ensemble Kalman method for PDE-constrained inverse problems”. In: Inverse Problems 32.2 (2016), p. 025002. DOI: 10.1088/0266-5611/32/2/025002.DOI
[31]D. Sanz-Alonso, A. Stuart, and A. Taeb. Inverse problems and data assimilation. Vol. 107. Cambridge University Press, 2023. DOI: 10.1017/9781009414319.
[32]D. Z. Huang et al. “Efficient derivative-free Bayesian inference for large-scale inverse problems”. In: Inverse Problems 38.12 (2022), p. 125006. DOI: 10.1088/1361-6420/ac99fa.DOI
[33]D. Straub and I. Papaioannou. “Bayesian updating with structural reliability methods”. In: Journal of Engineering Mechanics 141.3 (2015), p. 04014134. DOI: 10.1061/(asce)em.1943-7889.0000839.DOI
[34]S.-K. Au and J. L. Beck. “Estimation of small failure probabilities in high dimensions by subset simulation”. In: Probabilistic engineering mechanics 16.4 (2001), pp. 263–277. DOI: 10.1016/S0266-8920(01)00019-4.DOI
[35]W. Betz et al. “Bayesian inference with subset simulation: strategies and improvements”. In: Computer Methods in Applied Mechanics and Engineering 331 (2018), pp. 72–93. DOI: 10.1016/j.cma.2017.11.021.DOI
[36]P. Simon et al. “Parallelized adaptive Bayesian updating with structural reliability methods for inference of large engineering models”. In: Advances in Structural Engineering 29.1 (2026), pp. 150–175. DOI: 10.1177/13694332251346848.
[37]F. R. Crucinio and S. Pathiraja. “Sequential Monte Carlo approximations of Wasserstein–Fisher–Rao gradient flows”. In: arXiv preprint arXiv:2506.05905 (2025). DOI: 10.48550/arXiv.2506.05905.DOI
[38]Y. Yan, K. Wang, and P. Rigolet. “Learning Gaussian mixtures using the Wasserstein–Fisher–Rao gradient flow”. In: The Annals of Statistics 52.4 (2024), pp. 1774–1795. DOI: 10.1214/24-aos2416.DOI
[39]C. Bodnar et al. “A Foundation Model for the Earth System”. In: Nature (May 2025). DOI: 10.1038/s41586-025-09005-y.DOI
[40]J. Xiang et al. Native and Compact Structured Latents for 3D Generation. 2025. arXiv: 2512.14692 [cs.CV].arxiv.org/abs/2512.14692
[41]B. Guillard et al. “DeepMesh: Differentiable Iso-Surface Extraction”. In: IEEE Transactions on Pattern Analysis and Machine Intelligence 46.11 (2024), pp. 7072–7087. DOI: 10.1109/TPAMI.2024.3392291.DOI
[42]Y. Chen et al. “Efficient, multimodal, and derivative-free bayesian inference with Fisher–Rao gradient flows”. In: Inverse Problems 40.12 (2024), p. 125001. DOI: 10.1088/1361-6420/ad847b.DOI
[43]B. Z. I. Botev and D. P. Kroese. “Efficient Monte Carlo simulation via the generalized splitting method”. In: Statistics and Computing 22.1 (2012), pp. 1–16. DOI: 10.1007/s11222-010-9201-4.
[44]F. Cèrou et al. “Sequential Monte Carlo for rare event estimation”. In: Statistics and computing 22.3 (2012), pp. 795–808. DOI: 10.1007/s11222-011-9231-6.DOI
[45]P. Del Moral, A. Doucet, and A. Jasra. “Sequential monte carlo samplers”. In: Journal of the Royal Statistical Society Series B: Statistical Methodology 68.3 (2006), pp. 411–436. DOI: 10.1111/j.1467-9868.2006.00553.x.DOI
[46]I. Paipanomoun et al. “MCMC algorithms for subset simulation”. In: Probabilistic Engineering Mechanics 41 (2015), pp. 89–103. DOI: 10.1016/j.probengmech.2015.06.006.
[47]J.-J. Zhu. Inclusive KL Gradient Flows: Otto-Wasserstein, Fisher-Rao-Gaussian, and Local-Estimator Dynamics. 2026. DOI: 10.48550/arXiv.2411.00214. arXiv: 2411.00214 [stat.ML].DOI
[48]L. Goubergrits et al. “CT-based analysis of left ventricular hemodynamics using statistical shape modeling and computational fluid dynamics”. In: Frontiers in Cardiovascular Medicine 9 (2022), p. 901902. DOI: 10.3389/fcvm.2022.901902.DOI
[49]B. Thamsen et al. “Unsupervised learning and statistical shape modeling of the morphometry and hemodynamics of coarctation of the aorta”. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer. 2020, pp. 776–785. DOI: 10.1007/978-3-030-59719-1_75.DOI
[50]B. Thamsen et al. “Synthetic database of aortic morphometry and hemodynamics: overcoming medical imaging data availability”. In: IEEE Transactions on Medical Imaging 40.5 (2021), pp. 1438–1449. DOI: 10.1109/tmi.2021.3057496.DOI
[51]J. Versnjak et al. “Deep learning based assessment of hemodynamics in the coarctation of the aorta: comparison of bidirectional recurrent and convolutional neural networks”. In: Frontiers in Physiology 15 (2024), p. 1288339. DOI: 10.3389/fphys.2024.1288339.DOI
[52]P. Yevtushenko et al. “Deep learning based centerline-aggregated aortic hemodynamics: an efficient alternative to numerical modeling of hemodynamics”. In: IEEE Journal of Biomedical and Health Informatics 26.4 (2021), pp. 1815–1825. DOI: 10.1109/jbhi.2021.3116764.DOI
[53]F. Nicoud et al. “Using singular values to build a subgrid-scale model for large eddy simulations”. In: Physics of Fluids 23.8 (2011). DOI: 10.1063/1.3623274.DOI
[54]D. Arndt et al. “ExaDG: High-Order Discontinuous Galerkin for the Exa-Scale”. In: Software for Exascale Computing – SPPEXA 2016-2019. Ed. by H.-J. Bungartz et al. Cham: Springer International Publishing, 2020, pp. 189–224. DOI: 10.1007/978-3-030-47956-5_8.DOI
[55]G. Kohl, L.-W. Chen, and N. Thüerey. “Benchmarking autoregressive conditional diffusion models for turbulent flow simulation”. In: Neural Networks 199 (2026), p. 108641. DOI: https://doi.org/10.1016/j.neunet.2026.108641.
[56]T. Karras et al. “Elucidating the design space of diffusion-based generative models”. In: Advances in neural information processing systems 35 (2022), pp. 26565–26577. DOI: 10.5220/068431–1926.arxiv.org/abs/2206.00364

Inverse domain uncertainty quantification in hemodynamics: sampling virtual patient cohorts under constraints

Abstract

Introduction

Inverse domain uncertainty quantification

Problem setting

Ensemble Adjusted Kalman Inversion (EAKI)

Adaptive Bayesian Updating with Structural reliability (aBUS)

Wasserstein–Fisher–Rao sequential Monte Carlo (WFR)

Forward problem

Geometry dataset and shape registration

Parameter space

Parameter-to-observable map

Generative model for surface mesh

Mesh

Numerical solution

Observations

Geometry-conditioned autoregressive LDDMM stochastic interpolant

Numerical results

Centerline pressure drop

Scalar pressure drop

aBUS

Comparisons: aBUS vs. WFR

Conclusions

Acknowledgements

References

References

Publication

Taxonomy

Authors

AI contributions

More in Mathematics

More in Mathematics