Active Subspace Methods for Data-Intensive Inverse Problems (Final Report)

The Bayesian approach to inverse problems, in which the posterior probability distribution on an unknown field is sampled for the purposes of computing posterior expectations of quantities of interest, is starting to become computationally feasible for partial differential equation (PDE) inverse problems. Balancing the sources of error arising from finite-dimensional approximation of the unknown field, the PDE forward solution map and the sampling of the probability space under the posterior distribution are essential for the design of efficient computational Bayesian methods for PDE inverse problems. We study Bayesian inversion for a model elliptic PDE with an unknown diffusion coefficient. We provide complexity analyses of several Markov chain Monte Carlo (MCMC) methods for the efficient numerical evaluation of expectations under the Bayesian posterior distribution, given data δ. Particular attention is given to bounds on the overall work required to achieve a prescribed error level ε. Specifically, we first bound the computational complexity of ‘plain’ MCMC, based on combining MCMC sampling with linear complexity multi-level solvers for elliptic PDE. Our (new) work versus accuracy bounds show that the complexity of this approach can be quite prohibitive. Two strategies for reducing the computational complexity are then proposed and analyzed: first, a sparse, parametric and deterministic generalized polynomial chaos (gpc) ‘surrogate’ representation of the forward response map of the PDE over the entire parameter space, and, second, a novel multi-level Markov chain Monte Carlo strategy which utilizes sampling from a multi-level discretization of the posterior and the forward PDE. For both of these strategies, we derive asymptotic bounds on work versus accuracy, and hence asymptotic bounds on the computational complexity of the algorithms. In particular, we provide sufficient conditions on the regularity of the unknown coefficients of the PDE and on the approximation methods used, in order for the accelerations of MCMC resulting from these strategies to lead to complexity reductions over ‘plain’ MCMC algorithms for the Bayesian inversion of PDEs.

In this paper, we develop a non‐iterative way to solve biharmonic boundary value problems by using a radial basis meshless method. This is an original application of meshless method to solving inverse problems without any iteration, since traditional numerical methods for inverse boundary value problems mainly are iterative and hence very time‐consuming. Numerical examples are presented for inverse biharmonic boundary value problems and corresponding direct problems, since solving direct problems is a preliminary step for inverse problems. All our examples of direct and inverse problems are solved within seconds in CPU time on a standard PC, which makes our proposed technique a great potential candidate for wide‐spread applications to other inverse problems. Copyright © 2004 John Wiley & Sons, Ltd.

We study the convergence rates of a non-stationary asymptotical regularization method for linear statistical inverse problems from indirect noisy measurements. The convergence rates are usually derived under some range-type source conditions, which are difficult to be verified rigorously in general, especially for the models governed by some partial differential equations. We shall propose a new strategy associated with conditional stability estimates to obtain the error bounds in terms of mean-squared error. To the best of our knowledge, it is the first time to apply conditional stability estimates for convergence rates in statistical inverse problems. As an application, an inverse source problem in a wave equation is considered and its conditional stability estimate is verified rigorously, and therefore the convergence rates are derived immediately. Several numerical experiments for the inverse source problem are also provided to demonstrate the accuracy and efficiency of the non-stationary asymptotical regularization method.

<p>Multiobservable thermochemical tomography (MTT) is a recent computational approach to obtain estimates of the physical state (e.g. temperature distribution, compositional structure and rock properties) for the upper mantle [1]. It allows to jointly invert multiple independent datasets (e.g. gravity, seismic, magnetotelluric) within a thermodynamically-constrained and fully probabilistic framework. Evaluation of the plausibility of different physical states of the mantle with Markov Chain Monte Carlo (MCMC) simulations requires the solution of complex forward problems (e.g. Stokes flow, Maxwell’s equations, etc.) millions of times, making MTT computationally demanding for large-scale inverse problems. Furthermore, the number of parameters in a global study can easily reach several millions, making it increasingly difficult to 1) locate the regions of high probability and 2) sample these regions appropriately.</p><p>In order to overcome these limitations, we have combined and implemented a number of techniques, such as reduced-order modelling and efficient parallelization of both the forward problems and the MCMC algorithms, which dramatically accelerate the solution of the forward problems. Our software, LitMod1D_4INV and LitMod3D_4INV, allow to compute a proposal in less than 1 second, even when solving multiple complex forward problems together. We develop a multi-level parallel MPI driver for a collection of advanced MCMC sampling strategies to locate and sample high-probability regions efficiently. The massive amounts of data generated by large-scale MTT inversions need to be managed efficiently. We output results to open-source freeware formats, such as HDF5, TileDB, designed for big data problems. We emphasize that our methods and approaches are not only useful for MTT, but for any demanding inverse problem.</p><p>In this contribution, we will present applications of our software to complex, large-scale MTT problems and discuss its benefits, limitations and future improvements.</p><p>[1] J.C. Afonso, N. Rawlinson, Y. Yang, D. L. Schutt, A. G. Jones, J. Fullea, W. L. Griffin, 3‐D multiobservable probabilistic inversion for the compositional and thermal structure of the lithosphere and upper mantle: III. Thermochemical tomography in the Western‐Central U.S., Journal of Geophysical Research, 121, doi:10.1002/2016JB013049, 2016</p>

PurposeThe application of statistical inversion theory provides a powerful approach for solving estimation problems including the ability for uncertainty quantification (UQ) by means of Markov chain Monte Carlo (MCMC) methods and Monte Carlo integration. This paper aims to analyze the application of a state reduction technique within different MCMC techniques to improve the computational efficiency and the tuning process of these algorithms.Design/methodology/approachA reduced state representation is constructed from a general prior distribution. For sampling the Metropolis Hastings (MH) Algorithm and the Gibbs sampler are used. Efficient proposal generation techniques and techniques for conditional sampling are proposed and evaluated for an exemplary inverse problem.FindingsFor the MH-algorithm, high acceptance rates can be obtained with a simple proposal kernel. For the Gibbs sampler, an efficient technique for conditional sampling was found. The state reduction scheme stabilizes the ill-posed inverse problem, allowing a solution without a dedicated prior distribution. The state reduction is suitable to represent general material distributions.Practical implicationsThe state reduction scheme and the MCMC techniques can be applied in different imaging problems. The stabilizing nature of the state reduction improves the solution of ill-posed problems. The tuning of the MCMC methods is simplified.Originality/valueThe paper presents a method to improve the solution process of inverse problems within the Bayesian framework. The stabilization of the inverse problem due to the state reduction improves the solution. The approach simplifies the tuning of MCMC methods.

An introduction of Markov chain Monte Carlo method to geochemical inverse problems: Reading melting parameters from REE abundances in abyssal peridotites

Markov chain Monte Carlo (MCMC) methods form one of the algorithmic foundations of Bayesian inverse problems. The recent development of likelihood-informed subspace (LIS) methods offers a viable route to designing efficient MCMC methods for exploring high-dimensional posterior distributions via exploiting the intrinsic low-dimensional structure of the underlying inverse problem. However, existing LIS methods and the associated performance analysis often assume that the prior distribution is Gaussian. This assumption is limited for inverse problems aiming to promote sparsity in the parameter estimation, as heavy-tailed priors, e.g., Laplace distribution or the elastic net commonly used in Bayesian LASSO, are often needed in this case. To overcome this limitation, we consider a prior normalization technique that transforms any non-Gaussian (e.g. heavy-tailed) priors into standard Gaussian distributions, which makes it possible to implement LIS methods to accelerate MCMC sampling via such transformations. We also rigorously investigate the integration of such transformations with several MCMC methods for high-dimensional problems. Finally, we demonstrate various aspects of our theoretical claims on two nonlinear inverse problems.

This article aims at the acceleration of an inverse heat transfer problem solution within the Bayesian framework. The physical problem involves a spatially varying heat flux, which can reach very large magnitudes in small regions, such as in the heating imposed by high-power lasers. The inverse problem of estimating the imposed heat flux is solved by using the Markov chain Monte Carlo method, with simulated transient temperature measurements. The solution of the inverse problem is based on a reduced model, which consists of an improved lumped formulation of a linearized version of the original nonlinear problem. Two different priors are considered for the sought heat flux, including a total variation density and a Gaussian density. The Gaussian prior is based on the physics of the heat conduction problem. Parameters appearing in both priors are also estimated as part of the inference problem in hyperprior models. The Delayed Acceptance Metropolis-Hastings (DAMH) Algorithm and the Enhanced Approximation Error Model (AEM) are applied with the objective to improve the accuracy of the inverse problem solution.

Markov chain Monte Carlo methods such as the Gibbs sampler and the Metropolis algorithm can be used to sample solutions to non-linear inverse problems. In principle, these methods allow incorporation of prior information of arbitrary complexity. If an analytical closed form description of the prior is available, which is the case when the prior can be described by a multidimensional Gaussian distribution, such prior information can easily be considered. In reality, prior information is often more complex than can be described by the Gaussian model, and no closed form expression of the prior can be given. We propose an algorithm, called sequential Gibbs sampling, allowing the Metropolis algorithm to efficiently incorporate complex priors into the solution of an inverse problem, also for the case where no closed form description of the prior exists. First, we lay out the theoretical background for applying the sequential Gibbs sampler and illustrate how it works. Through two case studies, we demonstrate the application of the method to a linear image restoration problem and to a non-linear cross-borehole inversion problem. We demonstrate how prior information can reduce the complexity of an inverse problem and that a prior with little information leads to a hard inverse problem, practically unsolvable except when the number of model parameters is very small. Considering more complex and realistic prior information thus not only makes realizations from the posterior look more realistic but it can also reduce the computation time for the inversion dramatically. The method works for any statistical model for which sequential simulation can be used to generate realizations. This applies to most algorithms developed in the geostatistical community.

In this paper, we use the theory of symmetric Dirichlet forms to give a probabilistic interpretation of Calderón's inverse conductivity problem in terms of reflecting diffusion processes and their corresponding boundary trace processes. This probabilistic interpretation comes in three equivalent formulations which open up novel perspectives on the classical question of unique determinability of conductivities from boundary data. We aim to make this work accessible to both readers with a background in stochastic process theory as well as researchers working on deterministic methods in inverse problems.

For strongly non-linear and high-dimensional inverse problems, Markov chain Monte Carlo (MCMC) methods may fail to properly explore the posterior probability density function (PDF) given a realistic computational budget and are generally poorly amenable to parallelization. Particle methods approximate the posterior PDF using the states and weights of a population of evolving particles and they are very well suited to parallelization. We focus on adaptive sequential Monte Carlo (ASMC), an extension of annealed importance sampling (AIS). In AIS and ASMC, importance sampling is performed over a sequence of intermediate distributions, known as power posteriors, linking the prior to the posterior PDF. The AIS and ASMC algorithms also provide estimates of the evidence (marginal likelihood) as needed for Bayesian model selection, at basically no additional cost. ASMC performs better than AIS as it adaptively tunes the tempering schedule and performs resampling of particles when the variance of the particle weights becomes too large. We consider a challenging synthetic groundwater transport inverse problem with a categorical channelized 2D hydraulic conductivity field defined such that the posterior facies distribution includes two distinct modes. The model proposals are obtained by iteratively re-simulating a fraction of the current model using conditional multiple-point statistics (MPS) simulations. We examine how ASMC explores the posterior PDF and compare with results obtained with parallel tempering (PT), a state-of-the-art MCMC inversion approach that runs multiple interacting chains targeting different power posteriors. For a similar computational budget, ASMC outperforms PT as the ASMC-derived models fit the data better and recover the reference likelihood. Moreover, we show that ASMC partly retrieves both posterior modes, while none of them is recovered by PT. Lastly, we demonstrate how the power posteriors obtained by ASMC can be used to assess the influence of the assumed data errors on the posterior means and variances, as well as on the evidence. We suggest that ASMC can advantageously replace MCMC for solving many challenging inverse problems arising in the field of water resources.

Context.The characterization of the interior of an exoplanet is an inverse problem. The solution requires statistical methods such as Bayesian inference. Current methods employ Markov chain Monte Carlo (MCMC) sampling to infer the posterior probability of the planetary structure parameters for a given exoplanet. These methods are time-consuming because they require the evaluation of a planetary structure model ~105times.Aims.To speed up the inference process when characterizing an exoplanet, we propose to use conditional invertible neural networks to calculate the posterior probability of the planetary structure parameters.Methods.Conditional invertible neural networks (cINNs) are a special type of neural network that excels at solving inverse problems. We constructed a cINN following the framework for easily invertible architectures (FreIA). This neural network was then trained on a database of 5.6 × 106internal structure models to recover the inverse mapping between internal structure parameters and observable features (i.e., planetary mass, planetary radius, and elemental composition of the host star). We also show how observational uncertainties can be accounted for.Results.The cINN method was compared to a commonly used Metropolis-Hastings MCMC. To do this, we repeated the characterization of the exoplanet K2-111 b, using both the MCMC method and the trained cINN. We show that the inferred posterior probability distributions of the internal structure parameters from both methods are very similar; the largest differences are seen in the exoplanet water content. Thus, cINNs are a possible alternative to the standard time-consuming sampling methods. cINNs allow infering the composition of an exoplanet that is orders of magnitude faster than what is possible using an MCMC method. The computation of a large database of internal structures to train the neural network is still required, however. Because this database is only computed once, we found that using an invertible neural network is more efficient than an MCMC when more than ten exoplanets are characterized using the same neural network.

Traditional full-waveform inversion (FWI) methods only render a “best-fit” model that cannot account for uncertainties of the ill-posed inverse problem. Additionally, local optimization-based FWI methods cannot always converge to a geologically meaningful solution unless the inversion starts with an accurate background model. We seek the solution for FWI in the Bayesian inference framework to address those two issues. In Bayesian inference, the model space is directly probed by sampling methods such that we obtain a reliable uncertainty appraisal, determine optimal models, and avoid entrapment in a small local region of the model space. The solution of such a statistical inverse method is completely described by the posterior distribution, which quantifies the distributions for parameters and inversion uncertainties. To efficiently sample the posterior distribution, we introduce a sampling algorithm in which the proposal distribution is constructed by the local gradient and the diagonal approximate Hessian of the local log posterior. Our algorithm is called the gradient-based Markov chain Monte Carlo (GMCMC) method. The GMCMC FWI method can quantify inversion uncertainties with estimated posterior distribution given sufficiently long Markov chains. By directly sampling the posterior distribution, we obtain a global view of the model space. Theoretically speaking, statistical assessments do not depend on starting models. Our method is applied to the 2D Marmousi model with the frequency-domain FWI setting. Numerical results suggest that our method can be readily applied to 2D cases with affordable computational efforts.

Bayesian model updating has received considerable attention and has been extensively used in structural damage detection. It provides a rigorous statistical framework for realizing structural system identification and characterizing uncertainties associated with modeling and measurements. The Markov Chain Monte Carlo (MCMC) is a promising tool for inferring the posterior distribution of model parameters to avoid the intractable evaluation of multi-dimensional integration. However, the efficacy of most MCMC techniques suffers from the curse of parameter dimension, which restricts the application of Bayesian model updating to the damage detection of large-scale systems. In addition, there are several MCMC techniques that require users to properly choose application-specific models, based on the understanding of algorithm mechanisms and limitations. As seen in the literature, there is a lack of comprehensive work that investigates the performances of various MCMC algorithms in their application of structural damage detection. In this study, the Differential Evolutionary Adaptive Metropolis (DREAM), a multi-chain MCMC, is explored and adapted to Bayesian model updating. This paper illustrates how DREAM is used for model updating with many uncertainty parameters (i.e., 40 parameters). Furthermore, the study provides a tutorial to users who may be less experienced with Bayesian model updating and MCMC. Two advanced single-chain MCMC algorithms, namely, the Delayed Rejection Adaptive Metropolis (DRAM) and Transitional Markov Chain Monte Carlo (TMCMC), and DREAM are elaborately introduced to allow practitioners to understand better the concepts and practical implementations. Their performances in model updating and damage detection are compared through three different engineering applications with increased complexity, e.g., a forty-story shear building, a two-span continuous steel beam, and a large-scale steel pedestrian bridge.

Current deep learning-based solutions for image analysis tasks are commonly incapable of handling problems to which multiple different plausible solutions exist. In response, posterior-based methods such as conditional Diffusion Models and Invertible Neural Networks have emerged; however, their translation is hampered by a lack of research on adequate validation. In other words, the way progress is measured often does not reflect the needs of the driving practical application. Closing this gap in the literature, we present the first systematic framework for the application-driven validation of posterior-based methods in inverse problems. As a methodological novelty, it adopts key principles from the field of object detection validation, which has a long history of addressing the question of how to locate and match multiple object instances in an image. Treating modes as instances enables us to perform mode-centric validation, using well-interpretable metrics from the application perspective. We demonstrate the value of our framework through instantiations for a synthetic toy example and two medical vision use cases: pose estimation in surgery and imaging-based quantification of functional tissue parameters for diagnostics. Our framework offers key advantages over common approaches to posterior validation in all three examples and could thus revolutionize performance assessment in inverse problems.