PaperlayerView exact revision on arXiv

Source-linked AI summary

A computational framework for infinite-dimensional Bayesian inverse problems. Part I: The linearized case, with application to global seismic inversion

Tan Bui-Thanh, Omar Ghattas, James Martin, Georg Stadler

arXiv:1308.1313v1math.NAmath.OCstat.COstat.ME

TL;DR

The paper addresses uncertainty quantification for linearized infinite-dimensional inverse problems, where prior choice, convergent discretization, and high-dimensional posterior computation are challenging. It develops a Bayesian framework using Gaussian priors, consistent discretization, and low-rank matrix-free covariance approximations, and applies it to 3D global seismic inversion. The resulting approach requires a dimension-independent number of forward PDE solves for local covariance approximation and handles problems with up to 431,000 unknowns.

  • Problem

    Infinite-dimensional Bayesian inverse problems require well-posed priors, convergent discretizations, and tractable posterior-statistics computation despite expensive forward solves and very high parameter dimensions.

  • Method

    The framework uses Gaussian random-field priors, consistent finite-dimensional discretization, MAP-based linearization, and matrix-free Lanczos low-rank approximation of the data-informed posterior covariance.

  • Results

    The method requires a dimension-independent number of forward PDE solves for local covariance approximation and is demonstrated on global seismic problems with up to 431,000 unknown parameters.

  • Takeaways & Limitations

    Uncertainty quantification for the linearized inverse problem reduces to solving a fixed number of forward and adjoint PDEs while retaining a scalable posterior representation.

  • Takeaways & Limitations

    The framework restricts the prior to Gaussian random fields and uses a linearized Gaussian posterior approximation, which is reasonable when the parameter-to-observable map is nearly linear.

Abstract

from arXiv · show

We present a computational framework for estimating the uncertainty in the numerical solution of linearized infinite-dimensional statistical inverse problems. We adopt the Bayesian inference formulation: given observational data and their uncertainty, the governing forward problem and its uncertainty, and a prior probability distribution describing uncertainty in the parameter field, find the posterior probability distribution over the parameter field. The prior must be chosen appropriately in order to guarantee well-posedness of the infinite-dimensional inverse problem and facilitate computation of the posterior. Furthermore, straightforward discretizations may not lead to convergent approximations of the infinite-dimensional problem. And finally, solution of the discretized inverse problem via explicit construction of the covariance matrix is prohibitive due to the need to solve the forward problem as many times as there are parameters. Our computational framework builds on the infinite-dimensional formulation proposed by Stuart (A. M. Stuart, Inverse problems: A Bayesian perspective, Acta Numerica, 19 (2010), pp. 451-559), and incorporates a number of components aimed at ensuring a convergent discretization of the underlying infinite-dimensional inverse problem. The framework additionally incorporates algorithms for manipulating the prior, constructing a low rank approximation of the data-informed component of the posterior covariance operator, and exploring the posterior that together ensure scalability of the entire framework to very high parameter dimensions. We demonstrate this computational framework on the Bayesian solution of an inverse problem in 3D global seismic wave propagation with hundreds of thousands of parameters.

1. Introduction.

The paper develops a scalable Bayesian framework for uncertainty quantification in large-scale, infinite-dimensional inverse problems. It combines consistent discretization, posterior covariance approximation, and scalable computation, then applies the framework to global seismic inversion.

  • Bayesian inference provides a complete statistical description of parameters consistent with noisy data, model uncertainty, and prior information, beyond a single best-fit estimate.
  • Infinite-dimensional inverse problems require priors that ensure well-posedness, discretizations that converge, and scalable posterior covariance algorithms.
  • Posterior statistics are difficult because high-dimensional MCMC requires many expensive forward PDE solves, motivating the paper’s linearized inverse-problem setting.
  • The framework approximates the prior-preconditioned data-misfit Hessian with matrix-free Lanczos iterations and uses Sherman-Morrison-Woodbury to represent posterior covariance.
  • Up to 431,000 unknown parameters are treated in realistic 3D global-seismology problems, with related work extending the approach beyond one million parameters.

2. Bayesian framework for infinite-dimensional inverse problems.

The Bayesian formulation places the inverse problem on a function space and uses a Gaussian prior designed to support well-posed inference and computation. A MAP estimate and local linearization then produce the paper’s Gaussian posterior approximation.

  • The parameter field is defined on an open, bounded, sufficiently regular subset of R3, with observables generated by a PDE-based parameter-to-observable map.
  • The prior is a Gaussian random field with covariance C0 = A^-2, where A is Laplacian-like and its inverse provides smoothing needed for well-posedness.
  • The inverse problem seeks a posterior measure over parameter fields by combining a likelihood for observations with a prior measure through Bayes’ formula.
  • The prior construction yields bounded variance and almost surely continuous samples because A^-2 is trace class on L2(Ω).
  • The MAP estimate is obtained by solving an optimization problem, after which the parameter-to-observable map is linearized to obtain a Gaussian posterior approximation.
  • The linearized approximation is considered reasonable when the map is nearly linear and may also be useful in small-noise or many-observation regimes.

3. Discretization of the Bayesian inverse problem.

The discretization preserves the infinite-dimensional structure through finite-element spaces and mass-weighted inner products. It then represents priors, posteriors, samples, and variance fields while avoiding explicit dense covariance construction.

  • Careful discretization is necessary because naive finite-dimensional approximations may not converge to the desired infinite-dimensional solution.
  • A mass-matrix-weighted inner product better matches the infinite-dimensional L2 structure, although it introduces distinctions between transpose and adjoint operators.
  • The framework computes posterior samples and pointwise variance fields using finite-dimensional Gaussian representations and the low-rank posterior covariance form.
  • Finite-element discretization represents the parameter field in a subspace Vh, with its nodal coefficient vector serving as the finite-dimensional parameter.
  • The discretized prior uses stiffness and mass matrices, yielding Γprior = A^-2 and enabling scalable prior operations through elliptic-operator solves.
  • Explicit posterior covariance construction is prohibitive because forming the Jacobian generally requires n forward PDE solves when n is large.

4. Finding the MAP point.

The MAP point is computed by scalable inexact Newton–CG optimization without explicitly constructing the Hessian. Its PDE-based Hessian-vector products and mesh-independent iterations support large problems.

  • The MAP problem is solved using an inexact matrix-free Newton–CG method requiring only Hessian-vector products.Each product uses linearized forward-like and adjoint-like PDE solves, so the Hessian matrix is never formed explicitly.
  • PDE-like system solves give the method scalability with respect to parameter dimension.
  • Outer Newton and inner CG iteration counts can remain independent of mesh size.This follows from Newton optimization, compactness of the data-misfit Hessian, and prior-based preconditioning that makes the preconditioned Hessian a compact perturbation of identity.

5. Low rank approximation of the Hessian matrix.

The framework approximates the posterior covariance by compressing the prior-preconditioned data-misfit Hessian into a low-rank representation. Matrix-free eigensolvers and prior-aware formulas then support posterior variance computation and sampling at scalable cost.

  • Explicit Hessian construction is prohibitive because it requires one forward-like solve per uncertain parameter, while the data-informed subspace is typically low dimensional.The data-misfit Hessian behaves like a compact operator because sparse observations and smooth forward maps suppress many parameter modes.
  • The method uses Lanczos iterations to construct a low-rank approximation of the Gauss–Newton data-misfit Hessian and applies Sherman–Morrison–Woodbury to approximate posterior covariance.The same approach can target the full Hessian when noisy data make the Gauss–Newton approximation less adequate.
  • Low-rank posterior representations support posterior sampling and efficient pointwise variance computation.Sampling uses a covariance factorization, while variance evaluation uses the low-rank factors and prior covariance operations.
  • Rapid eigenvalue decay permits retaining only the dominant eigenpairs, while eigenvalues small relative to 1 can be neglected with controlled truncation error.
  • The posterior covariance equals prior uncertainty reduced by data information filtered through the prior.The covariance reduction can be used to compute pointwise variance fields and generate posterior samples.
  • Matrix-free Hessian actions require one incremental forward solve, one incremental adjoint solve, and scalable prior square-root operations, regardless of parameter dimension.Lanczos therefore avoids explicitly constructing the parameter-to-observable map.
  • The number of Lanczos iterations is independent of discretized parameter dimension when the misfit Hessian is compact and the prior covariance operator is continuous.
  • Posterior covariance estimation requires a constant number of forward and adjoint PDE solves independent of parameters, observations, and state variables.Overall scalability follows when the forward, adjoint, and prior elliptic solvers are scalable.

6. Application to global seismic statistical inversion.

The framework is applied to global seismic inversion, reconstructing heterogeneous acoustic wavespeed from synthetic noisy seismograms using a PREM-based prior and a structured, smooth parameterization. The prior incorporates anisotropic correlations and boundary effects relevant to Earth modeling.

  • 6. Application to global seismic statistical inversion: S20RTS supplies the synthetic ground-truth Earth model, while PREM provides prior knowledge and the starting point for MAP estimation.The inversion seeks to reconstruct S20RTS from noisy synthetic data and quantify the resulting uncertainty.
  • 6.1. Parameter space for seismic inversion: The seismic parameter field is the wavespeed anomaly relative to PREM, so its prior mean is zero.
  • 6.1. Parameter space for seismic inversion: The anomaly is modeled with continuous trilinear finite elements on an octree-based hexahedral mesh aligned with major Earth interfaces.The parameter mesh is locally refined to resolve seismic wavelengths and is shared with the wave-equation solver.
  • 6.2. The choice of prior: The prior covariance is chosen to produce smooth, almost surely continuous deviations with bounded, physically meaningful variance.Its precision operator is elliptic, and the prior is designed to support well-posedness of the infinite-dimensional inverse problem.
  • 6.2. The choice of prior: Green’s functions for the prior precision operator show anisotropy that is strongest for points closer to the surface.Larger Green’s-function values are shown as brighter gray shades and correspond directly to the covariance function.
  • 6.2. The choice of prior: The prior uses α = 1.5 · 10^-2 and θ = 4 · 10^-2, with θ introducing longer tangential than radial correlation lengths.The anisotropy decreases smoothly toward the sphere’s center.
  • 6.2. The choice of prior: The prior standard deviation is larger near the boundary, mainly because of the homogeneous Neumann boundary condition used to construct the prior square root.

6.3. The likelihood.

The likelihood is defined by solving the acoustic wave equation for a candidate wave speed, recording receiver velocities, and comparing truncated Fourier coefficients with observations under prescribed noise.

  • The acoustic wave initial-boundary value problem models seismic-wave propagation using velocity and strain dilatation, with specified initial and boundary conditions.The boundary condition sets strain dilatation to zero on the earth's boundary during the observation interval.
  • The parameter-to-observable map solves the wave equation for a given wave speed and records velocity at finitely many receivers over time.The recorded waveforms are subsequently transformed into frequency-domain observations.
  • The noise model uses a diagonal covariance matrix with constant variance.

6.4. Gradient and Hessian of the negative log posterior.

The paper computes derivatives of the negative log posterior through forward, adjoint, and incremental wave equations. A positive Gauss–Newton Hessian approximation omits selected nonlinear terms.

  • The adjoint wave equations run backward in time, use the data misfit as a source term, and otherwise resemble the forward equations.
  • The Hessian operator acts on a wave-speed variation through incremental forward and incremental adjoint wave propagation problems.
  • The gradient computation requires one forward and one adjoint wave solve, while a Hessian action additionally requires incremental-forward and incremental-adjoint solves.
  • The computations use a Gauss–Newton Hessian approximation that is guaranteed positive by neglecting terms containing ∇·w and a term involving d.

6.5. Discretization of the wave equation and implementation details.

The implementation combines continuous finite elements for the parameter with high-order discontinuous Galerkin wave solves and scalable solvers for derivative and covariance computations.

  • The parameter uses trilinear finite elements, while forward, adjoint, incremental forward, and incremental adjoint waves use high-order discontinuous Galerkin discretization.The same hexahedral mesh is used for the parameter and wave solution.
  • The discretized gradient and Hessian action agree with derivatives of the discretized posterior only in the limit as mesh size approaches zero.The stated inconsistency arises from additional numerical-flux jump terms at element interfaces.
  • Parameter-to-wave transfer prolongates the continuous parameter into the wave-solution space, while computed derivatives are restricted back to parameter space.
  • Parallel algebraic multigrid solves repeatedly apply A^-1, with elliptic-solve cost reported as negligible relative to time-dependent seismic-wave solves.
  • The implementation addresses large time-history storage requirements for backward adjoint calculations and Hessian-vector applications.
  • The seismic inverse problems require 1200–4096 processor cores for 10–20 hours, with most runtime spent on forward, adjoint, and incremental wave solves.

6.6. Setup of model problems.

The experiments use synthetic global-seismic observations generated from S20RTS, frequency-truncated receiver data, and two source–receiver configurations with different problem sizes.

  • The observations retain the first 101 Fourier modes, with coefficients between 10^-5 and 10^-1 and noise standard deviation 0.002.
  • Problem I uses one North Pole source, one receiver 45° south of the equator, 1800 s of propagation, and 78,558 parameter nodes.The forward discretization has about 21.4 million spatial wave unknowns and 2100 Runge–Kutta time steps.

6.7. Low rank approximation of the prior-preconditioned misfit Hes-

The prior-preconditioned misfit Hessian has a rapidly decaying spectrum whose approximation is effectively mesh-independent, enabling low-rank posterior computations. Dominant eigenvectors identify the best-observed earth modes, with smoother modes associated with larger eigenvalues.

  • Problem II requires at least 700 eigenvalues for a sensible low-rank Hessian approximation, compared with about 40 for Problem I.The greater requirement reflects the additional information provided by three simultaneous sources and 130 receivers in Problem II.
  • The three Problem II discretizations produce essentially overlapping eigenvalue spectra, indicating that the dominant spectrum is mesh-independent.This supports sufficiently fine meshes for resolving the dominant eigenvectors and approximating the infinite-dimensional problem.
  • Figure 6.6 visualizes the first, 3rd, 5th, 8th, and 13th eigenvectors for Problem I using a slice through source and receiver locations.
  • The low-rank approximation requires a number of wave propagation solutions that does not depend on the parameter dimension in this example.The number of Lanczos iterations is independent of the number of discrete parameters.
  • Dominant Hessian eigenvectors represent earth modes most observable from the source-and-receiver configuration.The largest eigenvalues correspond to the smoothest modes, while decreasing eigenvalues are associated with increasingly oscillatory modes.

6.8. Interpretation of the uncertainty in the solution of the inverse

The experiments show that sparse seismic data reduce uncertainty and recover wave speed most effectively in regions traversed by informative source–receiver paths, while less-informed fine scales retain prior-like variability.

  • Problem I’s dominant misfit-Hessian modes are smooth, consistent with sparse data reconstructing only limited information from the truth solution.
  • The posterior variance decreases relative to the prior, with lowest uncertainty at the surface halfway between the single source and receiver.The data also inform the core–mantle boundary through stronger reflected energy caused by its strong material contrast.
  • Problem II accurately recovers wave speed in the Northern Hemisphere region covered by sources and receivers.The comparison is made between MAP estimates and the ground-truth S20RTS model at depths of 67km, 670km, and 1340km.
  • Uncertainty reduction is greatest along source–receiver wave paths, particularly in the quarter of the Northern Hemisphere surface containing receivers.At 67km, the largest variance reduction occurs near sources and receivers, while the smallest occurs on the opposite side of Earth.
  • Posterior samples share large-scale features where data are informative, while fine-scale features remain qualitatively prior-like and vary between samples.The samples compare prior and posterior distributions; the MAP estimate is shown separately for comparison.

7. Conclusions.

The conclusions present a scalable framework for posterior sampling and uncertainty quantification in very high-dimensional inverse problems. Its implementation uses covariance factorization and Gaussian sampling constructions to explore the posterior without relying on dense covariance calculations.

  • The framework combines prior manipulation, low-rank posterior-covariance approximation, and posterior exploration for very high parameter dimensions.
  • Posterior Gaussian samples are generated by factoring the posterior covariance and applying the resulting linear map to a standard Gaussian sample.
  • The method constructs the required covariance factor through the spectral decomposition of B and an identity yielding the desired matrix.
  • The prior and posterior sampling derivation is formulated in Rn, with M −1/2 treated as a linear map from Rn to Rn.
  • The construction rewrites the sampling operator using the dimension n and a transformed map L = ˜LM −1/2.
Loading 1308.1313v1…