Data-Driven Model Reduction for the Bayesian Solution of Inverse Problems
One of the major challenges in the Bayesian solution of inverse problems governed by partial differential equations (PDEs) is the computational cost of repeatedly evaluating numerical PDE models, as required by Markov chain Monte Carlo (MCMC) methods for posterior sampling. This paper proposes a data-driven projection-based model reduction technique to reduce this computational cost. The proposed technique has two distinctive features. First, the model reduction strategy is tailored to inverse problems: the snapshots used to construct the reduced-order model are computed adaptively from the posterior distribution. Posterior exploration and model reduction are thus pursued simultaneously. Second, to avoid repeated evaluations of the full-scale numerical model as in a standard MCMC method, we couple the full-scale model and the reduced-order model together in the MCMC algorithm. This maintains accurate inference while reducing its overall computational cost. In numerical experiments considering steady-state flow in a porous medium, the data-driven reduced-order model achieves better accuracy than a reduced-order model constructed using the classical approach. It also improves posterior sampling efficiency by several orders of magnitude compared to a standard MCMC method.
1Introduction and motivation
An important and challenging task in computational modeling is the solution of inverse problems, which convert noisy and indirect observational data into useful characterizations of the unknown parameters of a numerical model. In this process, statistical methods—Bayesian methods in particular—play a fundamental role in modeling various information sources and quantifying the uncertainty of the model parameters . In the Bayesian framework, the unknown parameters are modeled as random variables and hence can be characterized by their posterior distribution. Markov chain Monte Carlo (MCMC) methods  provide a powerful and flexible approach for sampling from posterior distributions. The Bayesian framework has been applied to inverse problems in various fields, for example, geothermal reservoir modeling , groundwater modeling , ocean dynamics , remote sensing , and seismic inversion .
To generate a sufficient number of samples from the posterior distribution, MCMC methods require sequential evaluations of the posterior probability density at many different points in the parameter space. Each evaluation of the posterior density involves a solution of the forward model used to define the likelihood function, which typically is a computationally intensive undertaking (e.g., the solution of a system of PDEs). In this many-query situation, one way to address the computational burden of evaluating the forward model is to replace it with a computationally efficient surrogate. Surrogate models have been applied to inverse problems in several settings; for example,  employed generalized polynomial chaos expansions,  employed Gaussian process models, and  used projection-based reduced-order models. In this work, we also focus on projection-based reduced-order models (although the data-driven strategy underlying our approach should be applicable to other types of surrogate models). A projection-based reduced-order model reduces the computational complexity of the original or “full” forward model by solving a projection of the full model onto a reduced subspace. For the model reduction methods we consider, the construction of the reduced subspace requires evaluating the full model at representative samples drawn from the parameter space. The solutions of the full model at these samples are referred to as snapshots . Their span defines the reduced subspace, represented via an orthogonal basis.
The quality of the reduced-order model relies crucially on the choice of the samples for computing the snapshots. To construct reduced-order models targeting the Bayesian solution of the inverse problem, we employ existing projection-based model reduction techniques. Our innovation is in a data-driven approach that adaptively selects samples from the posterior distribution for the snapshot evaluations. This approach has two distinctive features:
We integrate the reduced-order model construction process into an adaptive MCMC algorithm that simultaneously samples the posterior and selects posterior samples for computing the snapshots. During the sampling process, the numerical accuracy of the reduced-order model is adaptively improved.
The approximate posterior distribution defined by the reduced-order model is used to increase the efficiency of MCMC sampling. We either couple the reduced-order model and the full model together to accelerate the sampling of the full posterior distribution
1, or directly explore the approximate posterior distribution induced by the reduced-order model. In the latter case, sampling the approximate distribution yields a biased Monte Carlo estimator, but the bias can be controlled using error indicators or estimators.
Compared to the classical offline approaches that build the reduced-order model before using it in the many-query situation, the motivation for collecting snapshots during posterior exploration is to build a reduced-order model that focuses on a more concentrated region in the parameter space. Because the solution of the inverse problem is unknown until the data are available, reduced-order models built offline have to retain a level of numerical accuracy over a rather large region in the parameter space, which covers the support of the posterior distributions for all the possible observed data sets. For example, samples used for computing the snapshots are typically drawn from the prior distribution . In comparison, our data-driven approach focuses only on the posterior distribution resulting from a particular data set. As the observed data necessarily increase the information divergence of the prior distribution from the posterior distribution , the support of the posterior distribution is more compact than that of the prior distribution.
Figure 1 illustrates a simple two-dimensional inference problem, where the prior distribution and the posterior distribution are represented by the blue and red contours, respectively. The left plot of Figure 1 shows 50 randomly drawn prior samples for computing the snapshots, each of which has a low probability of being in the support of the posterior. In comparison, the samples selected by our data-driven approach, as shown in the right plot of Figure 1, are scattered within the region of high posterior probability.
By retaining numerical accuracy only in a more concentrated region, the data-driven reduced-order model requires a basis of lower dimension to achieve the same level of accuracy compared to the reduced-order model built offline. For the same reason, the data-driven model reduction technique can potentially have better scalability with parameter dimension than the offline approach.
We note that goal-oriented model reduction approaches have been developed in the context of PDE-constrained optimization , in which the reduced-order model is simultaneously constructed during the optimization process. In these methods, the snapshots for constructing the reduced basis are only evaluated at points in the parameter space that are close to the trajectory of the optimization algorithm.
The remainder of this paper is organized as follows. In Section 2, we outline the Bayesian framework for solving inverse problems and discuss the sampling efficiency and Monte Carlo error of MCMC. In Section 3, we introduce the data-driven model reduction approach, and construct the data-driven reduced-order model within an adaptive delayed acceptance algorithm to speed up MCMC sampling. We also provide results on the ergodicity of the algorithm. Section 4 analyzes some properties of the data-driven reduced-order model. In Section 5, we discuss a modified framework that adaptively constructs the reduced-order model and simultaneously explores the induced approximate posterior distribution. We also provide an analysis of the mean square error of the resulting Monte Carlo estimator. In Section 6, we demonstrate and discuss various aspects of our methods through numerical examples. Section 7 offers concluding remarks.
2Sample-based inference for inverse problems
The first part of this section provides a brief overview of the Bayesian framework for inverse problems. Further details can be found in . The second part of this section discusses the efficiency of MCMC sampling for computationally intensive inverse problems.
2.1Posterior formulation and sampling
Given a physical system, let denote the -dimensional unknown parameter, and denote the -dimensional observed data. The forward model maps a given parameter to the observable model outputs .
In a Bayesian setting the first task is to construct the prior models and the likelihood function as probability distributions. The prior density is a stochastic model representing knowledge of the unknown before any measurements and is denoted by . The likelihood function specifies the probability density of the observation for a given set of parameters , denoted by . We assume that the data and the model parameters follow the stochastic relationship
where the random vector captures the measurement noise and other uncertainties in the observation-model relationship, including model errors. Without additional knowledge of the measurement process and model errors, is modeled as a zero mean Gaussian, , where is the covariance. Let
denote the data-misfit function. The resulting likelihood function is proportional to By Bayes’ formula, the posterior probability density is
is the normalizing constant.
In the general framework set up by , we can explore the posterior distribution given by (Equation 2) using MCMC methods such as the Metropolis-Hastings algorithm . The Metropolis-Hastings algorithm uses a proposal distribution and an acceptance/rejection step to construct the transition kernel of a Markov chain that converges to the desired target distribution.
2.2Sampling efficiency and Monte Carlo error
Once we draw a sufficient number of samples from the posterior distribution, the expectations of functions over the posterior distribution can be estimated by Monte Carlo integration. Suppose we wish to estimate the expectation of a function over the posterior distribution
by posterior samples, . The resulting Monte Carlo estimator of is
which is an unbiased estimator with the mean square error
Because the samples drawn by an MCMC algorithm are correlated, the variance of is dependent on the effective sample size
is the integrated autocorrelation of ; see  for a detailed derivation and further references.
In practice, we wish to improve the computational efficiency of the MCMC algorithm, which is defined by the effective sample size for a given budget of CPU time. There are two ways to achieve this:
Improve the statistical efficiency. To increase the effective sample size for a given number of MCMC iterations, one can reduce the sample correlation by designing proposals that traverse the parameter space more efficiently. For example, adaptive MCMC methods  use adaptation to learn the covariance structure of the posterior distribution online, and use this information to design proposal distributions. Algorithms such as stochastic Newton  and Riemannian manifold MCMC  use local derivative information to construct efficient proposals that adapt to the local geometry of the posterior distribution.
Increase the number of MCMC steps for a given amount of CPU time. Because simulating the forward model in the data-misfit function (Equation 1) is CPU-intensive, computing the posterior density at every iteration is the major bottleneck of the standard MCMC algorithm. By using a fast approximation to the forward model, the amount of computing time can be reduced. This effort complements the use of MCMC algorithms that require local derivatives, since fast approximations of the forward model also enable fast evaluations of its gradient and even higher-order derivatives.
One important remark is that by drawing samples directly from an approximate posterior distribution (i.e., the one induced by our approximation of the forward model), the resulting Monte Carlo estimator (Equation 5) is biased. Algorithms such as surrogate transition  or delayed acceptance  MCMC can be used to couple the posterior distribution and its approximation together to ensure the samples are drawn from the full posterior distribution. On the other hand, if the accuracy of the posterior approximation can be tightly controlled, then some bias may be acceptable if it enables significant variance reduction, and thus an overall reduction in the mean squared error of an estimate of a posterior expectation, for a given computational effort. We explore both options in this work.
3Data-driven reduced-order model and full target algorithm
This section introduces our data-driven model reduction approach, and the adaptive sampling framework for simultaneously constructing the reduced-order model and exploring the posterior distribution.
Suppose the system of interest can be described by a system of nonlinear steady partial differential equations (PDEs), with a finite element or finite difference discretization that results in a discretized system in the form of
In Equation (Equation 7), represents the discretized state of the system, is the dimension of the system (the number of unknowns), is a discretized linear operator, represents the discretized nonlinear terms of the governing PDE, and denotes the forcing terms. All of , , and can be parameterized by the unknown parameter . An observation operator maps the state of the system to the observable model outputs, i.e.,
Our data-driven model reduction approach selects a set of posterior samples to compute the snapshots by solving Equation (Equation 7) for each sample . By constructing the reduced basis , the state can be approximated by a linear combination of the reduced basis vectors, i.e., . Then Equation (Equation 7) can be approximated by applying the Galerkin projection:
and the associated model outputs are
If , the dimension of the unknown state in (Equation 9) and (Equation 10) is greatly reduced compared to that of the original full system (Equation 7) and (Equation 8). Thus (Equation 9) and (Equation 10) define a reduced-order model that maps a realization of the unknown parameter to an approximation of the observable model outputs .
Care must be taken to ensure efficient solution of the reduced-order model, since in the presence of general parametric dependence, the equations (Equation 9) and (Equation 10) have low state dimension but are not necessarily fast to solve. This is because for each new parameter , solution of the reduced-order model requires evaluating the full scale system matrices or residual, projecting those matrices or residual onto the reduced basis, and then solving the resulting low dimensional system. Since many elements of these computations depend on , the dimension of the original system, this process is typically computationally expensive (unless there is special structure to be exploited, such as affine parametric dependence). To reduce the computational cost of this process, methods such as the missing point estimation , the empirical interpolation  and its discrete variant  approximate the nonlinear term in the reduced-order model by selective spatial sampling.
We note that for systems that exhibit a wide range of behaviors, recently developed localization approaches such as  adaptively construct multiple local reduced bases, each tailored to a particular system behavior that is associated with a subdomain of the state space or the parameter space. Our use of adaptation is different, focusing on the adaptive selection of posterior samples for evaluating the snapshots. In this paper we consider only the case of a global reduced basis; however, future work could combine our adaptive posterior sampling approach with a localization strategy to build multiple local reduced bases, each adapted to the local structure of the posterior.
The data-misfit function (Equation 1) can be approximated by replacing the forward model with the reduced-order model , which gives
Then the resulting approximate posterior distribution has the form
where is the normalizing constant.
Our data-driven approach adaptively selects the sample for evaluating the next snapshot such that the scaled error of the reduced-order model outputs,
for the current reduced basis , is above a user-specified threshold. The whitening transformation computes the relative error of the reduced-order model compared with noise level of the observed data, which has standard deviation 1 after the transformation. In case that we wish to bypass the forward model evaluations, so that the error cannot be computed directly, an error indicator or an error estimator, , can be used. In this work, we use the dual weighted residual technique  to compute an error indicator.
3.2Full target algorithm
To achieve simultaneous model reduction and posterior exploration, we employ the adaptive delayed acceptance algorithm of . Suppose we have a reduced-order model constructed from an initial reduced basis. At each iteration of the MCMC sampling, we first sample the approximate posterior distribution based on the reduced-order model for a certain number of steps using a standard Metropolis-Hastings algorithm. This first-stage subchain simulation should have sufficient number of steps, so that its initial state and last state are uncorrelated. Then the last state of the subchain is used as a proposal candidate in the second-stage of the algorithm, and the acceptance probability for this candidate is computed based on the ratio of the full posterior density value to the approximate posterior density value. This delayed acceptance algorithm employs the first-stage subchain to decrease the sample correlation by sampling the computationally fast approximate posterior, then uses the acceptance/rejection in the second-stage to ensure that the algorithm correctly samples the full posterior distribution.
The statistical efficiency of the delayed acceptance algorithm relies on the accuracy of the approximate posterior. An approximate posterior induced by an inaccurate reduced-order model can potentially result in a higher second-stage rejection rate, which decreases the statistical efficiency of the delayed acceptance algorithm. This is because duplicated MCMC samples generated by rejections increase the sample correlation (we refer to  for a formal justification). To maintain statistical efficiency, we aim to construct a sufficiently accurate reduced-order model, so that the resulting second-stage acceptance probability is close to . To achieve this, we introduce adaptive reduced basis enrichment into the delayed acceptance algorithm.
After each full posterior density evaluation, the state of the associated forward model evaluation can be used as a potential new snapshot. We compute the scaled error (Equation 13) of the reduced-order model outputs at each new posterior sample, and the reduced basis is updated with the new snapshot when the error exceeds a user-given threshold . By choosing an appropriate threshold , we can control the maximum allowable amount of error of the reduced-order model.
The resulting reduced-order model is data-driven, since it uses the information provided by the observed data (in the form of the posterior distribution) to select samples for computing the snapshots. It is also an online model reduction approach, since the reduced-order model is built concurrently with posterior sampling.
Since the adaptive sampling procedure samples from the full posterior distribution, hereafter we refer to it as the “full target algorithm.” Details of the full target algorithm are given in Algorithm ?.
Lines 1–13 of Algorithm ? simulate a Markov chain with invariant distribution for steps. In Lines 14–19 of Algorithm ?, to ensure that the algorithm samples from the full posterior distribution, we post-process the last state of the subchain using a second-stage acceptance/rejection based on the full posterior distribution.
The second-stage acceptance probability is controlled by the accuracy of the reduced-order model and the subchain length . Ideally, we want to have large to give uncorrelated samples. However, if the second-stage acceptance probability is low, the effort spent on simulating the subchain will be wasted because the proposal is more likely to be rejected in the second step. To avoid this situation, we dynamically monitor the accuracy of the reduced-order model by evaluating its error indicator at each state of the subchain (Lines 10–12). The subchain is terminated if the L-infinity norm of the scaled error indicator exceeds a threshold .
Lines 20–22 of Algorithm ? describe the adaptation of the reduced basis, which is controlled by the scaled error (Equation 13) of the reduced-order model and a given threshold . Three criteria must be satisfied for the reduced basis to be updated: (1) the scaled error at must exceed the threshold ; (2) the dimensionality of the reduced basis must not exceed the maximum allowable dimensionality ; and (3) the finite adaptation criterion (Definition ?) should not yet be satisfied. The finite adaptation criterion is defined precisely as follows.
The finite adaptation criterion is a threshold for how infrequent updates to the reduced-order model should become before model adaptation is stopped entirely. When more and more MCMC steps are used between updates, the reduced-order model has satisfied the accuracy threshold over larger and larger regions of the posterior. The threshold of at least steps between updates can thus be understood as a measure of “sufficient coverage” by the reduced-order model. Further discussion and justification of the finite adaptation criterion is deferred to Section 4.
Algorithm ? uses the adaptation in Lines 20–22 to drive the online construction of the reduced-order model. Once the adaptation stops, the algorithm runs as a standard delayed acceptance algorithm, with an -step subchain in the first stage.
To prevent the situation where the reduced basis dimension is large so that the reduced-order model becomes computationally inefficient, we limit the reduced basis dimension to be less than a user-given threshold . When the reduced basis dimension reaches but the finite adaptation criterion is not satisfied, one strategy is to stop the reduced basis enrichment, and continue simulating Algorithm ?. This does not affect the ergodicity of Algorithm ?, as will be discussed in Section ?. However, the error of the resulting reduced-order model can be large in the subregion of the support of the posterior that remains unexplored. Consequently, this large reduced-order model error may decrease the second-stage acceptance rate and the statistical efficiency of the algorithm.
A wide range of proposal distributions can be used in Algorithm ?. In this work, we use the grouped-components adaptive Metropolis , which is a variant of adaptive Metropolis . More advanced algorithms such as stochastic Newton  or manifold MCMC  can also be used within Algorithm ?. The computational cost of evaluating derivative information for these algorithms can also be reduced by using the reduced-order models.
Ergodicity of the full target algorithm.
Throughout this work, the following assumption on the forward model and the reduced-order model is used:
Since the data-misfit functions and are quadratic, we have and . Assumption ? implies that there exists a constant such that We also note that Lipschitz continuity of implies the Lipschitz continuity of . Similarly, is Lipschitz continuous.
We first establish the ergodicity of a non-adaptive version of Algorithm ?, where the reduced-order model is assumed to be given and fixed.
The detailed balance condition and aperiodicity condition are satisfied by Algorithm ?, see . Since for all by Assumption ?, we have for all , thus the irreducibility condition is satisfied, see . Ergodicity follows from the detailed balance, aperiodicity, and irreducibility.
The adaptation used in Algorithm ? is different from that of standard adaptive MCMC algorithms such as  and the original adaptive delayed acceptance algorithm . These algorithms carry out an infinite number of adaptations because modifications to the proposal distribution or to the approximation continue throughout MCMC sampling. The number of adaptations in Algorithm ?, on the other hand, cannot exceed the maximum allowable dimension of the reduced basis, and hence it is finite. Furthermore, the adaptation stops after finite time because of the finite adaptation criterion. Since Algorithm ? is ergodic for each of the reduced-order models constructed by Lemma ?, Proposition 2 in  ensures the ergodicity of Algorithm ?. Lemma ? also reveals that Algorithm ? always converges to the full posterior distribution, regardless of the accuracy of the reduced-order model.
4Properties of the approximation
In Algorithm ?, the accuracy of the reduced-order model is adaptively improved during MCMC sampling, and thus it is important to analyze the error of the approximate posterior induced by the reduced-order model, compared to the full posterior. Since bounds on the bias of the resulting Monte Carlo estimator can be derived from the error of the approximate posterior, this analysis is particularly useful in situations where we want to utilize the approximate posterior directly for further CPU time reduction. This error analysis also justifies the finite adaptation criterion (see Definition ?) used in Algorithm ? to terminate the enrichment of the reduced basis.
We provide here an analysis of the Hellinger distance between the full posterior distribution and the approximate posterior distribution
The Hellinger distance translates directly into bounds on expectations, and hence we use it as a metric to quantify the error of the approximate posterior distribution.
The framework set by  is adapted here to analyze the Hellinger distance between the full posterior distribution and its approximation induced by the constructed in Algorithm ?, with respect to a given threshold . Given a set of samples where the snapshots are computed, and the associated reduced basis , we define the -feasible set and the associated posterior measure:
On any measurable subset of the -feasible set, the error of the reduced-order model is bounded above by . Then the Hellinger distance (Equation 14) can be bounded by the user-specified and the posterior measure of , which is the region that has not been well approximated by the reduced-order model. We formalize this notion though the following propositions and theorems.
This result directly follows from the definition of the -feasible set and Assumption ?.
From the definition of the -feasible set, we have a bound on the difference of the normalizing constants:
Since and are Lipschitz continuous and greater than zero, we have
for some constant . Then by Proposition ?, there exists a constant such that
There also exists a constant such that
Thus, we have .
Following Theorem 4.6 of  we have
Following the same derivation as Proposition ? we have
Thus there exist constants such that
Applying the bound
and Proposition ?, we have
for some constants .
Combining the above results we have
Since and , the above inequality can be rearranged as
for some constants and .
Under certain technical conditions , the pointwise error of the reduced-order model decreases asymptotically as the dimensionality of the reduced basis increases. Thus the -feasible set asymptotically grows with the reduced basis enrichment, and hence asymptotically decays. If the posterior distribution is sufficiently well-sampled such that
then the Hellinger distance (Equation 14) is characterized entirely by , as shown in Theorem ?. Thus, by adaptively updating the data-driven reduced-order model until the condition (Equation 16) is satisfied, we can build an approximate posterior distribution whose error is proportional to the user-specified error threshold .
In practice, it is not feasible to check the condition (Equation 16) within MCMC sampling, so we use heuristics to motivate the finite adaptation criterion in Definition ?. Consider a stationary Markov chain that has invariant distribution . We assume that the probability of visiting is proportional to its posterior measure. Suppose we have , then the probability of the Markov chain visiting is about , for some . In this case, the average number of MCMC steps needed for the next reduced basis refinement is about . Since the posterior measure of decays asymptotically with refinement of the reduced basis, the average number of steps needed for the basis refinement asymptotically increases. Thus we can treat the condition (Equation 16) as holding if the average number of steps used for reduced basis refinement exceeds , and thus terminate the adaptation. We recommend to choose a small value, for example , to delay the stopping time of the adaptation. In this way, the adaptive construction process can search the parameter space more thoroughly to increase the likelihood that the condition is satisfied.
5Extension to approximate Bayesian inference
The full target algorithm introduced in Section 3.2 has to evaluate the forward model after each subchain simulation to preserve ergodicity. The resulting Monte Carlo estimator (Equation 5) is unbiased. However, the effective sample size (ESS) for a given budget of computing time is still characterized by the number of full model evaluations. To increase the ESS, we also consider an alternative approach that directly samples from the approximate posterior distribution when the reduced-order model has sufficient accuracy.
Suppose that the scaled error indicator provides a reliable estimate of the scaled true error of the reduced-order model. Then the reliability and the refinement of the reduced-order model can be controlled by the error indicator. During MCMC sampling, we only evaluate the full model to correct a Metropolis acceptance and to update the reduced-order model when the error indicator exceeds the threshold . When a reduced-order model evaluation has error indicator less than , we treat the reduced-order model as a sufficient approximation to the full model. In this case, decisions in the MCMC algorithm are based on the approximate posterior distribution.
Compared to the full target algorithm that draws samples from the full posterior distribution, this approach only samples from an approximation that mimics the full posterior up to the error threshold . We refer to the proposed approach as the “-approximate algorithm.” Even though sampling the approximate posterior distribution will introduce bias to the Monte Carlo estimator (Equation 5), this bias may be acceptable if the resulting Monte Carlo estimator has a smaller mean squared error (MSE) for a given amount of computational effort, compared to the standard Metropolis-Hasting algorithm. In the remainder of this section, we provide details on the algorithm and analyze the bias of the resulting Monte Carlo estimator.
Algorithm ? details the -approximate algorithm. During the adaptation, for a proposed candidate , we discard the reduced-order model evaluation if its error indicator exceeds the new upper threshold . We set , which means that Algorithm ? does not use the information from the reduced-order model if its estimated error is greater than one standard deviation of the measurement noise. In this case (Lines 3–6), we run the full model directly to evaluate the posterior distribution and the acceptance probability. If the error indicator of the reduced-order model is between the lower threshold and the upper threshold (Lines 7–15), the reduced-order model is considered to be a reasonable approximation, and the delayed acceptance scheme is used to make the correction. If the error indicator is less than the lower threshold , or if the adaptation is stopped (Lines 17–19), the reduced-order model is considered to be a sufficiently accurate approximation to the full model, and is used to accept/reject the proposal directly.
The adaptation criterion used in Lines 2 of Algorithm ? has two conditions: the dimension of reduced basis should not exceed a specified threshold , and the finite adaptation criterion (Definition ?) should not yet be triggered. The reduced basis is updated using the full model evaluations at proposals accepted by MCMC.
When the reduced basis dimension reaches but the finite adaptation criterion is not satisfied, it is not appropriate to use the -approximate algorithm for the prescribed error threshold. This is because the large reduced-order model error can potentially result in unbounded bias in the Monte Carlo estimator. In this situation, we should instead use the full target algorithm, for which convergence is guaranteed.
5.2Monte Carlo error of the -approximate algorithm
To analyze the performance of the -approximate algorithm, we compare the MSE of the resulting Monte Carlo estimator with that of a standard single-stage MCMC algorithm that samples the full posterior distribution.
We wish to compute the expectation of a function over the posterior distribution , i.e.,
where the first and second moments of are assumed to be finite. Suppose a single-stage MCMC algorithm can sample the full posterior distribution for steps in a fixed amount of CPU time, and effective samples are produced. The resulting Monte Carlo estimator
which is characterized by the ESS and the variance of over .
By sampling the approximate posterior distribution, the expectation can be approximated by
Suppose we can sample the approximate posterior for steps in the same amount of CPU time as sampling the full posterior, and that these samples have effective sample size
where is the speedup factor that depends on the computational expense of the reduced-order model . The Monte Carlo estimator
has the MSE
where the bias is defined by
We are interested in the situation for which sampling the approximation leads to a smaller MSE than sampling the full posterior distribution, i.e.,
Rearranging the above inequality gives
Equation (Equation 17) reveals that when our target ESS for drawing from the full posterior distribution does not exceed the bound , sampling the approximate posterior will produce a smaller MSE. This suggests that the -approximate algorithm will be more accurate for a fixed computational cost than the single-stage MCMC when the target ESS of the single-stage MCMC satisfies (Equation 17). In such cases, the MSE of the Monte Carlo estimator is dominated by the variance rather than the bias.
The bias is characterized by the Hellinger distance between the full posterior distribution and its approximation (Lemma 6.37 of ), i.e.,
We assume that there exists an and a set of samples such that the resulting reduced-order model satisfies the condition (Equation 16), i.e., . Applying Theorem ?, the ratio of variance to squared bias can be simplified to
for some constant , and hence we have
For problems that have reliable error indicators or estimators, the -approximate algorithm provides a viable way to select the set of samples for computing the snapshots. However it is computationally infeasible to verify that the condition (Equation 16) holds in practice. We thus employ the finite adaptation criterion (Definition ?) to perform a heuristic check on the condition (Equation 16), as discussed in Section 4.
The bound is characterized by the user-given error threshold and the speedup factor , where is a problem-specific factor that is governed by the convergence rate and computational complexity of the reduced-order model. For a reduced-order model such that , there exists a so that the MSE of sampling the approximate posterior for steps will be less than the MSE of sampling the full posterior for the same amount of CPU time. In the regime where the reduced-order models have sufficiently large speedup factors, the bound is dominated by , and hence decreasing results in a higher bound . However there is a trade-off between the numerical accuracy and speedup factors. We should avoid choosing a very small value, because this can potentially lead to a high-dimensional reduced basis and a correspondingly expensive reduced-order model such that , where the ratio should be close to one for such an accurate reduced-order model. In this case, sampling the approximate posterior can be less efficient than sampling the full posterior.
6Numerical results and discussion
To benchmark the proposed algorithms, we use a model of isothermal steady flow in porous media, which is a classical test case for inverse problems.
Let be the problem domain, be the boundary of the domain, and denote the spatial coordinate. Let be the unknown permeability field, be the pressure head, and represent the source/sink. The pressure head for a given realization of the permeability field is governed by
where the source/sink term is defined by the superposition of four weighted Gaussian plumes with standard deviation , centered at , , , , and with weights . A zero-flux Neumann boundary condition
is prescribed, where is the outward normal vector on the boundary. To make the forward problem well-posed, we impose the extra boundary condition
In Section 6.2, we use a nine-dimensional example to carry out numerical experiments to benchmark various aspects of our algorithms. In this example, the spatially distributed permeability field is projected onto a set of radial basis functions, and hence inference is carried out on the weights associated with each of the radial basis functions. In Section 6.3, we apply our algorithms to a higher dimensional problem, where the parameters are defined on the computational grid and endowed with a Gaussian process prior. Both examples use fixed Gaussian proposal distributions, where the covariances of the proposals are estimated from a short run of the -approximate algorithm. In Section 6.4, we offer additional remarks on the performance of the -approximate algorithm.
6.2The 9D inverse problem
The permeability field is defined by radial basis functions:
where are the centers of the radial basis functions. These radial basis functions are shown in Figure 2. The prior distributions on each of the weights , are independent and log-normal, and hence we have
where and = 9. The true permeability field used to generate the test data, and the corresponding pressure head are shown in Figure 3. The measurement sensors are evenly distributed over with grid spacing , and the signal-to-noise ratio of the observed data is .
Numerical experiments for various choices of are carried out to test the computational efficiency of both algorithms and the dimensionality of reduced basis in the reduced-order model. For , we run the full target algorithm for iterations, with subchain length . To make a fair comparison in terms of the number of posterior evaluations, we run the -approximate algorithm for iterations, also with . For both algorithms, the reduced-order model construction process is started at the beginning of the MCMC simulation. We set in the finite adaptation criterion ?. As a reference, we run a single-stage MCMC algorithm for iterations using the same proposal distribution. The reference algorithm only uses the full posterior distribution. The first samples of the simulations generated by the full target algorithm are discarded as burn-in samples. Similarly, the first samples (to match the number of reduced-order model evaluations in the full target algorithm) are discarded for the simulations using the -approximate algorithm and the reference algorithm.
In the remainder of this subsection, we provide various benchmarks of our data-driven model reduction approach and the two algorithms. These include:
A comparison of the full target algorithm and the -approximate algorithm with the reference algorithm.
A comparison of the data-driven reduced-order model with the reduced-order model built with respect to the prior distribution.
A demonstration of the impact of observed data on our data-driven reduced-order model.
Table 1 summarizes the number of full model evaluations, the dimensionality of reduced basis, the CPU time, the ESS, and the speedup factor, comparing the full target algorithm and the -approximate algorithm with the reference algorithm. For the reduced-order models generated by the adaptive construction process, we provide estimates of the posterior measure of the complement of the -feasible set, . We also provide a summary of the average second-stage acceptance probability, , for the full target algorithm.
For the full target algorithm, the average second-stage acceptance probabilities for all three values are greater than in this test case. This shows that the reduced-order models produced by all three values are reasonably accurate compared to the full model, and hence simulating the approximate posterior distribution in the first stage usually yields the same Metropolis acceptance decision as simulating the full posterior distribution. As we enhance the accuracy of the reduced-order model by decreasing the value of , the dimensionality of the resulting reduced basis increases, and thus the reduced-order model takes longer to evaluate. Because the full target algorithm evaluates the full model for every reduced-order model evaluations, its computational cost is dominated by the number of full model evaluations. Thus the speedup factors for all three choices of are similar (approximately 40). Since all three reduced-order models are reasonably accurate here, the efficiency gain of using a small value is not significant. In this situation, one could consider simulating the subchain in the first stage for more iterations (by increasing the subchain length ) when the value of is small.
The -approximate algorithm produces speedup factors that are to times higher than the speedup factor of the full target algorithm in this test case. A larger value produces a larger speedup factor, because the dimension of the associated reduced basis is smaller.
To assess the sampling accuracy of the -approximate algorithm, Figure 4 provides a visual inspection of the marginal distributions of each component of the parameter , and the contours of the marginal distributions of each pair of components. The black lines represent the results generated by the reference algorithm, the blue lines represent results of the full target algorithm with , and red lines represent results of the -approximate algorithm with . The results from the more accurate simulations that use smaller values are not shown, as they are visually close to the case . The plots in Figure 4 suggest that all the algorithms generate similar marginal distributions in this test case. We note that both the reference algorithm and the full target algorithm sample from the full posterior distribution, and thus the small differences in the contours produced by various algorithms are likely caused by Monte Carlo error.
An alternative way to assess the computational efficiency and sampling accuracy of the -approximate algorithm is to compare the number of effective samples generated by the -approximate algorithm and the reference algorithm for a fixed amount of CPU time. As shown in Table 1, the -approximate algorithm with generates effective samples in seconds; the reference algorithm can only generate about effective samples in the same amount of CPU time. In the situation where the desired number of effective samples is at least an order of magnitude larger than the speedup factor, using the -approximate algorithm is clearly advantageous to using the reference algorithm.
For both the full target algorithm and the -approximate algorithm, and for each choice of , we use samples generated by the reference algorithm to compute the Monte Carlo estimator of the posterior measure of the complement of the -feasible set for the final reduced-order model produced by the adaptive construction process. As shown in Table 1, for all the reduced-order models, we have estimated . This suggests that the Hellinger distances between the full posterior distribution and its approximation can be characterized by the values in all three cases, and thus the finite adaptation criterion (Definition ?) with provides a useful indicator for terminating adaptation.
For , we note that the dimensions of the reduced bases produced by the full target algorithm and the -approximate algorithm are different. This is because the snapshots are evaluated at selected samples that are randomly drawn from the posterior. The spread of the sample set slightly affects the accuracy of the reduced-order model. Nonetheless, both reduced-order models achieve the desirable level of accuracy since the estimated posterior measures are less than in this case.
Numerical experiments with and in the finite adaptation criterion ? are also conducted. For both algorithms, choosing these smaller values leads only to one or two additional basis vectors being added in all the test cases, compared to the case . The resulting marginal distributions generated by using and are similar to the case . For brevity, the sampling results for these experiments are not reported. We consistently observe that the number of MCMC steps between adjacent basis enrichments increases as the adaptive construction progresses in these experiments. This is expected since the posterior measure asymptotically decreases with reduced basis enrichment. In this situation, choosing a smaller value leads only to minor increases in both of the numerical accuracy and the computational cost of the reduced-order model. Thus the sampling accuracy and the overall computational load of both sampling algorithms are not sensitive to the smaller values in this case.
Comparison with a reduced-order model built from the prior.
Now we compare the accuracy of the data-driven reduced-order model built with to that of a reduced-order model constructed with respect to the prior distribution (Equation 22). To construct the reduced-order model with respect to the prior, we use proper orthogonal decomposition (POD). random prior samples are drawn to compute the snapshots. The POD eigenspectrum is shown in the left plot of Figure 5. The eigendecomposition is truncated when it captures all but energy (relative 2-norm of the POD eigenvalues), leading to 110 reduced basis vectors being retained in the POD basis.
By using posterior samples generated from a separate reference algorithm, we compute the expectation of the norm of the scaled true error (Equation 13) over the full posterior distribution. The norm of the scaled true error gives the worst-case sensor error; its expectation over the the posterior quantifies the average numerical accuracy of the resulting reduced-order model. The right plot of Figure 5 shows this expectation with respect to the dimension of the reduced basis. For this test case, the data-driven reduced-order model undergoes a significant accuracy improvement once it includes at least 10 reduced basis vectors. The figure shows that the data-driven reduced-order model has a better convergence rate compared to the reduced-order model built from the prior.
The influence of posterior concentration.
The amount of information carried in the data affects the dimension of the data-driven reduced-order model, and hence has an impact on its computational efficiency. By adjusting the signal-to-noise ratio in the observed data, we examine the influence of the posterior concentration on the dimension of the reduced basis. We gradually increase the signal-to-noise ratio from 10 to 100, and record the number of reduced basis vectors in the reduced-order models. To quantify the posterior concentration, we use the “tightness” of the posterior distribution defined by
where is the standard deviation of the posterior marginal of , and is the standard deviation of the corresponding prior marginal. In Figure 6, we observe that the dimension of the reduced basis decreases as the signal-to-noise ratio increases. For this test problem, the larger amount of information in the data results in a lower dimensional reduced basis because our approach exploits the increasing concentration of the posterior.
6.3The high dimensional inverse problem
In the high-dimensional example, a log-normal distribution is employed to model the permeabilities as a random field. Let denote the coordinates of the grid points. Let be the permeability field defined at each grid point, then the latent field follows a Gaussian process prior
where is used to provide sufficient spatial variability. After applying the eigendecomposition of the prior covariance, the parameters are defined on eigenvectors that preserve energy of the prior distribution. To avoid an inverse crime, we use a “true” permeability field that is not directly drawn from the prior distribution. Figure 7 shows the true permeability field and the simulated pressure head. The setup of the measurement sensors is the same as the 9D example in Section 6.2.
Using the same setting as the 9D case, we simulate the full target algorithm with for iterations, with subchain length . For these full target MCMC simulations, the first samples are discarded as burn-in samples. We simulate the -approximate algorithm with , for iterations. The single-stage MCMC method is simulated for iterations as the reference. For all the -approximate MCMC simulations and the reference MCMC simulation, the first samples are discarded as burn-in samples.
Table 2 summarizes the number of full model evaluations, the number of reduced basis vectors, the CPU time, ESS, and speedup factor. The speedup factor of the full target algorithm is about 67. In comparison, the speedup factors of the -approximate algorithm range from to . The speedup factor increases as the error threshold increases. Figure 8 shows the mean and standard deviation at each spatial location of the permeability field, generated from the reference algorithm, the full target algorithm, and the least accurate setting () of the -approximate algorithm. We observe that all algorithms produce similar estimates of mean and standard deviation in this test case.
The values estimated from samples generated by the reference algorithm for all three values are also recorded in Table 2. In this test example, we have for all three values, and thus the Monte Carlo estimator provided by the -approximate algorithm can be characterized by the values. We note that some of the estimated posterior measures have zero values in Table 2, but these values do not necessarily mean that the posterior measures are exactly zero, because these are Monte Carlo estimates.
6.4Remarks on the -approximate algorithm
In the first case study, the -approximate algorithm offers some speedup compared to the full target algorithm (range from to ). In the second case study, the speedup factor of the -approximate algorithm compared to the full target algorithm drops to at most (with ) and it performs slightly worse than the full target algorithm for . The speedup factor of the -approximate algorithm decreases with in both cases. This result is to be expected, since the computational cost of the reduced-order model depends on the the dimensionality of the reduced basis, which grows as decreases. For in the second test case, the reduced-order model becomes computationally too expensive relative to the full model, and thus we lose the efficiency gain of the -approximate algorithm over the full target algorithm.
For problems that require only a limited number of effective samples, using the -approximate algorithm can be more advantageous. This is because we can use a relatively large value to keep the computational cost of the reduced-order model low, while the resulting MSE is still dominated by the variance of the estimator rather than the bias. If the goal is to obtain an accurate Monte Carlo estimator, in which the variance of the estimator is small compared to the bias resulting from sampling the approximate posterior distribution, we should use the full target algorithm. We also note that the accuracy and efficiency of the -approximate algorithm depend on reliable error indicators or error estimators, while the full target algorithm always samples the full posterior distribution, regardless of the error indicator.
We have introduced a new data-driven model reduction approach for solving statistical inverse problems. Our approach constructs the reduced-order model using adaptively-selected posterior samples to compute the snapshots. The reduced-order model construction process is integrated into the posterior sampling, to achieve simultaneous posterior exploration and model reduction.
Based on the data-driven reduced-order model, we have also developed two MCMC algorithms to sample the posterior distribution more efficiently than standard full-model MCMC algorithms. The full target algorithm aims to accelerate sampling of the full posterior distribution by coupling the full and approximate posterior distributions together. The -approximate algorithm samples an approximate posterior distribution, and attempts to reduce the mean squared error of the resulting Monte Carlo estimator, compared to a standard MCMC algorithm. Both algorithms adaptively construct the reduced-order model online through the MCMC sampling. The full target algorithm preserves ergodicity with respect to the true posterior. The -approximate algorithm does not sample the full posterior, but can provide further speedups for some problems.
In the case studies, we have demonstrated that both algorithms are able to accelerate MCMC sampling of computationally expensive posterior distributions by up to two orders of magnitude, and that the sampling accuracy of the -approximate algorithm is comparable to that of a reference full-model MCMC. We have also used the first case study to show the numerical accuracy of the data-driven reduced-order model, compared to a reduced-order model that is built offline with respect to the prior distribution. In this example, for the same number of reduced basis vectors, the posterior-averaged output error of the data-driven reduced-order model is several orders of magnitude smaller than that of the reduced-order model built with respect to the prior. Furthermore, we have demonstrated the impact of the amount of information carried in the observed data on the dimensionality of the reduced basis.
For solving statistical inverse problems, these results suggest that a data-driven reduced-order model is preferable to a reduced-order model built with respect to the prior, especially when the data are informative. Even though our approach is designed for constructing projection-based reduced-order models, the concept of building posterior-oriented surrogate models can be generalized to other approximation approaches such as Gaussian process regression and generalized polynomial chaos.
The authors thank Florian Augustin, Tan Bui-Thanh, Omar Ghattas, and Jinglai Li for many helpful comments and discussions. This work was supported by the United States Department of Energy Applied Mathematics Program, Awards DE-FG02-08ER2585 and DE-SC0009297, as part of the DiaMonD Multifaceted Mathematics Integrated Capability Center.
- The term “full posterior” refers to the posterior distribution induced by the original or full forward model.
- Using approximate derivatives does not impact the ergodicity of the MCMC algorithm, since the bias caused by the approximate derivatives is corrected by the acceptance/rejection step in the Metropolis-Hastings algorithm.
- Kaipio JP, Somersalo E. Statistical and computational inverse problems, vol. 160. Springer: New York, 2004.
- Tarantola A. Inverse Problem Theory and Methods for Model Parameter Estimation. Society for Industrial Mathematics: Philadelphia, 2005.
- Liu JS. Monte Carlo strategies in scientific computing. Springer: New York, 2001.
- Cui T, Fox C, O’Sullivan MJ. Bayesian calibration of a large-scale geothermal reservoir model by a new adaptive delayed acceptance Metropolis Hastings algorithm. Water Resource Research 2011; 47:W10 521.
- Higdon D, Lee H, Holloman C. Markov chain monte carlo-based approaches for inference in computationally intensive inverse problems. Bayesian Statistics 7, Bernardo JM, Bayarri MJ, Berger JO, Dawid AP, Heckerman D, Smith AFM, West M (eds.). Oxford University Press, Oxford 2003; 181–197.
- McKeague IW, Nicholls GK, Speer K, Herbei R. Statistical inversion of south atlantic circulation in an abyssal neutral density layer. Journal of Marine Research 2005; 63:683–704.
- Haario H, Laine M, Lehtinen M, Saksman E, Tamminen J. Markov chain Monte Carlo methods for high dimensional inversion in remote sensing. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 2004; 66:591–608.
- Martin J, Wilcox LC, Burstedde C, Ghattas O. A stochastic newton MCMC method for large-scale statistical inverse problems with application to seismic inversion. SIAM Journal on Scientific Computing 2012; 34(3):A1460–A1487.
- Marzouk YM, Najm HN, Rahn LA. Stochastic spectral methods for efficient bayesian solution of inverse problems. Journal of Computational Physics 2007; 224:560–586.
- Marzouk YM, Najm HN. Dimensionality reduction and polynomial chaos acceleration of Bayesian inference in inverse problems. Journal of Computational Physics 2009; 228:1862–1902.
- Bayarri MJ, Berger J, Kennedy M, Kottas A, Paulo R, Sacks J, Cafeo J, Lin C, Tu J. Predicting vehicle crashworthiness: Validation of computer models for functional and hierarchical data. Journal of the American Statistical Association 2009; 104:929–943.
- Galbally D, Fidkowski K, Willcox KE, Ghattas O. Nonlinear model reduction for uncertainty quantification in large scale inverse problems. International journal for numerical methods in engineering 2008; 81(12):1581–1608.
- Lipponen A, Seppänen A, Kaipio JP. Electrical impedance tomography imaging with reduced-order model based on proper orthogonal decomposition. Journal of Electronic Imaging 2013; 22:023 008.
- Lieberman C, Willcox KE, Ghattas O. Parameter and state model reduction for large-scale statistical inverse problems. SIAM Journal on Scientific Computing 2010; 32(5):2523–2542.
- Wang J, Zabaras N. Using bayesian statistics in the estimation of heat source in radiation. International Journal of Heat and Mass Transfer 2005; 48:15–29.
- Sirovich L. Turbulence and the dynamics of coherent structures. part 1: Coherent structures. Quarterly of Applied Mathematics 1987; 45:561–571.
- Cover TM, Thomas JA. Elements of information theory. Wiley-interscience, New York 2006.
- Arian E, Fahl M, Sachs EW. Trust-region proper orthogonal decomposition for flow control. Technical Report ICASE-2000-25, Institute for computer applications in science and engineering, Hampton, VA 2000.
- Ravindran SS. Adaptive reduced-order controllers for a thermal flow system using proper orthogonal decomposition. SIAM Journal on Scientific Computing 2002; 23(6):1924–1942.
- Kunisch K, Volkwein S. Proper orthogonal decomposition for optimality systems. ESAIM: Mathematical Modelling and Numerical Analysis 2008; 42(1):1–23.
- Carlberg K, Farhat C. A compact proper orthogonal decomposition basis for optimization-oriented reduced-order models. AIAA Paper 2008; 5964:10–12.
- Stuart AM. Inverse problems: a Bayesian perspective. Acta Numerica 2010; 19:451–559.
- Grenander U, Miller MI. Representations of knowledge in complex systems. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 1994; 56(4):549–603.
- Hastings W. Monte Carlo sampling using Markov chains and their applications. Biometrika 1970; 57:97–109.
- Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E. Equation of state calculations by fast computing machines. Journal of chemical physics 1953; 21:1087–1092.
- Haario H, Saksman E, Tamminen J. An adaptive Metropolis algorithm. Bernoulli 2001; 7(2):223–242.
- Roberts GO, Rosenthal JS. Coupling and ergodicity of adaptive Markov chain Monte Carlo algorithms. Journal of applied probability 2007; 44(2):458–475.
- Girolami M, Calderhead B. Riemann manifold Langevin and Hamiltonian Monte Carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 2011; 73(2):123–214.
- Liu JS, Chen R. Sequential Monte Carlo methods for dynamic systems. Journal of the American statistical association 1998; 93(443):1032–1044.
- Christen JA, Fox C. MCMC using an approximation. Journal of Computational and Graphical statistics 2005; 14(4):795–810.
- Astrid P, Weiland S, Willcox KE, Backx T. Missing point estimation in models described by proper orthogonal decomposition. IEEE Transactions on Automatic Control 2008; 53(10):2237–2251.
- Barrault M, Maday Y, Nguyen NC, Patera AT. An “empirical interpolation” method: application to efficient reduced-basis discretization of partial differential equations. Comptes Rendus Mathematique 2004; 339(9):667–672.
- Chaturantabut S, Sorensen DC. Nonlinear model reduction via discrete empirical interpolation. SIAM Journal on Scientific Computing 2010; 32(5):2737–2764.
- Haasdonk B, Dihlmann M, , Ohlberger M. A training set and multiple bases generation approach for parameterized model reduction based on adaptive grids in parameter space. Math. Comput. Model. Dyn. Syst. 2011; 17:423–442.
- Amsallem D, Zahr M, Farhat C. Nonlinear model order reduction based on local reduced-order bases. International Journal Numerical Methods Engineering 2012; 92:891–916.
- Eftang JL, Stamm B. Parameter multi?domain ‘hp’ empirical interpolation. International Journal for Numerical Methods in Engineering 2012; 90(4):412–428.
- Peherstorfer B, Butnaru D, Willcox KE, Bungartz HJ. Localized discrete empirical interpolation method. SIAM Journal on Scientific Computing 36; 2014(1):A168–A192.
- Meyer M, Matthies HG. Efficient model reduction in non-linear dynamics using the karhunen-loeve expansion and dual-weighted-residual methods. Computational Mechanics 2003; 31(1):179–191.
- Cui T. Bayesian calibration of geothermal reservoir models via markov chain monte carlo. PhD Thesis, The University of Auckland, Auckland 2010.
- Bui-Thanh T, Willcox KE, Ghattas O. Model reduction for large-scale systems with high-dimensional parametric input space. SIAM Journal on Scientific Computing 2008; 30(6):3270–3288.
- Patera AT, Rozza G. Reduced basis approximation and a posteriori error estimation for parametrized partial differential equations. MIT Pappalardo monographs in mechanical engineering (to appear), Copyright MIT (2006–2007), 2007.