An official website of the United States government
Bayesian inference provides a systematic framework for integration of data with mathematical models to quantify the uncertainty in the solution of the inverse problem. However, the solution of Bayesian inverse problems governed by complex forward models described by partial differential equations (PDEs) remains prohibitive with black-box Markov chain Monte Carlo (MCMC) methods. We present hIPPYlib-MUQ, an extensible and scalable software framework that contains implementations of state-of-the art algorithms aimed to overcome the challenges of high-dimensional, PDE-constrained Bayesian inverse problems. These algorithms accelerate MCMC sampling by exploiting the geometry and intrinsic low-dimensionality of parameter space via derivative information and low rank approximation. The software integrates two complementary open-source software packages, hIPPYlib and MUQ. hIPPYlib solves PDE-constrained inverse problems using automatically-generated adjoint-based derivatives, but it lacks full Bayesian capabilities. MUQ provides a spectrum of powerful Bayesian inversion models and algorithms, but expects forward models to come equipped with gradients and Hessians to permit large-scale solution. By combining these two complementary libraries, we created a robust, scalable, and efficient software framework that realizes the benefits of each and allows us to tackle complex large-scale Bayesian inverse problems across a broad spectrum of scientific and engineering disciplines. To illustrate the capabilities of hIPPYlib-MUQ, we present a comparison of a number of MCMC methods available in the integrated software on several high-dimensional Bayesian inverse problems. These include problems characterized by both linear and nonlinear PDEs, various noise models, and different parameter dimensions. The results demonstrate that large (∼ 50×) speedups over conventional black box and gradient-based MCMC algorithms can be obtained by exploiting Hessian information (from the log-posterior), underscoring the power of the integrated hIPPYlib-MUQ framework.
more »
« less
A mixed, unified forward/inverse framework for earthquake problems: fault implementation and coseismic slip estimate
SUMMARY We introduce a new finite-element (FE) based computational framework to solve forward and inverse elastic deformation problems for earthquake faulting via the adjoint method. Based on two advanced computational libraries, FEniCS and hIPPYlib for the forward and inverse problems, respectively, this framework is flexible, transparent and easily extensible. We represent a fault discontinuity through a mixed FE elasticity formulation, which approximates the stress with higher order accuracy and exposes the prescribed slip explicitly in the variational form without using conventional split node and decomposition discrete approaches. This also allows the first order optimality condition, that is the vanishing of the gradient, to be expressed in continuous form, which leads to consistent discretizations of all field variables, including the slip. We show comparisons with the standard, pure displacement formulation and a model containing an in-plane mode II crack, whose slip is prescribed via the split node technique. We demonstrate the potential of this new computational framework by performing a linear coseismic slip inversion through adjoint-based optimization methods, without requiring computation of elastic Green’s functions. Specifically, we consider a penalized least squares formulation, which in a Bayesian setting—under the assumption of Gaussian noise and prior—reflects the negative log of the posterior distribution. The comparison of the inversion results with a standard, linear inverse theory approach based on Okada’s solutions shows analogous results. Preliminary uncertainties are estimated via eigenvalue analysis of the Hessian of the penalized least squares objective function. Our implementation is fully open-source and Jupyter notebooks to reproduce our results are provided. The extension to a fully Bayesian framework for detailed uncertainty quantification and non-linear inversions, including for heterogeneous media earthquake problems, will be analysed in a forthcoming paper.
more »
« less
- Award ID(s):
- 2121666
- PAR ID:
- 10395792
- Publisher / Repository:
- Oxford University Press
- Date Published:
- Journal Name:
- Geophysical Journal International
- Volume:
- 230
- Issue:
- 2
- ISSN:
- 0956-540X
- Page Range / eLocation ID:
- p. 733-758
- Format(s):
- Medium: X
- Sponsoring Org:
- National Science Foundation
More Like this
-
-
We present an adjoint-based optimization method to invert for stress and frictional parameters used in earthquake modeling. The forward problem is linear elastodynamics with nonlinear rate-and-state frictional faults. The misfit functional quantifies the difference between simulated and measured particle displacements or velocities at receiver locations. The misfit may include windowing or filtering operators. We derive the corresponding adjoint problem, which is linear elasticity with linearized rate-and-state friction and, for forward problems involving fault normal stress changes, nonzero fault opening, with time-dependent coefficients derived from the forward solution. The gradient of the misfit is efficiently computed by convolving forward and adjoint variables on the fault. The method thus extends the framework of full-waveform inversion to include frictional faults with rate-and-state friction. In addition, we present a space-time dual-consistent discretization of a dynamic rupture problem with a rough fault in antiplane shear, using high-order accurate summation-by-parts finite differences in combination with explicit Runge–Kutta time integration. The dual consistency of the discretization ensures that the discrete adjoint-based gradient is the exact gradient of the discrete misfit functional as well as a consistent approximation of the continuous gradient. Our theoretical results are corroborated by inversions with synthetic data. We anticipate that adjoint-based inversion of seismic and/or geodetic data will be a powerful tool for studying earthquake source processes; it can also be used to interpret laboratory friction experiments.more » « less
-
SUMMARY Analysis of tectonic and earthquake-cycle associated deformation of the crust can provide valuable insights into the underlying deformation processes including fault slip. How those processes are expressed at the surface depends on the lateral and depth variations of rock properties. The effect of such variations is often tested by forward models based on a priori geological or geophysical information. Here, we first develop a novel technique based on an open-source finite-element computational framework to invert geodetic constraints directly for heterogeneous media properties. We focus on the elastic, coseismic problem and seek to constrain variations in shear modulus and Poisson’s ratio, proxies for the effects of lithology and/or temperature and porous flow, respectively. The corresponding nonlinear inversion is implemented using adjoint-based optimization that efficiently reduces the cost function that includes the misfit between the calculated and observed displacements and a penalty term. We then extend our theoretical and numerical framework to simultaneously infer both heterogeneous Earth’s structure and fault slip from surface deformation. Based on a range of 2-D synthetic cases, we find that both model parameters can be satisfactorily estimated for the megathrust setting-inspired test problems considered. Within limits, this is the case even in the presence of noise and if the fault geometry is not perfectly known. Our method lays the foundation for a future reassessment of the information contained in increasingly data-rich settings, for example, geodetic GNSS constraints for large earthquakes such as the 2011 Tohoku-oki M9 event, or distributed deformation along plate boundaries as constrained from InSAR.more » « less
-
Abstract The inverse problem for radiative transfer is important in many applications, such as optical tomography and remote sensing. Major challenges include large memory requirements and computational expense, which arise from high-dimensionality and the need for iterations in solving the inverse problem. Here, to alleviate these issues, we propose adaptive-mesh inversion: a goal-orientedhp-adaptive mesh refinement method for solving inverse radiative transfer problems. One novel aspect here is that the two optimizations (one for inversion, and one for mesh adaptivity) are treated simultaneously and blended together. By exploiting the connection between duality-based mesh adaptivity and adjoint-based inversion techniques, we propose a goal-oriented error estimator, which is cheap to compute, and can efficiently guide the mesh-refinement to numerically solve the inverse problem. We use discontinuous Galerkin spectral element methods to discretize the forward and the adjoint problems. Then, based on the goal-oriented error estimator, we propose anhp-adaptive algorithm to refine the meshes. Numerical experiments are presented at the end and show convergence speed-up and reduced memory occupation by the goal-oriented mesh adaptive method.more » « less
-
Abstract We consider Bayesian inference for large-scale inverse problems, where computational challenges arise from the need for repeated evaluations of an expensive forward model. This renders most Markov chain Monte Carlo approaches infeasible, since they typically require O ( 1 0 4 ) model runs, or more. Moreover, the forward model is often given as a black box or is impractical to differentiate. Therefore derivative-free algorithms are highly desirable. We propose a framework, which is built on Kalman methodology, to efficiently perform Bayesian inference in such inverse problems. The basic method is based on an approximation of the filtering distribution of a novel mean-field dynamical system, into which the inverse problem is embedded as an observation operator. Theoretical properties are established for linear inverse problems, demonstrating that the desired Bayesian posterior is given by the steady state of the law of the filtering distribution of the mean-field dynamical system, and proving exponential convergence to it. This suggests that, for nonlinear problems which are close to Gaussian, sequentially computing this law provides the basis for efficient iterative methods to approximate the Bayesian posterior. Ensemble methods are applied to obtain interacting particle system approximations of the filtering distribution of the mean-field model; and practical strategies to further reduce the computational and memory cost of the methodology are presented, including low-rank approximation and a bi-fidelity approach. The effectiveness of the framework is demonstrated in several numerical experiments, including proof-of-concept linear/nonlinear examples and two large-scale applications: learning of permeability parameters in subsurface flow; and learning subgrid-scale parameters in a global climate model. Moreover, the stochastic ensemble Kalman filter and various ensemble square-root Kalman filters are all employed and are compared numerically. The results demonstrate that the proposed method, based on exponential convergence to the filtering distribution of a mean-field dynamical system, is competitive with pre-existing Kalman-based methods for inverse problems.more » « less
