Seminars (UNQ)

From both mathematical and statistical perspectives, the fundamental goal of Uncertainty Quantification (UQ) is to ascertain uncertainties inherent to parameters, initial and boundary conditions, experimental data, and models themselves to make predictions with improved and quantified accuracy.  Some factors that motivate recent developments in mathematical UQ analysis include the following.  The first is the goal of quantifying uncertainties for models and applications whose complexity precludes sole reliance on sampling-based methods.  This includes simulation codes for discretized partial differential equation (PDE) models, which can require hours to days to run.  Secondly, models are typically nonlinearly parameterized thus requiring nonlinear statistical analysis.  Finally, there is often emphasis on extrapolatory or out-of-data predictions; e.g., using time-dependent models to predict future events.  This requires embedding statistical models within physical laws, such as conservation relations, to provide the structure required for extrapolatory predictions.  Within this context, the discussion will focus on techniques to isolate subsets and subspaces of inputs that are uniquely determined by data.  We will also discuss the use of stochastic collocation and Gaussian process techniques to construct and verify surrogate models, which can be used for Bayesian inference and subsequent uncertainty propagation to construct prediction intervals for statistical quantities of interest.  The presentation will conclude with discussion pertaining to the quantification of model discrepancies in a manner that preserves physical structures.

In this talk, I will present some recent works of EDF R&D about the uncertainty management of cpu-expensive computer codes (also called numerical models). I will focus on the following topics: visualization tools in uncertainty and sensitivity analysis, functional risk curve estimation via a metamodel-based methodology, sensitivity analysis with dependent inputs using Shapley effects and robustness analysis of quantile estimation by the way of a recent technique based on perturbation of inputs' probability distributions.

In this talk I present recent work in philosophy and decision theory on the classification, measurement and management on uncertainty.

We consider the determination of statistical information about outputs of interest that depend on the solution of a partial differential equation having random inputs, e.g., coefficients, boundary data, source term, etc. Monte Carlo methods are the most used approach used for this purpose. We discuss other approaches that, in some settings, incur far less computational costs. These include quasi-Monte Carlo, polynomial chaos, stochastic collocation, compressed sensing, reduced-order modeling, and multi-level and multi-fidelity methods for all of which we also discuss their relative strengths and weaknesses.

Co-authors: Ken Carslaw (University of Leeds), Carly Reddington (University of Leeds), Kirsty Pringle (University of Leeds), Graham Mann (University of Leeds), Oliver Wild (University of Lancaster), Edmund Ryan (University of Lancaster), Philip Stier (University of Oxford), Duncan Watson-Parris (University of Oxford)

I will introduce some of the applications of UQ in earth sciences and the challenges remaining that could be addressed during the programme. Earth science models are 3-d dynamic models whose CPU demands and data storage often limits the sample size for UQ. We often choose to use averages of the data and dimension reduction to carry out UQ but it is not always clear that the uncertainty quantified is the most useful for uncertainty reduction or increasing confidence in prediction. I will ask whether we should be applying the same techniques to understand and improve the model as those used to reduce uncertainty in predictions showing some examples where the end goal is different. I will look at UQ when constraint or calibration is the goal and how we incorporate uncertainty and use ‘real’ data. This will also raise the question of identifiability in our uncertainty quantification and how to deal with and accurately quantify irreducible uncertainty. Finally, I would like to discuss how we validate our methods in a meaningful way.

Demography, with its over 350 years of history, is renowned for its empiricism, firm links with statistics, and the lack of a strong theoretical background. The uncertainty of demographic processes is well acknowledged, and over the past three decades, methods have been designed for quantifying the uncertainty in population estimates and forecasts. In parallel, developments in model-based demographic simulation methods, such as agent-based modelling, have recently offered a promise of shedding some light on the complexity of population processes. However, the existing approaches are still far from fulfilling this promise, and are themselves fraught with epistemological pitfalls. Crucially, complex problems are not analytically tractable with the use of traditional methods. In this talk, I will discuss the potential of uncertainty quantification in bringing together the empirical data, statistical inference and computer simulations, with insights into behavioural and social theory and knowledge of social mechanisms. The discussion will be illustrated by an example of a new research programme on model-based migration studies.

In this talk, numerical computation - such as numerical solution of a PDE – will be treated as an inverse problem in its own right. The popular Bayesian approach to inversion is considered, wherein a posterior distribution is induced over the object of interest by conditioning a prior distribution on the same finite information that would be used in the classical numerical context. This approach to numerical computation opens the door to application of statistical techniques, and we discuss the relevance of decision theory and probabilistic graphical models in particular detail. The concepts will be briefly illustrated through numerical experiment.

The first part of this talk gives an introduction to low-rank tensors as a general tool for approximating high-dimensional functions. The second part deals with their application to partial differential equations with many parameters, as they arise in particular in uncertainty quantification problems. Here the focus is on error and complexity guarantees, as well as on low-rank approximability of solutions.

Normal and regular beating of the human heart is essential to maintain life. In each beat, as wave of electrical excitation arises in the heart's natural pacemaker, and spreads throughout the rest of the heart. This wave acts as a signal to initialise and synchronise mechanical contraction of the heart tissue, which in turn generates pressure in the chambers of the heart and acts to propel blood around the body. Models have been developed for the electrical and mechanical behaviour of the heart, as well for blood flow. In this talk I will concentrate on models of electrical activation because failures in the initiation and normal propagation of electrical activation can result in a disruption of normal mechanical behaviour, and underlie a range of common heart problems. Models of electrical activation in both single cells and in tissue are stiff, nonlinear, and have a large number of parameters. Until recently there has been little interest in how uncertainty in model parameters and other inputs influences model behaviour. However, the prospect of using these models for critical applications including drug safety testing and guiding interventions in patients has begun to stimulate activity in this area.

Approximate Bayesian computation (ABC) methods are widely used in some scientific disciplines for fitting stochastic simulators to data. They are primarily used in situations where the likelihood function of the simulator is unknown, but where it is possible to easily sample from the simulator. Methodological development of ABC methods has primarily focused on computational efficiency and tractability, rather than on careful uncertainty modelling. In this talk I'll briefly introduce ABC and its various extensions, and then interpret it from a UQ perspective and suggest how it may best be modified.

Careful assessment of model discrepancy is a crucial aspect of uncertainty quantification. We will discuss the different ways in which emulation may be used to support such assessment, illustrating with practical examples.

We address the numerical analysis of domain uncertainty in UQ for partial differential and integral equations. For small amplitude shape variation, a first order, kth moment perturbation analysis and sparse tensor discretization produces approximate k-point correlations at near optimal order: work and memory scale log-linearly w.r. to N, the number of degrees of freedom for approximating one instance of the nominal (mean-field) problem [1,3]. For large domain variations, the notion of shape holomorphy of the solution is introduced. It implies (the `usual') sparsity and dimension-independent convergence rates of gpc approximations (e.g., anisotropic stochastic collocation, least squares, CS, ...) of parametric domain-to-solution maps in forward UQ. This property holds for a broad class of smooth elliptic and parabolic boundary value problems. Shape holomorphy also implies sparsity of gpc expansions of certain posteriors in Bayesian inverse UQ [7], [->WS4]. We discuss consequences of gpc sparsity on some surrogate forward models, to be used e.g. in optimization under domain uncertainty [8,9]. We also report on dimension independent convergence rates of Smolyak and higher order Quasi-Monte Carlo integration [5,6,7]. Examples include the usual (anisotropic) diffusion problems, Navier-Stokes [2] and time harmonic Maxwell PDEs [4], and forward UQ for fractional PDEs. Joint work with Jakob Zech (ETH), Albert Cohen (Univ. P. et M. Curie), Carlos Jerez-Hanckes (PUC, Santiago, Chile). Work supported in part by the Swiss National Science Foundation. References: [1] A. Chernov and Ch. Schwab: First order k-th moment finite element analysis of nonlinear operator equations with stochastic data, Mathematics of Computation, 82 (2013), pp. 1859-1888. [2] A. Cohen and Ch. Schwab and J. Zech: Shape Holomorphy of the stationary Navier-Stokes Equations, accepted (2018), SIAM J. Math. Analysis, SAM Report 2016-45. [3] H. Harbrecht and R. Schneider and Ch. Schwab: Sparse Second Moment Analysis for Elliptic Problems in Stochastic Domains, Numerische Mathematik, 109/3 (2008), pp. 385-414. [4] C. Jerez-Hanckes and Ch. Schwab and J. Zech: Electromagnetic Wave Scattering by Random Surfaces: Shape Holomorphy, Math. Mod. Meth. Appl. Sci., 27/12 (2017), pp. 2229-2259. [5] J. Dick and Q. T. Le Gia and Ch. Schwab: Higher order Quasi Monte Carlo integration for holomorphic, parametric operator equations. SIAM Journ. Uncertainty Quantification, 4/1 (2016), pp. 48-79. [6] J. Zech and Ch. Schwab: Convergence rates of high dimensional Smolyak quadrature. In review, SAM Report 2017-27. [7] J. Dick and R. N. Gantner and Q. T. Le Gia and Ch. Schwab: Multilevel higher-order quasi-Monte Carlo Bayesian estimation. Math. Mod. Meth. Appl. Sci., 27/5 (2017), pp. 953-995. [8] P. Chen and Ch. Schwab: Sparse-grid, reduced-basis Bayesian inversion: Nonaffine-parametric nonlinear equations. Journal of Computational Physics, 316 (2016), pp. 470-503. [9] Ch. Schwab and J. Zech: Deep Learning in High Dimension. In review, SAM Report 2017-57.

Drawing meaningful conclusions on the way complex real life phenomena work and being able to predict the behavior of systems of interest require developing accurate and highly interpretable mathematical models whose parameters need to be estimated from observations. In modern applications, however, we are often challenged with the lack of such models, and even when these are available they are too computational demanding to be suitable for standard parameter optimization/inference methods. While probabilistic models based on Deep Gaussian Processes (DGPs) offer attractive tools to tackle these challenges in a principled way and to allow for a sound quantification of uncertainty, carrying out inference for these models poses huge computational challenges that arguably hinder their wide adoption. In this talk, I will present our contribution to the development of practical and scalable inference for DGPs, which can exploit distributed and GPU computing. In particular, I will introduce a formulation of DGPs based on random features that we infer using stochastic variational inference. Through a series of experiments, I will illustrate how our proposal enables scalable deep probabilistic nonparametric modeling and significantly advances the state-of-the-art on inference methods for DGPs.

Co-authors: Chris Pope (University of Leeds), Jill Johnson (University of Leeds), Stuart Barber (University of Leeds), Paul Blackwell (University of Sheffield)

Computationally expensive, complex computer programs are often used to model and predict real-world phenomena. The standard Gaussian process model has a drawback in that the computer code output is assumed to be homogeneous over the input space. Computer codes can behave very differently in various regions of the input space. Here, we introduce a piecewise Gaussian process model to deal with this problem where the input space is divided into separate regions using Voronoi tessellations (also known as Dirichlet tessellations, Thiessen polygons or the dual of the Delaunay triangulation). We demonstrate our method’s utility with an application in climate science.

Co-authors: Max Pfeffer (MPI MIS Leipzig), Manuel Marschall (WIAS Berlin), Reinhold Schneider (TU Berlin)

The Stochastic Galerkin Finite Element Method (SGFEM) is a common approach to numerically solve random PDEs with the aim to obtain a functional representation of the stochastic solution. As with any spectral method, the curse of dimensionality renders the approach challenging when the randomness depends on a large or countable infinite set of parameters. This makes function space adaptation and model reduction strategies a necessity. We review adaptive SGFEM based on reliable a posteriori error estimators for affine and non-affine parametric representations. Based on this, an adaptive explicit sampling-free Bayesian inversion in hierarchical tensor formats can be derived. As an outlook onto current research, a statistical learning viewpoint is presented, which connects concepts of UQ and machine learning from a Variational Monte Carlo perspective.

Co-author: Antoine Espinet (Cornell University)

This talk will describe general-purpose algorithms for global optimization These algorithms can be used to estimate model parameters to fit complex simulation models to data, to select among alternative options for design or management, or to quantify model uncertainty. In general the numerical results indicate these algorithms do very well in comparison to alternatives, including Gaussian Process based approaches.. Prof. Shoemaker’s group has developed open source (free) PySOT optimization software that is available online (18,000 downloads) . The algorithms can be run in serial or parallel. The focus of the talk will be on SOARS, an Uncertainty Quantification method for using optimization-based sampling to build a surrogate likelihood function followed by additional sampling The algorithms builds a surrogate approximation of the likelihood function based on simulations done during the optimization search. Then MCMC is performed by evaluating the surrogate likelihood function rather than the original expensive-to-evaluate function. Numerical results indicate the SOARS algorithm is very accurate when compared to the posterior densities computed when using the expensive exact likelihood function. I also discuss an application to a model of the underground movement of a plume of geologically sequestered carbon dioxide. The uncertainty in the parameter values obtained from the MCMC analysis on the surrogate likelihood function can be used to assess alternative strategies for identifying a cost-effective plan that will most efficiently give a reliable forecast of a carbon dioxide underground plume. This includes joint work with David Ruppert, Antoine Espinet, Nikolay Bliznyuk, and Yilun Wang.

Related Links

• https://www.isem.nus.edu.sg - NUS ISEM Department website
• https://sites.google.com/site/shoemakernusgroup/home - Shoemaker group page

Co-author: Dave Woods (University of Southampton)

Data collected from correlated processes arise in many diverse application areas including both computer and physical experiments, and studies in environmental science. Often, such data are used for prediction and optimisation of the process under study. For example, we may wish to construct an emulator of a computationally expensive computer model, or simulator, and then use this emulator to find settings of the controllable variables that maximise the predicted response. The design of the experiment from which the data are collected may strongly influence the quality of the model fit and hence the precision and accuracy of subsequent predictions and decisions. We consider Gaussian process models that are typically defined by a correlation structure that may depend upon unknown parameters. This parametric uncertainty may affect the choice of design points, and ideally should be taken into account when choosing a design. We consider a decision-theoretic Bayesian design for Gaussian process models which is usually computationally challenging as it requires the optimization of an analytically intractable expected loss function over high-dimensional design space. We use a new approximation to the expected loss to find decision-theoretic optimal designs. The resulting designs are illustrated through a number of simple examples.

Routine diagnostic checking of stationary Gaussian processes fitted to the output of complex computer codes often reveals nonstationary behaviour. There have been a number of approaches, both traditional and more recent, to modelling or accounting for this nonstationarity via the fitted process. These have included the fitting of complex mean functions to attempt to leave a stationary residual process (an idea that is often very difficult to get right in practice), using regression trees or other techniques to partition the input space into regions where different stationary processes are fitted (leading to arbitrary discontinuities in the modelling of the overall process), and other approaches which can be considered to live in one of these camps. In this work we allow the fitted process to be continuous by modelling the covariance kernel as a finite mixture of stationary covariance kernels and allowing the mixture weights to vary appropriately across parameter space. We introduce our method and compare its performance with the leading approaches in the literature for a variety of standard test functions and the cloud parameterisation of the French climate model. This is work led by my final-year PhD student, Victoria Volodina.

We investigate the merits of replication, and provide methods that search for optimal designs (including replicates), in the context of noisy computer simulation experiments. We first show that replication offers the potential to be beneficial from both design and computational perspectives, in the context of Gaussian process surrogate modeling. We then develop a lookahead based sequential design scheme that can determine if a new run should be at an existing input location (i.e., replicate) or at a new one (explore). When paired with a newly developed heteroskedastic Gaussian process model, our dynamic design scheme facilitates learning of signal and noise relationships which can vary throughout the input space. We show that it does so efficiently, on both computational and statistical grounds. In addition to illustrative synthetic examples, we demonstrate performance on two challenging real-data simulation experiments, from inventory management and epidemiology.

This presentation is joint work of:   Ivo Colombo, Dipartimento di Ingegneria Civile e Ambientale, Politecnico di Milano, Italy Fabio Nobile, CSQI-MATHICSE, Ecole Polytechnique Fédérale de Lausanne, Switzerland Giovanni Porta, Dipartimento di Ingegneria Civile e Ambientale, Politecnico di Milano, Italy Anna Scotti, MOX, Dipartimento di Matematica, Politecnico di Milano, Italy Lorenzo Tamellini, CNR - Istituto di Matematica Applicata e Tecnologie Informatiche “E. Magenes”, Pavia, Italy     In this work we propose an Uncertainty Quantification methodology for the evolution of sedimentary basins undergoing mechanical and geochemical compaction processes, which we model as a coupled, time-dependent, non-linear, monodimensional (depth-only) system of PDEs with uncertain parameters.   Specifically, we consider multi-layered basins, in which each layer is characterized by a different material. The multi-layered structure gives rise to discontinuities in the dependence of the state variables on the uncertain parameters. Because of these discontinuites, an appropriate treatment is needed for surrogate modeling techniques such as sparse grids to be effective.   To this end, we propose a two-steps methodology which relies on a change of coordinate system to align the discontinuities of the target function within the random parameter space. Once this alignement has been computed, a standard sparse grid approximation of the state variables can be performed. The effectiveness of this procedure is due to the fact that the physical locations of the interfaces among layers feature a smooth dependence on the random parameters and are therefore amenable to sparse grid polynomial approximations.   We showcase the capabilities of our numerical methodologies through some synthetic test cases.  

In this talk I will emphasize on the use of adaptive strategies in numerical algorithms to solve systems of ordinary and partial differential equations more efficiently and reliably. After a brief introduction to local and global error control for time integrators general approaches to combine adaptivity in space and time are discussed. Finally, I will speak about recent developments in using adaptive multilevel strategies for PDE-constrained optimization and uncertainty quantification. Throughout my talk I will present numerical results for academic as well as real-life applications including chemical reaction-diffusion systems, regional hyperthermia, electro-cardiology, magneto-quasistatics, glass cooling and complex turbulent flows.

Uncertainty Quantification (UQ) is established as an aspect of model-supported inference in scientific and engineering systems. With this expectation comes the desire for addressing increasingly complex applications that challenge the ability of UQ to scale. This talk will describe the view of UQ as a full-system modeling and analysis framework for scientific and engineering models, motivated by the example application of Carbon Capture technology development, and some of the challenging issues that have come up in this multi-level component to full-system modeling effort are discussed, as well as strategies for addressing some of these challenges. Another challenge area is dealing with multivariate responses, and some implications and challenges in this area will be discussed.

History Matching is a calibration technique that systematically reduces the input space in a numerical model. At every iteration, an implausibility measure discards combinations of input values that are unlikely to provide a match between model output and experimental observations. As the input space reduces, sampling becomes increasingly challenging due to the size of the relative volume of the non-implausible space and the fact that it can exhibit a complex, disconnected geometry. Since realistic numerical models are computationally expensive, surrogate models and dimensionality reduction are commonly employed. In this talk we will explore how Subset Simulation, a Markov chain Monte Carlo technique from engineering reliability analysis, can solve the sampling problem in History Matching. We will also explore alternative implausibility measures that can guide the selection of regions of the non-implausible in order to balance sampling exploration and exploitation.

In order to predict the behaviour of complex engineering systems (nuclear power plants, aircraft, infrastructure networks, etc.), analysts nowadays develop high-fidelity computational models that try and capture detailed physics. A single run of such simulators can take minutes to hours even on the most advanced HPC architectures. In the context of uncertainty quantification, methods based on Monte Carlo simulation are simply not affordable. This has led to the rapid development of surrogate modelling techniques in the last decade, e.g. polynomial chaos expansions, low-rank tensor representations, Kriging (a.k.a Gaussian process models) among others. Surrogate models have proven remarkable efficiency in the case of moderate dimensionality (e.g. tens to a hundred of inputs). In the case of high-dimensional problems (hundreds to thousands of inputs), or when the input is cast as time series, 2D maps, etc., the classical set-up of surrogate modelling does not apply straightforwardly. Usually, a pre-processing of the data is carried out to reduce this dimensionality, before a surrogate is constructed. In this talk, we show that the sequential use of compression algorithms (for dimensionality reduction (DR), e.g. kernel principal component analysis) and surrogate modelling (SM) is suboptimal. Instead, we propose a new general-purpose framework that cast the two sub-problems into a single DRSM optimization. In this set-up, the parameters of the DR step are selected so that as to maximize the quality of the subsequent surrogate model. The framework is versatile in the sense that the techniques used for DR and for SM can be freely selected and combined. Moreover, the method is purely data-driven. The proposed approach is illustrated on different engineering problems including 1D and 2D elliptical SPDEs and earthquake engineering applications.

Uncertainty quantification (UQ) is a fast growing research area which deals with the impact of parameter, data and model uncertainties in complex systems. We focus on models which are based on partial differential equations (PDEs) with random inputs. For deterministic PDEs there are many classical analytical results and numerical tools available. The treatment of PDEs with random inputs, however, requires novel ideas and tools. We illustrate the mathematical and algorithmic challenges of UQ for Bayesian inverse problems arising from geotechnical engineering and medicine, and an optimal control problem with uncertain PDE constraints.

We investigate the problem of computing a nested expectation of the form P[E[X|Y] >= 0] = E[H(E[X|Y])] where H is the Heaviside function. This nested expectation appears, for example, when estimating the probability of a large loss from a financial portfolio. We present a method that combines the idea of using Multilevel Monte Carlo (MLMC) for nested expectations with the idea of adaptively selecting the number of samples in the approximation of the inner expectation, as proposed by (Broadie et al., 2011). We propose and analyse an algorithm that adaptively selects the number of inner samples on each MLMC level and prove that the resulting MLMC method with adaptive sampling has an order e^-2|log(e)|^2 complexity to achieve a root mean-squared error e. The theoretical analysis is verified by numerical experiments on a simple model problem. Joint work with: Michael B. Giles (University of Oxford)

Co-authors: Sebastian Krämer, Christian Löbbert, Dieter Moser (RWTH Aachen) We introduce the concept of hierarchical low rank decompositions and approximations by use of the hierarchical Tucker format. In order to relate it to several other existing low rank formats we highlight differences, similarities as well as bottlenecks of these. One particularly difficult question is whether or not a tensor or multivariate function allows a priori a low rank representation or approximation. This question can be related to simple matrix decompositions or approximations, but still the question is not easy to answer, cf. the talks of Wolfgang Dahmen, Sergey Dolgov, Martin Stoll and Anthony Nouy. We provide numerical evidence for a model problem that the approximation can be efficient in terms of a small rank. In order to find such a decomposition or approximation we consider black box (cross) type non-intrusive sampling approaches. A special emphasis will be on postprocessing of the tensors, e.g. finding extremal points efficiently. This is of special interest in the context of model reduction and reliability analysis.

Co-authors: Robert Scheichl (University of Bath)

We consider the approximate solution of parametric PDEs using the low-rank Tensor Train (TT) decomposition. Such parametric PDEs arise for example in uncertainty quantification problems in engineering applications. We propose an algorithm that is a hybrid of the alternating least squares and the TT cross methods. It computes a TT approximation of the whole solution, which is particularly beneficial when multiple quantities of interest are sought. The new algorithm exploits and preserves the block diagonal structure of the discretized operator in stochastic collocation schemes. This disentangles computations of the spatial and parametric degrees of freedom in the TT representation. In particular, it only requires solving independent PDEs at a few parameter values, thus allowing the use of existing high performance PDE solvers. We benchmark the new algorithm on the stochastic diffusion equation against quasi-Monte Carlo and dimension-adaptive sparse grids methods. For sufficiently smooth random fields the new approach is orders of magnitude faster.

In this talk we discuss a weighted approach for the reduction of parametrized partial differential equations with random input, focusing in particular on weighted approaches based on reduced basis (wRB) [Chen et al., SIAM Journal on Numerical Analysis (2013); D. Torlo et al., submitted (2017)] and proper orthogonal decomposition (wPOD) [L. Venturi et al., submitted (2017)]. We will first present the wPOD approach. A first topic of discussion is related to the choice of samples and respective weights according to a quadrature formula. Moreover, to reduce the computational effort in the offline stage of wPOD, we will employ Smolyak quadrature rules. We will then introduce the wRB method for advection diffusion problems with dominant convection with random input parameters. The issue of the stabilisation of the resulting reduced order model will be discussed in detail. This work is in collaboration with Francesco Ballarin (SISSA, Trieste), Davide Torlo (University of Zurich), and Luca Venturi (Courant Institute for Mathematical Sciences, NYC)

Co-authors: Peter BINEV (University of South Carolina), Albert COHEN (University Pierre et Marie Curie), James NICHOLS (University Pierre et Marie Curie)

Parametric PDEs of the general form
$$\mathcal{P} (u,a) = 0$$
are commonly used to describe many physical processes, where $\cal P$ is a differential operator, $a$ is a high-dimensional vector of parameters and $u$ is the unknown solution belonging to some Hilbert space $V$. A typical scenario in state estimation is the following: for an unknown parameter $a$, one observes $m$ independent linear measurements of $u(a)$ of the form $\ell_i(u) = (w_i, u), i = 1, ..., m$, where $\ell_i \in V'$ and $w_i$ are the Riesz representers, and we write $W_m = \text{span}\{w_1,...,w_m\}$. The goal is to recover an approximation $u^*$ of $u$ from the measurements. Due to the dependence on a the solutions of the PDE lie in a manifold and the particular PDE structure often allows to derive good approximations of it by linear spaces Vn of moderate dimension n. In this setting, the observed measurements and Vn can be combined to produce an approximation $u^*$ of $u$ up to accuracy
$$
\Vert u -u^* \Vert \leq \beta(V_n, W_m) \text{dist}(u, V_n)
$$
where
$$
\beta(V_n, W_m) := \inf_{v\in V_n} \frac{\Vert P_{W_m} v \Vert}{\Vert v \Vert}
$$
plays the role of a stability constant. For a given $V_n$, one relevant objective is to guarantee that $\beta(V_n, W_m) \geq \gamma >0$ with a number of measurements $m \geq n$ as small as possible. We present results in this direction when the measurement functionals $\ell_i$ belong to a complete dictionary. If time permits, we will also briefly explain ongoing research on how to adapt the reconstruction technique to noisy measurements.

Related Links

• https://hal.archives-ouvertes.fr/hal-01638177/document - Preprint

We present two reduced-order model based approaches  for the efficient and accurate evaluation of the Conditional-Value-at-Risk  (CVaR) of quantities of interest (QoI) in engineering systems with uncertain parameters.  CVaR is used to model objective or constraint functions in risk-averse engineering design and optimization applications under uncertainty.  Estimating the CVaR of the QoI is expensive. While the distribution of the uncertain system parameters is known, the resulting QoI is a random variable that is implicitly determined via the state of the system. Evaluating the CVaR of the QoI requires  sampling in the tail of the QoI distribution and typically requires  many solutions of an expensive full-order model of the engineering system. Our reduced-order model approaches substantially reduce this computational expense.

Ice sheet models are currently unable to reproduce the retreat of the North American ice sheet through the last deglaciation, due to the large uncertainty in the boundary conditions. To successfully calibrate such a model, it is important to vary both the input parameters and the boundary conditions. These boundary conditions are derived from global climate model simulations, and hence the biases from the output of these models are carried through to the ice sheet output, restricting the range of ice sheet output that is possible. Due to the expense of running global climate models for the required 21,000 years, there are only a small number of such runs available; hence it is difficult to quantify the boundary condition uncertainty. We develop a methodology for generating a range of plausible boundary conditions, using a low-dimensional basis representation for the spatio-temporal input required. We derive this basis by combining key patterns, extracted from a small climate model ensemble of runs through the deglaciation, with sparse spatio-temporal observations. Varying the coefficients for the chosen basis vectors and ice sheet parameters simultaneously, we run ensembles of the ice sheet model. By emulating the ice sheet output, we history match iteratively and rule out combinations of the ice sheet parameters and boundary condition coefficients that lead to implausible deglaciations, reducing the uncertainty due to the boundary conditions.

In the applied mathematics community, reduced basis methods are typically used to reduce the computational cost of applying sampling methods to parameter-dependent partial differential equations (PDEs). When dealing with PDE models in particular, repeatedly running computer models (eg finite element solvers) for many choices of the input parameters, is computationally infeasible. The cost of obtaining each sample of the numerical solution is instead sought by projecting the so-called high fidelity problem into a reduced (lower-dimensional) space. The choice of reduced space is crucial in balancing cost and overall accuracy. In this talk, we do not consider sampling methods. Rather, we consider stochastic Galerkin finite element methods (SGFEMs) for parameter-dependent PDEs. Here, the idea is to approximate the solution to the PDE model as a function of the input parameters. We combine finite element approximation in physical space, with global polynomial approximation on the parameter domain. In the statistics community, the term intrusive polynomial chaos approximation is often used. Unlike samping methods, which require the solution of many deterministic problems, SGFEMs yield a single very large linear system of equations with coefficient matrices that have a characteristic Kronecker product structure. By reformulating the systems as multiterm linear matrix equations, we have developed [see: C.E. Powell, D. Silvester, V.Simoncini, An efficient reduced basis solver for stochastic Galerkin matrix equations, SIAM J. Comp. Sci. 39(1), pp A141-A163 (2017)] a memory-efficient solution algorithm which generalizes ideas from rational Krylov subspace approximation (which are known in the linear algebra community). The new approach determines a low-rank approximation to the solution matrix by performing a projection onto a reduced space that is iteratively augmented with problem-specific basis vectors. Crucially, it requires far less memory than standard iterative methods applied to the Kronecker formulation of the linear systems. For test problems consisting of elliptic PDEs, and indefinite problems with saddle point structure, we are able to solve systems of billions of equations on a standard desktop computer quickly and efficiently.

We present an extension of principal component analysis for functions of multiple random variables and an associated algorithm for the approximation of such functions using tree-based low-rank formats (tree tensor networks). A multivariate function is here considered as an element of a Hilbert tensor space of functions defined on a product set equipped with a probability measure. The algorithm only requires evaluations of functions on a structured set of points which is constructed adaptively. The algorithm constructs a hierarchy of subspaces associated with the different nodes of a dimension partition tree and a corresponding hierarchy of projection operators, based on interpolation or least-squares projection. Optimal subspaces are estimated using empirical principal component analysis of interpolations of partial random evaluations of the function. The algorithm is able to provide an approximation in any tree-based format with either a prescribed rank or a prescribed relative error, with a number of evaluations of the order of the storage complexity of the approximation format.

Bayesian optimal design is considered for experiments where the response distribution depends on the solution to a system of non-linear ordinary differential equations. The motivation is an experiment to estimate parameters in the equations governing the transport of amino acids through cell membranes in human placentas. Decision-theoretic Bayesian design of experiments for such nonlinear models is conceptually very attractive, allowing the formal incorporation of prior knowledge to overcome the parameter dependence of frequentist design and being less reliant on asymptotic approximations. However, the necessary approximation and maximization of the, typically analytically intractable, expected utility results in a computationally challenging problem. These issues are further exacerbated if the solution to the differential equations is not available in closed-form. This paper proposes a new combination of a probabilistic solution to the equations embedded within a Monte Carlo approximation to the expected utility with cyclic descent of a smooth approximation to find the optimal design. A novel precomputation algorithm reduces the computational burden, making the search for an optimal design feasible for bigger problems. The methods are demonstrated by finding new designs for a number of common models derived from differential equations, and by providing optimal designs for the placenta experiment.

Joint work with Antony Overstall and Ben Parker (University of Southampton)

In 1960 Rudolph Kalman published what is arguably the first paper to develop a systematic, principled approach to the use of data to improve the predictive capability of  the mathematical models developed to understand the world around us. As our ability to gather data grows at an enormous rate, the importance of this work continues to grow too. The lecture will describe this paper and developments that have stemmed from it, revolutionizing fields such space-craft control, weather prediction, oceanography, oil recovery, medical imaging and artificial intelligence. Some mathematical details will be also provided, but limited to simple concepts such as optimization and iteration; the talk is designed to be broadly accessible to anyone with an interest in quantitative science.

Many methods in uncertainty quantification suffer from the curse of dimensionality. I will discuss several approaches for identifying exploitable low-dimensional structure---e.g., active subspaces or likelihood-informed subspaces---that enable otherwise infeasible forward and inverse uncertainty quantification.

Related Links

• http://activesubspaces.org - Active subspaces

We discuss recent results in mathematical statistics that provide objective statistical guarantees for Bayesian algorithms in (possibly non-linear) noisy inverse problems. We focus in particular on the justification of Bayesian credible sets as proper frequentist confidence sets in the small noise limit via so-called `Bernstein - von Mises theorems', which provide Gaussian approximations to the posterior distribution, and introduce notions of such theorems in the infinite-dimensional settings relevant for inverse problems. We discuss in detail such a Bernstein-von Mises result for Bayesian inference on the unknown potential in the Schroedinger equation from an observation of the solution of that PDE corrupted by additive Gaussian white noise. See https://arxiv.org/abs/1707.01764 and also https://arxiv.org/abs/1708.06332

I will discuss two complementary efforts in Bayesian optimal experimental design for inverse problems. The first focuses on evaluating an experimental design objective: we describe a new computational approach for ``focused'' optimal Bayesian experimental design with nonlinear models, with the goal of maximizing expected information gain in targeted subsets of model parameters. Our approach considers uncertainty in the full set of model parameters, but employs a design objective that can exploit learning trade-offs among different parameter subsets. We introduce a layered multiple importance sampling scheme that provides consistent estimates of expected information gain in this focused setting, with significant reductions in estimator bias and variance for a given computational effort. The second effort focuses on optimization of information theoretic design objectives---in particular, from the combinatorial perspective of observation selection. Given many potential experiments, one may wish to choose a most informative subset thereof. Even if the data have in principle been collected, practical constraints on storage, communication, and computational costs may limit the number of observations that one wishes to employ. We introduce methods for selecting near-optimal subsets of the data under cardinality constraints. Our methods exploit the structure of linear inverse problems in the Bayesian setting, and can be efficiently implemented using low-rank approximations and greedy strategies based on modular bounds. This is joint work with Chi Feng and Jayanth Jagalur-Mohan.

Smooth supersaturated polynomials (Bates et al., 2014) have been proposed for emulating computer experiments. These models are simple to interpret and have spline-like properties inherited from the Sobolev-type of smoothing which is at the core of the method. An implementation of these models is available in the R package ssm. This talk aims to describe the method as well as discuss designs that could be performed with the help of smooth models. To illustrate the methodology, we use data from the fan blade assembly. This is joint work with H Wynn, R Bates and P Curtis.

Many problems in machine learning require the classification of high dimensional data. One methodology to approach such problems is to construct a graph whose vertices are identified with data points, with edges weighted according to some measure of affinity between the data points. Algorithms such as spectral clustering, probit classification and the Bayesian level set method can all be applied in this setting. The goal of the talk is to describe these algorithms for classification, and analyze them in the limit of large data sets. Doing so leads to interesting problems in the calculus of variations, Bayesian inverse problems and in Monte Carlo Markov Chain, all of which will be highlighted in the talk. These limiting problems give insight into the structure of the classification problem, and algorithms for it.    

Collaboration with:  
Andrea Bertozzi (UCLA)
Michael Luo (UCLA)
Kostas Zygalakis (Edinburgh)
https://arxiv.org/abs/1703.08816  
and  
Matt Dunlop (Caltech)
Dejan Slepcev (CMU)
Matt Thorpe (Cambridge)
(forthcoming paper)

Computer experiments using deterministic simulators are often used to supplement physical system experiments. A common problem is that a computer simulator may provide biased output for the physical process due to the simplified physics or biology used in the mathematical model. However, when physical observations are available, it may be possible to use these data to align the simulator output to be close to the true mean response by constructing a bias-corrected predictor (a process called calibration). This talk looks at two aspects of experimental design for prediction of physical system means using a Bayesian calibrated predictor. First, the empirical prediction accuracy over the output space of several different types of combined physical and simulator designs is discussed. In particular, designs constructed using the integrated mean squared prediction error seem to perform well. Secondly, a sequential design methodology for optimizing a physical manufacturing process when there are multiple, competing product objectives is described. The goal is to identify a set of manufacturing conditions each of which leads to outputs on the Pareto Front of the product objectives, i.e. identify manufacturing conditions which cannot be modified to improve all the product objectives simultaneously. A sequential design methodology which maximizes the posterior expected minimax fitness function is used to add data from either the simulator or the manufacturing process. The method is illustrated with an example from an injection molding study. The presentation is based on joint work with Thomas Santner, Erin Leatherman, and Po-Hsu Allen Chen.

Gaussian Process models have been used in a number of problems where an objective function f needs to be studied based on a drastically limited number of evaluations. Global optimization algorithms based on Gaussian Process models have been investigated for several decades, and have become quite popular notably in design of computer experiments. Also, further classes of problems involving the estimation of sets implicitly defined by f, e.g. sets of excursion above a given threshold, have inspired multiple research developments. In this talk, we will give an overview of recent results and challenges pertaining to the estimation of sets under Gaussian Process priors, with a particular interest for to the quantification and the sequential reduction of associated uncertainties. Based on a series of joint works primarily with Dario Azzimonti, François Bachoc, Julien Bect, Mickaël Binois, Clément Chevalier, Ilya Molchanov, Victor Picheny, Yann Richet and Emmanuel Vazquez.

Bayesian calibration and history matching are both well established tools for solving inverse problems: finding model inputs to make model outputs match observed data as closely as possible. I will discuss and compare both, within the context of decision-making. I will discuss the sometimes contentious issue of model discrepancy: how and whether we might account for an imperfect or misspecified model within the inference procedure. I will also present some work on history matching of a high dimensional individual based HIV transmission model (joint work with I. Andrianakis, N. McCreesh, I. Vernon, T. J. McKinley, R. N. Nsubuga, M. Goldstein and R. J. White).

Co-authors: Hoang Tran (Oak Ridge National Laboratory) & Clayton Webster (University of Tennessee & Oak Ridge National Laboratory) We present and analyze a novel sparse polynomial approximation method for the solution of PDEs with stochastic and parametric inputs. Our approach treats the parameterized problem as a problem of joint-sparse signal recovery, i.e., simultaneous reconstruction of a set of sparse signals, sharing a common sparsity pattern, from a countable, possibly infinite, set of measurements. In this setting, the support set of the signal is assumed to be unknown and the measurements may be corrupted by noise. We propose the solution of a linear inverse problem via convex sparse regularization for an approximation to the true signal. Our approach allows for global approximations of the solution over both physical and parametric domains. In addition, we show that the method enjoys the minimal sample complexity requirements common to compressed sensing-based approaches. We then perform extensive numerical experiments on several high-dimensional parameterized elliptic PDE models to demonstrate the recovery properties of the proposed approach.

Co-author: Mike Grosskopf (Los Alamos National Lab)

Model calibration uses outputs from a simulator and field data to build a predictive model for the physical system and to estimate unknown inputs. The conventional approach to model calibration assumes that the observations are continuous outcomes. In many applications this is not the case. The methodology proposed was motivated by an application in modeling photon counts at the Center for Exascale Radiation Transport. There, high performance computing is used for simulating the flow of neutrons through various materials. In this talk, new Bayesian methodology for computer model calibration to handle the count structure of our observed data allows closer fidelity to the experimental system and provides flexibility for identifying different forms of model discrepancy between the simulator and experiment.

Modeling complex environmental phenomena such as rainfall patterns has proven challenging due to the difficulty in capturing heavy-tailed behavior, such as extreme weather, in a meaningful way. Recently, a novel approach to this task has taken the form of so-called stochastic weather generators, which use statistical formulations to emulate the distributional patterns of an environmental process. However, while sampling from such models is usually feasible, they typically do not possess closed-form likelihood functions, rendering the usual approaches to model fitting infeasible. Furthermore, some of these stochastic weather generators are now becoming so complex that even simulating from them can be computationally expensive. We propose and compare two approaches to fitting computationally expensive stochastic weather generators motivated by Approximate Bayesian Computation and Stochastic Simulator Calibration methodologies. The methods are then demonstrated by estimating important parameters of a recent stochastic weather generator model applied to rainfall data from the continental USA.

Bayesian calibration of computer models tunes unknown input parameters by comparing outputs to observations. For model outputs distributed over space, this becomes computationally expensive due to the output size. To overcome this challenge, we employ a basis representations of the model outputs and observations: we match these decompositions to efficiently carry out the calibration. In a second step, we incorporate the nonstationary behavior, in terms of spatial variations of both variance and correlations, into the calibration. We insert two INLA-SPDE parameters into the calibration. A synthetic example and a climate model illustration highlight the benefits of our approach.

In this informal seminar, I will describe attempts to force emulators to have monotonic or convex outputs with respect to some of the input parameters. Various prior set-ups will be discussed including piecewise linear approximations, truncated Gaussian processes and non-linear transformations of Gaussian processes alongside computational methods such as Bayes linear updating and approximate Bayesian computations.

In the first part of the talk, we will introduce spatial Gaussian processes. Spatial Gaussian processes are widely studied from a statistical point of view, and have found applications in many fields, including geostatistics, climate science and computer experiments. Exact inference can be conducted for Gaussian processes, thanks to the Gaussian conditioning theorem. Furthermore, covariance parameters can be estimated, for instance by Maximum Likelihood. In the second part of the talk, we will introduce a class of iterative sampling strategies for Gaussian processes, called 'stepwise uncertainty reduction' (SUR). We will give examples of SUR strategies which are widely applied to computer experiments, for instance for optimization or detection of failure domains. We will provide a general consistency result for SUR strategies, together with applications to the most standard examples.

This talk will describe a flexible Bayesian model that can be used to predict the output of a deterministic simulator code. The model assumes that the output can be described as the sum of a smooth global trend plus deviations from the global trend. The global trend and the local deviations are modeled as draws from independent GPs with separable correlation functions subject to appropriate constraints to enforce smoothness of the global process compared with the local deviation process. The accuracy and limitations of predictions made using this model are demonstrated in a series of examples. The model is used to perform variable selection by identifying the most active inputs to the simulator. Inputs having ``smaller'' posterior distributions of the model's correlation parameters are judged to be more active. A reference inactive input is added to the data to judge the size of the correlation parameter for inactive inputs. Joint work with Casey Davis and Christopher Hans

When solving linear stochastic differential equations numerically, usually a high order spatial discretisation is used. Balanced truncation (BT) is a well-known projection technique in the deterministic framework which reduces the order of a control system and hence reduces computational complexity. We give an introduction to model order reduction (MOR) by BT and then consider a differential equation where the control is replaced by a noise term. We provide theoretical tools such as stochastic concepts for reachability and observability, which are necessary for balancing related MOR of linear stochastic differential equations with additive L'evy noise. Moreover, we derive error bounds for BT and provide numerical results for a specific example which support the theory. This is joint work with Martin Redmann (WIAS Berlin).

During this talk, we will present a numerically efficient method to bound the error that is made when approximating the output of a nonlinear problem depending on an unknown parameter (described by a probability distribution). The class of nonlinear problems under consideration includes high-dimensional nonlinear problems with a nonlinear output function. A goal-oriented probabilistic bound is computed by considering two phases. An offline phase dedicated to the computation of a reduced model during which the full nonlinear problem needs to be solved only a small number of times. The second phase is an online phase which approximates the output. This approach is applied to a toy model and to a nonlinear partial differential equation, more precisely the Burgers equation with unknown initial condition given by two probabilistic parameters. The savings in computational cost are evaluated and presented.

Selecting and performing experiments that produce the most useful data is extremely valuable in engineering and science applications where experiments are costly and resources are limited. Simulation-based experimental design thus provides a rigorous mathematical framework to systematically quantify and maximize the value of experiments while leveraging the existing knowledge and predictive capability of an available model.  We are particularly interested in design settings that accommodate nonlinear and computationally intensive models, such as those governed by ordinary and partial differential equations. Employing principles from Bayesian statistics to characterize and quantify uncertainty, we seek experiments that maximize the expected information gain. Computing these optimal designs using conventional approaches, however, is generally intractable. Major challenges include high dimensional parameter spaces, expensive model simulations, and numerical approximation and optimization of the expected information gain. We thus describe practical numerical methods to help overcome these obstacles, including global sensitivity analysis, surrogate modeling via polynomial chaos, and stochastic optimization.  The overall methodology is demonstrated through the design of combustion experiments for optimal learning of chemical rate parameters, and of configurations for a supersonic jet engine to obtain measurements most informative on turbulent flow parameters.