Search

In a typical inverse problem, a spatially distributed parameter in a physical model is estimated from measurements of model output. Since measurements are stochastic in nature, so is any parameter estimate. Moreover, in the Bayesian setting, the choice of regularization corresponds to the definition of the prior probability density function, which in turn is an uncertainty model for the unknown parameters. For both of these reasons, significant uncertainties exist in the solution of an inverse problem. Thus to fully understand the solution, quantifying these uncertainties is important. When the physical model is linear and the error model and prior are Gaussian, the posterior density function is Gaussian with a known mean and covariance matrix. However, the electrical impedance tomography inverse problem is nonlinear, and hence no closed form expression exists for the posterior density. The typical approach for such problems is to sample from the posterior and then use the samples to compute statistics (s...

High-dimensional inverse problems present a challenge for Markov chain Monte Carlo (MCMC)-type sampling schemes. Typically, they rely on finding an efficient proposal distribution, which can be dif...

In recent years, adaptive Markov Chain Monte Carlo (MCMC) methods have become a standard tool for Bayesian parameter estimation. In adaptive MCMC, the past iterations are used to tune the proposal distribution of the algorithm. The same adaptation mechanisms can be used in Simulated Annealing (SA), a popular optimization method based on MCMC. The difficulty in using adaptation directly in SA is that the target function changes along the iterations in the annealing process, and the adaptation should keep up with the annealing. In this paper, a mechanism for automatically tuning the proposal distribution in SA is proposed. The approach is based on the Adaptive Metropolis algorithm of Haario et al. (Bernoulli 7(2):223---242, 2001), combined with a weighting mechanism to account for the cooling target. The proposed adaptation mechanism does not add any computational complexity to the problem in terms of objective function evaluations. The effect of adaptation is demonstrated using two benchmark problems, showing that the proposed adaptation mechanism can significantly improve optimization results compared to non-adaptive SA. The approach is presented for continuous optimization problems and generalization to integer and mixed-integer problems is a topic of future research.

The ensemble Kalman filter (EnKF) is a technique for dynamic state estimation. EnKF approximates the standard extended Kalman filter (EKF) by creating an ensemble of model states whose mean and empirical covariance are then used within the EKF formulas. The technique has a number of advantages for large-scale, nonlinear problems. First, large-scale covariance matrices required within EKF are replaced by low rank and low storage approximations, making implementation of EnKF more efficient. Moreover, for a nonlinear state space model, implementation of EKF requires the associated tangent linear and adjoint codes, while implementation of EnKF does not. However, for EnKF to be effective, the choice of the ensemble members is extremely important. In this paper, we show how to use the conjugate gradient (CG) method, and the recently introduced CG sampler, to create the ensemble members at each filtering step. This requires the use of a variational formulation of EKF. The effectiveness of the method is demonstrated on both a large-scale linear, and a small-scale, non-linear, chaotic problem. In our examples, the CG-EnKF performs better than the standard EnKF, especially when the ensemble size is small.

Prior distributions for Bayesian inference that rely on the $l_1$-norm of the parameters are of considerable interest, in part because they promote parameter fields with less regularity than Gaussian priors (e.g., discontinuities and blockiness). These $l_1$-type priors include the total variation (TV) prior and the Besov $B^s_{1,1}$ space prior, and in general yield non-Gaussian posterior distributions. Sampling from these posteriors is challenging, particularly in the inverse problem setting where the parameter space is high-dimensional and the forward problem may be nonlinear. This paper extends the randomize-then-optimize (RTO) method, an optimization-based sampling algorithm developed for Bayesian inverse problems with Gaussian priors, to inverse problems with $l_1$-type priors. We use a variable transformation to convert an $l_1$-type prior to a standard Gaussian prior, such that the posterior distribution of the transformed parameters is amenable to Metropolized sampling via RTO. We demonstrate this approach on several deconvolution problems and an elliptic PDE inverse problem, using TV or Besov $B^s_{1,1}$ space priors. Our results show that the transformed RTO algorithm characterizes the correct posterior distribution and can be more efficient than other sampling algorithms. The variable transformation can also be extended to other non-Gaussian priors., Comment: Preprint 24 pages, 13 figures. v1 submitted to SIAM Journal on Scientific Computing on July 5, 2016

In model-based process optimization one uses a mathematical model to optimize a certain criterion, for example the product yield of a chemical process. Models often contain parameters that have to be estimated from data. Typically, a point estimate (e.g. the least squares estimate) is used to fix the model for the optimization stage. However, parameter estimates are uncertain due to incomplete and noisy data. In this article, it is shown how parameter uncertainty can be taken into account in process optimization. To quantify the uncertainty, Markov Chain Monte Carlo (MCMC) sampling, an emerging standard approach in Bayesian estimation, is used. In the Bayesian approach, the solution to the parameter estimation problem is given as a distribution, and the optimization criteria are functions of that distribution. The formulation and implementation of the optimization is studied, and numerical examples are used to show that parameter uncertainty can have a large effect in optimization results.

Classical optimal design criteria suffer from two major flaws when applied to nonlinear problems. First, they are based on linearizing the model around a point estimate of the unknown parameter and therefore depend on the uncertain value of that parameter. Second, classical design methods are unavailable in ill-posed estimation situations, where previous data lack the information needed to properly construct the design criteria. Bayesian optimal design can, in principle, solve these problems. However, Bayesian design methods are not widely applied, mainly due to the fact that standard implementations for efficient and robust routine use are not available. In this article, we point out a concrete recipe for implementing Bayesian optimal design, based on the concept of simulation-based design introduced by Muller, Sanso, and De Iorio (2004). We develop further a predictive variance criterion and introduce an importance weighting mechanism for efficient computation of the variances. The simulation-based appr...

SUMMARY The Kalman filter is a technique for estimating a time-varying state given a dynamical model for and indirect measurements of the state. It is used, for example, on the control problems associated with a variety of navigation systems. Even in the case of nonlinear state and/or measurement models, standard implementations require only linear algebra. However, for sufficiently large-scale problems, such as arise in weather forecasting and oceanography, the matrix inversion and storage requirements of the Kalman filter are prohibitive, and hence, approximations must be made. In this paper, we describe how the conjugate gradient iteration can be used within the Kalman filter for quadratic minimization, as well as for obtaining low-rank approximations of the covariance and inverse-covariance matrices required for its implementation. The approach requires that we exploit the connection between the conjugate gradient and Lanczos iterations. Copyright © 2011 John Wiley & Sons, Ltd.

Climate models contain closure parameters to which the model climate is sensitive. These parameters appear in physical parameterization schemes where some unresolved variables are expressed by predefined parameters rather than being explicitly modeled. Currently, best expert knowledge is used to define the optimal closure parameter values, based on observations, process studies, large eddy simulations, etc. Here, parameter estimation, based on the adaptive Markov chain Monte Carlo (MCMC) method, is applied for estimation of joint posterior probability density of a small number (n=4) of closure parameters appearing in the ECHAM5 climate model. The parameters considered are related to clouds and precipitation and they are sampled by an adaptive random walk process of the MCMC. The parameter probability densities are estimated simultaneously for all parameters, subject to an objective function. Five alternative formulations of the objective function are tested, all related to the net radiative flux at the top of the atmosphere. Conclusions of the closure parameter estimation tests with a low-resolution ECHAM5 climate model indicate that (i) adaptive MCMC is a viable option for parameter estimation in large-scale computational models, and (ii) choice of the objective function is crucial for the identifiability of the parameter distributions.

In this work, we consider the Bayesian optimization (BO) approach for tuning parameters of complex chaotic systems. Such problems arise, for instance, in tuning the sub-grid scale parameterizations in weather and climate models. For such problems, the tuning procedure is generally based on a performance metric which measures how well the tuned model fits the data. This tuning is often a computationally expensive task. We show that BO, as a tool for finding the extrema of computationally expensive objective functions, is suitable for such tuning tasks. In the experiments, we consider tuning parameters of two systems: a simplified atmospheric model and a low-dimensional chaotic system. We show that BO is able to tune parameters of both the systems with a low number of objective function evaluations and without the need of any gradient information.

The emergence of Markov chain Monte Carlo (MCMC) methods has opened a way for Bayesian analysis of complex models. Running MCMC samplers typically requires thousands of model evaluations, which can exceed available computer resources when this evaluation is computationally intensive. We will discuss two generally applicable techniques to improve the efficiency of MCMC. First, we consider a parallel version of the adaptive MCMC algorithm of Haario et al. (2001), implementing the idea of inter-chain adaptation introduced by Craiu et al. (2009). Second, we present an early rejection (ER) approach, where model simulation is stopped as soon as one can conclude that the proposed parameter value will be rejected by the MCMC algorithm. ¶ This work is motivated by practical needs in estimating parameters of climate and Earth system models. These computationally intensive models involve non-linear expressions of the geophysical and biogeochemical processes of the Earth system. Modeling of these processes, especially those operating in scales smaller than the model grid, involves a number of specified parameters, or ‘tunables’. MCMC methods are applicable for estimation of these parameters, but they are computationally very demanding. Efficient MCMC variants are thus needed to obtain reliable results in reasonable time. Here we evaluate the computational gains attainable through parallel adaptive MCMC and Early Rejection using both simple examples and a realistic climate model.