We describe a stochastic method to characterize uncertainties in numerical models that have a single real-valued output (observable). The rationale estimates the probability density function (PDF) of the observable by inverting the characteristic function (CF), which is computed by quadrature. This approach is seldom used due to the highly oscillatory integrals that define the CF, which cause traditional quadrature algorithms to fail. Instead, we overcome these high oscillations by using an oscillatory Levin quadrature algorithm. We specifically choose the implementation proposed by Webster and Evans for the one-dimensional case and generalize it to the multivariate case. Using the Webster-Evans/Levin rule (WEL), we build a probabilistic adjoint operator (PrAdO).
This work aims to enhance the small amount of data obtained from simulations using statistical methods, in order to improve the analysis of the underlying physical properties of plasma behavior, particularly electron density and electron temperature in the boundary region of the tokamak. Electron density and electron temperature are important observables for characterizing plasma behavior in the boundary region of the tokamak, as they influence transport properties and confinement regimes. While the simulations used in this work are fully deterministic, we aim to extract statistical information by interpreting numerical simulations as samples from underlying random fields. This approach allows for the construction of marginal and joint probability density functions (PDFs) that provide physical insight beyond standard deterministic interpretation and capture key structural features of the plasma behavior in the boundary region of the tokamak, including the role of triangularity, the impact of the diffusion coefficient, and the implications for plasma behavior. Numerical simulations are performed with the Gkeyll gyrokinetic code. Despite their resolution, the computational cost limits the number of available simulation data, suggesting a role for advanced analysis techniques capable of extracting meaningful physical insights from a small dataset. To this end, we explore a statistical framework based on probabilistic learning on manifolds (PLoM), a nonparametric method designed for small-sample inference.
Inverse uncertainty quantification is essential in applied mechanics, particularly for characterizing the uncertainty in model parameters based on uncertain observations of a model's output. Approaches via Bayesian inference often struggle with the intractability of the likelihood due to the complexities of real-world simulators and experiments. This paper addresses these challenges by proposing a novel approach to tractable inverse uncertainty quantification using probabilistic circuits. The method circumvents the need for an explicit likelihood by directly learning the joint distribution of system inputs and observations through density estimation. Additionally, this approach allows for various probabilistic queries, including marginals, posterior distributions, and sensitivity analysis on a single model, while ensuring computational tractability. Through numerical experiments on representative problems in applied mechanics, we demonstrate the effectiveness and flexibility of our method to real-world problems.
Shape and topology optimization with uncertainties typically requires to solve repeatedly the forward problems corresponding to many random samples during each optimization step, which occupies a large amount of computational cost. In this paper, we utilize the multimodes Monte Carlo method to develop eficient shape and topology optimization algorithms to solve exterior Bernoulli free boundary problems with random diffusion coeficients. The shape optimization with the shape functional of Kohn-Vogelius type and Neumann data-tracking type is considered to address the overdetermined problem. We construct the shape optimization algorithm to solve the Kohn-Vogelius problem by the shape gradient flow method. To tackle the Neumann data-tracking problem, the phase-ield method of Allen-Cahn type is proposed for topology optimization to track the evolution of the interface. The optimization process can be accelerated to some extent by using the multimodes Monte Carlo method to solve the governing stochastic equation. Numerical examples in 2D and 3D are presented to show the effectiveness and eficiency of the proposed algorithms.
Accurately estimating rare event probabilities for systems subject to uncertainty and randomness is an inherently difficult task that calls for dedicated tools and methods. One way to improve estimation efficiency on difficult rare event estimation problems is to leverage gradients of the computational model representing the system under consideration to explore the rare event faster and more reliably. We present a novel approach for estimating rare event probabilities using such model gradients by drawing on a technique to generate samples from non-normalized posterior distributions in Bayesian inference-the Stein variational gradient descent. We propagate samples generated from a tractable input distribution towards a near-optimal rare event importance sampling distribution by exploiting a similarity of the latter with Bayesian posterior distributions. Sample propagation takes the shape of passing samples through a sequence of invertible transforms such that their densities can be tracked and used to construct an unbiased importance sampling estimate of the rare event probability-the Stein variational rare event estimator. We discuss settings and parametric choices of the algorithm and suggest a method for balancing convergence speed with stability by choosing the step width or base learning rate adaptively. We analyze the method's performance on several analytical test functions and two engineering examples in low to high stochastic dimensions (d = 2-1500) and find that it consistently outperforms other state-of-the-art gradient-based rare event simulation methods.
A foundational challenge in uncertainty quantification involves estimating a probability measure on the space of uncertain parameters such that its push-forward through a computational model matches an observed probability measure on the output data associated with quantities of interest (QoI). When multiple, distinct sets of observational data are available, the desired parameter measure should simultaneously satisfy multiple push-forward constraints associated with various subsets of the QoI. In this work, we present a convergent measure-theoretic framework for solving this problem based on an iterative application of Data-Consistent Inversion (DCI). We first rigorously establish the theoretical optimality of the DCI solution to the standard problem, proving that it minimizes the f-divergence over the space of all possible pullback measures that satisfy the push-forward constraint. This optimality property provides the foundation for our iterative DCI scheme, which is shown to converge to a solution of the multiple push-forward constraint problem. This iterative solution minimizes the cumulative f-divergence across all constraints and, under uniform initializations, represents the maximal entropy solution (the I-projection) onto the intersection of the solution sets. We provide a rigorous convergence analysis for the proposed method and demonstrate its practical utility through numerical examples, including a high-dimensional parameter space governed by partial differential equations, where the iterative approach robustly avoids the complexities associated with approximating high-dimensional joint observed measures.
Predicting fuel assembly bow in pressurized water reactors requires solving tightly coupled fluid-structure interaction problems, whose direct simulations can be computationally prohibitive, making large-scale uncertainty quantification (UQ) very challenging. This work introduces a general mathematical framework for coupling Gaussian process (GP) surrogate models representing distinct physical solvers, aimed at enabling rigorous UQ in coupled multiphysics systems. A theoretical analysis establishes that the predictive variance of the coupled GP system remains bounded under mild regularity and stability assumptions, ensuring that uncertainty does not grow uncontrollably through the iterative coupling process. The methodology is then applied to the coupled hydraulic-structural simulation of fuel assembly bow, enabling global sensitivity analysis and full UQ at a fraction of the computational cost of direct code coupling. The results demonstrate accurate uncertainty propagation and stable predictions, establishing a solid mathematical basis for surrogate-based coupling in large-scale multiphysics simulations.
Deep Gaussian Processes (DGPs) are powerful surrogate models known for their flexibility and ability to capture complex functions. However, extending them to multi-output settings remains challenging due to the need for efficient dependency modeling. We propose the Deep Intrinsic Coregionalization Multi-Output Gaussian Process (deepICMGP) surrogate for computer simulation experiments involving multiple outputs, which extends the Intrinsic Coregionalization Model (ICM) by introducing hierarchical coregionalization structures across layers. This enables deepICMGP to effectively model nonlinear and structured dependencies between multiple outputs, addressing key limitations of traditional multi-output GPs. We benchmark deepICMGP against state-of-the-art models, demonstrating its competitive performance. Furthermore, we incorporate active learning strategies into deepICMGP to optimize sequential design tasks, enhancing its ability to efficiently select informative input locations for multi-output systems.
Gaussian process regression techniques have been used in fluid mechanics for the reconstruction of flow fields from a reduction-of-dimension perspective. A main ingredient in this setting is the construction of adapted covariance functions, or kernels, to obtain such estimates. In this paper, we derive physics-informed kernels for simulating two-dimensional velocity fields of an incompressible (divergence-free) flow around aerodynamic profiles. These kernels allow to define Gaussian process priors satisfying the incompressibility condition and the prescribed boundary conditions along the profile in a continuous manner. Such physical and boundary constraints can be applied to any pre-defined scalar kernel in the proposed methodology, which is very general and can be implemented with high flexibility for a broad range of engineering applications. Its relevance and performances are illustrated by numerical simulations of flows around a cylinder and a NACA 0412 airfoil profile, for which no observation at the boundary is needed at all.
The ability to design effective experiments is crucial for obtaining data that can substantially reduce the uncertainty in the predictions made using computational models. An optimal experimental design (OED) refers to the choice of a particular experiment that optimizes a particular design criteria, e.g., maximizing a utility function, which measures the information content of the data. However, traditional approaches for optimal experimental design typically require solving a large number of computationally intensive inverse problems to find the data that maximizes the utility function. Here, we introduce two novel OED criteria that are specifically crafted for the data consistent inversion (DCI) framework, but do not require solving inverse problems. DCI is a specific approach for solving a class of stochastic inverse problems by constructing a pullback measure on uncertain parameters from an observed probability measure on the outputs of a quantity of interest (QoI) map. While expected information gain (EIG) has been used for both DCI and Bayesian based OED, the characteristics and properties of DCI solutions differ from those of solutions to Bayesian inverse problems which should be reflected in the OED criteria. The new design criteria developed in this study, called the expected scaling effect and the expected skewness effect, leverage the geometric structure of pre-images associated with observable data sets, allowing for an intuitive and computationally efficient approach to OED. These criteria utilize singular value computations derived from sampled and approximated Jacobians of the experimental designs. We present both simultaneous and sequential (greedy) formulations of OED based on these innovative criteria. Numerical results demonstrate the effectiveness in our approach for solving stochastic inverse problems.
Serial ensemble filters implement triangular probability transport maps to reduce high-dimensional inference problems to sequences of state-by-state univariate inference problems. The univariate inference problems are solved by sampling posterior probability densities obtained by combining constructed prior densities with observational likelihoods according to Bayes' rule. Many serial filters in the literature focus on representing the marginal posterior densities of each state. However, rigorously capturing the conditional dependencies between the different univariate inferences is crucial to correctly sampling multidimensional posteriors. This work proposes a new serial ensemble filter, called the copula rank histogram filter (CoRHF), that seeks to capture the conditional dependency structure between variables via empirical copula estimates; these estimates are used to rigorously implement the triangular (state-by-state univariate) Bayesian inference. The success of the CoRHF is demonstrated on two-dimensional examples and the Lorenz '63 problem. A practical extension to the high-dimensional setting is developed by localizing the empirical copula estimation, and is demonstrated on the Lorenz '96 problem.
Hyper-differential sensitivity analysis with respect to model discrepancy was recently developed to enable uncertainty quantification for optimization problems. The approach consists of two primary steps: (i) Bayesian calibration of the discrepancy between high- and low-fidelity models, and (ii) propagating the model discrepancy uncertainty through the optimization problem. When high-fidelity model evaluations are limited, as is common in practice, the prior discrepancy distribution plays a crucial role in the uncertainty analysis. However, specification of this prior is challenging due to its mathematical complexity and many hyper-parameters. This article presents a novel approach to specify the prior distribution. Our approach consists of two parts: (1) an algorithmic initialization of the prior hyper-parameters that uses existing data to initialize a hyper-parameter estimate, and (2) a visualization framework to systematically explore properties of the prior and guide tuning of the hyper-parameters to ensure that the prior captures the appropriate range of uncertainty. We provide detailed mathematical analysis and a collection of numerical examples that elucidate properties of the prior that are crucial to ensure uncertainty quantification.
In this paper we introduce a novel approach for studying the influence of the choice of the prior distribution on Bayesian inference results. We define perturbed-law-based sensitivity indices (PLI) for Bayesian inference in order to provide a quantitative description of the impact of a lack of knowledge about the prior on Bayesian inference results. Based on recent work in the field of robustness analysis, these indices rely on perturbations of a reference prior, which are based on the concept of Fisher distance taken from information geometry. We also show that the proposed PLI can be reformulated as the relative variation of probabilities of rare events, which facilitates their practical computation. The proposed approach is showcased through several application examples involving Bayesian inverse problems with varying complexity. Results emphasize that the proposed approach enables the identification of parameters whose prior distribution choice has a significant impact on the inference results. Furthermore, it remains feasible in the case of Bayesian inverse problems with nonlinear forward models and possibly high-dimensional inputs, while allowing an arbitrary perturbation level for the prior.
In this paper, we analyze the numerical approximation of the Navier-Stokes problem over a bounded polygonal domain in ℝ^2, where the initial condition is modeled by a log-normal random field. This problem usually arises in the area of uncertainty quantification. We aim to compute the expectation value of linear functionals of the solution to the Navier-Stokes equations and perform a rigorous error analysis for the problem. In particular, our method includes the finite element, fully-discrete discretizations, truncated Karhunen-Loéve expansion for the realizations of the initial condition, and lattice-based quasi-Monte Carlo (QMC) method to estimate the expected values over the parameter space. Our QMC analysis is based on randomly-shifted lattice rules for the integration over the domain in high-dimensional space, which guarantees the error decays with 𝒪(N^-1+δ), where N is the number of sampling points, δ>0 is an arbitrary small number, and the constant in the decay estimate is independent of the dimension of integration.
In this paper, we introduce an advanced surrogate modeling approach using a cut high-dimensional model representation (cut-HDMR) framework, enhanced with clustering based multiple anchor points. The accuracy of cut-HDMR models is dependent upon their spatial proximity to these anchors. To optimize this dependency, we employ centroidal Voronoi tessellation (CVT) for efficient clustering, which systematically reduces the sum of errors between each sample within a cluster and its centroid, thereby setting the centroids as optimal anchor points. This setup significantly reduces the distance between system inputs and the nearest anchor, enabling precise cut-HDMR expansion tailored to each anchor. The selection of anchors for new input samples is directed by the nearest centroid principle, and accurate computation of the output via cut-HDMR is then expected. A thorough error analysis of our CVT-based multianchor HDMR is provided. Simulation and numerical experiment results involving high-dimensional integrals and elliptic stochastic partial differential equations indicate that the CVT-based multiple anchor points selection strategy not only mitigates the drawbacks of single, improperly placed anchor points but also markedly improves accuracy beyond that achieved by averaging multiple cut-HDMR expansions.
Computer models are widely used for the prediction of complex physical phenomena. Based on observations of these physical phenomena, it is possible to calibrate the model parameters. In most cases, such computer models are misspecified, and the calibration process must be improved by including a model error term. The model error hyperparameters are, however, rarely learned jointly with the model parameters to reduce the dimensionality of the problem. Sequential and nonsequential approaches have been introduced to estimate the hyperparameters. The former, such as the Kennedy and O'Hagan (KOH) framework, estimates the model error hyperparameters before calibrating the model parameters. The latter, such as the full maximum a posteriori (FMP), introduces a functional dependence between the model parameters and the model error hyperparameters. Despite being more reliable in some cases (bimodality, e.g.), the FMP method still fails to estimate correctly the posterior distribution shape. This work proposes a new methodology for treating the model error term in computer code calibration. It builds upon the KOH and FMP framework. Called the complete maximum a posteriori (CMP) method, it provides a closed-form expression for the marginalization integral over the model error hyperparameters, significantly reducing the dimensionality of the calibration problem. Such expression relies on a set of assumptions that are more general and less stringent than the ones usually employed. The CMP method is applied to four examples of increasing complexity, from elementary to real fluid dynamics problems, including or not bimodality. Compared to the true reference solution and unlike the KOH and FMP, the CMP method correctly captures the shape of the posterior distribution, including all modes and their weights. Moreover, it provides an accurate estimate of the distribution tails.
We consider an initial value problem of a linear fractional differential equation (FDE) with a Caputo derivative. The right-hand side of the FDE includes a coefficient, which depends on a random variable to model a variability in uncertainty quantification. Consequently, the solution of the FDE becomes a random process. We expand the random process into the generalized polynomial chaos, where a series consists of unknown time-dependent coefficients and predetermined basis polynomials. A stochastic Galerkin approach yields a larger deterministic linear system of FDEs, whose solution represents an approximation of the coefficient functions. We show that the approximations of the stochastic Galerkin systems converge to the exact random process, provided that some common assumptions are satisfied. Furthermore, we present results of numerical computations, where initial value problems of the stochastic Galerkin systems are solved. These results verify the deduced convergence property.
