PaperlayerView exact revision on arXiv

Source-linked AI summary

A computational framework for infinite-dimensional Bayesian inverse problems: Part II. Stochastic Newton MCMC with application to ice sheet flow inverse problems

Noemi Petra, James Martin, Georg Stadler, Omar Ghattas

arXiv:1308.6221v2stat.MEmath.NAmath.OCmath.STstat.CO

TL;DR

The paper addresses computationally difficult sampling for nonlinear, infinite-dimensional Bayesian inverse problems governed by PDEs. It introduces a stochastic Newton MCMC variant that reuses a low-rank Hessian computed at the MAP point, and reports best convergence in an ice sheet flow inverse problem while enabling posterior interpretation through prior- and data-informed directions.

  • Problem

    High-dimensional PDE-governed inverse problems make posterior sampling computationally difficult, especially when nonlinearities produce non-Gaussian posteriors.

  • Method

    The paper uses stochastic Newton MCMC with a low-rank Hessian computed at the MAP point and reused at every MCMC step, alongside Hessian-based posterior analysis.

  • Results

    The MAP-based Hessian stochastic Newton method achieves the best convergence in both sample count and PDE solves among the compared methods for the ice sheet flow inverse problem.

  • Takeaways & Limitations

    Spectral structure from the MAP Hessian and prior covariance makes high-dimensional posterior visualization and classification of prior- and observation-informed directions tractable.

  • Takeaways & Limitations

    Predictions for directions substantially influenced by both observations and prior may be too optimistic.

Abstract

from arXiv · show

We address the numerical solution of infinite-dimensional inverse problems in the framework of Bayesian inference. In the Part I companion to this paper (arXiv.org:1308.1313), we considered the linearized infinite-dimensional inverse problem. Here in Part II, we relax the linearization assumption and consider the fully nonlinear infinite-dimensional inverse problem using a Markov chain Monte Carlo (MCMC) sampling method. To address the challenges of sampling high-dimensional pdfs arising from Bayesian inverse problems governed by PDEs, we build on the stochastic Newton MCMC method. This method exploits problem structure by taking as a proposal density a local Gaussian approximation of the posterior pdf, whose construction is made tractable by invoking a low-rank approximation of its data misfit component of the Hessian. Here we introduce an approximation of the stochastic Newton proposal in which we compute the low-rank-based Hessian at just the MAP point, and then reuse this Hessian at each MCMC step. We compare the performance of the proposed method to the original stochastic Newton MCMC method and to an independence sampler. The comparison of the three methods is conducted on a synthetic ice sheet inverse problem. For this problem, the stochastic Newton MCMC method with a MAP-based Hessian converges at least as rapidly as the original stochastic Newton MCMC method, but is far cheaper since it avoids recomputing the Hessian at each step. On the other hand, it is more expensive per sample than the independence sampler; however, its convergence is significantly more rapid, and thus overall it is much cheaper. Finally, we present extensive analysis and interpretation of the posterior distribution, and classify directions in parameter space based on the extent to which they are informed by the prior or the observations.

1. Introduction and background.

The paper extends stochastic Newton MCMC to nonlinear infinite-dimensional Bayesian inverse problems and proposes reusing a MAP-based Hessian to reduce computational cost. It evaluates this approach on an ice sheet flow inverse problem and analyzes how prior and data inform posterior directions.

  • Motivation: 14 The high dimensionality of uncertain parameters and expensive forward models make conventional MCMC approaches intractable for large-scale inverse problems.These difficulties arise after discretizing infinite-dimensional inverse problems.
  • Existing approach: Stochastic Newton MCMC uses local Hessian-based Gaussian proposals to capture the highly stretched posterior contours typical of ill-posed inverse problems.The proposals exploit local curvature and guide sampling toward regions with higher acceptance probability.
  • Proposed method: The proposed method computes a low-rank Hessian at the MAP point and reuses it at every sample point, while retaining local gradient information.This avoids dynamically recomputing the Hessian during MCMC.
  • Evaluation: The method is compared with dynamically changing-Hessian stochastic Newton MCMC and a MAP-centered Gaussian independence sampler on a nonlinear Stokes ice sheet flow inverse problem.The experiment infers a basal boundary condition parameter field from surface velocity observations.
  • Contributions: The paper extends an infinite-dimensional computational framework from linearized to nonlinear inverse problems, requiring MCMC sampling and careful prior, likelihood, discretization, and inner-product treatment.The nonlinear setting produces non-Gaussian posteriors.
  • Posterior analysis: Posterior directions are interpreted using the spectral structure of the MAP-point Hessian and prior covariance, including directions jointly influenced by observations and prior.The paper notes that predictions for these mixed directions may be too optimistic.

2. Background on the infinite-dimensional Bayesian inverse problem, its consistent discretization, and characterization of the posterior.

The paper formulates Bayesian inference for infinite-dimensional parameter fields, develops a consistent finite-dimensional discretization, and characterizes the resulting posterior. It emphasizes the computational difficulty of posterior exploration when the forward model is expensive and nonlinear.

  • Problem formulation: Infinite-dimensional Bayesian inverse problems require addressing the absence of a Lebesgue-based pdf, appropriate prior-measure selection, consistent discretization, and difficult posterior exploration.These are distinct formulation and computational challenges.
  • Infinite-dimensional Bayesian model: The formulation uses a Radon-Nikodym derivative with respect to a Gaussian prior whose covariance is the inverse of a sufficiently high-order elliptic differential operator.The operator choice is intended to support a well-posed inverse problem and computation of the posterior.
  • Prior model: The prior covariance can encode spatially varying or anisotropic correlations through spatially dependent or tensor-valued coefficients.The coefficients also control correlation length and variance.
  • Likelihood model: Observational and model uncertainties are combined into a centered, additive Gaussian noise model for the likelihood.The noise covariance matrix represents uncertainty in the finite-dimensional observations.
  • Consistent discretization: Finite-dimensional approximations use a finite-element subspace, with mass-matrix-weighted inner products to inherit the infinite-dimensional L2 structure.The discretization also requires distinguishing matrix transposes from adjoints under the weighted inner product.
  • Posterior characterization: 49 Posterior evaluation becomes extremely challenging because each evaluation requires solving the forward model, and high-dimensional exploration needs very many evaluations.Methods that do not exploit posterior structure are therefore impractical for expensive forward models and high-dimensional parameter spaces.
  • Posterior characterization: For nonlinear parameter-to-observable maps, the posterior is generally non-Gaussian and cannot be represented solely by its mean and covariance.This necessitates alternative characterization methods such as MCMC.

3. A modified stochastic Newton MCMC method.

The paper modifies stochastic Newton MCMC for large-scale nonlinear inverse problems by computing a low-rank Hessian at the MAP point and reusing it, reducing repeated PDE costs while retaining local gradient information.

  • SNMAP: SNMAP computes a low-rank Hessian once at the MAP point and reuses it for every proposal, while recomputing the gradient locally.The proposal uses a global MAP-based Hessian and a locally evaluated gradient.
  • Original stochastic Newton MCMC: The original stochastic Newton method constructs a local Gaussian posterior approximation using the gradient and Hessian evaluated at the current sample point.Its forward proposal is drawn from this local quadratic approximation and uses a Metropolis-Hastings accept/reject step.
  • Computational cost: SNMAP avoids Hessian recomputation, leaving a pair of forward/adjoint PDE solves per sample for gradient evaluation after MAP initialization.The dominant savings arise because Hessian construction otherwise requires forward/adjoint PDE solves related to the effective rank of the data-misfit Hessian.
  • Tradeoff: SNMAP may reduce acceptance if posterior curvature changes rapidly, although its more stable proposal can also improve acceptance and chain convergence.The paper presents this as a tradeoff to be assessed numerically.
  • Independence sampler: ISMAP uses a fixed Gaussian proposal centered at the MAP point with the MAP Hessian as inverse covariance and omits gradient computation.Once the MAP Hessian is available, its per-sample cost is one forward PDE solve, and samples can be precomputed or generated in parallel.
  • Low-rank approximation: Low-rank Hessian approximations reduce the required Hessian operations after construction to inner products and vector sums, with dominant construction cost O(r) forward/adjoint PDE solves.These operations avoid further PDE solves during Hessian application and related linear algebra.

4. Application to the inversion of basal boundary conditions in ice flow problems.

The application formulates a Bayesian inverse problem for uncertain basal sliding in a two-dimensional ice-flow model, using surface velocity observations and adjoint-based derivatives for Hessian-based sampling.

  • Inverse problem: The inverse problem estimates the uncertain basal sliding coefficient field from pointwise surface velocity observations in an ice-flow model.The parameter-to-observable map restricts the modeled flow solution to ten measurement points on the upper surface.
  • Ice-flow model: Ice is modeled as a non-Newtonian, viscous, incompressible, isothermal fluid governed by mass and momentum balance with Glen’s flow law.The model uses a basal sliding boundary condition and appropriate conditions on the remaining domain boundary.
  • Basal boundary condition: The basal Robin coefficient exp(β) links tangential velocity to tangential traction and subsumes friction, bedrock roughness, and hydrological effects.The paper emphasizes that this coefficient is highly uncertain and is not itself a direct physical parameter.
  • Test geometry: The synthetic test problem uses a two-dimensional ISMIP-HOM geometry based on Haut Glacier d’Arolla, with a 5 km flowline and 100 m grid spacing.The displayed flow field is obtained by solving the ice-flow equations with the specified reference basal sliding field.
  • Observations: Synthetic observations are generated from the reference sliding field by applying the observation operator to the modeled flow and adding independent Gaussian noise.Two noise settings are considered: horizontal/vertical standard deviations of 62/10 and 18/3.
  • Prior: The Gaussian prior is defined through a surface Laplacian-based elliptic operator, with parameters controlling variance and correlation length.For the model problems, a = 10^-2 and b = 10^2, giving a correlation length roughly 5% of the glacier’s total length.
  • Derivative computation: Gradients and Hessian-vector products are computed efficiently with forward, adjoint, and incremental Stokes problems.The Hessian action combines the prior term with terms involving the forward, adjoint, incremental forward, and incremental adjoint fields.

5. Performance of algorithms.

The study computes a MAP estimate, assesses low-rank Hessian structure, and compares three MCMC algorithms using convergence and sampling-efficiency diagnostics.

  • 5.1. Computation of the MAP point.: The MAP point follows the prior mean in the glacier’s unobserved upper part because observations are available only in the lower half.Figure 5.1 compares the truth and MAP sliding coefficient fields and their corresponding surface velocities.
  • 5.1. Computation of the MAP point.: 8 outer Newton iterations reduced the nonlinear residual by a factor of 10^5 to compute the MAP point.Each outer iteration used an inner Newton solve for the nonlinear Stokes equation.
  • 5.1. Computation of the MAP point.: 208 linearized Stokes solves were required for the MAP computation, including 32 Hessian-vector products.The 32 conjugate-gradient iterations across eight outer Newton iterations each required a Hessian-vector product.
  • 5.3. Performance of proposed stochastic Newton MCMC method with MAP-based Hessian.: All 21 prior-preconditioned data misfit Hessian spectra decay rapidly, supporting an accurate low-rank inverse-Hessian approximation.The spectra are examined at the MAP point and at 21 points distributed over the posterior support, with λ = 1 used as a truncation reference.
  • 5.3. Performance of proposed stochastic Newton MCMC method with MAP-based Hessian.: The sampling runs use r = 20 eigenvectors for the low-rank prior-preconditioned data misfit Hessian approximation.The approximation costs 2(r+l) Stokes solves, with l = 5 Lanczos iterations used in the reported setup.
  • 5.3. Performance of proposed stochastic Newton MCMC method with MAP-based Hessian.: SNMAP achieves the best MPSRF values and largest effective sample size, indicating fastest convergence with respect to sample count.For Problem 2, SNMAP is about 10 times more efficient than SN and 5 times more efficient than ISMAP.

6. Analysis and interpretation of the solution of the Bayesian inverse problem.

The paper visualizes the posterior in physical coordinates and covariance-eigenvector coordinates, using these views to relate uncertainty reduction and parameter directions to observations and prior information.

  • 6. Analysis and interpretation of the solution of the Bayesian inverse problem.: Posterior visualization combines physical-coordinate marginals, covariance eigenvectors classified by data and prior influence, and one- and two-dimensional eigenvector marginals.The posterior approximation uses SNMAP samples and kernel density estimation for the marginal densities.
  • 6.1. Point marginals and samples from the posterior.: Point marginals are computed independently at each spatial point, so overlaid prior and posterior samples are needed to indicate spatial correlation.The independently computed marginal densities do not preserve relationships between neighboring points.
  • 6.1. Point marginals and samples from the posterior.: In the unobserved upper glacier, posterior point marginals resemble the prior, whereas observed regions show significantly reduced variance.The posterior therefore changes uncertainty more strongly where surface observations are available.
  • 6.2. Classification of posterior covariance eigenvectors.: The eigenvector analysis compares data-misfit and inverse-prior Rayleigh quotients with eigenvector norms in the observed lower and unobserved upper glacier.This combines information content with spatial localization across the two glacier regions.
  • 6.1. Point marginals and samples from the posterior.: The average width of spatial variations appears unchanged from prior to posterior, interpreted as insufficient evidence to update beliefs about that width.This observation applies in both the observed and unobserved parts of the glacier.
  • 6.2. Classification of posterior covariance eigenvectors.: Eigenvectors are ordered by the difference of squared Rayleigh coefficients so data-informed and prior-informed directions group together.Positive values indicate data dominance, while negative values indicate prior dominance; four groups are distinguished using these features and spatial concentration.

Data-informed eigenvectors.

Data-informed eigenvectors correspond to directions where observations dominate prior information and posterior variance is substantially reduced, especially in the observed glacier region.

  • Data-informed eigenvectors.: Data-informed eigenvectors are concentrated primarily in the lower glacier, where observations are available.Their spatial concentration matches the region in which the data provide information about the sliding coefficient field.
  • Data-informed eigenvectors.: These eigenvectors are relatively smooth and resemble the first nine Fourier modes in the observed region.The first and third eigenvectors illustrate this behavior.
  • Data-informed eigenvectors.: Features of the MAP point lying in the span of the first nine Fourier modes therefore receive support as features of the true basal sliding coefficient field.The paper presents this correspondence as evidence-based confidence about those MAP features.

Shadowed eigenvectors.

Shadowed eigenvectors identify parameter directions where the parameter-to-observable map is insensitive, especially in poorly observed glacier regions, so the posterior is largely prior-controlled.

  • Shadowed eigenvectors.: Shadowed eigenvectors generally concentrate where the parameter-to-observable map is insensitive to the basal sliding coefficient.The upper glacier is far from observation points, and local influence of the coefficient limits sensitivity there.
  • Shadowed eigenvectors.: The upper glacier contains most shadowed eigenvectors because surface-velocity observations provide little information about its basal sliding coefficient.These eigenvectors also resemble Fourier modes in the upper region.
  • Shadowed eigenvectors.: One shadowed eigenvector is concentrated at the glacier’s bottom edge, where boundary-determined flow limits the coefficient’s influence despite available observations.
  • Shadowed eigenvectors.: Mixed eigenvectors differ from shadowed directions because both observations and the prior significantly influence them, making posterior behavior harder to predict.

Mixed eigenvectors.

Mixed eigenvectors are directions in which observational data and prior information both significantly influence the posterior.

  • Mixed eigenvectors.: Mixed eigenvectors are directions for which observations and the prior both have a significant influence.

Prior-tail eigenvectors.

Prior-tail directions are governed mainly by a highly certain prior because observations provide insufficient evidence to change it; this Gaussian-dominated behavior has interpretive limits.

  • Prior-tail eigenvectors.: Prior-tail eigenvectors have small prior variance, while observations provide insufficient evidence to contradict or reinforce the prior.
  • Prior-tail eigenvectors.: The continuous inverse problem contains infinitely many prior-tail directions, each behaving similarly to its prior counterpart.
  • Prior-tail eigenvectors.: A representative prior-tail eigenvector qualitatively resembles a high-frequency Fourier mode.
  • Prior-tail eigenvectors.: Posterior analysis based only on mean and covariance obscures non-Gaussian behavior and may misrepresent behavior away from the MAP point.
  • Prior-tail eigenvectors.: Prior-tail directions have small posterior variance and remain essentially unchanged from the prior.
  • Prior-tail eigenvectors.: Posterior marginals appear near-Gaussian in directions dominated by the prior, while data-informed directions can also be near-Gaussian under small-noise conditions.
  • Prior-tail eigenvectors.: In selected directions, posterior marginals are close to the MAP Gaussian approximation, but such agreement does not necessarily imply a Gaussian posterior.
  • Prior-tail eigenvectors.: The paper develops efficient MCMC methods for infinite-dimensional inverse problems with expensive forward models.

7. Concluding remarks.

The paper extends stochastic Newton MCMC to nonlinear infinite-dimensional inverse problems and evaluates a MAP-based Hessian variant for an ice-sheet application. It also interprets high-dimensional posterior structure through eigenvector classifications and marginal visualizations.

  • 7. Concluding remarks.: Figure 6.4 displays posterior marginals and corresponding sample variances along covariance eigenvectors sorted by qualitative feature groups.
  • 7. Concluding remarks.: The MAP-based stochastic Newton MCMC variant reuses a Hessian evaluated at the MAP point to avoid recomputing it at every sample.
  • 7. Concluding remarks.: For the nonlinear Stokes ice-sheet inverse problem, the MAP-based method achieves the best convergence in both sample count and PDE solves.
  • 7. Concluding remarks.: Posterior visualizations include point marginals and eigenvector classifications based on the relative influence of observations and the prior.
  • 7. Concluding remarks.: The classification identifies directions dominated by a Gaussian prior and directions where data-informed nonlinearity can produce non-Gaussian behavior.
  • 7. Concluding remarks.: Figure 6.5 compares selected one- and two-dimensional posterior marginals with the Gaussian approximation at the MAP point.
Loading 1308.6221v2…