Computational Uncertainty Quantification Research Articles

CUQIpy: II. Computational uncertainty quantification for PDE-based inverse problems in Python

Inverse problems, particularly those governed by Partial Differential Equations (PDEs), are prevalent in various scientific and engineering applications, and uncertainty quantification (UQ) of solutions to these problems is essential for informed decision-making. This second part of a two-paper series builds upon the foundation set by the first part, which introduced CUQIpy, a Python software package for computational UQ in inverse problems using a Bayesian framework. In this paper, we extend CUQIpy’s capabilities to solve PDE-based Bayesian inverse problems through a general framework that allows the integration of PDEs in CUQIpy, whether expressed natively or using third-party libraries such as FEniCS. CUQIpy offers concise syntax that closely matches mathematical expressions, streamlining the modeling process and enhancing the user experience. The versatility and applicability of CUQIpy to PDE-based Bayesian inverse problems are demonstrated on examples covering parabolic, elliptic and hyperbolic PDEs. This includes problems involving the heat and Poisson equations and application case studies in electrical impedance tomography and photo-acoustic tomography, showcasing the software’s efficiency, consistency, and intuitive interface. This comprehensive approach to UQ in PDE-based inverse problems provides accessibility for non-experts and advanced features for experts.

CUQIpy: I. Computational uncertainty quantification for inverse problems in Python

This paper introduces CUQIpy, a versatile open-source Python package for computational uncertainty quantification (UQ) in inverse problems, presented as Part I of a two-part series. CUQIpy employs a Bayesian framework, integrating prior knowledge with observed data to produce posterior probability distributions that characterize the uncertainty in computed solutions to inverse problems. The package offers a high-level modeling framework with concise syntax, allowing users to easily specify their inverse problems, prior information, and statistical assumptions. CUQIpy supports a range of efficient sampling strategies and is designed to handle large-scale problems. Notably, the automatic sampler selection feature analyzes the problem structure and chooses a suitable sampler without user intervention, streamlining the process. With a selection of probability distributions, test problems, computational methods, and visualization tools, CUQIpy serves as a powerful, flexible, and adaptable tool for UQ in a wide selection of inverse problems. Part II of the series focuses on the use of CUQIpy for UQ in inverse problems with partial differential equations.

An Adaptive Sampling and Domain Learning Strategy for Multivariate Function Approximation on Unknown Domains

Many problems arising in computational science and engineering can be described in terms of approximating a smooth function of variables, defined over an unknown domain of interest , from sample data. Here both the underlying dimensionality of the problem (in the case ) as well as the lack of domain knowledge—with potentially irregular and/or disconnected—are confounding factors for sampling-based methods. Naïve approaches for such problems often lead to wasted samples and inefficient approximation schemes. For example, uniform sampling can result in upward of 20% wasted samples in some problems considered herein. In applications such as surrogate model construction in computational uncertainty quantification, the high cost of computing samples necessitates a more efficient sampling procedure. Over the last several years methods for computing such approximations from sample data have been studied in the case of irregular domains, and the advantages of computing sampling measures depending on an approximation space of have been shown. More specifically, such approaches confer advantages such as stability and well-conditioning, with an asymptotically optimal sample complexity scaling . The recently proposed adaptive sampling for general domains (ASGD) strategy is one such technique to construct these sampling measures. The main contribution of this paper is a procedure to improve upon the ASGD approach by adaptively updating the sampling measure in the case of unknown domains. We achieve this by first introducing a general domain adaptivity strategy, which computes an approximation of the function and domain of interest from the sample points. Second, we propose an adaptive sampling strategy, termed adaptive sampling for unknown domains (ASUD), which generates sampling measures over a domain that may not be known in advance, based on the ideas introduced in the ASGD approach. We then derive (weighted) least squares and augmented least squares techniques for polynomial approximation on unknown domains. We present numerical experiments demonstrating the efficacy of the adaptive sampling techniques with least squares–based polynomial approximation schemes. Our results show that the ASUD approach consistently achieves errors the same as or smaller than uniform sampling, but using fewer, and often significantly fewer, function evaluations.

CAS4DL: Christoffel adaptive sampling for function approximation via deep learning

The problem of approximating smooth, multivariate functions from sample points arises in many applications in scientific computing, e.g., in computational Uncertainty Quantification (UQ) for science and engineering. In these applications, the target function may represent a desired quantity of interest of a parameterized Partial Differential Equation (PDE). Due to the large cost of solving such problems, where each sample is computed by solving a PDE, sample efficiency is a key concern in these applications. Recently, there has been increasing focus on the use of Deep Neural Networks (DNN) and Deep Learning (DL) for learning such functions from data. In this work, we propose an adaptive sampling strategy, CAS4DL (Christoffel Adaptive Sampling for Deep Learning) to increase the sample efficiency of DL for multivariate function approximation. Our novel approach is based on interpreting the second to last layer of a DNN as a dictionary of functions defined by the nodes on that layer. With this viewpoint, we then define an adaptive sampling strategy motivated by adaptive sampling schemes recently proposed for linear approximation schemes, wherein samples are drawn randomly with respect to the Christoffel function of the subspace spanned by this dictionary. We present numerical experiments comparing CAS4DL with standard Monte Carlo (MC) sampling. Our results demonstrate that CAS4DL often yields substantial savings in the number of samples required to achieve a given accuracy, particularly in the case of smooth activation functions, and it shows a better stability in comparison to MC. These results therefore are a promising step towards fully adapting DL towards scientific computing applications.

Extrapolated Polynomial Lattice Rule Integration in Computational Uncertainty Quantification

Extrapolated Polynomial Lattice Rule Integration in Computational Uncertainty Quantification

Efficient Uncertainty Quantification and Variance-Based Sensitivity Analysis in Epidemic Modelling Using Polynomial Chaos

The use of epidemic modelling in connection with spread of diseases plays an important role in understanding dynamics and providing forecasts for informed analysis and decision-making. In this regard, it is crucial to quantify the effects of uncertainty in the modelling and in model-based predictions to trustfully communicate results and limitations. We propose to do efficient uncertainty quantification in compartmental epidemic models using the generalized Polynomial Chaos (gPC) framework. This framework uses a suitable polynomial basis that can be tailored to the underlying distribution for the parameter uncertainty to do forward propagation through efficient sampling via a mathematical model to quantify the effect on the output. By evaluating the model in a small number of selected points, gPC provides illuminating statistics and sensitivity analysis at a low computational cost. Through two particular case studies based on Danish data for the spread of Covid-19, we demonstrate the applicability of the technique. The test cases consider epidemic peak time estimation and the dynamics between superspreading and partial lockdown measures. The computational results show the efficiency and feasibility of the uncertainty quantification techniques based on gPC, and highlight the relevance of computational uncertainty quantification in epidemic modelling.

The Wasserstein distances between pushed-forward measures with applications to uncertainty quantification

In the study of dynamical and physical systems, the input parameters are often uncertain or randomly distributed according to a measure $\varrho$. The system's response $f$ pushes forward $\varrho$ to a new measure $f\circ \varrho$ which we would like to study. However, we might not have access to $f$ but only to its approximation $g$. We thus arrive at a fundamental question -- if $f$ and $g$ are close in $L^q$, does $g\circ \varrho$ approximate $f\circ \varrho$ well, and in what sense? Previously, we demonstrated that the answer to this question might be negative in terms of the $L^p$ distance between probability density functions (PDF). Here we show that the Wasserstein metric is the proper framework for this question. For any $p\geq 1$, we bound the Wasserstein distance $W_p (f\circ \varrho , g\circ \varrho) $ from above by $\|f-g\|_{q}$. Furthermore, we provide lower bounds for the cases of $p=1,2$. Finally, we apply our theory to the analysis of common numerical methods in the field of computational uncertainty quantification.

Combining Polynomial Chaos Expansions and the Random Variable Transformation Technique to Approximate the Density Function of Stochastic Problems, Including Some Epidemiological Models

In this paper, we deal with computational uncertainty quantification for stochastic models with one random input parameter. The goal of the paper is twofold: First, to approximate the set of probability density functions of the solution stochastic process, and second, to show the capability of our theoretical findings to deal with some important epidemiological models. The approximations are constructed in terms of a polynomial evaluated at the random input parameter, by means of generalized polynomial chaos expansions and the stochastic Galerkin projection technique. The probability density function of the aforementioned univariate polynomial is computed via the random variable transformation method, by taking into account the domains where the polynomial is strictly monotone. The algebraic/exponential convergence of the Galerkin projections gives rapid convergence of these density functions. The examples are based on fundamental epidemiological models formulated via linear and nonlinear differential and difference equations, where one of the input parameters is assumed to be a random variable.

Improved Efficiency of a Multi-Index FEM for Computational Uncertainty Quantification

We propose a multi-index algorithm for the Monte Carlo (MC) discretization of a linear, elliptic PDE with affine-parametric input. We prove an error versus work analysis which allows a multilevel f...

Uncertainty Quantification for Spectral Fractional Diffusion: Sparsity Analysis of Parametric Solutions

In bounded, polygonal domains $D \subset \Bbb{R}^d$, we analyze solution regularity and sparsity for computational uncertainty quantification for spectral fractional diffusion. Several types of unc...

Computational uncertainty quantification for random non-autonomous second order linear differential equations via adapted gPC: a comparative case study with random Fröbenius method and Monte Carlo simulation

AbstractThis paper presents a methodology to quantify computationally the uncertainty in a class of differential equations often met in Mathematical Physics, namely random non-autonomous second-order linear differential equations, via adaptive generalized Polynomial Chaos (gPC) and the stochastic Galerkin projection technique. Unlike the random Fröbenius method, which can only deal with particular random linear differential equations and needs the random inputs (coefficients and forcing term) to be analytic, adaptive gPC allows approximating the expectation and covariance of the solution stochastic process to general random second-order linear differential equations. The random inputs are allowed to functionally depend on random variables that may be independent or dependent, both absolutely continuous or discrete with infinitely many point masses. These hypotheses include a wide variety of particular differential equations, which might not be solvable via the random Fröbenius method, in which the random input coefficients may be expressed via a Karhunen-Loève expansion.

Computational uncertainty quantification for some strongly degenerate parabolic convection–diffusion equations

Strongly degenerate parabolic convection–diffusion equations arise as governing equations in a number of applications such as traffic flow with driver reaction and anticipation distance and sedimentation of solid–liquid suspensions in mineral processing and wastewater treatment. In these applications several parameters that define the convective flux function and the degenerating diffusion coefficient are subject to stochastic variability. A method to evaluate the variability of the solution of the governing partial differential equation in response to that of the parameters is presented. To this end, a generalized polynomial chaos (gPC) expansion of the solution is approximated by its projection onto a finite-dimensional space of piecewise polynomial functions defined on a suitable discretization of the stochastic domain, according to the basic principle of the hybrid stochastic Galerkin (HSG) approach. This approach is combined with a finite volume (FV) method, resulting in a so-called FV–HSG method, to compute the sought deterministic coefficient functions of the truncated polynomial-chaos-based expansion of the solution. Since the stochastic parameter space is now spanned by piecewise polynomial functions, one may employ the numerical result to compute the reconstruction of the numerical solution for arbitrary values of the random variables. The expectation, the variance or other stochastic quantities of the solution (as functions of time and position) can also be computed from these coefficient functions. The method is illustrated by a number of numerical examples.

Reproducibility of Experiments With Swirling Flow: Numerical Prediction With Polynomial Chaos

Achieving good reproducibility in fluid flow experiments can be challenging, in particular in scenarios where the experimental boundary conditions are obscure. We use computational uncertainty quantification (UQ) to evaluate the influence of uncertain inflow conditions on the reproducibility of experiments with swirling flow. Using a nonintrusive polynomial chaos method in combination with a computational fluid dynamics (CFD) code, we obtain the expectation and variance of the velocity fields downstream from symmetric and asymmetric swirl disturbance generators. Our results suggest that the flow patterns downstream from the asymmetric swirl disturbance generator are more reproducible than the flow patterns downstream from the symmetric swirl disturbance generator. This confirms that the inherent breaking of symmetry eliminates instability mechanisms in the wake of the disturber, thereby creating more stable swirling patterns that make the experiments more reproducible.

Electromagnetic wave scattering by random surfaces: Shape holomorphy

For time-harmonic electromagnetic waves scattered by either perfectly conducting or dielectric bounded obstacles, we show that the fields depend holomorphically on the shape of the scatterer. In the presence of random geometrical perturbations, our results imply strong measurability of the fields, in weighted spaces in the exterior of the scatterer. These findings are key to prove dimension-independent convergence rates of sparse approximation techniques of polynomial chaos type for forward and inverse computational uncertainty quantification. Also, our shape-holomorphy results imply parsimonious approximate representations of the corresponding parametric solution families, which are produced, for example, by greedy strategies such as model order reduction or reduced basis approximations. Finally, the presently proved shape holomorphy results imply convergence of shape Taylor expansions of far-field patterns for fixed amplitude domain perturbations in a vicinity of the nominal domain, thereby extending the widely used asymptotic linearizations employed in first-order, second moment domain uncertainty quantification.

Computational uncertainty quantification for a clarifier-thickener model with several random perturbations: A hybrid stochastic Galerkin approach

Continuous sedimentation processes in a clarifier-thickener unit can be described by a scalar nonlinear conservation law whose flux density function is discontinuous with respect to the spatial position. In the applications of this model, which include mineral processing and wastewater treatment, the rate and composition of the feed flow cannot be given deterministically. Efficient numerical simulation is required to quantify the effect of uncertainty in these control parameters in terms of the response of the clarifier-thickener system. Thus, the problem at hand is one of uncertainty quantification for nonlinear hyperbolic problems with several random perturbations. The presented hybrid stochastic Galerkin method is devised so as to extend the polynomial chaos approximation by multiresolution discretization in the stochastic space. This approach leads to a deterministic hyperbolic system, which is partially decoupled and therefore suitable for efficient parallelisation. Stochastic adaptivity reduces the computational effort. Several numerical experiments are presented.

Numerical Solution of Scalar Conservation Laws with Random Flux Functions

We consider scalar hyperbolic conservation laws in several space dimensions, with a class of random (and parametric) flux functions. We propose a Karhunen--Loève expansion on the state space of the random flux. For random flux functions which are continuously differentiable with respect to the state variable $u$, we prove the existence of a unique random entropy solution. Using a Karhunen--Loève spectral decomposition of the random flux into principal components with respect to the state variables, we introduce a family of parametric, deterministic entropy solutions on high-dimensional parameter spaces. We prove bounds on the sensitivity of the parametric and of the random entropy solutions on the Karhunen--Loève parameters. We also outline the convergence analysis for two classes of discretization schemes, the multilevel Monte Carlo finite volume method (MLMCFVM), developed in [S. Mishra and C. Schwab, Math. Comp., 81 (2012), pp. 1979--2018], [S. Mishra, C. Schwab, and J. Šukys, J. Comput. Phys., 231 (2012), pp. 3365--3388], and [S. Mishra, C. Schwab, and J. Šukys, Multi-level Monte Carlo finite volume methods for uncertainty quantification in nonlinear systems of balance laws, in Uncertainty Quantification in Computational Fluid Dynamics, Lecture Notes in Comput. Sci. Eng. 92, Springer, Heidelberg, 2013, pp. 225--294], and the stochastic collocation finite volume method (SCFVM) of [S. Tokareva, Stochastic Finite Volume Methods for Computational Uncertainty Quantification in Hyperbolic Conservation Laws, Ph.D. dissertation, ETH Diss. Nr. 21498, ETH Zürich, Zürich, Switzerland]. (A corrected version is attached.)

Adaptive finite difference solutions of Liouville equations in computational uncertainty quantification

In the context of computational uncertainty quantification, we present a finite difference based solution of a Liouville equation that describes the evolution of a conditional probability density function depending on state variables and model input parameters. The efficiency of the numerical solution is enhanced through (i) a quadrature-based sampling of random variables corresponding to model input parameters and (ii) time-adaptive methods for determining the computational grid in the space of state or response variables. The proposed method, which allows for accurate and computationally efficient determination of the conditional density and resulting statistical moments, is applied to a decay problem and a problem which contains nonlinear deterministic dynamics with multiple equilibrium points. Through comparison with exact solutions and direct numerical simulation on a fixed time-invariant grid, it is demonstrated that the proposed methodology not only is able to accurately predict the large time behavior of the response density but also can significantly reduce the computational costs.

On convergence of moments in uncertainty quantification based on direct quadrature

Theoretical results for the convergence of statistical moments in numerical quadrature based polynomial chaos computational uncertainty quantification are presented in this work. This is accomplished by considering the computation of the moments through a direct numerical quadrature method, which is shown to be equivalent to stochastic collocation. For problems which involve output variables which have a polynomial dependence on the random input variables, lower bound expressions are derived for the number of quadrature points required for convergence of arbitrary order moments. In addition, an error expression is derived for when this lower bound is used for problems which have a higher degree of continuity than what was assumed when the bounds are computed. The theoretical results are demonstrated through a simple random algebraic problem and a nonlinear plate problem. The results presented in this work provide further insight into the widely used polynomial chaos expansion method of uncertainty quantification along with presenting simple expressions which can be used for uncertainty quantification code verification.