Matthias Morzfeld

Alexandre J. Chorin, Matthias Morzfeld, Xuemin Tu

Implicit particle filters for data assimilation update the particles by first choosing probabilities and then looking for particle locations that assume them, guiding the particles one by one to the high probability domain. We provide a detailed description of these filters, with illustrative examples, together with new, more general, methods for solving the algebraic equations and with a new algorithm for parameter identification.

Matthias Morzfeld, Xuemin Tu, Ethan Atkins et al.

Implicit particle filters for data assimilation generate high-probability samples by representing each particle location as a separate function of a common reference variable. This representation requires that a certain underdetermined equation be solved for each particle and at each time an observation becomes available. We present a new implementation of implicit filters in which we find the solution of the equation via a random map. As examples, we assimilate data for a stochastically driven Lorenz system with sparse observations and for a stochastic Kuramoto-Sivashinski equation with observations that are sparse in both space and time.

Matthias Morzfeld, Daniel Hodyss, Chris Snyder

The ensemble Kalman filter (EnKF) is a reliable data assimilation tool for high-dimensional meteorological problems. On the other hand, the EnKF can be interpreted as a particle filter, and particle filters collapse in high-dimensional problems. We explain that these seemingly contradictory statements offer insights about how particle filters function in certain high-dimensional problems, and in particular support recent efforts in meteorology to "localize" particle filters, i.e., to restrict the influence of an observation to its neighborhood.

Matthias Morzfeld, Marcus S. Day, Ray W. Grout et al.

In parameter estimation problems one computes a posterior distribution over uncertain parameters defined jointly by a prior distribution, a model, and noisy data. Markov Chain Monte Carlo (MCMC) is often used for the numerical solution of such problems. An alternative to MCMC is importance sampling, which can exhibit near perfect scaling with the number of cores on high performance computing systems because samples are drawn independently. However, finding a suitable proposal distribution is a challenging task. Several sampling algorithms have been proposed over the past years that take an iterative approach to constructing a proposal distribution. We investigate the applicability of such algorithms by applying them to two realistic and challenging test problems, one in subsurface flow, and one in combustion modeling. More specifically, we implement importance sampling algorithms that iterate over the mean and covariance matrix of Gaussian or multivariate t-proposal distributions. Our implementation leverages massively parallel computers, and we present strategies to initialize the iterations using "coarse" MCMC runs or Gaussian mixture models.

John Bell, Marcus Day, Jonathan Goodman et al.

First-principles Markov Chain Monte Carlo sampling is used to investigate uncertainty quantification and uncertainty propagation in parameters describing hydrogen kinetics. Specifically, we sample the posterior distribution of thirty-one parameters focusing on the H2O2 and HO2 reactions resulting from conditioning on ninety-one experiments. Established literature values are used for the remaining parameters in the mechanism. The samples are computed using an affine invariant sampler starting with broad, noninformative priors. Autocorrelation analysis shows that O(1M) samples are sufficient to obtain a reasonable sampling of the posterior. The resulting distribution identifies strong positive and negative correlations and several non-Gaussian characteristics. Using samples drawn from the posterior, we investigate the impact of parameter uncertainty on the prediction of two more complex flames: a 2D premixed flame kernel and the ignition of a hydrogen jet issuing into a heated chamber. The former represents a combustion regime similar to the target experiments used to calibrate the mechanism and the latter represents a different combustion regime. For the premixed flame, the net amount of product after a given time interval has a standard deviation of less than 2% whereas the standard deviation of the ignition time for the jet is more than 10%. The samples used for these studies are posted online. These results indicate the degree to which parameters consistent with the target experiments constrain predicted behavior in different combustion regimes. This process provides a framework for both identifying reactions for further study from candidate mechanisms as well as combining uncertainty quantification and propagation to, ultimately, tie uncertainty in laboratory flame experiments to uncertainty in end-use numerical predictions of more complicated scenarios.

Matthias Morzfeld, Alexandre J. Chorin

Implicit particle filtering is a sequential Monte Carlo method for data assim- ilation, designed to keep the number of particles manageable by focussing attention on regions of large probability. These regions are found by min- imizing, for each particle, a scalar function F of the state variables. Some previous implementations of the implicit filter rely on finding the Hessians of these functions. The calculation of the Hessians can be cumbersome if the state dimension is large or if the underlying physics are such that derivatives of F are difficult to calculate. This is the case in many geophysical applica- tions, in particular for models with partial noise, i.e. with a singular state covariance matrix. Examples of models with partial noise include stochastic partial differential equations driven by spatially smooth noise processes and models for which uncertain dynamic equations are supplemented by con- servation laws with zero uncertainty. We make the implicit particle filter applicable to such situations by combining gradient descent minimization with random maps and show that the filter is efficient, accurate and reliable because it operates in a subspace whose dimension is smaller than the state dimension. As an example, we assimilate data for a system of nonlinear partial differential equations that appears in models of geomagnetism.

Matthias Morzfeld, Xin T. Tong, Youssef M. Marzouk

We investigate how ideas from covariance localization in numerical weather prediction can be used in Markov chain Monte Carlo (MCMC) sampling of high-dimensional posterior distributions arising in Bayesian inverse problems. To localize an inverse problem is to enforce an anticipated "local" structure by (i) neglecting small off-diagonal elements of the prior precision and covariance matrices; and (ii) restricting the influence of observations to their neighborhood. For linear problems we can specify the conditions under which posterior moments of the localized problem are close to those of the original problem. We explain physical interpretations of our assumptions about local structure and discuss the notion of high dimensionality in local problems, which is different from the usual notion of high dimensionality in function space MCMC. The Gibbs sampler is a natural choice of MCMC algorithm for localized inverse problems and we demonstrate that its convergence rate is independent of dimension for localized linear problems. Nonlinear problems can also be tackled efficiently by localization and, as a simple illustration of these ideas, we present a localized Metropolis-within-Gibbs sampler. Several linear and nonlinear numerical examples illustrate localization in the context of MCMC samplers for inverse problems.

Andrew Leach, Kevin K. Lin, Matthias Morzfeld

We study a class of importance sampling methods for stochastic differential equations (SDEs). A small-noise analysis is performed, and the results suggest that a simple symmetrization procedure can significantly improve the performance of our importance sampling schemes when the noise is not too large. We demonstrate that this is indeed the case for a number of linear and nonlinear examples. Potential applications, e.g., data assimilation, are discussed.

Xu-Hui Zhou, Lorenzo Beronilla, Michael K. Sleeman et al.

Data assimilation (DA) for compressible flows with shocks is challenging because many classical DA methods generate spurious oscillations and nonphysical features near uncertain shocks. We focus here on the ensemble Kalman filter (EnKF). We show that the poor performance of the standard EnKF may be attributed to the bimodal forecast distribution that can arise in the vicinity of an uncertain shock location; this violates the assumptions underpinning the EnKF, which assume a forecast which is close to Gaussian. To address this issue we introduce the new neural EnKF. The basic idea is to systematically embed neural function approximations within ensemble DA by mapping the forecast ensemble of shocked flows to the parameter space (weights and biases) of a deep neural network (NN) and to subsequently perform DA in that space. The nonlinear mapping encodes sharp and smooth flow features in an ensemble of NN parameters. Neural EnKF updates are therefore well-behaved only if the NN parameters vary smoothly within the neural representation of the forecast ensemble. We show that such a smooth variation of network parameters can be enforced via physics-informed transfer learning, and demonstrate that in so-doing the neural EnKF avoids the spurious oscillations and nonphysical features that plague the standard EnKF. The applicability of the neural EnKF is demonstrated through a series of systematic numerical experiments with an inviscid Burgers' equation, Sod's shock tube, and a two-dimensional blast wave.

Matthias Morzfeld, Xuemin Tu, Jon Wilkening et al.

Implicit sampling is a weighted sampling method that is used in data assimilation, where one sequentially updates estimates of the state of a stochastic model based on a stream of noisy or incomplete data. Here we describe how to use implicit sampling in parameter estimation problems, where the goal is to find parameters of a numerical model, e.g.~a partial differential equation (PDE), such that the output of the numerical model is compatible with (noisy) data. We use the Bayesian approach to parameter estimation, in which a posterior probability density describes the probability of the parameter conditioned on data and compute an empirical estimate of this posterior with implicit sampling. Our approach generates independent samples, so that some of the practical difficulties one encounters with Markov Chain Monte Carlo methods, e.g.~burn-in time or correlations among dependent samples, are avoided. We describe a new implementation of implicit sampling for parameter estimation problems that makes use of multiple grids (coarse to fine) and BFGS optimization coupled to adjoint equations for the required gradient calculations. The implementation is "dimension independent", in the sense that a well-defined finite dimensional subspace is sampled as the mesh used for discretization of the PDE is refined. We illustrate the algorithm with an example where we estimate a diffusion coefficient in an elliptic equation from sparse and noisy pressure measurements. In the example, dimension\slash mesh-independence is achieved via Karhunen-Loève expansions.

Fei Lu, Matthias Morzfeld, Xuemin Tu et al.

Polynomial chaos expansions are used to reduce the computational cost in the Bayesian solutions of inverse problems by creating a surrogate posterior that can be evaluated inexpensively. We show, by analysis and example, that when the data contain significant information beyond what is assumed in the prior, the surrogate posterior can be very different from the posterior, and the resulting estimates become inaccurate. One can improve the accuracy by adaptively increasing the order of the polynomial chaos, but the cost may increase too fast for this to be cost effective compared to Monte Carlo sampling without a surrogate posterior.

Jonathan Goodman, Kevin K. Lin, Matthias Morzfeld

Implicit samplers are algorithms for producing independent, weighted samples from multi-variate probability distributions. These are often applied in Bayesian data assimilation algorithms. We use Laplace asymptotic expansions to analyze two implicit samplers in the small noise regime. Our analysis suggests a symmetrization of the algo- rithms that leads to improved (implicit) sampling schemes at a rel- atively small additional cost. Computational experiments confirm the theory and show that symmetrization is effective for small noise sampling problems.