Stochastic Engine Final Report: Applying Markov Chain Monte Carlo Methods with Importance Sampling to Large-Scale Data-Driven Simulation

We present a new technique for statistical static timing analysis (SSTA) based on Markov chain Monte Carlo (MCMC), that allows fast and accurate estimation of the right-hand tail of the delay distribution. A ¿naive¿ MCMC approach is inadequate for SSTA. Several modifications and enhancements, presented in this paper, enable application of MCMC to SSTA. Moreover, such an approach overcomes inherent limitations of techniques such as importance sampling and Quasi-Monte Carlo. Our results on open source designs, with an independent delay variation model, demonstrate that our technique can obtain more than an order of magnitude improvement in computation time over simple Monte Carlo, given an estimation accuracy target at a point in the tail. Our approach works by providing a large number of samples in the region of interest. Open problems include extension of algorithm applicability to a broader class of synthesis conditions, and handling of correlated delay variations. In a broader context, this work aims to show that MCMC and associated techniques can be useful in rare event analyses related to circuits, particularly for high-dimensional problems.

Probabilistic inverse techniques, like the Markov Chain Monte Carlo (MCMC) algorithm, have had recent success in combining disparate data types into a consistent model. The Stochastic Engine (SE) initiative was a technique that developed this method and applied it to a number of earth science and national security applications. For instance, while the method was originally developed to solve ground flow problems (Aines et al.), it has also been applied to atmospheric modeling and engineering problems. The investigators of this proposal have applied the SE to regional-scale lithospheric earth models, which have applications to hazard analysis and nuclear explosion monitoring. While this broad applicability is appealing, tailoring the method for each application is inefficient and time-consuming. Stochastic methods invert data by probabilistically sampling the model space and comparing observations predicted by the proposed model to observed data and preferentially accepting models that produce a good fit, generating a posterior distribution. In other words, the method ''inverts'' for a model or, more precisely, a distribution of models, by a series of forward calculations. While powerful, the technique is often challenging to implement, as the mapping from model space to data needs to be ''customized'' for each data type. For example, all proposed models might need to be transformed through sensitivity kernels from 3-D models to 2-D models in one step in order to compute path integrals, and transformed in a completely different manner in the next step. We seek technical enhancements that widen the applicability of the Stochastic Engine by generalizing some aspects of the method (i.e. model-to-data transformation types, configuration, model representation). Initially, we wish to generalize the transformations that are necessary to match the observations to proposed models. These transformations are sufficiently general not to pertain to any single application. This is a new and innovative approach to the problem, providing a framework to increase the efficiency of its implementation. The overall goal is to reduce response time and make the approach as ''plug-and-play'' as possible, and will result in the rapid accumulation of new data types for a host of both earth science and non-earth science problems.

<p>For strongly non-linear inverse problems, Markov chain Monte Carlo (MCMC) methods may fail to properly explore the posterior probability density function (PDF). Particle methods are very well suited for parallelization and offer an alternative approach whereby the posterior PDF is approximated using the states and weights of a population of evolving particles. In addition, it provides reliable estimates of the evidence (marginal likelihood) that is needed for Bayesian model selection at essentially no cost. We consider adaptive sequential Monte Carlo (ASMC), which is an extension of annealed importance sampling (AIS). In these methods, importance sampling is performed over a sequence of intermediate distributions, known as power posteriors, linking the prior to the posterior PDF. The main advantages of ASMC with respect to AIS are that it adaptively tunes the tempering between neighboring distributions and it performs resampling of particles when the variance of the particle weights becomes too large. We consider a challenging synthetic groundwater transport inverse problem with a categorical channelized 2D hydraulic conductivity field designed such that the posterior facies distribution includes two distinct modes with equal probability. The model proposals are obtained by iteratively re-simulating a fraction of the current model using conditional multi-point statistics (MPS) simulations. We focus here on the ability of ASMC to explore the posterior PDF and compare it with previously published results obtained with parallel tempering (PT), a state-of-the-art MCMC inversion approach that runs multiple interacting chains targeting different power posteriors. For a similar computational budget involving 24 particles for ASMC and 24 chains for PT, the ASMC implementation outperforms the results obtained by PT: the models fit the data better and the reference likelihood value is contained in the ASMC sampled likelihood range, while this is not the case for PT range. Moreover, we show that ASMC recovers both reference modes, while none of them is recovered by PT. However, with 24 particles there is one of the modes that has a higher weight than the other while the approximation is improved when moving to a larger number of particles. As a future development, we suggest that including fast surrogate modeling (e.g., polynomial chaos expansion) within ASMC for the MCMC steps used to evolve the particles in-between importance sampling steps would strongly reduce the computational cost while still ensuring results of similar quality as the importance sampling steps could still be performed using the regular more costly forward solver.</p>

A natural method for approximating out-of-sample predictive evaluation is leave-one-out cross-validation (LOOCV) --- we alternately hold out each case from a full data set and then train a Bayesian model using Markov chain Monte Carlo (MCMC) without the held-out; at last we evaluate the posterior predictive distribution of all cases with their actual observations. However, actual LOOCV is time-consuming. This paper introduces two methods, namely iIS and iWAIC, for approximating LOOCV with only Markov chain samples simulated from a posterior based on a \textit{full} data set. iIS and iWAIC aim at improving the approximations given by importance sampling (IS) and WAIC in Bayesian models with possibly correlated latent variables. In iIS and iWAIC, we first integrate the predictive density over the distribution of the latent variables associated with the held-out without reference to its observation, then apply IS and WAIC approximations to the integrated predictive density. We compare iIS and iWAIC with other approximation methods in three real data examples that respectively use mixture models, models with correlated spatial effects, and a random effect logistic model. Our empirical results show that iIS and iWAIC give substantially better approximates than non-integrated IS and WAIC and other methods.

Markov chain Monte Carlo (MCMC) algorithms have become powerful tools for Bayesian inference. However, they do not scale well to large-data problems. Divide-and-conquer strategies, which split the data into batches and, for each batch, run independent MCMC algorithms targeting the corresponding subposterior, can spread the computational burden across a number of separate computer cores. The challenge with such strategies is in recombining the subposteriors to approximate the full posterior. By creating a Gaussian-process approximation for each log-subposterior density we create a tractable approximation for the full posterior. This approximation is exploited through three methodologies: firstly a Hamiltonian Monte Carlo algorithm targeting the expectation of the posterior density provides a sample from an approximation to the posterior; secondly, evaluating the true posterior at the sampled points leads to an importance sampler that, asymptotically, targets the true posterior expectations; finally, an alternative importance sampler uses the full Gaussian-process distribution of the approximation to the log-posterior density to re-weight any initial sample and provide both an estimate of the posterior expectation and a measure of the uncertainty in it.

For strongly non-linear and high-dimensional inverse problems, Markov chain Monte Carlo (MCMC) methods may fail to properly explore the posterior probability density function (PDF) given a realistic computational budget and are generally poorly amenable to parallelization. Particle methods approximate the posterior PDF using the states and weights of a population of evolving particles and they are very well suited to parallelization. We focus on adaptive sequential Monte Carlo (ASMC), an extension of annealed importance sampling (AIS). In AIS and ASMC, importance sampling is performed over a sequence of intermediate distributions, known as power posteriors, linking the prior to the posterior PDF. The AIS and ASMC algorithms also provide estimates of the evidence (marginal likelihood) as needed for Bayesian model selection, at basically no additional cost. ASMC performs better than AIS as it adaptively tunes the tempering schedule and performs resampling of particles when the variance of the particle weights becomes too large. We consider a challenging synthetic groundwater transport inverse problem with a categorical channelized 2D hydraulic conductivity field defined such that the posterior facies distribution includes two distinct modes. The model proposals are obtained by iteratively re-simulating a fraction of the current model using conditional multiple-point statistics (MPS) simulations. We examine how ASMC explores the posterior PDF and compare with results obtained with parallel tempering (PT), a state-of-the-art MCMC inversion approach that runs multiple interacting chains targeting different power posteriors. For a similar computational budget, ASMC outperforms PT as the ASMC-derived models fit the data better and recover the reference likelihood. Moreover, we show that ASMC partly retrieves both posterior modes, while none of them is recovered by PT. Lastly, we demonstrate how the power posteriors obtained by ASMC can be used to assess the influence of the assumed data errors on the posterior means and variances, as well as on the evidence. We suggest that ASMC can advantageously replace MCMC for solving many challenging inverse problems arising in the field of water resources.

Group Importance Sampling for particle filtering and MCMC

We present a new Bayesian method for estimating demographic and phylogenetic history using population genomic data. Several key innovations are introduced that allow the study of diverse models within an Isolation-with-Migration framework. The new method implements a 2-step analysis, with an initial Markov chain Monte Carlo (MCMC) phase that samples simple coalescent trees, followed by the calculation of the joint posterior density for the parameters of a demographic model. In step 1, the MCMC sampling phase, the method uses a reduced state space, consisting of coalescent trees without migration paths, and a simple importance sampling distribution without the demography of interest. Once obtained, a single sample of trees can be used in step 2 to calculate the joint posterior density for model parameters under multiple diverse demographic models, without having to repeat MCMC runs. Because migration paths are not included in the state space of the MCMC phase, but rather are handled by analytic integration in step 2 of the analysis, the method is scalable to a large number of loci with excellent MCMC mixing properties. With an implementation of the new method in the computer program MIST, we demonstrate the method's accuracy, scalability, and other advantages using simulated data and DNA sequences of two common chimpanzee subspecies: Pan troglodytes (P. t.) troglodytes and P. t. verus.

Bayesian statistics has attracted people's interest again recently because of the application of Markov Chain Monte Carlo (MCMC) theory. In particle filtering, the diversification of particles disappears in the process of importance sampling. However, this problem can be solved using Metropolis-Hastings (MH) sampling usually used in MCMC theory. As a modification to added MCMC (AMCMC)—an improved MCMC particle filter that can track variable number of targets at the same time, a new approach to optimize those rejected samples in MH sampling process by Mean Shift algorithm is proposed in this paper. Because the operation rate of particles in AMCMC is increased, the circles of sampling needed for the convergence of Markov chain is reduced. It is shown by experiment that, the optimized algorithm has better tracking performance under the condition of fewer particles.

<p>Multiobservable thermochemical tomography (MTT) is a recent computational approach to obtain estimates of the physical state (e.g. temperature distribution, compositional structure and rock properties) for the upper mantle [1]. It allows to jointly invert multiple independent datasets (e.g. gravity, seismic, magnetotelluric) within a thermodynamically-constrained and fully probabilistic framework. Evaluation of the plausibility of different physical states of the mantle with Markov Chain Monte Carlo (MCMC) simulations requires the solution of complex forward problems (e.g. Stokes flow, Maxwell’s equations, etc.) millions of times, making MTT computationally demanding for large-scale inverse problems. Furthermore, the number of parameters in a global study can easily reach several millions, making it increasingly difficult to 1) locate the regions of high probability and 2) sample these regions appropriately.</p><p>In order to overcome these limitations, we have combined and implemented a number of techniques, such as reduced-order modelling and efficient parallelization of both the forward problems and the MCMC algorithms, which dramatically accelerate the solution of the forward problems. Our software, LitMod1D_4INV and LitMod3D_4INV, allow to compute a proposal in less than 1 second, even when solving multiple complex forward problems together. We develop a multi-level parallel MPI driver for a collection of advanced MCMC sampling strategies to locate and sample high-probability regions efficiently. The massive amounts of data generated by large-scale MTT inversions need to be managed efficiently. We output results to open-source freeware formats, such as HDF5, TileDB, designed for big data problems. We emphasize that our methods and approaches are not only useful for MTT, but for any demanding inverse problem.</p><p>In this contribution, we will present applications of our software to complex, large-scale MTT problems and discuss its benefits, limitations and future improvements.</p><p>[1] J.C. Afonso, N. Rawlinson, Y. Yang, D. L. Schutt, A. G. Jones, J. Fullea, W. L. Griffin, 3‐D multiobservable probabilistic inversion for the compositional and thermal structure of the lithosphere and upper mantle: III. Thermochemical tomography in the Western‐Central U.S., Journal of Geophysical Research, 121, doi:10.1002/2016JB013049, 2016</p>

Markov chain Monte Carlo methods, in particular, the Gibbs sampler, are widely used algorithms both in application and theoretical works in the classical and Bayesian paradigms. However, these algorithms are often computer intensive. Samawi et al. [Steady-state ranked Gibbs sampler. J. Stat. Comput. Simul. 2012;82(8), 1223–1238. doi:10.1080/00949655.2011.575378] demonstrate through theory and simulation that the dependent steady-state Gibbs sampler is more efficient and accurate in model parameter estimation than the original Gibbs sampler. This paper proposes the independent steady-state Gibbs sampler (ISSGS) approach to improve the original Gibbs sampler in multidimensional problems. It is demonstrated that ISSGS provides accuracy with unbiased estimation and improves the performance and convergence of the Gibbs sampler in multidimensional problems.

Although Markov chain Monte Carlo methods have been widely used in many disciplines, exact eigen analysis for such generated chains has been rare. In this paper, a special Metropolis-Hastings algorithm, Metropolized independent sampling, proposed first in Hastings (1970), is studied in full detail. The eigenvalues and eigenvectors of the corresponding Markov chain, as well as a sharp bound for the total variation distance between the nth updated distribution and the target distribution, are provided. Furthermore, the relationship between this scheme, rejection sampling, and importance sampling are studied with emphasis on their relative efficiencies. It is shown that Metropolized independent sampling is superior to rejection sampling in two respects: asymptotic efficiency and ease of computation.

This paper presents a new methodology for reliability-based inspection planning focusing on robust and accurate computational strategies for fatigue-reliability updating using inspection results. The core of the proposed strategy is a conditioned sampling-based method, implemented by a Fast Probability Analyzer (FPA) software where efficiency is achieved by using the importance sampling principal. For a single component or limit state, FPA first generates Markov-Chain Monte Carlo (MCMC) samples in the failure domain, then applies an adaptive stratified importance sampling (ASIS) method to compute probability of failure (PoF) with error control. Once the MCMC samples have been created, solving a reliability updating problem is fairly straightforward and computationally robust relative to the conventional system reliability methods that rely on linearization of the limit states. The new approach is demonstrated using examples including stiffened panels of a ship-shaped vessel where reliability is updated using inspection results from 100 panel connections.

Computation of failure probability via hierarchical clustering

Chapters 5, 13 and 14 presented methods for making inference about. infectious diseases from available data. This is of course one of the main motivations for modeling: learning about important features, such as R, the initial growth rate, potential outbreak sizes and what effect different control measures might have in the context of specific infections. The models considered in these chapters have all been simple enough to obtain more or less explicit estimates. of just a few relevant parameters. In more complicated and parameter-rich models, and/or when analyzing large data sets, it is usually impossible to estimate key model parameters explicitly. In such situations there are (at least) two ways to proceed. One uses Bayesian statistical inference by means of Markov chain Monte Carlo methods (MCMC), and the other uses large scale simulations along with numerical optimization to fit parameters to data. This chapter mainly describes Bayesian inference using MCMC and only briefly some large simulation methods.